Week 9: Evaluating Spatial Hypotheses I: Clustering, Repulsion, and Randomness

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

Jeff Jacobs

jj1088@georgetown.edu

Wednesday, October 21, 2026

Roadmap to the Midterm!

• This Week (Oct 21): Evaluating Hypotheses for Point Data
• Next Week (Oct 28): Evaluating Hypotheses for Areal Data
• In-Class Midterm (Nov 4): Basically a “mini-homework” on sf → ppp → spatial hypothesis testing!

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

So… What Do We Do With Moran’s I?

	Doctor	Midterm Practice	Midterm
	Your temperature is 40° C	Moran’s I of Protestant Churches is 0.58	Moran’s I of waste dump sites is 0.67
	One cause could be malaria, 40% likelihood	One cause could be urbanization, 40% likelihood	One cause could be density of stigmatized group, 60% likelihood
	One cause could be lupus, 75% likelihood	One cause could be distance from Luther, 75% likelihood	One cause could be distance from UNOCI HQ, 75% likelihood
	I judge lupus as more likely cause, we should treat it	I judge distance from Luther as more likely, we should study Luther more	I judge distance from UNOCI as more likely; handle in UN not ICC

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

A pattern arising from CSR!
(Realizations \omega_{\tiny\triangle} and \omega_{\tiny\square} from same D)
From Parable of the Polygons

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

Namibian “Fairy Circles”; Smithsonian Magazine: Namibia’s mysterious fairy circles may be competition 👀 between plants

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)

Spatial clustering of birdsongs; from Crates et al. (2021)

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 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: Modeling 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() +
  transparent_bg()
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() +
  transparent_bg()
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() +
  transparent_bg()
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() +
  transparent_bg()
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() +
  transparent_bg()
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() +
  transparent_bg()
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 in 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") +
  transparent_bg()
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") +
  transparent_bg()
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") +
  transparent_bg()
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!

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) +
  transparent_bg()

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)
  ) +
    transparent_bg()

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) +
  transparent_bg()

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) +
  transparent_bg()

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")) +
  transparent_bg()

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

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) +
  transparent_bg()
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(1.5,0,0,2), las=2, oma=c(0,0,0,0), cex=2, bg=NA)
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) +
  transparent_bg()
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) +
  transparent_bg()

par(
  mar=c(0,0,1,2),
  # mai=c(0,0,0.1,0),
  las=2,
  oma=c(0,0,1,0)
)
sq_density <- density(sq_ppp)
plot(sq_density, main=NULL, xaxs='i', yaxs='i', bg=NA)
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) +
  transparent_bg()

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) +
  transparent_bg()

par(
  mar=c(0,0,1,2),
  las=2,
  oma=c(0,0,1,0)
)
csr_density <- density(csr_ppp)
plot(csr_density, main=NULL, xaxs="i", yaxs="i", bg=NA)
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) +
  transparent_bg()

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 <- 20
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) +
  transparent_bg()

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) +
  transparent_bg()

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) +
  transparent_bg()
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) +
  transparent_bg()
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)
# 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

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)

library(mdthemes) |> suppressPackageStartupMessages()
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") +
  transparent_bg()

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)
  ) +
  transparent_bg()

References

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

Crates, Ross, Naomi Langmore, Louis Ranjard, Dejan Stojanovic, Laura Rayner, Dean Ingwersen, and Robert Heinsohn. 2021. “Loss of Vocal Culture and Fitness Costs in a Critically Endangered Songbird.” Proceedings of the Royal Society B: Biological Sciences 288 (1947): 20210225.

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.