Spatial Point Patterns: Methodology and Applications with R
by Adrian Baddeley, Ege Rubak and Rolf Turner


R code for chapter 16

Below is the R code used to generate results and figures in chapter 16. The code is in a rather raw format extracted from the book manuscript files – please read the instructions for use if you haven’t done so yet. You can download the script here.

### R code from vignette source '16replicated.Rnw'
## Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner

###################################################
### code chunk number 1: 16replicated.Rnw:10-12
###################################################
source("R/startup.R")
source("R/short.output.R")


###################################################
### code chunk number 2: 16replicated.Rnw:15-19
###################################################
Terse <- 2
Digits <- 3
spatstat.options(terse=Terse)
options(digits=Digits)


###################################################
### code chunk number 3: 16replicated.Rnw:22-25
###################################################
#requireversion(book, "3.0-2")
requireversion(spatstat, "1.42-2.022")
#requireversion(multippm, "4.1-1")


###################################################
### code chunk number 4: 16replicated.Rnw:227-228
###################################################
spatstat.options(terse=0)


###################################################
### code chunk number 5: 16replicated.Rnw:236-237
###################################################
waterstriders


###################################################
### code chunk number 6: 16replicated.Rnw:243-244
###################################################
spatstat.options(terse=Terse)


###################################################
### code chunk number 7: 16replicated.Rnw:255-256
###################################################
sapply(waterstriders, npoints)


###################################################
### code chunk number 8: 16replicated.Rnw:259-260
###################################################
wins <- lapply(waterstriders, Window)


###################################################
### code chunk number 9: 16replicated.Rnw:281-282
###################################################
solist(rpoispp(100), rpoispp(100))


###################################################
### code chunk number 10: 16replicated.Rnw:299-300
###################################################
Y <- solapply(c(10, 30, 100), rpoispp)


###################################################
### code chunk number 11: 16replicated.Rnw:309-310
###################################################
K <- anylist(Kest(cells), Kest(redwood))


###################################################
### code chunk number 12: 16replicated.Rnw:316-317
###################################################
K <- as.anylist(lapply(waterstriders, Kest))


###################################################
### code chunk number 13: 16replicated.Rnw:324-325
###################################################
K <- anylapply(waterstriders, Kest)


###################################################
### code chunk number 14: 16replicated.Rnw:344-345 (eval = FALSE)
###################################################
## plot(K, main="", main.panel=letters[1:3], legend=FALSE)


###################################################
### code chunk number 15: fv3.Rnw:3-5
###################################################
newplot(12.5, 1.0)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 16: 16replicated.Rnw:363-364
###################################################
plot(K, main="", main.panel=letters[1:3], legend=FALSE)


###################################################
### code chunk number 17: 16replicated.Rnw:384-385
###################################################
D <- solapply(split(amacrine), density)


###################################################
### code chunk number 18: Amacrine2E.Rnw:3-5
###################################################
newplot(17,1)
setmargins(0.1)


###################################################
### code chunk number 19: 16replicated.Rnw:390-391
###################################################
plot(D, equal.ribbon=TRUE, main="", mar.panel=c(0,0,1,0), hsep=1)


###################################################
### code chunk number 20: 16replicated.Rnw:404-406 (eval = FALSE)
###################################################
## D <- solapply(split(amacrine), density)
## plot(D, equal.ribbon=TRUE, main="")


###################################################
### code chunk number 21: 16replicated.Rnw:515-520
###################################################
G <- hyperframe(X=c(0.23, 1.76, 3.14), T=anylist(sin,cos,tan),
                Y=letters[1:3], Z=factor(letters[1:3]),
                W=list(rpoispp(100),rpoispp(100), rpoispp(100)),
                U=42, V=rpoispp(100), stringsAsFactors=FALSE)
G


###################################################
### code chunk number 22: 16replicated.Rnw:536-537
###################################################
WS <- hyperframe(Striders=waterstriders)


###################################################
### code chunk number 23: 16replicated.Rnw:560-561
###################################################
G$W


###################################################
### code chunk number 24: 16replicated.Rnw:569-572
###################################################
G$U <- letters[24:26]
G$animal <- c("horse","dog","deer")
G


###################################################
### code chunk number 25: 16replicated.Rnw:575-579
###################################################
WS <- hyperframe()
WS$larvae <- waterstriders
WS$experiment <- factor(1:3)
WS


###################################################
### code chunk number 26: 16replicated.Rnw:593-594
###################################################
spatstat.options(terse=0)


###################################################
### code chunk number 27: 16replicated.Rnw:596-597
###################################################
G[3,2]


###################################################
### code chunk number 28: 16replicated.Rnw:599-600
###################################################
spatstat.options(terse=Terse)


###################################################
### code chunk number 29: 16replicated.Rnw:605-606
###################################################
G[3,2,drop=TRUE]


