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


R code for chapter 14

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

###################################################
### code chunk number 1: 14multitype.Rnw:9-13
###################################################
source("R/startup.R")
requireversion(spatstat, "1.41-1.028")
require(RandomFields)
require(dixon)


###################################################
### code chunk number 2: 14multitype.Rnw:16-21
###################################################
## shorter output
Digits <- 4
dopt <- options(digits=Digits)
Terse <- 3
spatstat.options(terse=Terse)


###################################################
### code chunk number 3: MucosaL.Rnw:3-5
###################################################
newplot(5, 0.5)
setmargins(0)


###################################################
### code chunk number 4: 14multitype.Rnw:52-53
###################################################
plot(mucosa, main="", chars=c(16, 3), cex=c(1, 0.6))


###################################################
### code chunk number 5: Mucosa2.Rnw:3-5
###################################################
newplot(10, 1)
setmargins(0)


###################################################
### code chunk number 6: 14multitype.Rnw:383-386
###################################################
plot(split(mucosa), main="", main.panel="",
     mar.panel=0, hsep=1,
     panel.args=function(i) if(i == 1) list(pch=16) else list(pch=3, cex=0.7))


###################################################
### code chunk number 7: 14multitype.Rnw:679-682
###################################################
mydata <- data.frame(x=runif(10), y=runif(10), 
           m=factor(sample(letters[1:3], 10, replace=TRUE)))
X <- as.ppp(mydata, square(1))


###################################################
### code chunk number 8: 14multitype.Rnw:693-696 (eval = FALSE)
###################################################
## copyExampleFiles("amacrine")
## W <- owin(c(0,1060/662), c(0,1))
## X <- scanpp("amacrine.txt", factor.marks=TRUE, window=W)


###################################################
### code chunk number 9: 14multitype.Rnw:713-714 (eval = FALSE)
###################################################
## marks(X) <- sample(marks(X), npoints(X))


###################################################
### code chunk number 10: 14multitype.Rnw:768-770
###################################################
Y <- split(mucosa)
Y


###################################################
### code chunk number 11: 14multitype.Rnw:784-785
###################################################
superimpose(A=runifpoint(5), B=runifpoint(7))


###################################################
### code chunk number 12: 14multitype.Rnw:789-790
###################################################
superimpose(runifpoint(5), runifpoint(7))


###################################################
### code chunk number 13: 14multitype.Rnw:825-827
###################################################
X <- betacells
marks(X) <- marks(X)[,1]


###################################################
### code chunk number 14: 14multitype.Rnw:830-833
###################################################
XX <- do.call(superimpose, split(X))
identical(X, XX)
max(nncross(X, XX)$dist)


###################################################
### code chunk number 15: 14multitype.Rnw:841-845
###################################################
X <- amacrine
Y <- split(amacrine)
Y$on <- rjitter(Y$on, 0.1)
split(X) <- Y


###################################################
### code chunk number 16: UnitL.Rnw:3-5
###################################################
newplot(9, 0.7)
setmargins(0.1+ c(0,2,0,0))


###################################################
### code chunk number 17: 14multitype.Rnw:905-906
###################################################
setmargins(0.05+c(0,1.5, 0, 0))


###################################################
### code chunk number 18: 14multitype.Rnw:910-911
###################################################
plot(lansing, main="", cex=0.6, leg.args=list(cex=1))


###################################################
### code chunk number 19: Unit2x3title.Rnw:3-5
###################################################
newplot(6.5, 0.8)
setmargins(0)


###################################################
### code chunk number 20: 14multitype.Rnw:930-931
###################################################
plot(split(lansing), main="", pch="+", mar.panel=c(0,0,1,0), hsep=1, vsep=1)


###################################################
### code chunk number 21: 14multitype.Rnw:943-944 (eval = FALSE)
###################################################
## plot(nbfires, which.marks="fire.type")


###################################################
### code chunk number 22: 14multitype.Rnw:949-950 (eval = FALSE)
###################################################
## plot(subset(nbfires, year == "1996"), which.marks="fire.type")


###################################################
### code chunk number 23: 14multitype.Rnw:953-954 (eval = FALSE)
###################################################
## plot(subset(nbfires, year == "1996", select=fire.type))


###################################################
### code chunk number 24: 14multitype.Rnw:958-959
###################################################
nbf <- subset(nbfires, year == "1996")


###################################################
### code chunk number 25: nbfiresL.Rnw:3-5
###################################################
newplot(5.76, 0.6)
setmargins(0)


###################################################
### code chunk number 26: 14multitype.Rnw:964-965
###################################################
plot(nbf, which.marks="fire.type", main="")


