Week 8: Point Processes: Null Models, Clustering, and Regularity

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2025

Author

Affiliation

Jeff Jacobs

jj1088@georgetown.edu

Open slides in new tab →

[Naïve] Spatial Autocorrelation

• Where we left off: Location \(i\) has high value \(\implies\) locations near \(i\) more likely to have high values

source("../dsan-globals/_globals.r")
set.seed(6805)
library(tidyverse) |> suppressPackageStartupMessages()
library(sf) |> suppressPackageStartupMessages()
library(spatstat) |> suppressPackageStartupMessages()
library(mdthemes) |> suppressPackageStartupMessages()
library(mapview) |> suppressPackageStartupMessages()
library(leaflet) |> suppressPackageStartupMessages()

Moran’s \(I\)

\[ I = \underset{\text{Inverse of Variance}}{\boxed{\frac{n}{\sum_{i=1}^{n}(y_i - \overline{y})^2}}} \frac {\overbrace{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(y_i - \overline{y})(y_j - \overline{y})}^{\text{Weighted Covariance}}} {\underbrace{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}}_{\text{Normalize Weights}}} \]

• \(I\) is Large when:
  • \(y_i\) and \(y_j\) are neighbors: \(w_{ij}\), and
  • \(y_i\) and \(y_j\) large at the same time: \((y_i - \overline{y})(y_j - \overline{y})\)

Local Indicators of Spatial Autocorrelation (LISA)

• tldr: See how much a given location \(\mathbf{s}\) contributes to overall Moran’s \(I\)
• Local Moran’s \(I\):

\[ I_i = \frac{y_i - \overline{y}}{S_i^2}\sum_{j=1}^{n}w_{ij}(y_j - \overline{y}) \]

Agenda for Point Patterns: Something “Noteworthy” is Happening if Pattern Unlikely to Arise Randomly Under CSR!

Less clustering than expected

Roughly expected amount of clustering under CSR

More clustering than expected

Is this likely to arise under CSR?
But is it because of repulsion?
BBC, Aug 10 2025

Is this likely to arise under CSR?
Is it because of homophily?
(Realization \(\omega_1\) in red, \(\omega_2\) in blue)

…Or, [Non-Human] Nature Examples

“Fairy circles” in Namibia Smithsonian Magazine: Namibia’s mysterious fairy circles might actually be caused by competition 👀 between grasses

Termite mounds in Mozambique

Free stock photo by @natalinadmay

Three kestral (bird) colonies form a triangle straddling the Basilicata and Apulia regions of Italy; from Morinay et al. (2023)

But, a Challenge…

Typically depends on context and/or level of resolution and/or level of analysis

From Waller and Gotway (2004)

Naïve Clustering

Weeks 8-9 (My labels)	Naïve Clustering	→	Point Process Models	→	Autocorrelation in Point Data	→	Autocorrelation in Lattice Data
Waller and Gotway (2004)	(A bunch of stuff we already learned)	→	Ch 5: Analysis of Spatial Point Patterns	→	Ch 6: Point Data (Cases and Controls)	→	Ch 7: Regional Count Data (First mention of autocorrelation on pg. 227)
Schabenberger and Gotway (2004)	Ch 1 (pg. 14): Autocorrelation	→	(A bunch of more complicated stuff we’ll learn at the end)

Why “Naïve”?

• Regional data: easiest way to visualize / operationalize simple measures of clustering via spatial autocorrelation
• So, we develop intuitions for measures like Moran’s \(I\) by looking at lattice data, but…
• Naïve because we’re not modeling the lattice regions (cells) themselves!
• Non-naïve clustering: Knowing the clusters but also what process led to their formation! Which brings us to…

…Point Process Models!

• Our New Library: spatstat!
• Homepage: spatstat.org
• GitHub: github.com/spatstat
• Book: Baddeley, Rubak, and Turner (2015) [Companion website]
• PDF: here

Spatial Randomness