###################################################
### code chunk number 30: 16replicated.Rnw:629-630
###################################################
G[,2] <- anylist(sqrt,exp,log)


###################################################
### code chunk number 31: 16replicated.Rnw:653-654
###################################################
B <- subset(G, npoints(W) > 100)


###################################################
### code chunk number 32: 16replicated.Rnw:662-663
###################################################
B <- subset(G, minnndist(W) < 0.02, select=Z:V)


###################################################
### code chunk number 33: 16replicated.Rnw:674-675
###################################################
head(pyramidal)


###################################################
### code chunk number 34: 16replicated.Rnw:692-693
###################################################
pg <- split(pyramidal, pyramidal$group)


###################################################
### code chunk number 35: Unit5x2t.Rnw:4-5
###################################################
setmargins(0,0,1,0)


###################################################
### code chunk number 36: 16replicated.Rnw:714-715
###################################################
plot(simba, main="", nrows=2, mar.panel=c(0,0,1,0), hsep=1, vsep=1)


###################################################
### code chunk number 37: PromptOff.Rnw:1-2
###################################################
options(prompt="  ")


###################################################
### code chunk number 38: 16replicated.Rnw:736-737 (eval = FALSE)
###################################################
## plot(h, e)


###################################################
### code chunk number 39: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")


###################################################
### code chunk number 40: 16replicated.Rnw:743-745 (eval = FALSE)
###################################################
## plot(demohyper, 
##      quote({ plot(Image, main=""); plot(Points, add=TRUE) }))


###################################################
### code chunk number 41: Unit3R.Rnw:3-5
###################################################
newplot(22,0.9)
setmargins(0,0,0,1)


###################################################
### code chunk number 42: 16replicated.Rnw:753-756
###################################################
plot(demohyper, 
     quote({ plot(Image, main=""); plot(Points, add=TRUE) }),
     main="", parargs=list(mar=1.2 * c(1,1,1,2)))


###################################################
### code chunk number 43: 16replicated.Rnw:767-769 (eval = FALSE)
###################################################
## H <- hyperframe(Bugs=waterstriders)
## plot(H, quote(plot(Kest(Bugs))))


###################################################
### code chunk number 44: 16replicated.Rnw:787-789
###################################################
df <- data.frame(A=1:10, B=10:1)
with(df, A-B)


###################################################
### code chunk number 45: PromptOff.Rnw:1-2
###################################################
options(prompt="  ")


###################################################
### code chunk number 46: 16replicated.Rnw:800-801 (eval = FALSE)
###################################################
## with(h,e)


###################################################
### code chunk number 47: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")


###################################################
### code chunk number 48: 16replicated.Rnw:809-812
###################################################
H <- hyperframe(Bugs=waterstriders)
with(H, npoints(Bugs))
D <- with(H, distmap(Bugs))


###################################################
### code chunk number 49: 16replicated.Rnw:846-847
###################################################
simba$Dist <- with(simba, distmap(Points))


###################################################
### code chunk number 50: 16replicated.Rnw:856-857
###################################################
with(simba, npoints(Points))


###################################################
### code chunk number 51: 16replicated.Rnw:862-864
###################################################
H <- hyperframe(Gerris=waterstriders)
K <- with(H, Kest(Gerris))


###################################################
### code chunk number 52: 16replicated.Rnw:869-871
###################################################
H <- hyperframe(Gerris=waterstriders)
m <- with(H, nndist(Gerris))


###################################################
### code chunk number 53: 16replicated.Rnw:876-877
###################################################
with(H, min(nndist(Gerris)))


###################################################
### code chunk number 54: 16replicated.Rnw:882-886
###################################################
set.seed(190105)
lambda <- rexp(3, rate=1/50)
H <- hyperframe(lambda=lambda)
H$Points <- with(H, rpoispp(lambda))


###################################################
### code chunk number 55: 16replicated.Rnw:893-896 (eval = FALSE)
###################################################
## lambda <- rexp(3, rate=1/50)
## H <- hyperframe(lambda=lambda)
## H$Points <- with(H, rpoispp(lambda))


###################################################
### code chunk number 56: 16replicated.Rnw:899-900 (eval = FALSE)
###################################################
## plot(H, quote(plot(Points, main=lambda)))


###################################################
### code chunk number 57: 16replicated.Rnw:903-905 (eval = FALSE)
###################################################
## H$Title <- with(H, parse(text=paste("lambda==", signif(lambda, 3))))
## plot(H, quote(plot(Points, main=Title)))


###################################################
### code chunk number 58: Unit3t.Rnw:3-5
###################################################
newplot(16.1,0.9)
setmargins(0,0,1,0)


###################################################
### code chunk number 59: 16replicated.Rnw:914-915
###################################################
setmargins(0,0,2,0)