###################################################
### code chunk number 27: 14multitype.Rnw:976-978
###################################################
## print warnings in summary.ppp
spatstat.options(terse=2)


###################################################
### code chunk number 28: 14multitype.Rnw:983-985
###################################################
lansing
summary(lansing)


###################################################
### code chunk number 29: 14multitype.Rnw:992-994
###################################################
## back to terse printout
spatstat.options(terse=Terse)


###################################################
### code chunk number 30: 14multitype.Rnw:1000-1001
###################################################
intensity(rescale(lansing))


###################################################
### code chunk number 31: 14multitype.Rnw:1007-1008
###################################################
intensity(rescale(lansing, 5280/924))


###################################################
### code chunk number 32: 14multitype.Rnw:1014-1015
###################################################
head(as.data.frame(amacrine), 3)


###################################################
### code chunk number 33: 14multitype.Rnw:1108-1110
###################################################
sum(intensity(lansing))
intensity(unmark(lansing))


###################################################
### code chunk number 34: 14multitype.Rnw:1117-1118
###################################################
sqrt(intensity(lansing)/area(Window(lansing)))


###################################################
### code chunk number 35: 14multitype.Rnw:1124-1129
###################################################
Nx <- Ny <- 4
QC <- lapply(split(lansing), quadratcount, nx=Nx, ny=Ny)
QC <- lapply(QC, as.vector)
v <- sapply(QC, var)
(se <- sqrt(Nx * Ny * v)/area(Window(lansing)))


###################################################
### code chunk number 36: 14multitype.Rnw:1174-1175
###################################################
denlan <- density(split(lansing))


###################################################
### code chunk number 37: listof2x3.Rnw:3-5
###################################################
newplot(16.5, 1)
setmargins(0)


###################################################
### code chunk number 38: 14multitype.Rnw:1180-1181
###################################################
plot(denlan, main="", mar.panel=c(0,0,0,1), hsep=1)


###################################################
### code chunk number 39: 14multitype.Rnw:1198-1199 (eval = FALSE)
###################################################
## plot(density(split(lansing)))


###################################################
### code chunk number 40: 14multitype.Rnw:1218-1219 (eval = FALSE)
###################################################
## plot(density(split(lansing), se=TRUE)$SE)


###################################################
### code chunk number 41: 14multitype.Rnw:1252-1255
###################################################
lambda <- intensity(lansing)
probs <- lambda/sum(lambda)
probs


###################################################
### code chunk number 42: pLansing
###################################################
ProbU <- relrisk(lansing)
ProbD  <- relrisk(lansing, diggle=TRUE)


###################################################
### code chunk number 43: 14multitype.Rnw:1339-1342
###################################################
coco <- if(monochrome) grey(seq(1,0,length=128)) else NULL
plot(ProbD, main="", mar.panel=0.1+c(0,0,1,2), 
     equal.ribbon=TRUE, col=coco)


###################################################
### code chunk number 44: 14multitype.Rnw:1355-1358 (eval = FALSE)
###################################################
## lami <- density(split(lansing))
## lamdot <- Reduce("+", lami)
## probi <- solapply(lami, "/", lamdot)


###################################################
### code chunk number 45: 14multitype.Rnw:1373-1374
###################################################
ProbX <- relrisk(lansing, at="points")


###################################################
### code chunk number 46: 14multitype.Rnw:1376-1377
###################################################
ProbX[1:5]


###################################################
### code chunk number 47: 14multitype.Rnw:1435-1436 (eval = FALSE)
###################################################
## b <- bw.relrisk(lansing)


###################################################
### code chunk number 48: 14multitype.Rnw:1438-1440
###################################################
b <- bw.relrisk(lansing)
b2 <- bw.relrisk(lansing, hmin=0.025, hmax=0.1, nh=32) 


###################################################
### code chunk number 49: 14multitype.Rnw:1452-1453 (eval = FALSE)
###################################################
## Prob <- relrisk(lansing, sigma=b)


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


###################################################
### code chunk number 51: 14multitype.Rnw:1462-1467
###################################################
par(mfrow=c(1,2))
oa <- list(col=if(monochrome) 1 else 4, lty=2)
plot(b, main="", optargs=oa)
plot(b2, main="", optargs=oa)
par(mfrow=c(1,1))


###################################################
### code chunk number 52: 14multitype.Rnw:1490-1493
###################################################
dominant <- im.apply(ProbD, which.max)
species <- levels(marks(lansing))
dominant <- eval.im(factor(dominant, levels=1:6, labels=species))


###################################################
### code chunk number 53: 14multitype.Rnw:1502-1503
###################################################
setmargins(0.2 + c(0,0,0,3))


