Lecture 6

Spatial point processes

Sara Martino
sara.martino@ntnu.no

Dept. of Mathematical Science, NTNU

Janine Illian
Janine.Illian@glasgow.ac.uk

University of Glasgow

Jafet Belmont
Jafet.Belmont@glasgow.ac.uk

University of Glasgow

Point process data

Many of the ecological and environmental processes of interest can be represented by a spatial point process or can be viewed as an aggregation of one.

  • Many contemporary data sources collect georeferenced information about the location where an event has occur (e.g., species occurrence, wildfire, flood events).
  • This point-based information provides valuable insights into ecosystem dynamics.

Point patterns vs. geostatistical data

Point patterns:

  • Data format: x,y coordinates
  • Optional: marks
  • Aim: model locations as random

Geostatistical data:

  • Data format: x,y coordinates
  • Measurements mandatory
  • Aim: model continuous process at fixed locations

models of spatial patterns: modelling locations and properties (“marks”) of objects, events, individuals in space and time

aim: understanding mechanisms that generated the pattern

Defining a Point Process

  • Consider a fixed geographical region \(A\).

  • The set of locations at which events occur are denoted by \(\mathbf{s} = (\mathbf{s}_1, \ldots, \mathbf{s}_n)\).

  • We let \(N(A)\) be a random variable which represents the total number of events in every subset of region \(A\).

  • Our primary interest is in measuring where events occur, so the locations are our data.

  • Let \(\lambda\) the intensity or point density, i.e, the mean number of points per unit area

  • In some cases, the intensity will be constant over space (homogeneous), while in other cases it can vary by location (inhomogeneous).

  • If our intensity is homogeneous, we can define it as

    \[ \lambda(s) = \frac{\mathbb{E}[N(A)]}{|A|} = \frac{\lambda |A|}{|A|} = \lambda. \]

Complete spatial randomness

We can use the concept of intensity to help us define complete spatial randomness (CSR).

  • For any spatial region \(A\), CSR requires that:

  1. Uniformity and Independent scattering : Given the number of events \(N(A) = n\) in a region, the \(n\) events are independently and have an equal probability of occurring anywhere in the study area.

  2. Poisson distribution of point counts: The number of points in any set \(A_i\) follows a Poisson distribution with mean \(\lambda|A_i|\), that is \[N(A_i) \sim \text{Poisson}(\lambda \,|A_i|).\]

  • If these conditions are satisfied, we can describe our process as a homogeneous Poisson process.

Homogeneous Poisson Process

  • A simplest point process model is the homogeneous Poisson process (HPP).

  • The likelihood of a point pattern \(\mathbf{y} = \left[ \mathbf{s}_1, \ldots, \mathbf{s}_n \right]^\intercal\) distributed as a HPP with intensity \(\lambda\) and observation window \(\Omega\) is

    \[ p(\mathbf{y} | \lambda) \propto \lambda^n e^{ \left( - |\Omega| \lambda \right)} , \]

    • \(|\Omega|\) is the size of the observation window.

    • \(\lambda\) is the expected number of points per unit area.

    • \(|\Omega|\lambda\) the total expected number of points in the observation window.

  • A key property of a Poisson process is that the number of points within any subset \(A_i\) of region \(A\) is Poisson distributed with constant rate \(|A_i|\lambda\).

Inhomogeneous Poisson process

Let \(\mathbf{y} = s_1,\ldots,s_n\) the \(n\) number of observed events/points in an observation window \(\Omega\)

For an IPP with an intensity \(\lambda(s)\), the likelihood is given by:

\[ p(\mathbf{y} | \lambda) \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i). \]

  • If the case of an HPP the integral in the likelihood can easily be computed as \(\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} =|\Omega|\lambda\)

  • For an HPP with an intensity \(\lambda\), the log-likelihood is given by: \[ l(\beta;\mathbf{y}) = n\log(\lambda) -\lambda|\Omega|, \]

  • The maximum likelihood estimators is \(\hat{\lambda} = n/|\Omega|\).

  • For IPP, the integral in the likelihood has to be approximateda as a weighted sum.

Inhomogeneous Poisson process

This integral is approximated as \(\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \approx \sum_{j=1}^J w_j \lambda(\mathbf{s}_j)\)

  • \(w_j\) are the integration weights

  • \(\mathbf{s}_j\) are the quadrature locations.

This serves two purposes:

  1. Approximating the integral

  2. re-writing the inhomogeneous Poisson process likelihood as a regular Poisson likelihood.

Inhomogeneous Poisson process

The idea behind the trick to rewrite the approximate likelihood is to introduce a dummy vector \(\mathbf{z}\) and an integration weights vector \(\mathbf{w}\) of length \(J + n\)

\[\mathbf{z} = \left[\underbrace{0_1, \ldots,0_J}_\text{quadrature locations}, \underbrace{1_1, \ldots ,1_n}_{\text{data points}} \right]^\intercal\]