###################################################
### code chunk number 60: 16replicated.Rnw:920-926
###################################################
options(digits=6)
plot(H, quote(plot(Points, 
                   main=parse(text=paste("lambda==", signif(lambda, 3))))),
     marsize=1,
     main="")
options(digits=Digits)


###################################################
### code chunk number 61: 16replicated.Rnw:937-938
###################################################
H$X <- with(H, rpoispp(50))


###################################################
### code chunk number 62: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()


###################################################
### code chunk number 63: 16replicated.Rnw:1006-1009
###################################################
plot(bdspots, equal.scales=TRUE, 
     cex=0.3, pch=16, main="", main.panel="", 
     mar.panel=0, hsep=0.6)


###################################################
### code chunk number 64: OsteoFullT.Rnw:3-5
###################################################
newplot(9.5, 1)
setmargins(0)


###################################################
### code chunk number 65: 16replicated.Rnw:1055-1057
###################################################
plot(osteo$pts, main="", main.panel=as.character(osteo$shortid), 
     ncols=8, mar.panel=0.1+c(0,0,1,0), pch=21, bg='white')


###################################################
### code chunk number 66: 16replicated.Rnw:1074-1075
###################################################
head(osteo)


###################################################
### code chunk number 67: 16replicated.Rnw:1105-1107
###################################################
pyr <- pyramidal
levels(pyr$group) <- c("c", "af", "ph")


###################################################
### code chunk number 68: Unit4x8t.Rnw:3-5
###################################################
newplot(12, 1)
setmargins(0)


###################################################
### code chunk number 69: 16replicated.Rnw:1112-1115
###################################################
plot(pyr$Neurons, main.panel=pyr$group, main="", 
     ncols=8, mar.panel=0, hsep=1, vsep=2, equal.scales=TRUE,
     pch=16, cex=0.7, cex.main=0.8)


###################################################
### code chunk number 70: Unit2L.Rnw:3-5
###################################################
newplot(13, 0.9)
setmargins(0,4,0,0)


###################################################
### code chunk number 71: 16replicated.Rnw:1210-1218
###################################################
pa <- function(i){ if(i == 1) 
                     list(chars=c(16,3), cex=c(1,0.6), 
                          leg.args=list(cex=1)) else
                   list(chars=c(16,7)) }
plot(flu[c(12, 24), 1], 
     panel.args=pa,
     main="", main.panel="", equal.scales=TRUE,
     mar.panel=0, hsep=1)


###################################################
### code chunk number 72: 16replicated.Rnw:1269-1273
###################################################
py <- pyramidal
py$n <- with(py, npoints(Neurons))
py$area <- with(py, area(Neurons))
py <- as.data.frame(py, warn=FALSE)


###################################################
### code chunk number 73: fvSquat.Rnw:3-5
###################################################
newplot(6, 0.65)
setmargins(0.5+c(3,3,0,0))


###################################################
### code chunk number 74: 16replicated.Rnw:1281-1282
###################################################
setmargins(0.5+c(3,3,0,3))


###################################################
### code chunk number 75: 16replicated.Rnw:1286-1290
###################################################
plot(sqrt(n) ~ group, data=py)
v <- pretty(py$n)
axis(4, at=sqrt(v), labels=v)
mtext(side=4, line=2, expression(italic(n)))


###################################################
### code chunk number 76: 16replicated.Rnw:1314-1315
###################################################
sapply(split(py$n/py$area, py$group), mean)


###################################################
### code chunk number 77: 16replicated.Rnw:1319-1322
###################################################
ntot <- sapply(split(py$n, py$group), sum)
atot <- sapply(split(py$area, py$group), sum)
ntot/atot


###################################################
### code chunk number 78: 16replicated.Rnw:1326-1330
###################################################
fitn <- glm(n ~ offset(log(area)) + group, 
           family=poisson, data=py)
newd <- data.frame(area=1, group=levels(py$group))
predict(fitn, newdata=newd, type="response")


###################################################
### code chunk number 79: 16replicated.Rnw:1335-1336 (eval = FALSE)
###################################################
## anova(fitn, test="LRT")


###################################################
### code chunk number 80: 16replicated.Rnw:1338-1339
###################################################
skipblanklines(anova(fitn, test="LRT"))


###################################################
### code chunk number 81: 16replicated.Rnw:1352-1354
###################################################
py$ce <- with(pyramidal, clarkevans(Neurons, correction="Donnelly"))
sapply(split(py$ce, py$group), mean)


###################################################
### code chunk number 82: 16replicated.Rnw:1363-1365 (eval = FALSE)
###################################################
## py$z <- with(pyramidal, clarkevans.test(Neurons)$statistic)
## anova(lm(z ~ group, data=py))


###################################################
### code chunk number 83: 16replicated.Rnw:1367-1369
###################################################
py$z <- with(pyramidal, clarkevans.test(Neurons)$statistic)
skipblanklines(anova(lm(z ~ group, data=py)))