###################################################
### code chunk number 54: 14multitype.Rnw:1508-1509
###################################################
textureplot(dominant, main="")


###################################################
### code chunk number 55: LansingSegTest
###################################################
seglan <- segregation.test(lansing)


###################################################
### code chunk number 56: 14multitype.Rnw:1553-1554 (eval = FALSE)
###################################################
## segregation.test(lansing, nsim=19)


###################################################
### code chunk number 57: 14multitype.Rnw:1556-1557
###################################################
seglan


###################################################
### code chunk number 58: 14multitype.Rnw:1569-1570
###################################################
mucoRRnonparam <- relrisk(mucosa, casecontrol=FALSE)[[1]]


###################################################
### code chunk number 59: MucosaR.Rnw:3-5
###################################################
newplot(5, 0.5)
setmargins(0,0,0,1.2)


###################################################
### code chunk number 60: 14multitype.Rnw:1581-1583
###################################################
plot(mucoRRnonparam, main="")
contour(mucoRRnonparam, add=TRUE, col="white", drawlabels=FALSE)


###################################################
### code chunk number 61: 14multitype.Rnw:1696-1697 (eval = FALSE)
###################################################
## plot(relrisk(chorley, relative=TRUE, control="lung"))


###################################################
### code chunk number 62: 14multitype.Rnw:1704-1705 (eval = FALSE)
###################################################
## plot(relrisk(chorley, relative=TRUE, control="lung", se=TRUE)$SE)


###################################################
### code chunk number 63: 14multitype.Rnw:1709-1712
###################################################
rrcho <- relrisk(chorley,relative=TRUE,control="lung", 
                 hmin=10, hmax=20, nh=64,
                 se=TRUE)


###################################################
### code chunk number 64: 14multitype.Rnw:1714-1715
###################################################
spatstat.options(image.colfun=lighttodark)


###################################################
### code chunk number 65: Chorley2LR.Rnw:3-5
###################################################
newplot(18, 1)
setmargins(0, 2, 0, 2)


###################################################
### code chunk number 66: 14multitype.Rnw:1719-1720
###################################################
setmargins(0,0,0,2)


###################################################
### code chunk number 67: 14multitype.Rnw:1724-1731
###################################################
zz <- range(0.061, range(rrcho[[1]]))
pan <- function(i) { if(i == 1) list(zlim=zz) else list() }
pen <- function(i, y, ...) contour(y, ..., drawlabels=FALSE, col="white")
plot(as.solist(rrcho), 
     main="", main.panel="", mar.panel=0, hsep=2.8, 
     equal.scales=TRUE, panel.end=pen, panel.args=pan,
     ribscale=100)


###################################################
### code chunk number 68: 14multitype.Rnw:1824-1832
###################################################
## generate homogeneous and inhomogeneous multitype Poisson patterns.
set.seed(66888)
X1   <- rmpoispp(100, types=c("A","B"))
X2   <- rmpoispp(c(100,20), types=c("A","B"))
X3   <- rmpoispp(function(x,y,m){300*exp(-3*x)},types=c("A","B"))
lamb <- function(x,y,m) {ifelse(m=="A",
            300*exp(-4*x),300*exp(-4*(1-x))) }
X4   <- rmpoispp(lamb, types=c("A","B"))


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


###################################################
### code chunk number 70: 14multitype.Rnw:1838-1839
###################################################
plot(solist(X1, X2),main="",main.panel="",chars=c(16,3))


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


###################################################
### code chunk number 72: 14multitype.Rnw:1915-1916
###################################################
plot(solist(X3,X4),main="",main.panel="", chars=c(16,3))


###################################################
### code chunk number 73: 14multitype.Rnw:1994-1999 (eval = FALSE)
###################################################
## X1 <- rmpoispp(100, types=c("A","B"))
## X2 <- rmpoispp(c(100,20), types=c("A","B"))
## X3 <- rmpoispp(function(x,y,m){300*exp(-3*x)}, types=c("A","B"))
## ll <- function(x,y,m){ 300*exp(-4*ifelse(m=="A", x, (1-x))) }
## X4 <- rmpoispp(ll, types=c("A","B"))


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


###################################################
### code chunk number 75: 14multitype.Rnw:2067-2068 (eval = FALSE)
###################################################
## ppm(X ~ marks)


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


###################################################
### code chunk number 77: 14multitype.Rnw:2076-2077
###################################################
ppm(lansing ~ marks)


###################################################
### code chunk number 78: 14multitype.Rnw:2094-2095
###################################################
ppm(lansing ~ 1)


###################################################
### code chunk number 79: 14multitype.Rnw:2112-2113
###################################################
ppm(unmark(lansing) ~ 1)