library(tidyverse)
library(spatstat)
set.seed(6805)
N <- 60
r_core <- 0.05
obs_window <- square(1)
# Regularity via Inhibition
#reg_sims <- rMaternI(N, r=r_core, win=obs_window)
cond_reg_sims <- rSSI(r=r_core, N)
# CSR data
#csr_sims <- rpoispp(N, win=obs_window)
cond_sr_sims <- rpoint(N, win=obs_window)
### Clustered data
#clust_sims <- rMatClust(kappa=6, r=2.5*r_core, mu=10, win=obs_window)
#clust_sims <- rMatClust(mu=5, kappa=1, scale=0.1, win=obs_window, n.cond=N, w.cond=obs_window)
#clust_sims <- rclusterBKBC(clusters="MatClust", kappa=10, mu=10, scale=0.05, verbose=FALSE)
# Each cluster consist of 10 points in a disc of radius 0.2
nclust <- function(x0, y0, radius, n) {
    #print(n)
    return(runifdisc(10, radius, centre=c(x0, y0)))
}
cond_clust_sims <- rNeymanScott(kappa=5, expand=0.0, rclust=nclust, radius=2*r_core, n=10)
# And PLOT
plot_w <- 400
plot_h <- 400
plot_scale <- 2.25
cond_reg_plot <- cond_reg_sims |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  dsan_theme()
ggsave("images/cond_reg.png", cond_reg_plot, width=plot_w, height=plot_h, units="px", scale=plot_scale)
cond_sr_plot <- cond_sr_sims |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  dsan_theme()
ggsave("images/cond_sr.png", cond_sr_plot, width=plot_w, height=plot_h, units="px", scale=plot_scale)
cond_clust_plot <- cond_clust_sims |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  dsan_theme()
ggsave("images/cond_clust.png", cond_clust_plot, width=plot_w, height=plot_h, units="px", scale=plot_scale)

Autocorrelation	\(I = -1\)	←	\(I = 0\)	→	\(I = 1\)
Description	Negative Autocorr		No Autocorr		Positive Autocorr
Event at \(\mathbf{s} = (x,y)\) Implies	Less likely to find another point nearby		No information about nearby points		More likely to find another point nearby
Resulting Pattern	Regularity		Reg/Clustered Mix		Clustering
Process(es) Which Could Produce Pattern	1st Order: Random within even-spaced grid 2nd Order: Competition		1st Order: i.i.d. points 2nd Order: i.i.d. distances		1st Order: Tasty food at clust centers 2nd Order: Cooperation
Fixed \(N\)	60		60		60
Problem with fixed \(N\): Learning there are \(n_1\) points in region \(R_1\) \(\Rightarrow\) there are \(N - n_1\) points outside \(R_1\)

Complete Spatial Randomness (CSR)

library(tidyverse)
library(spatstat)
set.seed(6807)
lambda <- 60
r_core <- 0.05
obs_window <- square(1)
# Regularity via Inhibition
# Regularity via Inhibition
reg_sims <- rMaternI(lambda, r=r_core, win=obs_window)
# CSR data
csr_sims <- rpoispp(N, win=obs_window)
### Clustered data
clust_mu <- 10
clust_sims <- rMatClust(kappa=lambda / clust_mu, scale=2*r_core, mu=10, win=obs_window)
# And PLOT
plot_w <- 400
plot_h <- 400
plot_scale <- 2.25
reg_plot <- reg_sims |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  labs(title=paste0("N = ",reg_sims$n)) +
  dsan_theme()
ggsave("images/reg.png", reg_plot, width=plot_w, height=plot_h, units="px", scale=plot_scale)
csr_plot <- csr_sims |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  labs(title=paste0("N = ",csr_sims$n)) +
  dsan_theme()
ggsave("images/csr.png", csr_plot, width=plot_w, height=plot_h, units="px", scale=plot_scale)
clust_plot <- clust_sims |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  labs(title=paste0("N = ",clust_sims$n)) +
  dsan_theme()
ggsave("images/clust.png", clust_plot, width=plot_w, height=plot_h, units="px", scale=plot_scale)