###################################################
### code chunk number 84: 16replicated.Rnw:1372-1375 (eval = FALSE)
###################################################
## sdf <- as.data.frame(simba, warn=FALSE)
## sdf$z <- with(simba, clarkevans.test(Points)$statistic)
## anova(lm(z ~ group, data=sdf))


###################################################
### code chunk number 85: 16replicated.Rnw:1377-1380
###################################################
sdf <- as.data.frame(simba, warn=FALSE)
sdf$z <- with(simba, clarkevans.test(Points)$statistic)
skipblanklines(anova(lm(z ~ group, data=sdf)))


###################################################
### code chunk number 86: 16replicated.Rnw:1389-1391 (eval = FALSE)
###################################################
## plot(pyramidal, quote(plot(density(Neurons), main=group)))
## plot(with(pyramidal, density(Neurons)))


###################################################
### code chunk number 87: 16replicated.Rnw:1395-1396 (eval = FALSE)
###################################################
## plot(with(demohyper, rhohat(Points, Image)))


###################################################
### code chunk number 88: 16replicated.Rnw:1488-1490
###################################################
Keach <- lapply(waterstriders, Kest, ratio=TRUE)
Keach[[1]]


###################################################
### code chunk number 89: 16replicated.Rnw:1494-1495
###################################################
K <- pool(Keach[[1]], Keach[[2]], Keach[[3]])


###################################################
### code chunk number 90: 16replicated.Rnw:1498-1499
###################################################
K <- pool(as.anylist(Keach))


###################################################
### code chunk number 91: 16replicated.Rnw:1549-1551 (eval = FALSE)
###################################################
## plot(K, cbind(pooliso, pooltheo, loiso, hiiso) ~ r, 
##         shade=c("loiso", "hiiso"))


###################################################
### code chunk number 92: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))


###################################################
### code chunk number 93: 16replicated.Rnw:1563-1565
###################################################
plot(K, cbind(pooliso, pooltheo, loiso, hiiso) ~ r, 
     shade=c("loiso", "hiiso"), main="")


###################################################
### code chunk number 94: 16replicated.Rnw:1579-1582
###################################################
os <- osteo
os$K <- with(os, K3est(pts, ratio=TRUE))
Kanimal <- anylapply(split(os$K, os$id), pool)


###################################################
### code chunk number 95: fv3.Rnw:3-5
###################################################
newplot(12.5, 1.0)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 96: 16replicated.Rnw:1590-1593
###################################################
plot(Kanimal[-1], cbind(pooliso, pooltheo, loiso, hiiso)/1e4 ~ r, 
     shade=c("loiso", "hiiso"),
     main="", legend=FALSE, nrows=1)


###################################################
### code chunk number 97: 16replicated.Rnw:1614-1620
###################################################
pa <- pyramidal
pa$L <- with(pa, Lest(Neurons))
Lsplit <- split(pa$L, pa$group)
Leach <- anylapply(Lsplit, collapse.fv, 
                    same="theo", different="iso")
Lpool <- anylapply(Lsplit, pool)


###################################################
### code chunk number 98: fv3.Rnw:3-5
###################################################
newplot(12.5, 1.0)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 99: 16replicated.Rnw:1626-1630
###################################################
plot(Lpool, 
     cbind(pooliso,pooltheo, hiiso, loiso)-r~ r, 
     shade=c("hiiso", "loiso"), xlim=c(0, 0.2), 
     equal.scales=TRUE, legend=FALSE, main="")


###################################################
### code chunk number 100: 16replicated.Rnw:1643-1649 (eval = FALSE)
###################################################
## pa <- pyramidal
## pa$L <- with(pa, Lest(Neurons))
## Lsplit <- split(pa$L, pa$group)
## Lpool <- anylapply(Lsplit, pool)
## plot(Lpool, cbind(pooliso,pooltheo,hiiso,loiso) - r ~ r, 
##      shade=c("hiiso", "loiso"), xlim=c(0, 0.2), equal.scales=TRUE)


###################################################
### code chunk number 101: 16replicated.Rnw:1657-1660 (eval = FALSE)
###################################################
## Leach <- anylapply(Lsplit, collapse.fv, 
##                     same="theo", different="iso")
## plot(Leach, legend=FALSE, xlim=c(0, 0.2), ylim=c(0, 0.2))


###################################################
### code chunk number 102: fv3.Rnw:3-5
###################################################
newplot(12.5, 1.0)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 103: 16replicated.Rnw:1666-1669
###################################################
plot(Leach, xlim=c(0, 0.2), ylim=c(0, 0.2),
     equal.scales=TRUE, legend=FALSE, main="",
     col=1, lty=1)


###################################################
### code chunk number 104: 16replicated.Rnw:1719-1721
###################################################
set.seed(195505)
testpyramidal <- studpermu.test(pyramidal, Neurons ~ group)


