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


R code for chapter 15

Below is the R code used to generate results and figures in chapter 15. 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 '15multivariate.Rnw'
## Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner

###################################################
### code chunk number 1: 15multivariate.Rnw:9-10
###################################################
source("R/startup.R")


###################################################
### code chunk number 2: SprucesL.Rnw:3-5
###################################################
newplot(7.1, 0.9)
setmargins(0,0.5,0,0)


###################################################
### code chunk number 3: 15multivariate.Rnw:133-134
###################################################
plot(spruces, main="", markscale=4)


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


###################################################
### code chunk number 5: 15multivariate.Rnw:171-172
###################################################
setmargins(0,0.7,0,0)


###################################################
### code chunk number 6: 15multivariate.Rnw:176-177
###################################################
plot(finpines, main="", markscale=0.1, mar.panel=0, hsep=1.2)


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


###################################################
### code chunk number 8: 15multivariate.Rnw:216-220
###################################################
par(mfrow=c(1,2))
hist(marks(spruces), main="")
hist(marks(longleaf), main="")
par(mfrow=c(1,1))


###################################################
### code chunk number 9: 15multivariate.Rnw:234-236 (eval = FALSE)
###################################################
## Hp <- anylapply(marks(finpines), hist, plot=FALSE) 
## plot(Hp)


###################################################
### code chunk number 10: 15multivariate.Rnw:240-241 (eval = FALSE)
###################################################
## pairs(marks(finpines))


###################################################
### code chunk number 11: 15multivariate.Rnw:244-246 (eval = FALSE)
###################################################
## pairs(marks(finpines), diag.panel=panel.histogram)
## with(marks(finpines), plot(diameter ~ height))


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


###################################################
### code chunk number 13: 15multivariate.Rnw:257-258
###################################################
with(marks(finpines), plot(diameter ~ height))


###################################################
### code chunk number 14: 15multivariate.Rnw:268-270
###################################################
vols <- with(marks(finpines), (pi/12) * height * diameter^2)
X <- finpines %mark% vols


###################################################
### code chunk number 15: 15multivariate.Rnw:289-290 (eval = FALSE)
###################################################
## pairs(as.data.frame(longleaf))


###################################################
### code chunk number 16: 15multivariate.Rnw:311-312
###################################################
Y <- cut(longleaf, breaks=c(0, 30, Inf), labels=c("J","A"))


###################################################
### code chunk number 17: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)


###################################################
### code chunk number 18: 15multivariate.Rnw:325-327
###################################################
X <- cut(longleaf,breaks=c(0,30,Inf),right=FALSE,labels=c("J","A"))
plot(X,main="",chars=c(20,1))


###################################################
### code chunk number 19: 15multivariate.Rnw:351-352 (eval = FALSE)
###################################################
## NB <- cut(nbfires, "fnl.size", breaks=4)


###################################################
### code chunk number 20: 15multivariate.Rnw:360-361
###################################################
big.old <- subset(nbfires, fnl.size > 10 & dis.date < "2001-01-01")


###################################################
### code chunk number 21: 15multivariate.Rnw:368-370
###################################################
longsmooth <- Smooth(longleaf, bw.smoothppp)
longnear   <- nnmark(longleaf)


###################################################
### code chunk number 22: 15multivariate.Rnw:396-397 (eval = FALSE)
###################################################
## longsmooth <- Smooth(longleaf, bw.smoothppp)


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


###################################################
### code chunk number 24: 15multivariate.Rnw:408-411
###################################################
plot(solist(longsmooth, longnear),
     main="", main.panel="", 
     equal.ribbon=TRUE, mar.panel=0, hsep=1)


###################################################
### code chunk number 25: 15multivariate.Rnw:496-497
###################################################
msd <- sqrt(markvar(longleaf, bw.smoothppp))


