# NOT RUN {
if(require("spatstat", quietly=TRUE)) {
## Testing complete spatial randomness (CSR)
#-------------------------------------------
# }
# NOT RUN {
nsim <- 1999 # Number of simulations
# }
# NOT RUN {
# }
# NOT RUN {
pp <- unmark(spruces)
# Generate nsim simulations under CSR, calculate L-function for the data and simulations
env <- envelope(pp, fun="Lest", nsim=nsim,
savefuns=TRUE, # save the functions
correction="translate", # edge correction for L
simulate=expression(runifpoint(ex=pp))) # Simulate CSR
# The rank envelope test (extreme rank length (ERL) breaking of ties)
res <- global_envelope_test(env)
# Plot the result.
# - The central curve is now obtained from env[['theo']], which is the
# value of the L-function under the null hypothesis (L(r) = r).
# - Three different plot styles are provided:
# a) ggplot2 style (requires ggplot2)
plot(res, plot_style="ggplot2")
# b) spatstat's style (requires spatstat)
plot(res, plot_style="fv")
# c) a basic style
plot(res, plot_style="basic")
## Advanced use:
# Choose the interval of distances [r_min, r_max] (at the same time create a curve_set from 'env')
curve_set <- crop_curves(env, r_min=1, r_max=7)
# For better visualisation, take the L(r)-r function
curve_set <- residual(curve_set, use_theo=TRUE)
# Do the rank envelope test (erl)
res <- global_envelope_test(curve_set); plot(res, ylab=expression(italic(L(r)-r)))
# }
# NOT RUN {
## Random labeling test
#----------------------
mpp <- spruces
# 1) Perform simulations under the random labelling hypothesis and calculate
# the test function T(r) for the data pattern (mpp) and each simulation.
# The command below specifies that the test function is T(r) = \hat{L}_mm(r),
# which is an estimator of the mark-weighted L function, L_mm(r),
# with translational edge correction.
nsim <- 1999 # Number of simulations
env <- envelope(mpp, fun=Kmark, nsim = nsim, f=function(m1, m2) { m1*m2 },
correction="translate", returnL=TRUE,
simulate=expression(rlabel(mpp, permute=TRUE)), # Permute the marks
savefuns=TRUE) # Save the functions
# 2)
# Crop curves to desired r-interval
curve_set <- crop_curves(env, r_min=1.5, r_max=9.5)
# Center the functions, i.e. take \hat{L}_mm(r)-T_0(r).
# Below T_0(r) = \hat{L}(r) is the mean of simulated functions.
# (This is only for visualization, does not affect the test result.)
curve_set <- residual(curve_set)
# 3) Do the rank envelope test
res <- global_envelope_test(curve_set)
# 4) Plot the test result
plot(res, ylab=expression(italic(L[m](r)-L(r))))
## Goodness-of-fit test (typically conservative, see ?GET.composite for adjusted tests)
#-----------------------------------------------
pp <- unmark(spruces)
# Minimum distance between points in the pattern
min(nndist(pp))
# Fit a model
fittedmodel <- ppm(pp, interaction=Hardcore(hc=1)) # Hardcore process
# Simulating Gibbs process by 'envelope' is slow, because it uses the MCMC algorithm
#env <- envelope(fittedmodel, fun="Jest", nsim=999, savefuns=TRUE,
# correction="none", r=seq(0, 4, length=500))
# Using direct algorihm can be faster, because the perfect simulation is used here.
simulations <- NULL
nsim <- 999 # Number of simulations
for(j in 1:nsim) {
simulations[[j]] <- rHardcore(beta=exp(fittedmodel$coef[1]),
R=fittedmodel$interaction$par$hc,
W=pp$window)
if(j%%10==0) cat(j, "...", sep="")
}
env <- envelope(pp, simulate=simulations, fun="Jest", nsim=length(simulations),
savefuns=TRUE, correction="none", r=seq(0, 4, length=500))
curve_set <- crop_curves(env, r_min=1, r_max=3.5)
res <- global_envelope_test(curve_set, type="erl"); plot(res, ylab=expression(italic(J(r))))
# }
# NOT RUN {
# A combined global envelope test
#--------------------------------
# As an example test CSR of the saplings point pattern by means of
# L, F, G and J functions.
data(saplings)
X <- saplings
# }
# NOT RUN {
nsim <- 499 # Number of simulations
# }
# NOT RUN {
# }
# NOT RUN {
# Specify distances for different test functions
n <- 500 # the number of r-values
rmin <- 0; rmax <- 20; rstep <- (rmax-rmin)/n
rminJ <- 0; rmaxJ <- 8; rstepJ <- (rmaxJ-rminJ)/n
r <- seq(0, rmax, by=rstep) # r-distances for Lest
rJ <- seq(0, rmaxJ, by=rstepJ) # r-distances for Fest, Gest, Jest
# Perform simulations of CSR and calculate the L-functions
env_L <- envelope(X, nsim=nsim,
simulate=expression(runifpoint(ex=X)),
fun="Lest", correction="translate",
transform=expression(.-r), # Take the L(r)-r function instead of L(r)
r=r, # Specify the distance vector
savefuns=TRUE, # Save the estimated functions
savepatterns=TRUE) # Save the simulated patterns
# Take the simulations from the returned object
simulations <- attr(env_L, "simpatterns")
# Then calculate the other test functions F, G, J for each simulated pattern
env_F <- envelope(X, nsim=nsim,
simulate=simulations,
fun="Fest", correction="Kaplan", r=rJ,
savefuns=TRUE)
env_G <- envelope(X, nsim=nsim,
simulate=simulations,
fun="Gest", correction="km", r=rJ,
savefuns=TRUE)
env_J <- envelope(X, nsim=nsim,
simulate=simulations,
fun="Jest", correction="none", r=rJ,
savefuns=TRUE)
# Crop the curves to the desired r-interval I
curve_set_L <- crop_curves(env_L, r_min=rmin, r_max=rmax)
curve_set_F <- crop_curves(env_F, r_min=rminJ, r_max=rmaxJ)
curve_set_G <- crop_curves(env_G, r_min=rminJ, r_max=rmaxJ)
curve_set_J <- crop_curves(env_J, r_min=rminJ, r_max=rmaxJ)
res <- global_envelope_test(curve_sets=list(curve_set_L, curve_set_F,
curve_set_G, curve_set_J))
plot(res, labels=c("L(r)-r", "F(r)", "G(r)", "J(r)"))
}
## A test based on a low dimensional random vector
#-------------------------------------------------
# Let us generate some example data.
X <- matrix(c(-1.6,1.6),1,2) # data pattern X=(X_1,X_2)
if(requireNamespace("mvtnorm", quietly=TRUE)) {
Y <- mvtnorm::rmvnorm(200,c(0,0),matrix(c(1,0.5,0.5,1),2,2)) # simulations
plot(Y, xlim=c(min(X[,1],Y[,1]), max(X[,1],Y[,1])), ylim=c(min(X[,2],Y[,2]), max(X[,2],Y[,2])))
points(X, col=2)
# Test the null hypothesis is that X is from the distribution of Y's (or if it is an outlier).
# Case 1. The test vector is (X_1, X_2)
cset1 <- create_curve_set(list(r=1:2, obs=as.vector(X), sim_m=t(Y)))
res1 <- global_envelope_test(cset1)
plot(res1)
# Case 2. The test vector is (X_1, X_2, (X_1-mean(Y_1))*(X_2-mean(Y_2))).
t3 <- function(x, y) { (x[,1]-mean(y[,1]))*(x[,2]-mean(y[,2])) }
cset2 <- create_curve_set(list(r=1:3, obs=c(X[,1],X[,2],t3(X,Y)), sim_m=rbind(t(Y), t3(Y,Y))))
res2 <- global_envelope_test(cset2)
plot(res2)
}
# }
Run the code above in your browser using DataLab