###################################################
### code chunk number 105: 16replicated.Rnw:1723-1724 (eval = FALSE)
###################################################
## testpyramidal <- studpermu.test(pyramidal, Neurons ~ group)


###################################################
### code chunk number 106: 16replicated.Rnw:1726-1727
###################################################
testpyramidal


###################################################
### code chunk number 107: 16replicated.Rnw:1755-1761
###################################################
set.seed(12345)
sample1 <- rpoispp(100, nsim=10)
sample2 <- rMaternII(110, 0.02, nsim=10)
patterns <- list(Poisson = sample1, MaternII = sample2)
studpermu.test(patterns, summaryfunction = pcf, nperm=199, 
               interval=c(0,0.15))


###################################################
### code chunk number 108: 16replicated.Rnw:1792-1794
###################################################
set.seed(1809981)
stp <- studpermu.test(pyramidal, summaryfunction = Lest, use.Tbar = TRUE)


###################################################
### code chunk number 109: 16replicated.Rnw:1796-1797 (eval = FALSE)
###################################################
## studpermu.test(pyramidal, summaryfunction = Lest, use.Tbar = TRUE)


###################################################
### code chunk number 110: 16replicated.Rnw:1799-1800
###################################################
stp


###################################################
### code chunk number 111: 16replicated.Rnw:1815-1817
###################################################
waka2 <- split(cut(waka, breaks=c(0,20,200), 
                   labels = c("small", "large")))


###################################################
### code chunk number 112: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)


###################################################
### code chunk number 113: 16replicated.Rnw:1823-1824
###################################################
plot(waka2, main="", mar.panel=0, hsep=1, pch=16)


###################################################
### code chunk number 114: 16replicated.Rnw:1839-1841 (eval = FALSE)
###################################################
## waka2 <- split(cut(waka, breaks=c(0,20,200), 
##                    labels = c("small", "large")))


###################################################
### code chunk number 115: 16replicated.Rnw:1844-1845
###################################################
Lwak <- anylapply(waka2, Lest, correction="iso")


###################################################
### code chunk number 116: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))


###################################################
### code chunk number 117: 16replicated.Rnw:1855-1863
###################################################
plot(Lwak[[1]], main="",  . - r ~ r, 
  lty = c("dashed", "dotdash"), lwd=c(2,1), col =1,
  ylim=c(-.4, .8),  main="", legend = FALSE)
plot(Lwak[[2]], iso - r ~ r, lwd = 2, add=TRUE)
# legend("topright", 
#       c("big", "small", "Poisson"), 
#       lty=c("solid", "dashed", "dotdash"),
#       lwd=c(2,2,1), bty="n")


###################################################
### code chunk number 118: 16replicated.Rnw:1877-1879
###################################################
subwindows <- quadrats(waka, 3, 3)
treegroups <- lapply(waka2, split, f = subwindows)


###################################################
### code chunk number 119: 16replicated.Rnw:1890-1891
###################################################
set.seed(1989001)


###################################################
### code chunk number 120: 16replicated.Rnw:1893-1895
###################################################
wakatest <- studpermu.test(treegroups, summaryfunction = Lest, 
                           rinterval=c(5,15))


###################################################
### code chunk number 121: 16replicated.Rnw:1900-1901
###################################################
options(digits=5) # print.htest subtracts 3


###################################################
### code chunk number 122: 16replicated.Rnw:1903-1904
###################################################
wakatest


###################################################
### code chunk number 123: 16replicated.Rnw:1906-1907
###################################################
options(digits=Digits)


###################################################
### code chunk number 124: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))


###################################################
### code chunk number 125: 16replicated.Rnw:1937-1938
###################################################
plot(wakatest, . - r ~ r, lwd.mean=2)


###################################################
### code chunk number 126: cache
###################################################
sizes <- cut(marks(waka), breaks=c(0,20,200), 
            labels = c("small", "large"))
arbres <- nestsplit(waka, sizes, quadrats(waka, nx=3, ny=3))


###################################################
### code chunk number 127: 16replicated.Rnw:1972-1973 (eval = FALSE)
###################################################
## arbres <- nestsplit(waka, sizes, nx=3, ny=3)


###################################################
### code chunk number 128: 16replicated.Rnw:1976-1977
###################################################
set.seed(195507)


###################################################
### code chunk number 129: 16replicated.Rnw:1979-1980
###################################################
epreuveWaka <- studpermu.test(arbres, pts ~ f1)


###################################################
### code chunk number 130: 16replicated.Rnw:2023-2026
###################################################
bro <- unmark(bronzefilter)
mx <- median(coords(bro)$x)
halves <- chop.tess(Window(bro), infline(v=mx))


###################################################
### code chunk number 131: Bronze.Rnw:3-4
###################################################
setmargins(0)