\[\mathbf{w} = \left[ \underbrace{w_1, \ldots, w_J}_\text{quadrature locations}, \underbrace{0_1, \ldots, 0_n}_\text{data points} \right]^\intercal\]

Then the approximate likelihood can be written as

\[ \begin{aligned} p(\mathbf{z} | \lambda) &\propto \prod_{i=1}^{J + n} \eta_i^{z_i} \exp\left(-w_i \eta_i \right) \\ \eta_i &= \log\lambda(\mathbf{s}_i) = \mathbf{x}(s)'\beta \end{aligned} \]

  • This is similar to a product of Poisson distributions with means \(\eta_i\), exposures \(w_i\) and observations \(z_i\)

Limitations with IPP

  • IPP models assume that data points are conditionally independent given the covariates, meaning that any spatial variation is fully explained by environmental and sampling factors.
  • Unmeasured endogenous and exogenous factors can create spatial dependence.
  • Ignoring them can lead to bias in our conclusions.

The Log-Gaussian Cox Process

  • Log-Gaussian Cox processes (LGCP) extend the IPP by allowing the intensity function to vary spatially according to a structured spatial random effect, i.e. the intensity is random

\[ \log~\lambda(s)= \mathbf{x}(s)'\beta + \xi(s) \]

  • The events are then assumed to be independent given the covariates and \(\xi(s)\) - a GRF with Matérn covariance.

How do we model \(\xi(s)\) ?

  • We use an SPDE model!
  • The software inlabru has implemented some integration schemes that are especially well suited to integrating the intensity in models with an SPDE effect.
  • By default, inlabru uses the same points to define the SPDE approximation and to approximate the integral in the likelihood, but this can be changed.

Example: Forest fires in Castilla-La Mancha

  • In this example we model the location of forest fires in the Castilla-La Mancha region of Spain between 1998 and 2007.

  • We are now going to use the elevation as a covariate to explain the variability of the intensity \(\lambda(s)\) over the domain of interest and a spatially structured SPDE model.

\[ \log\lambda(s) = \beta_0 + \beta_1 \text{elevation}(s) + \xi(s) \]

pp  %>% print(n = 6)
Simple feature collection with 148 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 34.96293 ymin: 34.87499 xmax: 378.819 ymax: 357.528
CRS:           NA
First 6 features:
      cause burnt.area       date julian.date                  geometry
1 lightning       0.01 2004-01-08        2198 POINT (165.5546 74.87501)
2 lightning       0.04 2004-01-08        2198  POINT (215.7145 254.875)
3 lightning       0.80 2004-01-08        2198   POINT (263.326 284.166)
4 lightning       0.00 2004-01-08        2198  POINT (265.2747 284.875)
5 lightning       0.25 2004-01-08        2198   POINT (236.799 270.223)
6 lightning       0.01 2004-01-09        2199   POINT (297.892 211.545)

Example: Forest fires in Castilla-La Mancha

The IPP Model \[ \begin{aligned} p(\mathbf{y} | \lambda) & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) \\ \eta(s) & = \log ~\lambda(s) = \beta_0 + \beta_1 \,x(s) \end{aligned} \] The code

# define model component
cmp = ~ Intercept(1) + elev(elev_raster, model = "linear")

# define model predictor
eta  = geometry ~ Intercept +  elev

# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = pp,
              ips = ips)

# fit the model
fit = bru(cmp, lik)
n.int = 1000
ips <- st_sample(region, size = n.int, type = "regular") # May not be exactly n.int points
ips <- new_fm_int(
  ips,
  weight = st_area(region) / length(ips),
  name = "geometry"
)

# The mesh
mesh = fm_mesh_2d(boundary = region,
                  max.edge = c(5, 10),
                  cutoff = 4, crs = NA)

# build integration scheme
ips = fm_int(mesh,
             samplers = region)

Example: Forest fires in Castilla-La Mancha

The LGCP Model

\[ \begin{aligned} p(\mathbf{y} | \lambda) & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) \\ \eta(s) & = \log(\lambda(s)) = \color{#FF6B6B}{\boxed{\beta_0}} + \color{#FF6B6B}{\boxed{\beta_1 x(s)}} + \color{#FF6B6B}{\boxed{ \omega(s)}}\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + elev(elev_raster, model = "linear") +
  space(geometry, model = spde_model)


# define model predictor
eta  = geometry ~ Intercept +  elev + space


# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = pp,
              ips = ips)


# fit the model
fit = bru(cmp, lik)
# The mesh
mesh = fm_mesh_2d(boundary = region,
                  max.edge = c(5, 10),
                  cutoff = 4, crs = NA)

# The SPDE model
spde_model =  inla.spde2.pcmatern(mesh,
                                  prior.sigma = c(1, 0.5),
                                  prior.range = c(100, 0.5))

# build integration scheme
ips = fm_int(mesh,
             samplers = region)

Example: Forest fires in Castilla-La Mancha

The LGCP Model

