Spatial Point Processes (R)¶
spatstat — Ripley's K, inhomogeneous Poisson (ppm), and LGCP (kppm)¶
The R counterpart to pointproc_python.ipynb, using spatstat (Baddeley, Rubak & Turner) — the definitive point-pattern package. envelope+Lest test complete spatial randomness, ppm fits the inhomogeneous Poisson process on covariates, and kppm fits the log-Gaussian Cox process. Same bei rainforest data (3,604 trees + elevation/gradient).
options(repr.plot.width=7.5, repr.plot.height=3.8)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(spatstat))
data(bei)
cat(bei$n, "trees on a", diff(bei$window$xrange), "x", diff(bei$window$yrange), "m plot\n")
plot(bei, pch=".", cols="darkgreen", main="bei: rainforest tree locations")
Warning message: "package 'spatstat' was built under R version 4.6.1"
Warning message: "package 'spatstat.data' was built under R version 4.6.1"
Warning message: "package 'spatstat.univar' was built under R version 4.6.1"
Warning message: "package 'spatstat.geom' was built under R version 4.6.1"
Warning message: "package 'spatstat.random' was built under R version 4.6.1"
Warning message: "package 'spatstat.explore' was built under R version 4.6.1"
Warning message: "package 'spatstat.model' was built under R version 4.6.1"
Warning message: "package 'spatstat.linnet' was built under R version 4.6.1"
3604 trees on a 1000 x 500 m plot
1. The pattern and a covariate¶
The tree locations and the elevation surface — density visibly tracks the terrain.
plot(bei.extra$elev, main="Elevation (m)"); plot(bei, add=TRUE, pch=".", cols="white")
2. Ripley's K/L against CSR¶
envelope with Lest computes the centred L-function and a simulation envelope from complete spatial randomness. The observed curve above the envelope means clustering.
set.seed(1)
E <- envelope(bei, Lest, nsim=19, correction="translate", verbose=FALSE)
plot(E, . - r ~ r, main="Ripley's L (centred) with CSR envelope", legend=FALSE)
cat("The observed L-r rises far above the CSR envelope: the trees are strongly clustered at all scales.\n")
The observed L-r rises far above the CSR envelope: the trees are strongly clustered at all scales.
3. Inhomogeneous Poisson with ppm¶
ppm(bei ~ elev + grad) fits a Poisson point process whose log-intensity is linear in elevation and gradient — the covariate-driven intensity surface.
fit <- ppm(bei ~ elev + grad, covariates=bei.extra)
print(round(coef(summary(fit))[,1:2], 3))
# the Python notebook standardises elev and grad, so convert to the same scale for comparison
sd_e <- sd(bei.extra$elev[]); sd_g <- sd(bei.extra$grad[])
ce <- coef(fit)
cat(sprintf("\nin standardised units (x sd of each covariate, as the Python notebook reports them):\n"))
cat(sprintf(" elev %.3f x %.3f = %.2f grad %.3f x %.3f = %.2f\n",
ce[["elev"]], sd_e, ce[["elev"]]*sd_e, ce[["grad"]], sd_g, ce[["grad"]]*sd_g))
plot(predict(fit), main="Inhomogeneous Poisson: fitted intensity (covariates only)")
cat("Intensity rises with elevation and gradient; but ppm's standard errors assume independence given covariates,\n")
cat("which the Ripley test contradicts -- so the Cox process is needed for honest inference and residual clustering.\n")
Estimate S.E. (Intercept) -8.564 0.341 elev 0.021 0.002 grad 5.846 0.256
in standardised units (x sd of each covariate, as the Python notebook reports them):
elev 0.021 x 8.056 = 0.17 grad 5.846 x 0.059 = 0.34
Intensity rises with elevation and gradient; but ppm's standard errors assume independence given covariates,
which the Ripley test contradicts -- so the Cox process is needed for honest inference and residual clustering.
4. Log-Gaussian Cox process with kppm¶
kppm(..., "LGCP") fits a log-Gaussian Cox process: the same covariate trend plus a latent Gaussian random field for clustering beyond the covariates. Its covariate estimates come with correct (wider) uncertainty, and it reports the clustering scale.
set.seed(1)
fitl <- kppm(bei ~ elev + grad, "LGCP", covariates=bei.extra)
print(round(coef(fitl), 3))
cat("\nkppm returns the SAME trend coefficients as ppm -- both solve the same Poisson estimating\n")
cat("equations, so the clustering changes the UNCERTAINTY rather than the point estimates:\n")
se_p <- sqrt(diag(vcov(fit))); se_k <- tryCatch(sqrt(diag(vcov(fitl))), error=function(e) rep(NA_real_, length(se_p)))
cat(sprintf(" %-12s ppm se %.4f kppm se %.4f %.1fx wider\n", names(se_p), se_p,
se_k[seq_along(se_p)], se_k[seq_along(se_p)]/se_p))
cat(sprintf("LGCP latent field: variance (sigma2) %.2f, scale %.0f m -- substantial clustering beyond terrain\n",
fitl$clustpar["var"], fitl$clustpar["scale"]))
plot(predict(fitl), main="LGCP fitted intensity (covariates + latent field)")
(Intercept) elev grad
-8.564 0.021 5.846
kppm returns the SAME trend coefficients as ppm -- both solve the same Poisson estimating
equations, so the clustering changes the UNCERTAINTY rather than the point estimates:
(Intercept) ppm se 0.3411 kppm se 3.6358 10.7x wider elev ppm se 0.0023 kppm se 0.0246 10.8x wider grad ppm se 0.2558 kppm se 2.9079 11.4x wider
LGCP latent field: variance (sigma2) 1.58, scale 48 m -- substantial clustering beyond terrain
5. Summary¶
spatstat reproduces the from-scratch analysis. envelope/Lest confirm the bei trees are strongly clustered relative to CSR; ppm fits the inhomogeneous Poisson intensity from elevation and gradient; and kppm(..., "LGCP") fits the log-Gaussian Cox process, capturing the substantial residual clustering the covariates miss and giving honest uncertainty — matching pointproc_python.ipynb, where the LGCP is the areal Poisson-CAR (BYM) model on a grid.
spatstat is the field-standard point-pattern toolkit and the package mirror of the from-scratch Ripley/Poisson/LGCP engine. Point processes complete the spatial trio alongside the areal (CAR/SAR) and geostatistical (kriging) subsections. The arc's final piece is the spatiotemporal capstone.