Validation of spatial point process

validation_spp.R

library(spatstat)

#------- HOMOGENEOUS

# spatial point pattern from a homogeneous Poisson process with intensity lambda = 100 
# in the window A = [0,1]x[0k2] which has area |A| = 2

phom <- rpoispp(lambda = 100,
                win = owin(xrange = c(0, 1), yrange = c(0, 2)))
phom$n
plot(phom)

# test complete spatial randomness
envelope_K_csr <- envelope(phom, Kest, nsim = 99)
plot(envelope_K_csr, main = "CSR Test with Envelopes")

# black is the empirical K, red dashed line is theoretical CSR value
# if the black curve stays close to the red one, the pattern looks like CSR. 
# if it goes above, it suggests clustering; if below, inhibition.

# simulate a Thomas cluster process for comparison
clustered_pattern <- rThomas(kappa = 10, 
                             scale = 0.05, 
                             mu = 10, 
                             win = owin(xrange = c(0, 1), yrange = c(0, 2)))
plot(clustered_pattern, main = "Clustered Point Pattern (Thomas Process)")

# test
envelope_K_clustered <- envelope(clustered_pattern, Kest, nsim = 99)
plot(envelope_K_clustered, main = "Clustered Pattern with Envelopes")

# black line outside envelope, and above, therefore clustering

#-------- INHOMOGENEOUS

# we consider an intensity function lambda(x,y) = 10 + 100 * x + 200 * y 
# and an observation window equal to [0,1] x [0,2]

# intensity function
fnintensity <- function(x, y){return(10 + 100 * x + 200 * y)}

# grid
xseq <- seq(0, 1, length.out = 50)
yseq <- seq(0, 2, length.out = 100)
grid <- expand.grid(xseq, yseq)

# evaluation of the function on a grid
z <- mapply(fnintensity, grid$Var1, grid$Var2)

# plot
library(fields)
zmat <- matrix(z, 50, 100)
fields::image.plot(xseq, yseq, zmat, xlab = "x", ylab = "y",
                   main = "lambda(x, y)", asp = 1)

# now simulate
pinhom <- rpoispp(lambda = fnintensity,
                  win = owin(xrange = c(0, 1), yrange = c(0, 2)))
pinhom$n
plot(pinhom, main = "Inhomogeneous")

# naive test 
envelope_K_clustered <- envelope(pinhom, Kest, nsim = 99)
plot(envelope_K_clustered, main = "Naive K-function (assumes homogeneity)")
# we see the empirical K above the CSR line — but this is misleading

# need to estimate the intensity at each observed point
# one option is a kernel estimate using density():

# estimate intensity at each point location
sigma_val <- bw.diggle(pinhom)
lambda_hat <- density(pinhom, sigma = sigma_val, positive = TRUE)  # automatic bandwidth
env_inhom <- envelope(pinhom, Kinhom, nsim = 99,
                simulate = expression(rpoispp(lambda_hat)),
                sigma = sigma_val)
plot(env_inhom, main = "Inhomogeneous K-function with Envelopes")

# Now we're asking: 
# "Given the inhomogeneous intensity, are the points still more/less clustered than expected?"

# if from a fitted model, just replace lambda_hat in rpoispp(lambda_hat)
# by the estimated lambda