###################################################
### code chunk number 132: 16replicated.Rnw:2039-2041
###################################################
plot(halves, main="", lwd=2)
plot(bro, add=TRUE, pch=19, cex=0.5)


###################################################
### code chunk number 133: 16replicated.Rnw:2052-2053
###################################################
Ks <- anylapply(split(bro, halves), Kscaled)


###################################################
### code chunk number 134: fv2.Rnw:3-5
###################################################
newplot(12, 0.95)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 135: 16replicated.Rnw:2060-2063
###################################################
#plot(collapse.fv(left=Ks[[1]], right=Ks[[2]], same="theo", different="iso"),
#     main="")
plot(Ks, main="", main.panel=c("left", "right"), legend=FALSE, xlim=c(0, 1.5))


###################################################
### code chunk number 136: 16replicated.Rnw:2074-2076
###################################################
bronzerho <- rhohat(bro, "x", method="tr")
bronzelambda <- predict(bronzerho)


###################################################
### code chunk number 137: fvSquat.Rnw:3-5
###################################################
newplot(6, 0.65)
setmargins(0.5+c(3,3,0,0))


###################################################
### code chunk number 138: 16replicated.Rnw:2088-2089
###################################################
plot(bronzerho, main="", ylim=c(0, 14))


###################################################
### code chunk number 139: 16replicated.Rnw:2100-2101
###################################################
Ki <- anylapply(split(bro, halves), Kinhom, lambda=bronzelambda)


###################################################
### code chunk number 140: fv2.Rnw:3-5
###################################################
newplot(12, 0.95)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 141: 16replicated.Rnw:2108-2111
###################################################
#plot(collapse.fv(left=Ki[[1]], right=Ki[[2]], same="theo", different="iso"),
#     main="")
plot(Ki, main="", main.panel=c("left", "right"), legend=FALSE, xlim=c(0,1))


###################################################
### code chunk number 142: 16replicated.Rnw:2134-2135
###################################################
bronze6 <- quantess(bro, "x", 6)


###################################################
### code chunk number 143: Bronze.Rnw:3-4
###################################################
setmargins(0)


###################################################
### code chunk number 144: 16replicated.Rnw:2142-2144
###################################################
plot(bronze6, main="")
plot(bro, add=TRUE, pch=19, cex=0.5)


###################################################
### code chunk number 145: 16replicated.Rnw:2154-2155
###################################################
b63 <- nestsplit(bro, bronze6, ny=3)


###################################################
### code chunk number 146: 16replicated.Rnw:2161-2162
###################################################
head(b63)


###################################################
### code chunk number 147: 16replicated.Rnw:2166-2167
###################################################
b63$inten <- factor(as.integer(b63$f1) <= 3, labels=c("Hi","Lo"))


###################################################
### code chunk number 148: 16replicated.Rnw:2170-2171
###################################################
set.seed(4242401)


###################################################
### code chunk number 149: 16replicated.Rnw:2173-2177
###################################################
locTest <- studpermu.test(b63, pts ~ inten, summaryfunction=Kscaled,
                          rinterval=c(0, 1.5))
corrTest <- studpermu.test(b63, pts ~ inten, summaryfunction=Kinhom,
                          lambda=bronzelambda, rinterval = c(0, 0.7))


###################################################
### code chunk number 150: 16replicated.Rnw:2180-2182
###################################################
locTest
corrTest


###################################################
### code chunk number 151: 16replicated.Rnw:2215-2216
###################################################
mppm(Points ~ 1, simba)


###################################################
### code chunk number 152: 16replicated.Rnw:2238-2239
###################################################
mppm(Points ~ group, simba)


###################################################
### code chunk number 153: 16replicated.Rnw:2243-2244
###################################################
mppm(Points ~ id, simba)


###################################################
### code chunk number 154: 16replicated.Rnw:2249-2250
###################################################
simba2 <- simba[c(FALSE,TRUE), ]


###################################################
### code chunk number 155: 16replicated.Rnw:2258-2259
###################################################
mppm(Points ~ Image, data=demohyper)


###################################################
### code chunk number 156: 16replicated.Rnw:2275-2277
###################################################
mppm(Points ~ Group/Image, data=demohyper)
mppm(Points ~ (Group-1)/Image, data=demohyper)


###################################################
### code chunk number 157: 16replicated.Rnw:2331-2332 (eval = FALSE)
###################################################
## mppm(Points ~ 1, simba, random = ~ 1 | id)


###################################################
### code chunk number 158: 16replicated.Rnw:2355-2367
###################################################
H <- hyperframe(P=waterstriders)
mod <- mppm(P ~ 1, data=H, random=~1|id)
## extract fixed effect intercept 
mu <- fixef(mod)
## extract random effect mean and variance 
ss <- summary(mod)
sigma <- exp(as.numeric(ss$ranef$reStruct$id))
lambda <- exp(mu + sigma^2/2)
#
mu     <- signif(mu, Digits)
sigma  <- signif(sigma, Digits)
lambda <- signif(lambda, Digits)