###################################################
### code chunk number 80: 14multitype.Rnw:2132-2133 (eval = FALSE)
###################################################
## ppm(lansing ~ marks + x)


###################################################
### code chunk number 81: 14multitype.Rnw:2135-2138
###################################################
fit <- ppm(lansing ~ marks + x)
co <- coef(fit)
co <- signif(co, Digits)


###################################################
### code chunk number 82: 14multitype.Rnw:2140-2141
###################################################
fit


###################################################
### code chunk number 83: 14multitype.Rnw:2181-2182
###################################################
coX <- signif(coef(ppm(lansing ~ marks * x)), Digits)


###################################################
### code chunk number 84: 14multitype.Rnw:2186-2187
###################################################
ppm(lansing ~ marks * x)


###################################################
### code chunk number 85: 14multitype.Rnw:2224-2227
###################################################
mod3  <- ppm(lansing ~ polynom(x,y,3) * marks)
pred3 <- predict(mod3)
names(pred3) <- levels(marks(lansing))


###################################################
### code chunk number 86: Unit2x3tR.Rnw:4-6
###################################################
newplot(15, 1)
setmargins(0)


###################################################
### code chunk number 87: 14multitype.Rnw:2232-2238
###################################################
pen <- function(i, y, ...) contour(y, ..., drawlabels=FALSE, col="white")
plot(pred3, main="",mar.panel=c(0,0,0,1), hsep=1, panel.end=pen)
if(FALSE){
plot(pred3, equal.ribbon=TRUE,main="",mar.panel=rep(0.1,4), 
     ribmar=c(1,0,0,2))
}


###################################################
### code chunk number 88: 14multitype.Rnw:2260-2263 (eval = FALSE)
###################################################
## f0 <- ppm(lansing ~ polynom(x,y,3) + marks)
## f1 <- ppm(lansing ~ polynom(x,y,3) * marks)
## anova(f0, f1, test="LR")


###################################################
### code chunk number 89: 14multitype.Rnw:2278-2281
###################################################
Ants <- ants
levels(marks(Ants)) <- c("C", "M")
ppm(Ants ~ marks * (side - 1), data=ants.extra)


###################################################
### code chunk number 90: 14multitype.Rnw:2349-2350
###################################################
fit <- ppm(mucosa ~ marks * polynom(x,y,3))


###################################################
### code chunk number 91: cache
###################################################
probs <- relrisk(fit, casecontrol=FALSE)


###################################################
### code chunk number 92: 14multitype.Rnw:2361-2362 (eval = FALSE)
###################################################
## rusk <- relrisk(fit, relative=TRUE, control="other")


###################################################
### code chunk number 93: 14multitype.Rnw:2367-2369
###################################################
mucofit <- ppm(mucosa ~ marks * polynom(x,y,3))
mucoRRparam <- relrisk(mucofit, casecontrol=FALSE)[[1]]


###################################################
### code chunk number 94: MucosaR.Rnw:3-5
###################################################
newplot(5, 0.5)
setmargins(0,0,0,1.2)


###################################################
### code chunk number 95: 14multitype.Rnw:2375-2377
###################################################
plot(mucoRRparam, main="")
contour(mucoRRparam, add=TRUE, col="white", drawlabels=FALSE)


###################################################
### code chunk number 96: 14multitype.Rnw:2434-2436
###################################################
d <- nndist(amacrine, by=marks(amacrine))
head(d, 3)


###################################################
### code chunk number 97: 14multitype.Rnw:2447-2449
###################################################
d <- nndist(amacrine, by=marks(amacrine))
a <- aggregate(d, by=list(from=marks(amacrine)), min)


###################################################
### code chunk number 98: 14multitype.Rnw:2457-2461
###################################################
m <- marks(amacrine)
n <- nnwhich(amacrine)
m1 <- m[n]
table(m, m1)


###################################################
### code chunk number 99: 14multitype.Rnw:2514-2515
###################################################
nncorr(amacrine)


###################################################
### code chunk number 100: 14multitype.Rnw:2552-2553
###################################################
marktable(amacrine, N=1, collapse=TRUE)


###################################################
### code chunk number 101: 14multitype.Rnw:2560-2563
###################################################
dL <- dixon(as.data.frame(lansing))
dA <- dixon(as.data.frame(amacrine))
dM <- dixon(as.data.frame(mucosa))


###################################################
### code chunk number 102: 14multitype.Rnw:2565-2567 (eval = FALSE)
###################################################
## library(dixon)
## dixon(as.data.frame(lansing))$tablaC


###################################################
### code chunk number 103: 14multitype.Rnw:2569-2570
###################################################
dL$tablaC


