# NOT RUN {
# Goodness-of-fit testing for simple hypothesis
if(require("spatstat", quietly=TRUE)) {
# Testing complete spatial randomness (CSR)
#==========================================
X <- unmark(spruces)
# }
# NOT RUN {
nsim <- 1999 # Number of simulations
# }
# NOT RUN {
# }
# NOT RUN {
# Illustration of general workflow for simple hypotheses
#=======================================================
# First illustrate the general workflow for the test by this example
# of CSR test for a point pattern X using the empirical L-function.
# Define the argument values at which the functions are evaluated
obs.L <- Lest(X, correction = "translate")
r <- obs.L[['r']]
# The test function for the data
obs <- obs.L[['trans']] - r
# Prepare simulations and calculate test functions for them at same r as 'obs'
sim <- matrix(nrow = length(r), ncol = nsim)
for(i in 1:nsim) {
sim.X <- runifpoint(ex = X) # simulation under CSR
sim[, i] <- Lest(sim.X, correction = "translate", r = r)[['trans']] - r
}
# Create a curve_set containing argument values, observed and simulated functions
cset <- create_curve_set(list(r = r, obs = obs, sim_m = sim))
# Perform the test
res <- global_envelope_test(cset, type="erl")
plot(res, ylab = expression(italic(hat(L)(r)-r)))
# Simple hypothesis for a point pattern utilizing the spatstat package
#=====================================================================
# Generate nsim simulations under CSR, calculate L-function for the data and simulations
env <- envelope(X, fun="Lest", nsim=nsim,
savefuns=TRUE, # save the functions
correction="translate", # edge correction for L
transform = expression(.-r), # centering
simulate=expression(runifpoint(ex=X))) # Simulate CSR
# The rank envelope test (ERL)
res <- global_envelope_test(env, type="erl")
# Plot the result
plot(res)
## Advanced use:
# Choose the interval of distances [r_min, r_max] (at the same time create a curve_set from 'env')
cset <- crop_curves(env, r_min=1, r_max=7)
# Do the rank envelope test (erl)
res <- global_envelope_test(cset, type="erl")
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[mm](r)-L(r))))
# 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
nsim <- 499 # Number of simulations
nsim <- 19 # Number of simulations
# 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)"))
# }
# NOT RUN {
}
# 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