###################################################
### code chunk number 159: 16replicated.Rnw:2373-2375 (eval = FALSE)
###################################################
## H <- hyperframe(P=waterstriders)
## mppm(P ~ 1, data=H, random=~1|id)


###################################################
### code chunk number 160: 16replicated.Rnw:2377-2378
###################################################
mod


###################################################
### code chunk number 161: 16replicated.Rnw:2400-2401 (eval = FALSE)
###################################################
## mppm(Neurons ~ group, data=pyramidal, random=~1|id)


###################################################
### code chunk number 162: 16replicated.Rnw:2403-2404
###################################################
m <- mppm(Neurons ~ group, data=pyramidal, random=~1|id)


###################################################
### code chunk number 163: 16replicated.Rnw:2406-2408
###################################################
# Don't print all 31 random effects estimates..
short.output(m, 1-10, 36-60, excised="\n [...] \n ")


###################################################
### code chunk number 164: 16replicated.Rnw:2486-2488 (eval = FALSE)
###################################################
## fitpyr <-mppm(Neurons ~ group, data=pyramidal, random=~1|id) 
## anova(fitpyr)


###################################################
### code chunk number 165: 16replicated.Rnw:2490-2491
###################################################
anova(m)


###################################################
### code chunk number 166: 16replicated.Rnw:2524-2525
###################################################
mppm(Points ~ id, data=simba2, interaction=Strauss(0.07))


###################################################
### code chunk number 167: 16replicated.Rnw:2545-2548
###################################################
fiteach <- mppm(Points ~ id, data=simba2, interaction=Strauss(0.07),
                iformula = ~Interaction:id)
fiteach


###################################################
### code chunk number 168: 16replicated.Rnw:2558-2560
###################################################
fitgr <- mppm(Points ~ group, simba2, Strauss(0.07), 
              iformula = ~Interaction:group)


###################################################
### code chunk number 169: 16replicated.Rnw:2572-2573
###################################################
fitgr


###################################################
### code chunk number 170: 16replicated.Rnw:2575-2577
###################################################
co <- coef(fitgr)
si <- function(x) { signif(x, 4) }


###################################################
### code chunk number 171: 16replicated.Rnw:2579-2585
###################################################
# sanity check
expectednames <- c("(Intercept)", 
                   "grouptreatment",
                   "groupcontrol:Interaction",
                   "grouptreatment:Interaction")
stopifnot(identical(names(coef(fitgr)), expectednames))


###################################################
### code chunk number 172: 16replicated.Rnw:2589-2590
###################################################
coef(fitgr)


###################################################
### code chunk number 173: 16replicated.Rnw:2631-2632
###################################################
radii <- with(simba2, mean(nndist(Points)))


###################################################
### code chunk number 174: 16replicated.Rnw:2637-2639
###################################################
Rad <- hyperframe(R=radii)
Str <- with(Rad, Strauss(R))


###################################################
### code chunk number 175: 16replicated.Rnw:2642-2644
###################################################
Int <- hyperframe(str=Str)
mppm(Points ~ 1, simba2, interaction=Int)


###################################################
### code chunk number 176: 16replicated.Rnw:2658-2659 (eval = FALSE)
###################################################
## mppm(Points ~ id, simba2, interaction=Int, iformula = ~str:group)


###################################################
### code chunk number 177: 16replicated.Rnw:2669-2672
###################################################
h <- hyperframe(Y=waterstriders)
g <- hyperframe(po=Poisson(), str4 = Strauss(4), str7= Strauss(7))
mppm(Y ~ 1, data=h, interaction=g, iformula=~str4)


###################################################
### code chunk number 178: PromptOff.Rnw:1-2
###################################################
options(prompt="  ")


###################################################
### code chunk number 179: 16replicated.Rnw:2684-2685 (eval = FALSE)
###################################################
## interaction=hyperframe(po=Poisson(), str=Strauss(0.07))


###################################################
### code chunk number 180: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")


###################################################
### code chunk number 181: PromptOff.Rnw:1-2
###################################################
options(prompt="  ")


###################################################
### code chunk number 182: 16replicated.Rnw:2691-2692 (eval = FALSE)
###################################################
## iformula=~ifelse(group=="control", po, str)


###################################################
### code chunk number 183: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")


###################################################
### code chunk number 184: PromptOff.Rnw:1-2
###################################################
options(prompt="  ")


###################################################
### code chunk number 185: 16replicated.Rnw:2705-2706 (eval = FALSE)
###################################################
## iformula=~I((group=="control")*po) + I((group=="treatment") * str)


###################################################
### code chunk number 186: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")