###################################################
### code chunk number 104: 14multitype.Rnw:2574-2575 (eval = FALSE)
###################################################
## dixon(as.data.frame(amacrine))$tablaC


###################################################
### code chunk number 105: 14multitype.Rnw:2577-2578
###################################################
dA$tablaC


###################################################
### code chunk number 106: 14multitype.Rnw:2642-2643 (eval = FALSE)
###################################################
## K <- Kcross(amacrine, "on", "off")


###################################################
### code chunk number 107: 14multitype.Rnw:2657-2658
###################################################
Kall <- alltypes(rescale(amacrine), Kcross)


###################################################
### code chunk number 108: fasp2x2.Rnw:3-5
###################################################
newplot(12, 0.75)
setmargins(4, 4, 0, 0.5)


###################################################
### code chunk number 109: 14multitype.Rnw:2680-2681
###################################################
plot(Kall, lwd=2, banner=FALSE, xlim=c(0, 100))


###################################################
### code chunk number 110: 14multitype.Rnw:2828-2829 (eval = FALSE)
###################################################
## ma <- markconnect(amacrine, "on", "off")


###################################################
### code chunk number 111: 14multitype.Rnw:2846-2847
###################################################
AM <- alltypes(rescale(amacrine), markconnect)


###################################################
### code chunk number 112: fasp2x2.Rnw:3-5
###################################################
newplot(12, 0.75)
setmargins(4, 4, 0, 0.5)


###################################################
### code chunk number 113: 14multitype.Rnw:2852-2853
###################################################
plot(AM,banner=FALSE)


###################################################
### code chunk number 114: 14multitype.Rnw:2943-2945 (eval = FALSE)
###################################################
## AG <- Gcross(amacrine, "on", "off")
## AJ <- Jcross(amacrine, "on", "off")


###################################################
### code chunk number 115: 14multitype.Rnw:2957-2958
###################################################
AAG <- alltypes(rescale(amacrine), Gcross)


###################################################
### code chunk number 116: fasp2x2.Rnw:3-5
###################################################
newplot(12, 0.75)
setmargins(4, 4, 0, 0.5)


###################################################
### code chunk number 117: 14multitype.Rnw:2963-2964
###################################################
plot(AAG, lwd=2,banner=FALSE)


###################################################
### code chunk number 118: 14multitype.Rnw:2993-2995 (eval = FALSE)
###################################################
## plot(Gcross(amacrine, "on", "off"), km ~ r)
## plot(Fest(split(amacrine)$off), km ~ r, add=TRUE, lty=2)


###################################################
### code chunk number 119: 14multitype.Rnw:3103-3113
###################################################
X  <- rescale(amacrine)
AG <- Gdot(X, "on")
Jdif <- function(X, ..., i) { 
  Jidot <- Jdot(X, ..., i=i)
  Jii <- Jcross(X,i=i,j=i)
  dif <- eval.fv(Jidot-Jii)
  return(dif)
}
AJ <- Jdif(X,i="on")
PA <- function(i){switch(EXPR=i,list(ylim=c(0,1)),list(ylim=c(-0.3,0.3)))}


###################################################
### code chunk number 120: fv2Rolf.Rnw:3-5
###################################################
newplot(12, 1)
zeromargins()


###################################################
### code chunk number 121: 14multitype.Rnw:3118-3120
###################################################
plot(anylist(AG,AJ),main="",main.panel="",
     legend=FALSE,panel.args=PA,xlim=c(0,45))


###################################################
### code chunk number 122: 14multitype.Rnw:3147-3149
###################################################
aG <- alltypes(rescale(amacrine), "G")
Phi <- function(x) asin(sqrt(x))


###################################################
### code chunk number 123: 14multitype.Rnw:3151-3152 (eval = FALSE)
###################################################
## plot(aG, Phi(.) ~ Phi(theo))


###################################################
### code chunk number 124: fasp2x2.Rnw:3-5
###################################################
newplot(12, 0.75)
setmargins(4, 4, 0, 0.5)


###################################################
### code chunk number 125: 14multitype.Rnw:3159-3160
###################################################
plot(aG, Phi(.) ~ Phi(theo), banner=FALSE)


###################################################
### code chunk number 126: 14multitype.Rnw:3172-3177
###################################################
lansL <- alltypes(lansing, Lcross)
subL <- lansL[1:3, 1:3]
hickL <- lansL[2, ]
oaks <- c("redoak", "blackoak", "whiteoak")
oakL <- lansL[oaks, oaks]


###################################################
### code chunk number 127: 14multitype.Rnw:3181-3182
###################################################
Lmh <- lansL[3,2,drop=TRUE]