Autocorrelation	\(I = -1\)	←	\(I = 0\)	→	\(I = 1\)
Description	Negative Autocorr		No Autocorr		Positive Autocorr
Event at \(\mathbf{s} = (x,y)\) Implies	Less likely to find another point nearby		No information about nearby points		More likely to find another point nearby
Resulting Pattern	Regularity		Reg/Clustered Mix		Clustering
Process(es) Which Could Produce Pattern	1st Order: Random within even-spaced grid 2nd Order: Competition		1st Order: i.i.d. points 2nd Order: i.i.d. distances		1st Order: Tasty food at clust centers 2nd Order: Cooperation
Fixed Intensity \(\lambda\)	60		60		60
Random \(N\)

Caveat 1: Measures are Relative to Window of Observation

Same data can be spatially random within one window, clustered or regular in others!

N <- 60
obs_window <- square(1)
window_scale <- 3.5
csr_sims_square <- rpoispp(N, win=obs_window)
# Triangular window
obs_window_tri <- st_sfc(st_polygon(list(
    matrix(c(0.3,0.1,0.7,0.1,0.5,0.5,0.3,0.1), byrow=TRUE, nrow=4)
)))
obs_window_geom <- st_sfc(st_linestring(
    matrix(c(0,0,1,0,1,1,0,1,0,0), byrow=TRUE, nrow=5)
))
csr_sims_tri <- ppp(csr_sims_square$x, csr_sims_square$y, window=as.owin(obs_window_tri))
tri_plot <- csr_sims_tri |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  dsan_theme("quarter");
ggsave("images/window_tri.png", tri_plot, width=plot_w, height=plot_h, units="px", scale=window_scale)
# Square window
square_plot <- csr_sims_square |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  geom_sf(data=obs_window_tri |> sf::st_boundary()) +
  dsan_theme("quarter")
ggsave("images/window_square.png", square_plot, width=plot_w, height=plot_h, units="px", scale=window_scale)
# Circular window
obs_window_disc <- st_sfc(st_point(c(1, 0.5))) |> st_buffer(1.2)
csr_sims_circ <- ppp(csr_sims_square$x, csr_sims_square$y, window=as.owin(obs_window_disc))
circ_plot <- csr_sims_circ |> sf::st_as_sf() |>
  ggplot() +
  geom_sf() +
  geom_sf(data=obs_window_geom) +
  dsan_theme("quarter")
ggsave("images/window_circ.png", circ_plot, width=plot_w, height=plot_h, units="px", scale=window_scale)

Regular	CSR	Clustered

Caveat 2: Summary Statistics Like \(I\) are Not Models!

• Moran’s \(I\) is to GISers what a thermometer is to doctors
• Measures symptoms; many possible underlying causes!
• Need to ask why autocorrelation seems to be present!

[] Why Do Events Appear Where They Do?

center_l_function <- function(x, ...) {
  if (!spatstat.geom::is.ppp(x) && !spatstat.geom::is.fv(x)) {
    stop("Please provide either ppp or fv object.")
  }
  if (spatstat.geom::is.ppp(x)) {
    x <- spatstat.explore::Lest(x, ...)
  }
  r <- x$r
  l_centered <- spatstat.explore::eval.fv(x - r)
  return(l_centered)
}
cond_clust_sf <- cond_clust_sims |> sf::st_as_sf()
pines_plot <- cond_clust_sf |>
  ggplot() +
  geom_sf() +
  theme_classic(base_size=12)
ggsave("images/pines.png", pines_plot)
# density() calls density.ppp() if the argument is a ppp object
den <- density(cond_clust_sims, sigma = 0.1)
#summary(den)
png("images/intensity_plot.png")
plot(den, main = "Intensity λ(s)")
contour(den, add = TRUE) # contour plot
dev.off()
# And Gest / Kest / Lest
# saveRDS(cond_clust_sims, "cond_clust_sims.rds")
pcf_result <- spatstat.explore::pcf.ppp(
  cond_clust_sims,
  divisor="d",
  stoyan=0.25,
  bw=0.05,
  r=seq(from=0.0, to=1.0, by=0.001)
)
png("images/spatstat_pcf.png")
plot(pcf_result, xlim=c(0, 1), main="pcf")
dev.off()
# kest_result <- Kest(cond_clust_sims, rmax=0.5, correction="best")
# lest_result <- center_l_function(cond_clust_sims, rmax=0.5)
# png("images/lest.png")
# plot(lest_result, main="K(h)")
# dev.off()

	First-Order	Second-Order
	Events considered individually \(\Rightarrow\) Intensity function \(\lambda(\mathbf{s})\)	Events considered pairwise \(\Rightarrow\) Pairwise Correlation Function \(\textrm{pcf}(\vec{h})\)