\[ \begin{aligned} p(\mathbf{y} | \lambda) & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) \\ \color{#FF6B6B}{\boxed{\eta(s)}} & = \log(\lambda(s)) = \color{#FF6B6B}{\boxed{\beta_0 + \beta_1 x(s) + \omega(s)}}\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + elev(elev_raster, model = "linear") +
  space(geometry, model = spde_model)


# define model predictor
eta  = geometry ~ Intercept +  elev + space


# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = pp,
              ips = ips)


# fit the model
fit = bru(cmp, lik)
# The mesh
mesh = fm_mesh_2d(boundary = region,
                  max.edge = c(5, 10),
                  cutoff = 4, crs = NA)

# The SPDE model
spde_model =  inla.spde2.pcmatern(mesh,
                                  prior.sigma = c(1, 0.5),
                                  prior.range = c(100, 0.5))

# build integration scheme
ips = fm_int(mesh,
             samplers = region)

Example: Forest fires in Castilla-La Mancha

The LGCP Model

\[ \begin{aligned} \color{#FF6B6B}{\boxed{p(\mathbf{y} | \lambda)}} & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) \\ \eta(s) & = \log(\lambda(s)) = \beta_0 + \beta_1x(s) + \omega(s)\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + elev(elev_raster, model = "linear") +
  space(geometry, model = spde_model)


# define model predictor
eta  = geometry ~ Intercept +  elev + space


# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = pp,
              ips = ips)


# fit the model
fit = bru(cmp, lik)
# The mesh
mesh = fm_mesh_2d(boundary = region,
                  max.edge = c(5, 10),
                  cutoff = 4, crs = NA)

# The SPDE model
spde_model =  inla.spde2.pcmatern(mesh,
                                  prior.sigma = c(1, 0.5),
                                  prior.range = c(100, 0.5))

# build integration scheme
ips = fm_int(mesh,
             samplers = region)

Example: Forest fires in Castilla-La Mancha

The LGCP Model

\[ \begin{aligned} p(\mathbf{y} | \lambda) & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) \\ \eta(s) & = \log(\lambda(s)) = \beta_0 + \beta_1 x(s) + \omega(s)\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + elev(elev_raster, model = "linear") +
  space(geometry, model = spde_model)


# define model predictor
eta  = geometry ~ Intercept +  elev + space


# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = pp,
              ips = ips)


# fit the model
fit = bru(cmp, lik)
# The mesh
mesh = fm_mesh_2d(boundary = region,
                  max.edge = c(5, 10),
                  cutoff = 4, crs = NA)

# The SPDE model
spde_model =  inla.spde2.pcmatern(mesh,
                                  prior.sigma = c(1, 0.5),
                                  prior.range = c(100, 0.5))

# build integration scheme
ips = fm_int(mesh,
             samplers = region)

Example: Forest fires in Castilla-La Mancha

Model Predictions

We can use the predict() method to visualize the posterior mean intensity, log-intensity, std dev, and the SPDE.

Note

Note that the mean field is, by definition, smoother than any realization of the field. Instead we can draw simulations of the field from the posterior distribution using the generate() function.

Posterior quantities

We can estimate the total number of events over the whole region as follows:

\[ \mathbb{E}(N_\Omega) = \int_\Omega \exp(\lambda(s)) ds \]

To do so, we simulate possible realizations of \(N_\omega\) to include also the likelihood variability in our estimate using the generate() method

N_fires = generate(fit, ips,
                      formula = ~ {
                        lambda = sum(weight * exp(elev + Intercept + space))
                        rpois(1, lambda)},
                    n.samples = 2000)

Take-home messages

  • Point processes are a stochastic processes that describe the locations where events occur

  • Unlike geostatistical data where the locations are fixed, here the locations have a stochastic nature; the locations are our data!

  • IPPs allow the intensity of the point process to vary across space through spatially varying covariates.

  • Numerical integration schemes are required to estimate the parameters of an IPP

  • LGCP are a doubly stochastic process that extends IPP models by allowing the intensity function to vary spatially according to a structured spatial random effect

References

Rivera et al. (2025) Laxton et al. (2023) Palmı́-Perales et al. (2025) Simpson et al. (2016)

Laxton, Megan R, Óscar Rodrı́guez de Rivera, Andrea Soriano-Redondo, and Janine B Illian. 2023. “Balancing Structural Complexity with Ecological Insight in Spatio-Temporal Species Distribution Models.” Methods in Ecology and Evolution 14 (1): 162–72.
Palmı́-Perales, Francisco, Finn Lindgren, and Virgilio Gómez-Rubio. 2025. “Bayesian Inference for Case-Control Point Pattern Data in Spatial Epidemiology with the Inlabru r Package.” arXiv Preprint arXiv:2503.14954.
Rivera, Óscar Rodrı́guez de, Juncal Espinosa, Javier Madrigal, Marta Blangiardo, and Antonio López-Quı́lez. 2025. “Spatio-Temporal Marked Point Process Model to Understand Forest Fires in the Mediterranean Basin: Ór d. Rivera Et Al.” Journal of Agricultural, Biological and Environmental Statistics 30 (3): 700–729.
Simpson, Daniel, Janine Baerbel Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.