### 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 a 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")
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)

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



### 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.