The Tree-Grid Mystery

You’ve been hired as an archaeologist – congratulations! Your job: determine whether arrangement of trees formed:

• Naturally, via a process of resource competition, or
• Artificially, via an ancient civilization planting in a grid…

Two Possible Histories…

Hypothesis \(\mathcal{H}_{\textsf{Art}}\): Artificial Formation

ha_base <- 28
square_sf <- sf::st_as_sf(spatstat.geom::square(1))
grid_sf <- sf::st_as_sf(sf::st_make_grid(square_sf))
grid_buffer_sf <- grid_sf |> sf::st_buffer(dist=-0.01, singleSide = TRUE)
grid_buffer_sf |> ggplot() +
  geom_sf() +
  theme_classic(base_size=ha_base)

grid_points <- sf::st_sample(grid_buffer_sf, size=rep(1,100))
grid_buffer_sf |> ggplot(aes(shape='Cell')) +
  geom_sf() +
  geom_sf(data=grid_points) +
  scale_shape_manual("Shape", values=c('Cell'=19)) +
  theme_classic(base_size=ha_base) +
  theme(
    legend.title = element_blank(),
    # legend.text = element_text(size=18)
  )

grid_ppp <- as.ppp(grid_points, W=spatstat.geom::square(1))
grid_ppp |> sf::st_as_sf() |> ggplot() +
  geom_sf() +
  theme_classic(base_size=ha_base)

Hypothesis \(\mathcal{H}_{\textsf{Nat}}\): Natural Formation

hn_base <- 28
r <- 0.05
pois_ppp <- rpoispp(150)
pois_sf <- pois_ppp |> sf::st_as_sf()
pois_sf |> ggplot() +
  geom_sf() +
  theme_classic(base_size=hn_base)

age <- runif(npoints(pois_ppp))
pair_dists <- pairdist(pois_ppp)
close <- (pair_dists < r)
later <- outer(age, age, ">")
killed <- apply(close & later, 1, any)
killed_ppp <- pois_ppp[killed]
alive_ppp <- pois_ppp[!killed]
pois_window_sf <- pois_ppp |> sf::st_as_sf() |> filter(label=="window")
pois_killed_sf <- killed_ppp |> sf::st_as_sf() |> filter(label=="point")
pois_alive_sf <- alive_ppp |> sf::st_as_sf() |> filter(label=="point")
alive_buff_sf <- pois_alive_sf |> sf::st_buffer(r) |> sf::st_union() |> sf::st_intersection(pois_window_sf)
ggplot() +
  geom_sf(data=pois_window_sf) +
  geom_sf(data=alive_buff_sf, aes(color='Inhibition', shape='Inhibition'), linetype='dashed') +
  geom_sf(data=pois_killed_sf, aes(color='Dead', shape='Dead'), size=2, stroke=2) +
  geom_sf(data=pois_alive_sf, aes(color='Alive', shape='Alive'), size=1, stroke=1) +
  scale_shape_manual(name=NULL, values=c("Alive"=19, "Dead"=4, 'Inhibition'=21), labels=c("Alive", "Dead", "Inhibition")) +
  scale_color_manual(name=NULL, values=c("Alive"="black", "Dead"=cb_palette[1], "Inhibition"="black"), labels=c("Alive", "Dead", "Inhibition")) +
  guides(shape=guide_legend(override.aes=list(fill = "white"))) +
  theme_classic(base_size = hn_base) +
  theme(plot.margin = unit(c(0,0,0,0), "cm"))

alive_ppp |> sf::st_as_sf() |> ggplot() +
  geom_sf() +
  theme_classic(base_size=hn_base)

What Tools Do We Have for Distinguishing Between These Cases?

• “First-Order” measures vs. “Second-Order” measures

Why Do Events Appear Where They Do?