###################################################
### code chunk number 128: 14multitype.Rnw:3188-3189 (eval = FALSE)
###################################################
## aGfish <- eval.fasp(Phi(aG))


###################################################
### code chunk number 129: 14multitype.Rnw:3280-3283
###################################################
RA <- rescale(amacrine)
MCA <- markcorr(RA)
IA <- Iest(RA, r=seq(0, 100, length=512))


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


###################################################
### code chunk number 131: 14multitype.Rnw:3289-3290
###################################################
plot(anylist(MCA, IA), main="", main.panel="", legend=FALSE)


###################################################
### code chunk number 132: 14multitype.Rnw:3383-3389
###################################################
bD <- bw.diggle(mucosa)
bR <- bw.relrisk(mucosa)
LMD <- Lcross.inhom(mucosa, "ECL", "other", sigma=bD)
LMR <- Lcross.inhom(mucosa, "ECL", "other", sigma=bR)
bD3 <- signif(bD, 3)
bR3 <- signif(bR, 3)


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


###################################################
### code chunk number 134: 14multitype.Rnw:3395-3396
###################################################
plot(anylist(LMD, LMR), main="", main.panel="", legend=FALSE)


###################################################
### code chunk number 135: 14multitype.Rnw:3428-3430 (eval = FALSE)
###################################################
## LMD <- Lcross.inhom(mucosa, "ECL", "other", sigma=bw.diggle(mucosa))
## LMR <- Lcross.inhom(mucosa, "ECL", "other", sigma=bw.relrisk(mucosa))


###################################################
### code chunk number 136: 14multitype.Rnw:3594-3597 (eval = FALSE)
###################################################
## Amac <- rescale(amacrine)
## E <- envelope(Amac, Lcross, nsim=39, i="on", j="off",
##               simulate=expression(rshift(Amac, radius=150)))


###################################################
### code chunk number 137: 14multitype.Rnw:3614-3618
###################################################
set.seed(40063)
Amac <- rescale(amacrine)
E <- envelope(Amac, Lcross, nsim=39, i="on", j="off",
              simulate=expression(rshift(Amac, radius=150)))


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


###################################################
### code chunk number 139: 14multitype.Rnw:3625-3626
###################################################
plot(E, . - r ~ r, main="",legend=FALSE)


###################################################
### code chunk number 140: 14multitype.Rnw:3690-3696
###################################################
Jdif <- function(X, ..., i) {
  Jidot <- Jdot(X, ..., i=i)
  J <- Jest(X, ...)
  dif <- eval.fv(Jidot-J)
  return(dif)
}


###################################################
### code chunk number 141: 14multitype.Rnw:3709-3710
###################################################
set.seed(58742)


###################################################
### code chunk number 142: 14multitype.Rnw:3712-3714
###################################################
EPJ <- envelope(paracou, Jdif, nsim=39, i="juvenile",
                simulate=expression(rlabel(paracou)))


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


###################################################
### code chunk number 144: 14multitype.Rnw:3740-3741
###################################################
plot(EPJ, main="",legend=FALSE)


###################################################
### code chunk number 145: 14multitype.Rnw:3773-3775
###################################################
X <- rescale(amacrine)
aE <- alltypes(X, Lcross, nsim=39, envelope=TRUE, global=TRUE)


###################################################
### code chunk number 146: fasp2x2.Rnw:3-5
###################################################
newplot(12, 0.75)
setmargins(4, 4, 0, 0.5)


###################################################
### code chunk number 147: 14multitype.Rnw:3786-3787
###################################################
plot(aE, . - r  ~ r, banner=FALSE, samex=TRUE, samey=TRUE)


###################################################
### code chunk number 148: 14multitype.Rnw:3859-3860 (eval = FALSE)
###################################################
## ppm(amacrine ~ marks, MultiHard())


###################################################
### code chunk number 149: 14multitype.Rnw:4200-4201 (eval = FALSE)
###################################################
## ppm(amacrine ~ marks, interaction=MultiHard())


###################################################
### code chunk number 150: 14multitype.Rnw:4229-4230
###################################################
ppm(hyytiala ~ marks, Hardcore())


###################################################
### code chunk number 151: 14multitype.Rnw:4245-4246
###################################################
ppm(hamster ~ marks, Strauss(0.02))


###################################################
### code chunk number 152: 14multitype.Rnw:4321-4322
###################################################
ppm(amacrine ~ marks, MultiHard())


###################################################
### code chunk number 153: 14multitype.Rnw:4340-4342
###################################################
Beta <- betacells
marks(Beta) <- marks(Beta)[,"type"]


