### R code from vignette source '10tests.Rnw' ## Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner ################################################### ### code chunk number 1: 10tests.Rnw:9-10 ################################################### source("R/startup.R") ################################################### ### code chunk number 2: 10tests.Rnw:13-16 ################################################### # library(envel) # superseded # requireversion(envel, "1.6-4") requireversion(spatstat, "1.42-2.023") ################################################### ### code chunk number 3: 10tests.Rnw:19-22 ################################################### spatstat.options(terse=2) Digits <- 3 options(digits=Digits) ################################################### ### code chunk number 4: 10tests.Rnw:145-159 ################################################### xx <- seq(-3,3,length=1024) showdens <- function(x=xx) { plot(x, dnorm(x), type="l", main="", lwd=2, ylim=c(0, dnorm(0) * 1.1), xlab="t", ylab="probability density", axes=FALSE, xaxs='i', yaxs='i') box() } shadeup <- function(x0, ..., x=xx) { xup <- x[x >= x0] polygon(c(xup, rev(xup)), c(dnorm(xup), rep(0, length(xup))), ...) segments(x0, 0, x0, dnorm(x0), lwd=2) } ################################################### ### code chunk number 5: fvSquat.Rnw:3-5 ################################################### newplot(6, 0.65) setmargins(0.5+c(3,3,0,0)) ################################################### ### code chunk number 6: 10tests.Rnw:164-165 ################################################### setmargins(3,1.5,0,0) ################################################### ### code chunk number 7: 10tests.Rnw:169-173 ################################################### showdens() shadeup(1.2, border=NA, col="grey") axis(1, at=1.2, labels=expression(t[obs])) text(1.6, dnorm(1.6)/2, expression(p)) ################################################### ### code chunk number 8: 10tests.Rnw:175-180 ################################################### showdens() shadeup(0.6, border=NA, col="grey") axis(1, at=0.6, labels=expression(t[crit])) axis(1, at=1.2, labels=expression(t[obs])) text(1.25, 0.07, expression(alpha)) ################################################### ### code chunk number 9: 10tests.Rnw:297-299 ################################################### mur <- lapply(murchison, rescale, s=1000, unitname="km") mur$dfault <- with(mur, distfun(faults)) ################################################### ### code chunk number 10: 10tests.Rnw:303-305 ################################################### mfit0 <- ppm(gold ~ greenstone, data=mur) mfit1 <- ppm(gold ~ greenstone + dfault, data=mur) ################################################### ### code chunk number 11: 10tests.Rnw:308-309 ################################################### mfit1 <- update(mfit0, . ~ . + dfault) ################################################### ### code chunk number 12: 10tests.Rnw:330-333 ################################################### Dist <- with(copper, distfun(SouthLines)) cfit0 <- ppm(SouthPoints ~ 1, data=copper) cfit1 <- ppm(SouthPoints ~ Dist, data=copper) ################################################### ### code chunk number 13: 10tests.Rnw:374-375 ################################################### anova(mfit0, mfit1, test="LRT") ################################################### ### code chunk number 14: 10tests.Rnw:396-397 ################################################### anova(cfit0, cfit1, test="LRT") ################################################### ### code chunk number 15: chorleyload ################################################### if(draftversion) { reload.or.compute(datafilepath("chorleyCoarse.rda"), { lung <- split(chorley)$lung larynx <- split(chorley)$larynx smo <- density(lung, sigma=0.15, eps=0.1) smo <- eval.im(pmax(smo, 1e-10)) Q <- quadscheme(larynx, eps=0.5) }) draftmessage <- "using coarse quadrature scheme, eps=0.5, in draft version" } else { reload.or.compute(datafilepath("chorley.rda"), { lung <- split(chorley)$lung larynx <- split(chorley)$larynx smo <- density(lung, sigma=0.15, eps=0.1) smo <- eval.im(pmax(smo, 1e-10)) Q <- quadscheme(larynx, eps=0.1) }) draftmessage <- NULL } ################################################### ### code chunk number 16: 10tests.Rnw:425-428 (eval = FALSE) ################################################### ## lung <- split(chorley)$lung ## larynx <- split(chorley)$larynx ## Q <- quadscheme(larynx, eps=0.1) ################################################### ### code chunk number 17: chorleyDfitAgain ################################################### d2incin <- function(x, y, xincin=354.5, yincin=413.6) { (x - xincin)^2 + (y - yincin)^2 } raisin <- function(x,y, alpha, beta) { 1 + alpha * exp( - beta * d2incin(x,y)) } chorleyDfit <- ippm(Q ~ offset(log(smo) + log(raisin)), start=list(alpha=5, beta=1)) ################################################### ### code chunk number 18: 10tests.Rnw:447-448 ################################################### chorleyDfit ################################################### ### code chunk number 19: 10tests.Rnw:452-453 ################################################### chorley0fit <- update(chorleyDfit, . ~ offset(log(smo))) ################################################### ### code chunk number 20: 10tests.Rnw:468-469 ################################################### anova(chorley0fit, chorleyDfit, test="LRT") ################################################### ### code chunk number 21: 10tests.Rnw:513-514 ################################################### options(digits=3) ################################################### ### code chunk number 22: 10tests.Rnw:520-521 ################################################### coef(summary(mfit1)) ################################################### ### code chunk number 23: 10tests.Rnw:527-528 ################################################### options(digits=Digits) ################################################### ### code chunk number 24: 10tests.Rnw:539-541 ################################################### V <- coef(summary(mfit1))["dfault", "Zval"] pnorm(V, lower.tail=TRUE) ################################################### ### code chunk number 25: 10tests.Rnw:601-602 ################################################### anova(mfit0, mfit1, test="Rao") ################################################### ### code chunk number 26: 10tests.Rnw:667-668 ################################################### anova(cfit0, cfit1, test="score") ################################################### ### code chunk number 27: 10tests.Rnw:705-707 ################################################### dcop <- distfun(copper$SouthLines) berman.test(cfit0, dcop) ################################################### ### code chunk number 28: 10tests.Rnw:903-906 ################################################### fit2e <- ppm(bei ~ polynom(elev,2), data=bei.extra) M <- quadrat.test(fit2e, nx=4, ny=2) M ################################################### ### code chunk number 29: Bei.Rnw:3-5 ################################################### newplot(12,0.8) setmargins(0) ################################################### ### code chunk number 30: 10tests.Rnw:915-917 ################################################### plot(bei, main="", cols="lightgrey", cex=0.5) plot(M, add=TRUE, col=if(monochrome) "black" else "red") ################################################### ### code chunk number 31: 10tests.Rnw:926-927 ################################################### rr <- round(range(residuals(M))) ################################################### ### code chunk number 32: 10tests.Rnw:971-972 (eval = FALSE) ################################################### ## quadrat.test(fit2e, nx=4, ny=2, CR=0) ################################################### ### code chunk number 33: 10tests.Rnw:985-987 ################################################### V <- quantess(Window(bei), bei.extra$grad, 4) quadrat.test(fit2e, tess=V) ################################################### ### code chunk number 34: 10tests.Rnw:1018-1020 ################################################### X <- bdspots[[2]][square(c(-200,200))] qX <- quadratcount(X, 20, 20) ################################################### ### code chunk number 35: 10tests.Rnw:1027-1029 ################################################### tX <- table(factor(as.numeric(qX), levels=0:max(qX))) tX ################################################### ### code chunk number 36: 10tests.Rnw:1035-1036 ################################################### options(digits=2) ################################################### ### code chunk number 37: 10tests.Rnw:1038-1040 ################################################### eX <- 400 * dpois(0:max(qX), mean(qX)) eX ################################################### ### code chunk number 38: 10tests.Rnw:1044-1051 ################################################### fX <- factor(as.numeric(qX), levels=0:max(qX)) fX <- mergeLevels(fX, ">=5"=5:7) tX <- table(fX) tX p04 <- dpois(0:4, mean(qX)) eX <- 400 * c(p04, 1 - sum(p04)) eX ################################################### ### code chunk number 39: 10tests.Rnw:1053-1054 ################################################### options(digits=Digits) ################################################### ### code chunk number 40: 10tests.Rnw:1057-1059 ################################################### X2 <- sum((tX - eX)^2/eX) (pval <- pchisq(X2, df=length(tX)-1, lower.tail=FALSE)) ################################################### ### code chunk number 41: 10tests.Rnw:1262-1264 ################################################### grad <- bei.extra$grad cdf.test(fit2e, grad, test="ks") ################################################### ### code chunk number 42: 10tests.Rnw:1351-1363 ################################################### set.seed(19141917) ## make list of 20 point patterns. Plist <- list() ## First entry is 'cells' data Plist[[1]] <- cells ## other entries are randomised for(i in 2:20) Plist[[i]] <- runifpoint(npoints(cells), Window(cells)) Plist <- as.solist(Plist) ## Evaluate Clark-Evans statistic for 5 patterns clarkP <- unlist(lapply(Plist, clarkevans, correction="none")) clarkP <- round(clarkP, 3) ################################################### ### code chunk number 43: 10tests.Rnw:1366-1367 ################################################### zeromargins() ################################################### ### code chunk number 44: 10tests.Rnw:1371-1374 ################################################### plot(Plist, main="", main.panel="", equal.scales=TRUE, mar.panel=0, hsep=1/2, vsep=1/2, nrows=4, pch=16, cex=0.8) ################################################### ### code chunk number 45: 10tests.Rnw:1384-1385 ################################################### set.seed(19141917) ################################################### ### code chunk number 46: 10tests.Rnw:1409-1412 ################################################### opa <- options(width=80) sort(clarkP[-1]) par(opa) ################################################### ### code chunk number 47: 10tests.Rnw:1578-1579 ################################################### (Tobs <- clarkevans(cells, correction="none")) ################################################### ### code chunk number 48: 10tests.Rnw:1585-1590 ################################################### Tsim <- numeric(19) for(i in 1:19) { X <- runifpoint(npoints(cells), Window(cells)) Tsim[i] <- clarkevans(X, correction="none") } ################################################### ### code chunk number 49: 10tests.Rnw:1600-1603 ################################################### # check this! if(!(Tobs > max(Tsim))) stop("Unusual MC test result- pick another seed") ################################################### ### code chunk number 50: 10tests.Rnw:1608-1609 ################################################### (preg <- (1 + sum(Tsim > Tobs))/(1 + length(Tsim))) ################################################### ### code chunk number 51: 10tests.Rnw:1612-1613 ################################################### (pclus <- (1 + sum(Tsim < Tobs))/(1 + length(Tsim))) ################################################### ### code chunk number 52: 10tests.Rnw:1616-1617 ################################################### (peither <- 2 * min(pclus, preg)) ################################################### ### code chunk number 53: 10tests.Rnw:1623-1626 (eval = FALSE) ################################################### ## quadrat.test(redwood, nx=5, alternative="clustered", ## method="MonteCarlo", nsim=999) ## clarkevans.test(redwood, alternative="clustered", nsim=999) ################################################### ### code chunk number 54: 10tests.Rnw:1741-1747 ################################################### xxx <- seq(0, 10, length=1001) scal <- 4 * sqrt(2 * pi) set.seed(103) xobs <- qnorm(0.9, mean=5) xsim <- rnorm(10, mean=5) ysim <- xobs + (xsim-5) ################################################### ### code chunk number 55: 10tests.Rnw:1752-1753 ################################################### setmargins(0) ################################################### ### code chunk number 56: 10tests.Rnw:1757-1791 ################################################### topline <- FALSE if(topline) { plot(c(0,22),c(-2,19),type="n", axes=FALSE, xlab="", ylab="", main="") } else { plot(c(0,22),c(-2,13),type="n", axes=FALSE, xlab="", ylab="", main="") } lines(xxx, scal * dnorm(xxx, mean=5)) points(xobs, 0, pch=16) points(xsim, rep(0, length(xsim)), pch="|") segments(0,6,10,6) points(xobs, 6, pch=16) lines(xxx, 8 + scal * dnorm(xxx, mean=5)) if(topline) { for(i in 2:8) lines(xxx, 14+scal*dnorm(xxx, mean=i), lty=3) for(i in 2:8) lines(12+xxx, 14+scal*dnorm(xxx, mean=i), lty=3) } for(i in 5:9) lines(12+xxx, 8+scal*dnorm(xxx, mean=i), lty=if(i==5) 1 else 3) segments(12,6,22,6) points(12+xobs, 6, pch=16) points(12+xobs, 0, pch=16) lines(12+xxx, scal*dnorm(xxx, mean=xobs), lty=2) points(12+ysim, rep(0,length(ysim)), pch="|") if(topline) { text(5, 18.5, "Model") text(17, 18.5, "Model") } text(5, 12.5, "Simple null hypothesis") text(17, 12.5, "Composite null hypothesis") text(5.5, 6.7, "Data") text(17.5, 6.7, "Data") text(2.3,3, "Null model") text(14,3, "Fitted null model") text(5, -1.25, "Simulations from null model") text(17, -1.25, "Simulations from fitted model") ################################################### ### code chunk number 57: 10tests.Rnw:1852-1861 ################################################### dotest <- function(X, nsim=19) { Tobs <- clarkevans(X, correction="none") Tsim <- numeric(nsim) for(i in 1:nsim) { Xsim <- rpoispp(intensity(X), win=Window(X)) Tsim[i] <- clarkevans(Xsim, correction="none") } return(Tobs > max(Tsim)) } ################################################### ### code chunk number 58: 10tests.Rnw:1869-1877 (eval = FALSE) ################################################### ## N <- 10000 ## lambda <- 50 ## rejected <- logical(N) ## for(i in 1:N) { ## X <- rpoispp(lambda) ## rejected[i] <- dotest(X) ## } ## alpha <- mean(rejected) ################################################### ### code chunk number 59: 10tests.Rnw:2025-2027 ################################################### set.seed(11111918) KE <- envelope(cells, Kest, nsim=39) ################################################### ### code chunk number 60: fv.Rnw:3-5 ################################################### newplot(6, 0.5) setmargins(0.5+c(3,3,1,0)) ################################################### ### code chunk number 61: 10tests.Rnw:2032-2033 ################################################### plot(KE, main="", lwd=4, legend=FALSE) ################################################### ### code chunk number 62: 10tests.Rnw:2162-2164 ################################################### mad.test(cells, Lest, nsim=99, ginterval=c(0, 0.2), fix.n=TRUE, verbose=FALSE) ################################################### ### code chunk number 63: 10tests.Rnw:2186-2190 ################################################### set.seed(28283) envKglob19 <- envelope(cells, Kest, nsim=19, fix.n=TRUE, global=TRUE, savepatterns=TRUE) envLglob19 <- envelope(envKglob19, Lest, global=TRUE, savefuns=TRUE) ################################################### ### code chunk number 64: fv2.Rnw:3-5 ################################################### newplot(12, 0.95) setmargins(0.5+c(3,3,0,1)) ################################################### ### code chunk number 65: 10tests.Rnw:2195-2198 ################################################### par(mfrow=c(1,2)) plot(envKglob19, main="", lwd=c(4,1,1,1), legend=FALSE) plot(envLglob19, main="", lwd=c(4,1,1,1), legend=FALSE) ################################################### ### code chunk number 66: 10tests.Rnw:2300-2301 ################################################### dclf.test(cells, Lest, nsim=99, fix.n=TRUE, verbose=FALSE) ################################################### ### code chunk number 67: 10tests.Rnw:2377-2380 ################################################### set.seed(4242422) dc <- dclf.test(cells, Lest, nsim=19, fix.n=TRUE, alternative="less", verbose=FALSE) ################################################### ### code chunk number 68: 10tests.Rnw:2382-2383 (eval = FALSE) ################################################### ## dclf.test(cells, Lest, nsim=19, fix.n=TRUE, alternative="less") ################################################### ### code chunk number 69: 10tests.Rnw:2385-2386 ################################################### dc ################################################### ### code chunk number 70: 10tests.Rnw:2439-2441 (eval = FALSE) ################################################### ## plot(envelope(redwood, Lest, nsim=39)) ## plot(envelope(redwood, Lest, nsim=39, global=TRUE)) ################################################### ### code chunk number 71: 10tests.Rnw:2444-2446 ################################################### EP <- envelope(redwood, Lest, nsim=39, savefuns=TRUE) EG <- envelope(EP, global=TRUE) ################################################### ### code chunk number 72: fv2.Rnw:3-5 ################################################### newplot(12, 0.95) setmargins(0.5+c(3,3,0,1)) ################################################### ### code chunk number 73: 10tests.Rnw:2451-2454 ################################################### par(mfrow=c(1,2)) plot(EP, main="", legend=FALSE) plot(EG, main="", legend=FALSE) ################################################### ### code chunk number 74: 10tests.Rnw:2592-2596 (eval = FALSE) ################################################### ## numata <- residualspaper$Fig1 ## fit <- ppm(numata ~ polynom(x,y,3)) ## E <- envelope(fit, Lest, nsim=19, global=TRUE, correction="border") ## plot(E) ################################################### ### code chunk number 75: 10tests.Rnw:2603-2605 (eval = FALSE) ################################################### ## B <- bdspots[[1]] ## fit <- ppm(B ~ 1, Hardcore()) ################################################### ### code chunk number 76: 10tests.Rnw:2609-2610 (eval = FALSE) ################################################### ## EH <- envelope(fit, Lest, nsim=39, global=TRUE) ################################################### ### code chunk number 77: 10tests.Rnw:2629-2630 (eval = FALSE) ################################################### ## e <- expression(rpoispp(100)) ################################################### ### code chunk number 78: 10tests.Rnw:2644-2649 (eval = FALSE) ################################################### ## n <- npoints(B) ## W <- Window(B) ## r <- min(nndist(B))*n/(n+1) ## e <- expression(rSSI(r, n, W)) ## E <- envelope(B, Lest, nsim=39, global=TRUE, simulate=e) ################################################### ### code chunk number 79: 10tests.Rnw:2664-2667 (eval = FALSE) ################################################### ## e <- expression(rlabel(amacrine)) ## E <- envelope(amacrine, Lcross, nsim=19, global=TRUE, simulate=e) ## plot(E) ################################################### ### code chunk number 80: 10tests.Rnw:2682-2684 (eval = FALSE) ################################################### ## Xlist <- runifpoint(42, nsim=99) ## envelope(cells, Kest, nsim=99, simulate=Xlist) ################################################### ### code chunk number 81: 10tests.Rnw:2689-2691 (eval = FALSE) ################################################### ## EK <- envelope(cells, Kest, nsim=99, savepatterns=TRUE) ## Ep <- envelope(cells, pcf, nsim=99, simulate=EK) ################################################### ### code chunk number 82: 10tests.Rnw:2702-2704 (eval = FALSE) ################################################### ## EK <- envelope(cells, Kest, nsim=99, savepatterns=TRUE) ## Ep <- envelope(EK, pcf) ################################################### ### code chunk number 83: 10tests.Rnw:2711-2713 (eval = FALSE) ################################################### ## A <- mad.test(cells, Lest, nsim=99, savefuns=TRUE) ## B <- dclf.test(A) ################################################### ### code chunk number 84: 10tests.Rnw:2774-2776 (eval = FALSE) ################################################### ## D <- density(X) ## envelope(X, Kinhom, simulate=expression(rpoispp(D))) ################################################### ### code chunk number 85: 10tests.Rnw:2789-2791 (eval = FALSE) ################################################### ## D <- density(X) ## envelope(X, Kinhom, lambda=D, simulate=expression(rpoispp(D))) ################################################### ### code chunk number 86: 10tests.Rnw:2840-2841 (eval = FALSE) ################################################### ## envelope(redwood, Lest, global=TRUE) ################################################### ### code chunk number 87: 10tests.Rnw:2880-2883 (eval = FALSE) ################################################### ## fisher <- function(x) { asin(sqrt(x)) } ## envelope(redwood, Fest, global=TRUE, ## transform=expression(fisher(.))) ################################################### ### code chunk number 88: 10tests.Rnw:2902-2903 (eval = FALSE) ################################################### ## envelope(cells, Kest, nsim=100, VARIANCE=TRUE) ################################################### ### code chunk number 89: 10tests.Rnw:2921-2923 ################################################### a1 <- envelope(cells, Lest, nsim=19, alternative="less") a2 <- envelope(cells, Lest, nsim=19, alternative="greater") ################################################### ### code chunk number 90: fv2.Rnw:3-5 ################################################### newplot(12, 0.95) setmargins(0.5+c(3,3,0,1)) ################################################### ### code chunk number 91: 10tests.Rnw:2928-2931 ################################################### plot(anylist(a1, a2), main="", main.panel="", equal.scales=TRUE, mar.panel=c(4,4,1,1), hsep=1, legend=FALSE) ################################################### ### code chunk number 92: 10tests.Rnw:2962-2965 (eval = FALSE) ################################################### ## E1 <- envelope(redwood, Kest, savepatterns=TRUE) ## E2 <- envelope(E1, Gest, global=TRUE, ## transform=expression(fisher(.))) ################################################### ### code chunk number 93: 10tests.Rnw:2972-2975 (eval = FALSE) ################################################### ## A1 <- envelope(redwood, Kest, nsim=39, savefuns=TRUE) ## A2 <- envelope(A1, global=TRUE, nsim=19, ## transform=expression(sqrt(./pi))) ################################################### ### code chunk number 94: 10tests.Rnw:2991-2994 (eval = FALSE) ################################################### ## E1 <- envelope(cells, Kest, nsim=10, savefuns=TRUE) ## E2 <- envelope(cells, Kest, nsim=20, savefuns=TRUE) ## E <- pool(E1, E2) ################################################### ### code chunk number 95: 10tests.Rnw:3023-3024 (eval = FALSE) ################################################### ## dclf.sigtrace(cells, Lest, nsim=19) ################################################### ### code chunk number 96: 10tests.Rnw:3026-3027 ################################################### a <- dclf.sigtrace(cells, Lest, nsim=19) ################################################### ### code chunk number 97: fv.Rnw:3-5 ################################################### newplot(6, 0.5) setmargins(0.5+c(3,3,1,0)) ################################################### ### code chunk number 98: 10tests.Rnw:3044-3045 ################################################### plot(a, main="") ################################################### ### code chunk number 99: 10tests.Rnw:3085-3089 ################################################### set.seed(1234567) X <- rescale(swedishpines) prog19 <- dclf.progress(X, Lest, nsim=19) prog1999 <- dclf.progress(X, Lest, nsim=1999, nrank=100) ################################################### ### code chunk number 100: fv2.Rnw:3-5 ################################################### newplot(12, 0.95) setmargins(0.5+c(3,3,0,1)) ################################################### ### code chunk number 101: 10tests.Rnw:3096-3099 ################################################### par(mfrow=c(1,2)) plot(prog19, main="", legend=FALSE) plot(prog1999, main="", legend=FALSE) ################################################### ### code chunk number 102: 10tests.Rnw:3129-3131 (eval = FALSE) ################################################### ## X <- rescale(swedishpines) ## plot(dclf.progress(X, Lest, nsim=19)) ################################################### ### code chunk number 103: 10tests.Rnw:3220-3221 ################################################### D <- dg.test(cells, Lest, nsim=19) ################################################### ### code chunk number 104: 10tests.Rnw:3223-3224 (eval = FALSE) ################################################### ## dg.test(cells, Lest, nsim=19) ################################################### ### code chunk number 105: 10tests.Rnw:3226-3227 ################################################### D ################################################### ### code chunk number 106: 10tests.Rnw:3232-3233 (eval = FALSE) ################################################### ## dg.test(cells, Lest, nsim=19, exponent=Inf) ################################################### ### code chunk number 107: 10tests.Rnw:3279-3280 ################################################### set.seed(180981) ################################################### ### code chunk number 108: 10tests.Rnw:3282-3285 ################################################### swp <- rescale(swedishpines) E <- dg.envelope(swp, Lest) Eo <- dg.envelope(swp, Lest, alternative="less") ################################################### ### code chunk number 109: fv2.Rnw:3-5 ################################################### newplot(12, 0.95) setmargins(0.5+c(3,3,0,1)) ################################################### ### code chunk number 110: 10tests.Rnw:3291-3294 ################################################### plot(anylist(E, Eo), main="", main.panel="", equal.scales=TRUE, mar.panel=c(4,4,1,1), hsep=1, legend=FALSE)