library(tidyverse)
library(spatstat)
set.seed(6809)
N <- 60
r_core <- 0.05
obs_window <- square(1)
### Clustered data
clust_ppp <- rMatClust(
  kappa=6,
  scale=r_core,
  mu=10
)
clust_sf <- clust_ppp |> sf::st_as_sf()
clust_plot <- clust_sf |>
  ggplot() +
  geom_sf(size=2) +
  theme_classic(base_size=18)
ggsave("images/clust_ppp.png", clust_plot, width=3, height=3)
# Intensity fn
clust_intensity <- density(clust_ppp, sigma = 0.1)
png("images/clust_intensity.png")
par(mar=c(0,0,0,2), las=2, oma=c(0,0,0,0), cex=2)
plot(clust_intensity, main=NULL)
contour(clust_intensity, add = TRUE)
dev.off()
### PCF
clust_pcf <- spatstat.explore::pcf(
  clust_ppp, divisor="d",
  r=seq(from=0.00, to=0.50, by=0.01)
)
clust_pcf_plot <- clust_pcf |> ggplot(aes(x=r, y=iso)) +
  geom_hline(yintercept=1, linetype='dashed', linewidth=1) +
  geom_area(color='black', fill=cb_palette[1], alpha=0.75) +
  scale_x_continuous(breaks=seq(from=0.0, to=1.0, by=0.1)) +
  labs(x="Distance", y="Density") +
  theme_classic(base_size=14)
ggsave("images/clust_pcf.png", clust_pcf_plot, width=3, height=3)

Original Data	First-Order	Second-Order
\(N = 60\) Events	Events modeled individually \(\Rightarrow\) Intensity Function \(\lambda(\mathbf{s})\)	Events modeled pairwise \(\Rightarrow\) Pairwise-Corr Function \(\textrm{pcf}(\vec{h})\)

What Do These Functions “Detect”?

sq_base <- 16
sq_psize <- 2.5
obs_window <- square(1)
r0 <- 0.2
sq_df <- tibble::tribble(
  ~x, ~y,
  0.5-r0,0.5-r0,
  0.5+r0,0.5+r0,
  0.5-r0,0.5+r0,
  0.5+r0,0.5-r0
)
sq_sf <- sf::st_as_sf(
  sq_df,
  coords = c("x","y")
)
sq_ppp <- as.ppp(sq_sf, W=obs_window)
sq_ppp |> sf::st_as_sf() |> ggplot() +
  geom_sf(size=sq_psize) +
  theme_classic(base_size=sq_base)

par(mar=c(0,0,0,2), las=2, oma=c(0,0,1,0))
sq_density <- density(sq_ppp)
plot(sq_density, main=NULL, xaxs='i', yaxs='i')
contour(sq_density, xaxs='i', yaxs='i', add = TRUE)

### PCF
pcf_result <- spatstat.explore::pcf(
  sq_ppp,
  divisor="d",
  r=seq(from=0.00, to=0.8, by=0.01)
)
pcf_result |> ggplot(aes(x=r, y=iso)) +
  geom_hline(yintercept=1, linetype='dashed', linewidth=1.5) +
  geom_area(color='black', fill=cb_palette[1], alpha=0.75) +
  scale_x_continuous(breaks=seq(from=0.0, to=1.0, by=0.1)) +
  theme_classic(base_size=sq_base)

csr_ppp <- spatstat.random::rpoispp(
  lambda = 60,
  win=obs_window
)
csr_ppp |> sf::st_as_sf() |> ggplot() +
  geom_sf(size=sq_psize) +
  theme_classic(base_size=sq_base)

par(
  mar=c(0,0,0,2),
  las=2,
  oma=c(0,0,1,0)
)
csr_density <- density(csr_ppp)
plot(csr_density, main=NULL, xaxs="i", yaxs="i")
contour(csr_density, xaxs="i", yaxs="i", add = TRUE, lwd=1.5)

csr_pcf_result <- spatstat.explore::pcf(
  csr_ppp,
  divisor="d",
  r=seq(from=0.00, to=0.8, by=0.01)
)
csr_pcf_result |> ggplot(aes(x=r, y=iso)) +
  geom_hline(yintercept=1, linetype='dashed', linewidth=1.5) +
  geom_area(color='black', fill=cb_palette[1], alpha=0.75) +
  scale_x_continuous(breaks=seq(from=0.0, to=1.0, by=0.1)) +
  theme_classic(base_size=sq_base)

