# Goodness-of-fit testing for simple hypothesis
if(require("spatstat.explore", quietly=TRUE)) {
# Testing complete spatial randomness (CSR)
#==========================================
X <- unmark(spruces)
nsim <- 1999 # Number of simulations
nsim <- 19 # Number of simulations
# 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 <- curve_set(r=r, obs=obs, sim=sim)
# Perform the test
res <- global_envelope_test(cset, type="erl")
plot(res) + ggplot2::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) + ggplot2::ylab(expression(italic(L(r)-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 <- as.ppp(saplings, W=square(75))
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
r <- r[1:50]; rJ <- rJ[1:50]
# 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(L=curve_set_L, F=curve_set_F,
G=curve_set_G, J=curve_set_J))
plot(res)
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 <- curve_set(r=1:2, obs=as.vector(X), sim=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 <- curve_set(r=1:3, obs=c(X[,1],X[,2],t3(X,Y)), sim=rbind(t(Y), t3(Y,Y)))
res2 <- global_envelope_test(cset2)
plot(res2)
}
Run the code above in your browser using DataLab