###################################################
### code chunk number 26: 15multivariate.Rnw:503-511
###################################################
mfit <- Smooth(longleaf, bw.smoothppp, at="points")
res <- marks(longleaf) - mfit
stuff <- solist(layered(msd, 
                        plotargs=list(ribside="bottom")), 
                layered(setmarks(longleaf, res), 
                        plotargs=list(leg.side="bottom")),
                layered(idw(longleaf),
                        plotargs=list(ribside="bottom")))


###################################################
### code chunk number 27: 15multivariate.Rnw:516-517
###################################################
setmargins(0)


###################################################
### code chunk number 28: 15multivariate.Rnw:522-525
###################################################
plot(stuff,
     main="", main.panel="", equal.scales=TRUE, valign=TRUE,
     mar.panel=c(1,0,0,0), hsep=0.5)


###################################################
### code chunk number 29: 15multivariate.Rnw:541-541
###################################################



###################################################
### code chunk number 30: 15multivariate.Rnw:553-555 (eval = FALSE)
###################################################
## mfit <- Smooth(longleaf, bw.smoothppp, at="points")
## res <- marks(longleaf) - mfit


###################################################
### code chunk number 31: 15multivariate.Rnw:703-705 (eval = FALSE)
###################################################
## envelope(spruces, markcorr, nsim=39, 
##          simulate=expression(rlabel(spruces)))


###################################################
### code chunk number 32: 15multivariate.Rnw:713-716
###################################################
set.seed(27291)
emcs <- envelope(spruces, markcorr, nsim=39, 
                 simulate=expression(rlabel(spruces)))


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


###################################################
### code chunk number 34: 15multivariate.Rnw:722-724
###################################################
plot(emcs, main="", ylim.covers=c(0.5, 1.1), 
     legendpos="top", legendargs=list(cex=0.85))


###################################################
### code chunk number 35: 15multivariate.Rnw:787-792
###################################################
set.seed(6876109)
smki <- Kmark(spruces, returnL=TRUE)
smke <- envelope(spruces, Kmark, returnL=TRUE,
                 simulate=expression(rlabel(spruces)),
                 nsim=1999, nrank=50, savefuns=TRUE)


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


###################################################
### code chunk number 37: 15multivariate.Rnw:798-799 (eval = FALSE)
###################################################
## Kmark(X, f, ...)


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


###################################################
### code chunk number 39: 15multivariate.Rnw:808-812 (eval = FALSE)
###################################################
## smki <- Kmark(spruces, returnL=TRUE)
## smke <- envelope(spruces, Kmark, returnL=TRUE,
##                  simulate=expression(rlabel(spruces)),
##                  nsim=1999, nrank=50, savefuns=TRUE)


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


###################################################
### code chunk number 41: 15multivariate.Rnw:824-827
###################################################
plot(anylist(smki, smke), main.panel="", main="",
     panel.args=function(i) if(i == 1) list() else list(. - mmean ~ r),
     legend=FALSE)


###################################################
### code chunk number 42: 15multivariate.Rnw:847-849
###################################################
Lspruce <- Lest(spruces)
diffL <- eval.fv(smki - Lspruce)


###################################################
### code chunk number 43: 15multivariate.Rnw:902-903
###################################################
mv <- markvario(spruces)


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


###################################################
### code chunk number 45: 15multivariate.Rnw:908-909
###################################################
plot(mv, main="", legendpos="top")


###################################################
### code chunk number 46: 15multivariate.Rnw:951-953
###################################################
ES <- Emark(spruces)
VS <- Vmark(spruces)


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


###################################################
### code chunk number 48: 15multivariate.Rnw:958-960
###################################################
plot(anylist(ES, VS), main="", main.panel="",
     panel.args=function(i) list(ylim.covers=if(i==1) c(0,0.3) else 0))


###################################################
### code chunk number 49: 15multivariate.Rnw:984-985
###################################################
nnmean(spruces)


###################################################
### code chunk number 50: 15multivariate.Rnw:997-998
###################################################
nnvario(spruces)


###################################################
### code chunk number 51: 15multivariate.Rnw:1029-1031
###################################################
nncorr(spruces)
nncorr(finpines)