A Menagerie of Models

Open Interactive (WebR) Version

Poisson Point Processes (CSR)

spatstat.random::rpoispp(lambda, win)

• \(N \sim \text{Pois}(\lambda)\)
• For \(i \in \{1, \ldots, N\}\):
  • Generate \(X_i, Y_i \sim \mathcal{U}(\texttt{win})\)

sim_base <- 22
sim_psize <- 2
sim_xticks <- seq(from=0.0, to=1.0, by=0.2)
sim_yticks <- seq(from=0.0, to=1.0, by=0.2)
gen_pois_df <- function(num_sims=1) {
  pois_sims <- spatstat.random::rpoispp(
    lambda = 60, nsim=num_sims
  )
  return(tibble::as_tibble(pois_sims))
}
#pois_dfs <- gen_pois_df()
#pois_dfs |> head()
pois_sims <- spatstat.random::rpoispp(
  lambda = 60, nsim=3
)
to_sim_df <- function(cur_sim, sim_name) {
  cur_df <- tibble::as_tibble(cur_sim) |> mutate(sim=sim_name)
  return(cur_df)
}
combined_df <- imap(.x=pois_sims, .f=to_sim_df) |> bind_rows()
combined_df |> ggplot(aes(x=x, y=y)) +
  geom_point(size=sim_psize) +
  facet_wrap(vars(sim)) +
  coord_equal() +
  theme_classic(base_size=sim_base) +
  theme(panel.spacing.x = unit(2, "lines")) +
  scale_x_continuous(breaks=sim_xticks) +
  scale_y_continuous(breaks=sim_yticks)

Simple Sequential Inhibition (SSI)

spatstat.random::rSSI(r, n=Inf, giveup=1000, win)

• \(\mathbf{S} = \varnothing\)
• While not done:
  • Generate \(\mathbf{E} = (X, Y) \sim \mathcal{U}(\texttt{win})\)
  • Check if \(\mathbf{E}\) within r units of any existing point in \(\mathbf{S}\)
    • If it is, throw \(\mathbf{E}\) away. Otherwise, add \(\mathbf{E}\) to \(\mathbf{S}\)
  • done=TRUE if \(\mathbf{S}\) has n points OR has been the same for giveup steps

capture.output(ssi_sims <- spatstat.random::rSSI(
  r = 0.05, n=60, nsim=3
), file=nullfile())
combined_df <- imap(.x=ssi_sims, .f=to_sim_df) |> bind_rows()
combined_df |> ggplot(aes(x=x, y=y)) +
  geom_point(size=sim_psize) +
  facet_wrap(vars(sim)) +
  coord_equal() +
  theme_classic(base_size=24) +
  theme(panel.spacing.x = unit(3, "lines")) +
  scale_x_continuous(breaks=sim_xticks) +
  scale_y_continuous(breaks=sim_yticks)

Matérn Cluster Process

spatstat.random::rMatClust(kappa, scale, mu, win)

Docs

• Generate \(K(\kappa)\) parent points via Poisson Point Process with intensity \(\lambda = \kappa\)
• For each parent point \(\mathbf{s}_i \in \left\{\mathbf{s}_1, \ldots, \mathbf{s}_{K(\kappa)}\right\}\):
  • Generate \(N(\mu)\) offspring points via Poisson Point Process with intensity \(\lambda = \mu\), distributed uniformly within a circle of radius scale centered at \(\mathbf{s}_i\)
• Offspring points form the outcome (parent points are thrown away)

matclust_sims <- rMatClust(
  kappa = 6,
  scale = 0.075,
  mu = 10,
  nsim = 3
)
matclust_df <- imap(.x=matclust_sims, .f=to_sim_df) |> bind_rows()
matclust_plot <- matclust_df |> ggplot(aes(x=x, y=y)) +
  geom_point(size=sim_psize) +
  facet_wrap(vars(sim), nrow=1) +
  coord_equal() +
  theme_classic(base_size=sim_base) +
  theme(panel.spacing.x = unit(2, "lines")) +
  scale_x_continuous(breaks=sim_xticks) +
  scale_y_continuous(breaks=sim_yticks)
