# NOT RUN {
if(require("spatstat", quietly=TRUE)) {
data(saplings)
# First choose the r-distances
rmin <- 0.3; rmax <- 10; rstep <- (rmax-rmin)/500
r <- seq(0, rmax, by=rstep)
# Number of simulations
nsim <- 19 # Increase nsim for serious analysis!
# Option 1: Give the fitted model object to GET.composite
#--------------------------------------------------------
# This can be done and is preferable when the model is
# a point process model of spatstat.
# 1. Fit the Matern cluster process to the pattern
# (using minimum contrast estimation with the K-function)
M1 <- kppm(saplings~1, clusters = "MatClust", statistic="K")
summary(M1)
# 2. Make the adjusted global area rank envelope test using the L(r)-r function
adjenvL <- GET.composite(X = M1, nsim = nsim,
testfuns = list(L = list(fun="Lest", correction="translate",
transform = expression(.-r), r=r)), # passed to envelope
type = "area", r_min=rmin, r_max=rmax)
# Plot the test result
plot(adjenvL)
# Option 2: Generate the simulations "by yourself"
#-------------------------------------------------
# and provide them as curve_set or envelope objects
# Preferable when you want to have a control
# on the simulations yourself.
# 1. Fit the model
M1 <- kppm(saplings~1, clusters = "MatClust", statistic="K")
# 2. Generate nsim simulations by the given function using the fitted model
X <- envelope(M1, nsim=nsim, savefuns=TRUE,
fun="Lest", correction="translate",
transform = expression(.-r), r=r)
plot(X)
# 3. Create another set of simulations to be used to estimate
# the second-state p-value (as proposed by Baddeley et al., 2017).
simpatterns2 <- simulate(M1, nsim=nsim)
# 4. Calculate the functions for each pattern
simf <- function(rep) {
# Fit the model to the simulated pattern Xsims[[rep]]
sim_fit <- kppm(simpatterns2[[rep]], clusters = "MatClust", statistic="K")
# Generate nsim simulations from the fitted model
envelope(sim_fit, nsim=nsim, savefuns=TRUE,
fun="Lest", correction="translate",
transform = expression(.-r), r=r)
}
X.ls <- parallel::mclapply(X=1:nsim, FUN=simf, mc.cores=1) # list of envelope objects
# 5. Perform the adjusted test
res <- GET.composite(X=X, X.ls=X.ls, type="area",
r_min=rmin, r_max=rmax)
plot(res)
# Option 2: Generate the simulations "by yourself"
# Writing with for-loops (instead of envelope calls)
# Serves as an example for non-spatstat models.
#-------------------------------------------------
# 1. Fit the model
M1 <- kppm(saplings~1, clusters = "MatClust", statistic="K")
# Generate nsim simulations by the given function using the fitted model
# 2. Simulate a sample from the fitted null model and
# compute the test vectors for data and each simulation
nsim <- nsimsub <- 19 # Number of simulations
# Observed function
obs.L <- Lest(saplings, correction = "translate", r = r)
obs <- obs.L[['trans']] - r
# Simulated functions
sim <- matrix(nrow = length(r), ncol = nsimsub)
for(i in 1:nsimsub) {
sim.X <- simulate(M1)[[1]]
sim[, i] <- Lest(sim.X, correction = "translate", r = r)[['trans']] - r
}
# Create a curve_set
X <- create_curve_set(list(r = r, obs = obs, sim_m = sim))
# 3. Simulate another sample from the fitted null model.
sims2 <- simulate(M1, nsim=nsim) # list of simulations
# 4. Fit the null model to each of the patterns,
# simulate a sample from the null model,
# and compute the test vectors for all.
X.ls <- list()
for(rep in 1:nsim) {
M2 <- kppm(sims2[[rep]], clusters = "LGCP", statistic="K")
obs2 <- Lest(sims2[[rep]], correction = "translate", r = r)[['trans']] - r
sim2 <- matrix(nrow = length(r), ncol = nsimsub)
for(i in 1:nsimsub) {
sim2.X <- simulate(M2)[[1]]
sim2[, i] <- Lest(sim2.X, correction = "translate", r = r)[['trans']] - r
}
X.ls[[rep]] <- create_curve_set(list(r = r, obs = obs2, sim_m = sim2))
}
# 5. Perform the adjusted test
res <- GET.composite(X=X, X.ls=X.ls, type="area",
r_min=rmin, r_max=rmax)
plot(res)
# Option 3: Provide fitfun and simfun functions
#----------------------------------------------
# fitfun = a function to fit the model
# simfun = a function to make a simulation from the model
# calcfun = a function to calculate the test vector from data/simulation
fitf <- function(X) {
kppm(X, clusters = "MatClust", statistic="K")
}
simf <- function(M) {
simulate(M, nsim=1)[[1]]
}
calcf <- function(X) {
rmin <- 0.3; rmax <- 10; rstep <- (rmax-rmin)/500
r <- seq(0, rmax, by=rstep)
L <- Lest(X, correction = "translate", r = r)[['trans']] - r
L[r >= rmin & r <= rmax]
}
res <- GET.composite(X = saplings, nsim=nsim, fitfun = fitf, simfun=simf, calcfun=calcf,
type="area")
plot(res, xlab="Index")
# Note: r-values are not passed; r in the x-axis is from 1 to length of r.
}
# }
Run the code above in your browser using DataLab