###################################################
### code chunk number 187: 16replicated.Rnw:2716-2721
###################################################
g <- hyperframe(po=Poisson(), str=Strauss(0.07))
fit2 <- mppm(Points ~ 1, simba2, g, 
             iformula=~I((group=="control")*po) 
                     + I((group=="treatment") * str))
fit2


###################################################
### code chunk number 188: 16replicated.Rnw:2731-2735
###################################################
fit2a <- mppm(Points ~ 1, simba2, g, 
             iformula=~I((group=="treatment") * str))
fit2b <- mppm(Points ~ 1, simba2, Strauss(0.07), 
             iformula=~I((group=="treatment") * Interaction))


###################################################
### code chunk number 189: 16replicated.Rnw:2756-2758
###################################################
## temporarily suppress standard errors etc
op <- spatstat.options(print.ppm.SE="never")


###################################################
### code chunk number 190: 16replicated.Rnw:2766-2769
###################################################
H <- hyperframe(W=waterstriders)
fit <- mppm(W ~ 1, H)
subfits(fit)


###################################################
### code chunk number 191: 16replicated.Rnw:2794-2796
###################################################
fitI <- mppm(W ~ id, H)
subI <- subfits(fitI)


###################################################
### code chunk number 192: 16replicated.Rnw:2799-2800
###################################################
subII <- with(H, ppm(W ~ 1))


###################################################
### code chunk number 193: 16replicated.Rnw:2805-2807
###################################################
## reinstate standard errors
spatstat.options(op)


###################################################
### code chunk number 194: 16replicated.Rnw:2821-2823
###################################################
fit <- mppm(P ~ x, hyperframe(P=waterstriders))
res <- residuals(fit, type="Pearson")


###################################################
### code chunk number 195: 16replicated.Rnw:2832-2833
###################################################
smor <- with(hyperframe(res=res), Smooth(res, sigma=4))


###################################################
### code chunk number 196: Unit3R.Rnw:3-5
###################################################
newplot(22,0.9)
setmargins(0,0,0,1)


###################################################
### code chunk number 197: 16replicated.Rnw:2843-2844
###################################################
plot(smor, main="", main.panel="", equal.ribbon=TRUE)


###################################################
### code chunk number 198: 16replicated.Rnw:2857-2858
###################################################
sapply(res, integral.msr)


###################################################
### code chunk number 199: 16replicated.Rnw:2862-2866
###################################################
mod <- mppm(Neurons ~ group * x, data=pyramidal)
res <- residuals(mod, type="raw")
df <- as.data.frame(pyramidal, warn=FALSE)
df$resid <- sapply(res, integral.msr)


###################################################
### code chunk number 200: 16replicated.Rnw:2868-2869 (eval = FALSE)
###################################################
## plot(resid ~ group, df)


###################################################
### code chunk number 201: fvSquat.Rnw:3-5
###################################################
newplot(6, 0.65)
setmargins(0.5+c(3,3,0,0))


###################################################
### code chunk number 202: 16replicated.Rnw:2876-2877
###################################################
plot(resid~group, df)


###################################################
### code chunk number 203: 16replicated.Rnw:2891-2893
###################################################
fit <- mppm(P ~ 1, hyperframe(P=waterstriders))
subs <- hyperframe(Model=subfits(fit))


###################################################
### code chunk number 204: 16replicated.Rnw:2895-2896 (eval = FALSE)
###################################################
## plot(subs, quote(diagnose.ppm(Model)))


###################################################
### code chunk number 205: fv3.Rnw:3-5
###################################################
newplot(12.5, 1.0)
setmargins(0.5+c(3,3,0,1))


###################################################
### code chunk number 206: 16replicated.Rnw:2908-2911
###################################################
plot(subs, 
     quote(diagnose.ppm(Model, plot.neg="contour", xlab="", ylab="", rlab="")),
     main="")


###################################################
### code chunk number 207: 16replicated.Rnw:2927-2935
###################################################
H <- hyperframe(P = waterstriders)
fitall <- mppm(P ~ 1, H)
together <- subfits(fitall)
separate <- with(H, ppm(P))
Fits <- hyperframe(Together=together, Separate=separate)
dr <- with(Fits, unlist(coef(Separate)) - unlist(coef(Together)))
dr
exp(dr)


###################################################
### code chunk number 208: 16replicated.Rnw:2951-2958
###################################################
H <- hyperframe(X=waterstriders)
# Poisson with constant intensity for all patterns
fit1 <- mppm(X~1, H)
quadrat.test(fit1, nx=2)
# uniform Poisson with different intensity for each pattern
fit2 <- mppm(X ~ id, H)
quadrat.test(fit2, nx=2)


###################################################
### code chunk number 209: PromptOff.Rnw:1-2
###################################################
options(prompt="  ")


###################################################
### code chunk number 210: 16replicated.Rnw:2979-2980 (eval = FALSE)
###################################################
## cdf.test(model, covariate, test)


###################################################
### code chunk number 211: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")