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

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

Class Sessions
Author
Affiliation

Jeff Jacobs

Published

Wednesday, October 15, 2025

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

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

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!

Spatial Randomness

Code
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)

Code
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!

Code
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?

Code
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
Code (Step 1: Grid Creation)
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)

Code (Step 2: Point Generation)
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)
  )

Code (Step 3: Observed Result)
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
Code (Step 1: Tree Generation)
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)

Code (Step 2: Competition)
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"))

Code (Step 3: Observed Result)
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?

Code
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”?

Code (Fixed Points)
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)

Code
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)

Code
### 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)

Code (CSR Points)
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)

Code
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)

Code
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})\)
Code
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
Code
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)
  • 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)
Code
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
Code
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")
Code
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]
Code
# 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!

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:

Code
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:

Code
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
Code
attributes(pois_ppp)$names
[1] "window"     "n"          "x"          "y"          "markformat"
Code
pois_ppp$x |> head(4)
[1] 0.51779807 0.62213583 0.84467836 0.05354204

ppp \(\leftrightarrow\) sf Conversion

ppp to sf Conversion:

Code
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)
Code
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:

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