matclust_plot

Matérn Inhibition Process (I and II)

spatstat.random::rMaternI(kappa, r, win)
spatstat.random::rMaternII(kappa, r, win)

rMaternI() [Docs]	rMaternII() [Docs]
Generate events \(\mathbf{S} = \{\mathbf{s}_1, \ldots, \mathbf{s}_{N(\lambda)}\}\) via Poisson point process with \(\lambda = \kappa\)	Generate events \(\mathbf{S} = \{\mathbf{s}_1, \ldots, \mathbf{s}_{N(\lambda)}\}\) via Poisson point process with \(\lambda = \kappa\), plus timestamp \(t_i \sim \mathcal{U}(0,1)\) for each \(\mathbf{s}_i\)
Delete all pairs of points \(\mathbf{s}_i\), \(\mathbf{s}_j\) for which \(\textsf{dist}(\mathbf{s}_i, \mathbf{s}_j) < \texttt{r}\)	For each pair of points \(\mathbf{s}_i\), \(\mathbf{s}_j\): if \(\textsf{dist}(\mathbf{s}_i, \mathbf{s}_j) < \texttt{r}\), delete the later point

mI_sims <- rMaternI(
  kappa = 60, r = 0.075, nsim=2
)
mII_sims <- rMaternII(
  kappa = 60, r = 0.075, nsim=2
)
mI_combined_df <- imap(.x=mI_sims, .f=to_sim_df) |> bind_rows() |> mutate(sim=paste0("MI ",sim))
mII_combined_df <- imap(.x=mII_sims, .f=to_sim_df) |> bind_rows() |> mutate(sim=paste0("MII ",sim))
m_combined_df <- bind_rows(mI_combined_df, mII_combined_df)
m_plot <- m_combined_df |> ggplot(aes(x=x, y=y)) +
  geom_point(size=sim_psize) +
  facet_wrap(vars(sim), nrow=1) +
  coord_equal() +
  theme_classic(base_size=sim_base) +
  theme(panel.spacing.x = unit(2, "lines")) +
  scale_x_continuous(breaks=sim_xticks) +
  scale_y_continuous(breaks=sim_yticks)
m_plot

Cox Processes: Random Intensity → Random Events

spatstat.random::rLGCP(model, mu, param, win)
models=c("exponential", "gauss", "stable", "gencauchy", "matern")

Docs

cox_pcol <- "black"
cox_bg <- "grey90"
cox_pch <- 21
# inhomogeneous LGCP with Gaussian covariance function
m <- as.im(function(x, y){
  5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2
}, W=owin())
lgcp_sims <- rLGCP("gauss", m, var=0.15, scale =0.5, nsim=3)

Warning: 32750 out of 65536 terms (50%) in FFT calculation of matrix square
root were negative, and were set to zero. Range: [-2.96, 1910]

# lgcp_combined_df <- imap(.x=ssi_sims, .f=to_sim_df) |> bind_rows()
plot_lgcp <- function(lgcp_sim) {
  plot(attr(lgcp_sim, "Lambda"), main=NULL)
  points(lgcp_sim, col=cox_pcol, bg=cox_bg, pch=cox_pch)
}
par(mfrow=c(1,3), mar=c(0,0,0,2), oma=c(0,0,0,0), las=2)
nulls <- lapply(X=lgcp_sims, FUN=plot_lgcp)

HW3 Tips!

• ppp Objects
• Converting between ppp and sf
• Plotting

[From Last Week] Our New Library: spatstat!

• Homepage: spatstat.org
• GitHub: github.com/spatstat
• Book: Baddeley, Rubak, and Turner (2015) [Companion website]
• PDF: here

ppp Objects

• The main datatype used to represent Planar Point Patterns [spatstat book p. 41]
• Unlike sf objects, which contain data+geometries for any desired collection of \(N\) entities, ppp objects are required to have at least an observation window!

sf Creation:

tree_df <- tibble::tibble(lon=runif(100,0,1), lat=runif(100,0,1), age=runif(100,0,1))
tree_sf <- sf::st_as_sf(
  tree_df,
  coords = c('lon', 'lat')
)
tree_sf |> head(4)

age	geometry
0.8593935	POINT (0.5035446 0.05519595)
0.9770798	POINT (0.2606818 0.3874983)
0.9517643	POINT (0.5464588 0.8271631)
0.0592741	POINT (0.7999682 0.5852131)

ppp Creation:

pois_ppp <- spatstat.random::rpoispp(
  lambda=100, win=spatstat.geom::square(1)
)
pois_ppp

Planar point pattern: 93 points
window: rectangle = [0, 1] x [0, 1] units

attributes(pois_ppp)$names

[1] "window"     "n"          "x"          "y"          "markformat"

pois_ppp$x |> head(4)

[1] 0.51779807 0.62213583 0.84467836 0.05354204

ppp \(\leftrightarrow\) sf Conversion

ppp to sf Conversion:

pois_sf <- pois_ppp |> sf::st_as_sf()
pois_sf |> head(4)

label	geom
window	POLYGON ((0 0, 1 0, 1 1, 0 …
point	POINT (0.5177981 0.553718)
point	POINT (0.6221358 0.3933781)
point	POINT (0.8446784 0.8187611)

pois_sf |> ggplot() +
  geom_sf(data=pois_sf |> filter(label=="window"), aes(fill='grey')) +
  geom_sf(data=pois_sf |> filter(label != "window"), aes(color='black')) +
  md_theme_classic(base_size=26) +
  scale_fill_manual(name=NULL, values=c("gray90"), labels=c("<span style='font-family: mono'>label == 'window'</span>")) +
  scale_color_manual(name=NULL, values=c("black"), labels=c("<span style='font-family: mono'>label == 'point'</span>")) +
  labs(title = "<span style='font-family: mono'>ppp</span> &rarr; <span style='font-family: mono'>sf</span> Result")

sf to ppp Conversion:

square_sfc <- sf::st_polygon(list(
  matrix(c(0,0,1,0,1,1,0,1,0,0), nrow=5, byrow=TRUE)
)) |> sf::st_sfc()
tree_ppp <- as.ppp(
  sf::st_as_sfc(tree_sf),
  W=as.owin(square_sfc)
)
tree_ppp

Planar point pattern: 100 points
window: polygonal boundary
enclosing rectangle: [0, 1] x [0, 1] units

tree_ppp_sf <- tree_ppp |> sf::st_as_sf()
tree_ppp_sf |> ggplot() +
  geom_sf(aes(fill='gray90')) +
  geom_sf(data=tree_ppp_sf |> filter(label != "window"), aes(color='black')) +
  md_theme_classic(base_size=26) +
  scale_fill_manual(name=NULL, values=c("gray90"), labels=c("<span style='font-family: mono'>tree_ppp$window</span>")) +
  scale_color_manual(name=NULL, values=c("black"), labels=c("<span style='font-family: mono'>tree_ppp${x,y}</span>")) +
  labs(
    title = "<span style='font-family: mono'>sf</span> &rarr; <span style='font-family: mono'>ppp</span> Result",
    x="<span style='font-family: mono'>tree_ppp$x</span>",
    y="<span style='font-family: mono'>tree_ppp$y</span>"
  ) + 
  guides(fill = guide_legend(order = 1), 
              color = guide_legend(order = 2))

References

Baddeley, Adrian, Ege Rubak, and Rolf Turner. 2015. Spatial Point Patterns: Methodology and Applications with R. CRC Press.

Morinay, Jennifer, Louise Riotte-Lambert, Geert Aarts, Federico De Pascalis, Simona Imperio, Michelangelo Morganti, Carlo Catoni, et al. 2023. “Within-Colony Segregation of Foraging Areas: From Patterns to Processes.” Oikos 2023 (8): e09926.

Schabenberger, Oliver, and Carol A. Gotway. 2004. Statistical Methods for Spatial Data Analysis. CRC Press.

Waller, Lance A., and Carol A. Gotway. 2004. Applied Spatial Statistics for Public Health Data. John Wiley & Sons.