### NMST543 Spatial Statistics, exercise 4 ###
##############################################


library(spatstat)

# Estimates of the K-function and other characteristics are variable,
# even if we know the true model we never see a perfect alignment
# with the theoretical curve.

# Try plotting several times an estimate of the K-function
# for simulated realizations of a Poisson process to get
# some intuition about the variability of the estimates.

plot(Kest(rpoispp(50),correction="translate"),ylim=c(0,0.20))


### Deviation tests ###
#######################

# Integral deviation test (Diggle-Cressie-Loosmore-Ford):
dclf.test(cells, Lest, nsim=39)   # p-value = 0.0250, rank = 1
dclf.test(cells, Lest, nsim=99)   # p-value = 0.0100, rank = 1
dclf.test(cells, Lest, nsim=999)  # p-value = 0.0010, rank = 1
dclf.test(cells, Lest, nsim=9999) # p-value = 0.0001, rank = 1

# Supremum deviation test (MAD = Maximum Absolute Deviation):
m <- mad.test(cells, Lest, verbose=FALSE, rinterval=c(0, 0.1), nsim=19)
m

# Extract the p-value:
m$p.value 

# Note that we can also use a variant with "savefuns=TRUE" or
# "savepatterns=TRUE" if we plan to recycle the simulations.


### Envelope tests ###
######################

# Point-wise envelopes:
Ep <- envelope(cells, Kest, nsim = 39, rank = 1, do.pwrong=TRUE)
plot(Ep, main = "Pointwise envelopes")

# Check the significance level of the envelope test (the "pwrong" value):
Ep

# "Global" envelopes (fewer simulations, in fact one-sided scalar MAD test):
Eg <- envelope(cells, Kest, nsim = 19, rank = 1, global = TRUE)
plot(Eg, main = "Global envelopes")

# Stabilization of the variance (K-function => L-function) results
# in a more powerful test:
EL <- envelope(cells, Lest, nsim = 19, rank = 1, global = TRUE)
plot(EL, main = "global envelopes of L(r)")


### Refined envelope test ###
#############################

# How many simulated realizations are needed in order to obtain
# a "reasonable" significance level of the envelope test?
E <- envelope(cells, Kest, nsim = 999, rank = 1, do.pwrong=TRUE)
plot(E, main = "Pointwise envelopes")
E  # pwrong = 0.112 (may be different when repeating the simulations)

E <- envelope(cells, Kest, nsim = 1999, rank = 1, do.pwrong=TRUE)
plot(E, main = "Pointwise envelopes")
E  # pwrong = 0.0685 (may be different when repeating the simulations)


### Exact envelope test ###
###########################

# Package for the "Global Envelope Tests"
library(GET)

pp <- unmark(spruces) # forget about the marks
plot(pp)

# The classical approach (we save the estimated functions):
env0 <- envelope(pp, fun="Lest", nsim=39, savefuns=TRUE, correction="translate")
plot(env0)

# Exact test with the same (number of) realizations:
res0 <- global_envelope_test(env0, type="erl")
plot(res0)

# Generate "nsim" realizations of a Poisson process, estimate the L-function
# both for the data and for the simulated realizations:
env <- envelope(pp, fun="Lest", nsim=1999, savefuns=TRUE, correction="translate")
plot(env)
res <- global_envelope_test(env, type="erl")
plot(res)


### More advanced work with the GET package ###
###############################################

# Form a new set of curves by specifying the range of values [r_min, r_max]:
curve_set <- crop_curves(env, r_min = 1, r_max = 7)

# Use the transformation L(r)-r for a better visualisation:
curve_set <- residual(curve_set, use_theo = TRUE)
plot(curve_set)

# Perform the exact envelope test:
res <- global_envelope_test(curve_set, type="erl")
plot(res)


### Combining more test functions together ###
##############################################

data(saplings)
X <- as.ppp(saplings, W = square(75))
plot(X)
nsim <- 499 # 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 saving the simulated patterns and 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
                       verbose = FALSE)

# The simulations can be obtained from the returned object:
simulations <- attr(env_L, "simpatterns")

# The other test functions F, G, J should be calculated for each simulated pattern:
env_F <- envelope(X, nsim = nsim,
                     simulate = simulations,
                     fun = "Fest", correction = "Kaplan", r = rJ,
                     savefuns = TRUE, verbose = FALSE)
env_G <- envelope(X, nsim = nsim,
                     simulate = simulations,
                     fun = "Gest", correction = "km", r = rJ,
                     savefuns = TRUE, verbose = FALSE)
env_J <- envelope(X, nsim = nsim,
                     simulate = simulations,
                     fun = "Jest", correction = "none", r = rJ,
                     savefuns = TRUE, verbose = FALSE)

# All the curves can be cropped to the desired r-interval I, if needed:
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)

# The combined global envelope is calculated and plotted:
res_combined <- global_envelope_test(curve_sets = list(curve_set_L, curve_set_F, curve_set_G, curve_set_J))
plot(res_combined, labels = c("L(r)-r", "F(r)", "G(r)", "J(r)"))




### Individual task: ###
########################

# Using different approaches test the hypothesis that the "japanesepines" dataset
# is a realization of a homogeneous Poisson process. Choose different types of tests
# and perhaps different types of summary characteristics.