###################################################
### code chunk number 154: 14multitype.Rnw:4345-4346
###################################################
Beta <- subset(betacells, select=type)


###################################################
### code chunk number 155: 14multitype.Rnw:4354-4357
###################################################
R <- diag(c(60,80))
fit0 <- ppm(Beta ~ marks, MultiStrauss(R))
fit0


###################################################
### code chunk number 156: 14multitype.Rnw:4361-4364
###################################################
R2 <- matrix(c(60,60,60,80), nrow=2, ncol=2)
fit1 <- ppm(Beta ~ marks, MultiStrauss(R2))
fit1


###################################################
### code chunk number 157: 14multitype.Rnw:4372-4373
###################################################
anova(fit0, fit1, test="LR")


###################################################
### code chunk number 158: 14multitype.Rnw:4389-4391 (eval = FALSE)
###################################################
## fitP <- ppm(Beta ~ marks + polynom(x,y,3), MultiStrauss(R))
## fitI <- ppm(Beta ~ marks * polynom(x,y,3), MultiStrauss(R))


###################################################
### code chunk number 159: 14multitype.Rnw:4414-4415
###################################################
fitBeta <- ppm(Beta ~ marks, MultiStraussHard(iradii=R2))


###################################################
### code chunk number 160: 14multitype.Rnw:4417-4418 (eval = FALSE)
###################################################
## plot(fitin(fitBeta))


###################################################
### code chunk number 161: fasp2x2.Rnw:3-5
###################################################
newplot(12, 0.75)
setmargins(4, 4, 0, 0.5)


###################################################
### code chunk number 162: 14multitype.Rnw:4431-4432
###################################################
plot(fitin(fitBeta),banner=FALSE)


###################################################
### code chunk number 163: 14multitype.Rnw:4453-4455
###################################################
set.seed(24819)
Beta.sim <- rmh(fitBeta, verbose=FALSE)


###################################################
### code chunk number 164: 14multitype.Rnw:4457-4458 (eval = FALSE)
###################################################
## Beta.sim <- rmh(fitBeta)


###################################################
### code chunk number 165: Beta2L.Rnw:3-5
###################################################
newplot(12, 0.85)
setmargins(0)


###################################################
### code chunk number 166: 14multitype.Rnw:4466-4468
###################################################
plot(solist(Beta,Beta.sim),main="",main.panel="",
     mar.panel=0, hsep=1, equal.scales=TRUE)


###################################################
### code chunk number 167: 14multitype.Rnw:4483-4491 (eval = FALSE)
###################################################
## HC <- matrix(c(20,10,10,20),2,2)
## R  <- matrix(c(30,20,20,30),2,2)
## G  <- matrix(0.3,2,2)
## B  <- c(5,10)*1e-5
## M  <- rmhmodel(cif="straushm",
##                par=list(beta=B,gamma=G,iradii=R,hradii=HC),
##                types=c("A","B"), w=square(1000))
## X  <- rmh(M)


###################################################
### code chunk number 168: UnitL.Rnw:3-5
###################################################
newplot(9, 0.7)
setmargins(0.1+ c(0,2,0,0))


###################################################
### code chunk number 169: 14multitype.Rnw:4546-4547
###################################################
plot(bramblecanes, main="")


###################################################
### code chunk number 170: 14multitype.Rnw:4854-4858
###################################################
Ants <- ants 
levels(marks(Ants)) <- c("Cat", "Mess")
rmat <- matrix(c(65,45,45,30), nrow=2, ncol=2)
fit <- ppm(Ants ~ marks, HierStrauss(rmat, archy=c(2,1)))


###################################################
### code chunk number 171: 14multitype.Rnw:4860-4862
###################################################
fit
coef(summary(fit))


###################################################
### code chunk number 172: 14multitype.Rnw:4864-4866
###################################################
g <- summary(fitin(fit))$sensible$param$gammas
g <- signif(g, Digits)


###################################################
### code chunk number 173: 14multitype.Rnw:4870-4876 (eval = FALSE)
###################################################
## rseq   <- seq(20, 70, by=5)
## s   <- expand.grid(r11=rseq, r12=rseq, r22=rseq)
## pfn <- function(r11, r12, r22) {
##   HierStrauss(matrix(c(r11,r12,r12,r22),2,2), archy=2:1)
## }
## prf <- profilepl(s, pfn, ants ~ marks, correction="translate")


###################################################
### code chunk number 174: 14multitype.Rnw:4894-4896
###################################################
fit0 <- update(fit, HierStrauss(diag(diag(rmat))))
anova(fit0, fit, test="LR")


###################################################
### code chunk number 175: 14multitype.Rnw:4979-4980
###################################################
set.seed(873987)


