Tags :: Spatial Statistics

• Simulating in R
• References

Because \(\pi h^2\) is quadratic, it is more difficult to evaluate the K-function estimate visually. I.e. it does not center around 0.

Under CSR, \[ K(h) = \pi h^2 \rightarrow \sqrt{\frac{K(H)}{\pi}} = h \] therefore, as an alternative we define the L-function where \[ L(H) = \left( \frac{K(h)}{\pi} \right)^{\frac{1}{2}} \] which we estimate as \[ \hat{L}(H) = \left( \frac{\hat{K}(h)}{\pi} \right)^{\frac{1}{2}} \] An edge correction function, \(\hat{L}_{ec}(h)\) exists.

What is the expected L(h)?

Under CSR we would expect \[ \hat{L} = H \; \rightarrow \; \hat{L} - h = 0 \] if we plot \(\hat{L}(h) - h\) versus \(h\) and look for deviations about 0, we can make the following interpretations:

• Peaks in positive values at distance \(h\) imply clustering at \(h\)
• Troughs in negative values at distance \(h\) imply regularity at \(h\)

Simulating in R

We start with the same code we used to simulate the K-function with the spatial domain \(X = [-2, 2], Y=[0, 4]\) and an intensity of \(\lambda = 5\)

library(maps)
library(splancs)
set.seed(15)

area <- 4*4
lambda <- 5

N <- rpois(1, lambda * area) # simulate number of points to use
bbox <- as.points(
    c(-2, 2, 2, -2, -2),
    c(0, 0, 4, 4, 0)
)

pts.csr <- csr(bbox, N) # Generate from CSR process
polymap(bbox)
pointmap(pts.csr, add=TRUE)

h <- seq(0.1, 2.5, 0.1)
kpts <- splacs::khat(pts.csr, bbox, h) # calc K-function

plot(h, kpts, type="l")

we then transform our calculated \(\hat{K}\) into \(\hat{L}\) values and then plot.

lpts <- sqrt(kpts / pi) - h

plot(h, lpts, type="l")
abline(h=0, col=2)

Simulating envelope

We want approximate the the test distributions of interest and evaluate our envelope to allow inference on out L-function values. First we want to simulate several CSR processes on the same spatial domain with the same number of events.

# Simulate 999 iterations
nsim <- 999
sim <- matrix(0, nsim, length(h))

for (i in 1:nsim) {
    sim[i, ] <- sqrt(splancs::khat(
                splancs::csr(bbox, N), bbox, h) / pi
                 ) - h
}

From these simulations we can calcuate the envelope and perform inference on our L values.

env <- apply(sim, 2, range) # Calculate upper and lower bounds for envelope

# plot
matplot(h, t(env), type="l", lty=2, col=3, ylab="L")
lines(h, lpts)
abline(h=0, col=2)

Alternatively we can calculate the envelope using various quantiles instead of the min and max.

# 25% and 97.5% quantiles
env.quant <- apply(sim, 2, quantile, c(0.025, 0.975))

# plot
matplot(h,t(env),type="l",lty=2,col=3,ylab="L") matplot(h,t(env.quant),type="l",lty=3,col=4,add=TRUE)
lines(h,lpts)

References