###################################################
### code chunk number 52: 15multivariate.Rnw:1043-1045
###################################################
rr <- nncorr(spruces)[["correlation"]]
sqrt(2 * (1-rr^2)/(npoints(spruces)-4))


###################################################
### code chunk number 53: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)


###################################################
### code chunk number 54: 15multivariate.Rnw:1078-1080
###################################################
X <- osteo$pts[[36]]
plot(X, main="", pch=21, bg='white', cex=1.2)


###################################################
### code chunk number 55: 15multivariate.Rnw:1136-1137
###################################################
X <- osteo$pts[[36]]


###################################################
### code chunk number 56: 15multivariate.Rnw:1145-1149 (eval = FALSE)
###################################################
## copyExampleFiles("osteo")
## xyz <- read.table("osteo36.txt", header=TRUE)
## b <- box3(c(0,81),c(0,100),c(-100,0),unitname=c("micron", "microns"))
## X <- with(xyz, pp3(x, y, z, b))


###################################################
### code chunk number 57: 15multivariate.Rnw:1151-1153
###################################################
X
summary(X)


###################################################
### code chunk number 58: 15multivariate.Rnw:1199-1200 (eval = FALSE)
###################################################
## plot(X, main="", pch=21, bg='white', cex=1.2)


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


###################################################
### code chunk number 60: 15multivariate.Rnw:1270-1271 (eval = FALSE)
###################################################
## K3est(X, ..., rmax, nrval, correction)


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


###################################################
### code chunk number 62: 15multivariate.Rnw:1284-1286
###################################################
KX <- K3est(X)
EK <- envelope(X, K3est, nsim=1999, nrank=50, nrval=512)


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


###################################################
### code chunk number 64: 15multivariate.Rnw:1292-1297
###################################################
pan <- function(i) {
  if(i == 1) list() else list(fmla= sqrt(.) ~ r, legendargs=list(cex=0.85))
}
plot(anylist(KX, EK), main="", main.panel="", 
     mar.panel=c(3,3,0,0), hsep=2, panel.args=pan)


###################################################
### code chunk number 65: 15multivariate.Rnw:1340-1342 (eval = FALSE)
###################################################
## E <- envelope(X, K3est, nsim=1999, nrank=50, nrval=512)
## plot(E, sqrt(.) ~ r)


###################################################
### code chunk number 66: 15multivariate.Rnw:1380-1381
###################################################
p3 <- pcf3est(X)


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


###################################################
### code chunk number 68: 15multivariate.Rnw:1386-1387
###################################################
plot(p3, main="")


###################################################
### code chunk number 69: 15multivariate.Rnw:1425-1426
###################################################
G3 <- G3est(X)


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


###################################################
### code chunk number 71: 15multivariate.Rnw:1431-1432
###################################################
plot(G3, main="", legend=FALSE)


###################################################
### code chunk number 72: 15multivariate.Rnw:1449-1452
###################################################
v <- eroded.volumes(domain(X), c(0, 10, 20, 30))
f <- v/v[1]
round(f, 2)


###################################################
### code chunk number 73: 15multivariate.Rnw:1530-1531
###################################################
F3 <- F3est(X, sphere="digital")


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


###################################################
### code chunk number 75: 15multivariate.Rnw:1537-1538
###################################################
plot(F3, main="")


###################################################
### code chunk number 76: 15multivariate.Rnw:1615-1622
###################################################
df <- data.frame(x=runif(100, max=3), 
                 y=runif(100, max=3), 
                 z=runif(100, max=2),
                 t=runif(100))
bb <- boxx(c(0,3), c(0,3), c(0,2), c(0,1))
X <- ppx(data=df, domain=bb, coord.type=c("s","s", "s", "t"))
X


###################################################
### code chunk number 77: 15multivariate.Rnw:1631-1633
###################################################
marks(X) <- with(as.hyperframe(df), disc(centre=c(x,y)))
X