###################################################
### code chunk number 176: 14multitype.Rnw:4982-4991
###################################################
P <- dirichlet(runifpoint(10))
logLambda <- rMosaicField(P, rnorm, 
                          dimyx=512, rgenargs=list(mean=4, sd=1))
Lambda <- exp(logLambda)
X <- rpoispp(Lambda)
Xlinked <- X %mark% factor(sample(c("a","b"), npoints(X), 
                                  replace=TRUE, prob=c(2,3)/5))
Y <- rpoispp(100)
Xbalanced <- Y %mark% factor(ifelse(Lambda[Y] > exp(4), "a", "b"))


###################################################
### code chunk number 177: 14multitype.Rnw:4994-4995
###################################################
set.seed(676532)


###################################################
### code chunk number 178: 14multitype.Rnw:4997-5003
###################################################
X1 <- rLGCP("exp", 4, var=0.2, scale=.1)
Z1 <- log(attr(X1, "Lambda"))
Z2 <- log(attr(rLGCP("exp", 4, var=0.1, scale=.1), "Lambda"))
Lam1 <- exp(Z1)
Lam2 <- exp((Z1 + Z2)/2)
Xlgcp <- superimpose(a=rpoispp(Lam1), b=rpoispp(Lam2))


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


###################################################
### code chunk number 180: 14multitype.Rnw:5009-5013
###################################################
plot(solist(Xlinked, Xlgcp, Xbalanced),
     main="", main.panel="", nrows=1, 
     equal.scales=TRUE, mar.panel=0, hsep=1,
     legend=FALSE, chars=c(16,3))


###################################################
### code chunk number 181: 14multitype.Rnw:5056-5057
###################################################
set.seed(29872)


###################################################
### code chunk number 182: 14multitype.Rnw:5059-5066
###################################################
ftcf <- function(x0, y0, radius, mu) {
  n <- rpois(1, mu)
  Y   <- runifdisc(n, radius=radius, centre=c(x0, y0))
  Z   <- factor(sample(1:2, npoints(Y), replace=TRUE), levels=1:2)
  return(Y %mark% Z)
}
Xmc2 <- rNeymanScott(15,0.1,ftcf, radius=0.1, mu=5)


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


###################################################
### code chunk number 184: 14multitype.Rnw:5073-5074
###################################################
plot(Xmc2, chars=c(1,16),main="")


###################################################
### code chunk number 185: 14multitype.Rnw:5128-5139
###################################################
V <- rmpoispp(c(50,50), types=letters[1:2])
Va <- split(V)$a
Vb <- split(V)$b
oka <- nncross(Va,Vb, what="dist") >0.07
okb <- nncross(Vb,Va, what="dist") >0.07
U <- superimpose(a=Va[oka], b=Vb[okb])
A <- rpoispp(50)
B <- rpoispp(50)
ok <- nncross(B,A,what="dist") > 0.07
B <- B[ok]
AB <- superimpose(a=A,b=B,W=owin())


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


###################################################
### code chunk number 187: 14multitype.Rnw:5146-5148
###################################################
plot(solist(U,AB), chars=c(1,16), 
     main.panel="", main="", mar.panel=0, hsep=1, legend=FALSE)


###################################################
### code chunk number 188: 14multitype.Rnw:5187-5195
###################################################
m  <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=owin())
set.seed(31852)
RL <- attr(rLGCP("gauss", m, var=0.15, scale =0.5),"Lambda")
RL <- cut(RL,4,labels=letters[1:4])
X  <- rpoispp(100,win=Window(RL))
Y <- X %mark% RL[X]
layerplotargs(RL) <- list(ribside="left")
layerplotargs(Y) <- list(leg.side="right")


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


###################################################
### code chunk number 190: 14multitype.Rnw:5198-5199
###################################################
setmargins(0,1,0,1)


###################################################
### code chunk number 191: 14multitype.Rnw:5203-5207
###################################################
plot(solist(RL, X, Y),
     main="", main.panel="",
     equal.scales=TRUE, valign=TRUE,
     mar.panel=0, hsep=1, box=TRUE)


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


###################################################
### code chunk number 193: 14multitype.Rnw:5467-5476 (eval = FALSE)
###################################################
## nt <- length(levels(marks(X)))
## Mf <- function(rwithin, rbetween) { 
##   R <- matrix(rbetween, nt, nt)
##   diag(R) <- rwithin
##   return(MultiStrauss(R))
## }
## df <- expand.grid(rwithin=seq(rmin, rmax, length=5),
##                   rbetween=seq(rmin, rmax, length=5))
## profilepl(df, Mf, X ~ trend)


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