Lecture 5

Modelling areal data

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

Areal Data

Areal data are data which come from well-defined geographical units such as postcode areas, health boards, or pixels on a satellite image.
As with other types of spatial modelling, our goal is to observe and explain spatial variation in our data.
We are trying to understand and account for spatial dependence.
Generally, we aim to produce a smoothed map that summarises the spatial patterns observed in our data.

Overview of Areal processes

An areal process (or lattice process) is a stochastic process defined on a set of regions that form a partition of our region of interest \(D\).

Let \(B_1, \ldots B_m\) be our set of \(m\) distinct regions such that: \[\bigcup\limits_{i=1}^m \hspace{1mm}B_i = D.\]
Here we require that our regions are non-overlapping, with \[B_i \cap B_j = \emptyset.\]
Then our areal process is simply the stochastic process \[\{Z(B_i); i=1,\ldots,m\}.\]

Neighbourhood structures

Each of our regions \(B_i\) has a set of other regions nearby which can be considered neighbours
We might expect that areas have more in common with their neighbours than with regions further away.
Therefore, we can construct dependence structures based on the principle that neighbours are correlated and non-neighbours are uncorrelated.
First, we need to come up with a sensible way of defining what a neighbour is in this context.

Defining a Neighbourhood

There are many different ways to define a region’s neighbours.
The most common ones fall into two main categories - those based on borders, and those based on distance.

Common borders

Assume that regions which share a border on a map are neighbours.

Simple and easy to implement
Treats all borders the same, regardless of length, which can be unrealistic.
Areas very close together are not neighbours if there is even a small gap between them.

Distance-based

Assume that regions which are a within a certain distance of each other area neighbours.

What distance do you choose? How do you decide that?
Where do you measure from? (e.g., nearest border or a central point).

Neighbourhood matrix

Once we have identified a set of neighbours we construct a neighbourhood matrix (or proximity matrix), which defines how each of our \(m\) regions relate to each other.
Let \(W\) denote an \(m \times m\) matrix where the \((i,j)\)th entry, \(w_{ij}\) denotes the proximity between regions \(B_i\) and \(B_j\).

Note

The values of this matrix can be discrete (which regions are neighbours) or continuous (how far apart are the regions).

The spdep package in R can be used on an sf spatial polygon object to define \(W\) using the poly2nb() and nb2mat() functions.

Binary Neighbourhood matrix

The by far most common approach is to use a binary neighbourhood matrix, \(W\), denoted by

\[ \begin{aligned} w_{ij} &= 1 \hspace{2mm} \mbox{ if areas} (B_i, B_j) \mbox{ are neighbours.}\\ w_{ij} &= 0 \hspace{2mm} \mbox{ otherwise.} \end{aligned} \]

Dependence structures are described through this spatial weights matrix
Binary matrices are used for simplicity

How do we model this?

Imagine we have animal counts in each region. We can model them as Poisson counts

\[ y_i \sim \mathrm{Poisson}(e^{\eta_i}) \] How do we model the linear predictor \(\eta_i\)?

We could model the number of animals in each region independently

\[ \eta_i = \beta_0 + \mathbf{x}'\beta + u_i ~~~ u_i \overset{iid}{\sim} N(0,\sigma^2_i) \]

regional differences are accounted for through a random effect

How do we model this?

Imagine we have animal counts in each region. We can model them as Poisson counts

\[ y_i \sim \mathrm{Poisson}(e^{\eta_i}) \] How do we model the linear predictor \(\eta_i\)?

But… what if the distribution varies across space following a trend, i.e. is structured in space?
Do the covariates account for those structures?
If there’s an area where the animal is rare, we’ll get lots of zero counts…

How do we model this?

Imagine we have animal counts in each region. We can model them as Poisson counts.

\[ y_i \sim \mathrm{Poisson}(e^{\eta_i}) \] How do we model the linear predictor \(\eta_i\)?

We could model some dependence across regions:
- Nearby regions should have similar counts

\[ \eta_i = \beta_0 + \mathbf{x}'\beta + u_i ~~~ u_i \overset{iid}{\sim} N(0,\Sigma) \]

Now the random effect \(u_i \sim N(0, \Sigma)\) is correlated.

How do we do this?

Sparse Matrixes

\[ \mathcal{u}\sim\mathcal{N}(0,\Sigma) \]

Force the covariance matrix \(\Sigma\) to be sparse
Force the precision matrix \(\Sigma^{-1}\) to be sparse

Sparse Matrixes

\[ \mathcal{u}\sim\mathcal{N}(0,\Sigma) \]

Force the covariance matrix \(\Sigma\) to be sparse
- \(\Sigma_{ij} = \text{Cov}(u_i, u_j)\) : Covariance between \(u_i\) and \(u_j\)
- \(\Sigma_{ij} = 0\) \(\longrightarrow\) \(u_i\) and \(u_j\) are independent
- A sparse covariance matrix implies that many elements of \(\mathbf{u}\) are mutually independent…..is this desirable?
Force the precision matrix \(\mathbf{Q} = \Sigma^{-1}\) to be sparse

Sparse Matrixes

\[ \mathcal{u}\sim\mathcal{N}(0,\Sigma) \]

Force the covariance matrix \(\Sigma\) to be sparse
- \(\Sigma_{ij} = \text{Cov}(u_i, u_j)\) : Covariance between \(u_i\) and \(u_j\)
- \(\Sigma_{ij} = 0\) \(\longrightarrow\) \(u_i\) and \(u_j\) are independent
- A sparse covariance matrix implies that many elements of \(\mathbf{u}\) are mutually independent…..is this desirable?
Force the precision matrix \(\mathbf{Q} = \Sigma^{-1}\) to be sparse
- What does \(Q_{ij}\) represent?
- What does a sparse precision matrix imply?

Gaussian Markov Random Fields

Recall: The AR(1) proces

Definition

\[ \begin{aligned} \mathbf{i=1}&: x_1 \sim N\left(0, \frac{1}{1-\phi^2}\right)\\ \mathbf{i=2,\dots,T}&: x_i = \phi\ x_{i-1} +\epsilon_i,\ \epsilon_i\sim\mathcal{N}(0,1) \end{aligned} \]

Very common to model dependence in time.
The joint distribution of \(\mathbf{x}=x_1,x_2,\dots\) is Gaussian.
What do covariance and precision matrices look like?

Covariance and Precision Matrix for AR1

Covariance Matrix

\[ \Sigma = \frac{1}{1-\phi^2} \begin{bmatrix} 1& \phi & \phi^2 & \dots& \phi^N \\ \phi & 1& \phi & \dots& \phi^{N-1} \\ \phi^2 & \phi & 1 & \dots& \phi^{N-2} \\ \dots& \dots& \dots& \dots& \dots& \\ \phi^{N} & \phi^{N-1}& \phi^{N-2} & \dots& 1\\ \end{bmatrix} \]

This is a dense matrix.
All elements of the \(\mathbf{x}\) vector are dependent.

Covariance and Precision Matrix for AR1

Precision Matrix

\[ \mathbf{Q} = \Sigma^{-1} = \begin{bmatrix} 1& -\phi & 0 & 0 &\dots& 0 \\ -\phi & 1 + \phi^2& -\phi & 0 & \dots& 0 \\ 0 & -\phi & 1-\phi^2 &-\phi & \dots& 0 \\ 0 & 0 & -\phi &1-\phi^2 & \dots & \dots \\ \dots& \dots& \dots& \dots& \dots& \dots& \\ 0 &0 & 0 & \dots & -\phi& 1\\ \end{bmatrix} \]

This is a tridiagonal matrix, it is sparse.
The tridiagonal form of \(\mathbf{Q}\) can be exploited for quick calculations.

What is the key property of this example that causes \(\mathbf{Q}\) to be sparse?

Conditional independence

The key lies in the full conditionals

\[ x_t|\mathbf{x}_{-t}\sim\mathcal{N}\left(\frac{\phi}{1-\phi^2}(x_{t-1}+x_{t+1}), \frac{1}{1+\phi^2}\right) \]

The circles represent the values of \(x\) at individual time points
There is a line between them if they are conditionally dependent
Each point in time is only conditionally dependent on the two closest neighbours
The nonzero pattern in the precision matrix is given by the neighborhood structure of the process

Conditional autoregressive models

Markov in Space:

First order conditional autoregressive model or a CAR(1) model.

Every node is conditionally dependent on its four nearest neighbours

Models based on neighbourhood have a name in statistics: they are Markovian models
Markovian models are specified entirely through “neighbourhood structures”
Recall our first Model:

\[ \begin{aligned} y_i &\sim \mathrm{Poisson}(e^{\eta_i}) \\ \eta_i &= \beta_0 + \mathbf{x}'\beta + u_i ~~~ u_i \sim N(0,Q^{-1}) \end{aligned} \]
This mean \(u_i\) is independent of all the other parameters \(\mathbf{u}_{-i}\), given the set of its neighbors.
for any pair of elements (\(i,j\)) in \(\mathbf{u}\), \(u_i\perp u_j|\mathbf{u}_{-ij}\Longleftrightarrow Q{ij} =0\)
\(Q{ij} \neq 0\) only if \(j\in \{i,\mathcal{N}(i)\}\)

(Informal) definition of a GMRF

Gauss Markov random field (GMRF) is

a Gaussian distribution where the non-zero elements of the precision matrix are defined by a neighbourhood matrix (or graph structure)
each region conditionally has a Gaussian distribution with
mean equal to the average of the neighbours and
precision proportional to the number of neighbours

GMRFs are key to the many inferential approaches:

computationally efficient
other data structures can be approximated by a GMRF in a clever way

Modelling spatial similarity

An example of a GMRF is the Besag model a.k.a. the Intrinsic Conditional Autoregressive (ICAR) model. The conditional distribution for \(u_i\) is

\[ u_i|\mathbf{u}_{-i},\tau_u, \sim N\left(\frac{1}{d_i}\sum_{j\sim i}u_j,\frac{1}{d_i\tau_u}\right) \]

\(\mathbf{u}_{-i} = (u_i,\ldots,u_{i-1},u_{i+1},\ldots,u_n)^T\)
\(\tau_u\) is the precision parameter (inverse variance).
\(d_i\) is the number of neighbours
The mean of \(u_i\) is equivalent to the mean of the effects over all neighbours, and the precision is proportional to the number of neighbours (e.g., if an area has many neighbours then its variance will be smaller)

Modelling spatial similarity

The joint distribution is given by:

\[ \mathbf{u}|\tau_u \sim N\left(0,\frac{1}{\tau_u}Q^{-1}\right), \]

where \(Q\) denotes the precision matrix defined as

\[ Q_{i,j} = \begin{cases} d_i, & i = j \\ -1, & i \sim j \\ 0, &\text{otherwise} \end{cases} \]

This structure matrix directly defines the neighbourhood structure and is sparse.

What does it look like?

Larynx cancer relative risk map

What does it look like?

Connecting all the neighbouring areas yields the following graph

What does it look like?

Let us focus on one small part of the graph

What does it look like?

We apply an ICAR model where each region conditionally has a Gaussian distribution with mean equal to the average of the neighbours and a precision proportional to the number of neighbour

\[ x_9\mid\mathbf{x}_{-9}\sim N\left(\frac{1}{6}(x_7+x_{11}+x_{12}+x_{13}+x_{14}+x_{15}),\frac{1}{6\tau}\right) \]

What does it look like?

The sub graph leads to a precision matrix with 21.6% non-zero elements.

What does it look like?

The full graph leads to a precision matrix with 0.1% non-zero elements.

BYM model

The ICAR model accounts only for spatially structured variability and does not include a limiting case where no spatial structure is present.

We typically add an unstructured random effect \(z_i|\tau_z \sim N(0,\tau_{z}^{-1})\) to account for region-specific differences.
The resulting model \(v_i = u_i + z_i\) is known as the Besag-York-Mollié model (BYM)
The structured spatial effect is controlled by \(\tau_u\) which controls the degree of smoothing:
- Higher \(\tau_u\) values lead to stronger smoothing (less spatial variability).
- Lower \(\tau_u\) values allow for greater local variation.
In the following example we will illustrate how to fit this model using inlabru.

Respiratory hospitalisations in Glasgow

In this example we model the number of respiratory hospitalisations across Intermediate Zones (IZ) that make up the Greater Glasgow and Clyde health board in Scotland.

In epidemiology, disease risk is assessed using Standardized Mortality Ratios (SMR):

\[ SMR_i = \dfrac{Y_i}{E_i} \]

A value \(SMR > 1\) indicates a high risk area.
A value \(SMR<1\) suggests a low risk area.

Respiratory hospitalisation in Glasgow

Recall that a LGM consists of three elements:

Stage 1: We assume the responses are Poisson distributed: \[ \begin{aligned}y_i|\eta_i & \sim \text{Poisson}(E_i\lambda_i)\\\text{log}(\lambda_i) = \color{#FF6B6B}{\boxed{\eta_i}} & = \color{#FF6B6B}{\boxed{\beta_0 + \beta_1 \mathrm{pm10} + u_i + z_i} }\end{aligned} \]

Stage 2: \(\eta_i\) is a linear function of four components: an intercept, pm10 effect , a spatially structured effect \(u\) and an unstructured iid random effect \(z\):
\[ \eta_i = \beta_0 + \beta_1 \mathrm{pm10} + u_i + z_i \]

Stage 3: \(\{\tau_{z},\tau_u\}\): Precision parameters for the random effects

The latent field is \(\mathbf{u}= (\beta_0, \beta_1, u_1, u_2,\ldots, u_n,z_1,...)\), the hyperparameters are \(\boldsymbol{\theta} = (\tau_u,\tau_z)\), and must be given a prior.

Respiratory hospitalisation in Glasgow

The Model

\[ \begin{aligned} y_i|\eta_t & \sim \text{Poisson}(E_i\lambda_i)\\ \text{log}(\lambda_i) = \eta_i & = \color{#FF6B6B}{\boxed{\beta_0}} + \color{#FF6B6B}{\boxed{\beta_1 \mathrm{pm10}}} + \color{#FF6B6B}{\boxed{u_i}} + \color{#FF6B6B}{\boxed{z_i}} \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) +  pm10(pm10, model = "linear") +
  ui(region_id, model = "besag", graph = Q) + 
  zi(region_id, model = "iid") 

# define model predictor
eta =  observed ~ Intercept + ui + vi + pm10

# build the observation model
lik = bru_obs(formula = eta, 
              family = "poisson",
              E = expected,
              data = resp_cases)

# fit the model
fit = bru(cmp, lik)

neighbourhood structure

library(spdep)
W.nb <- poly2nb(GGHB.IZ,queen = TRUE)

Respiratory hospitalisation in Glasgow

The Model

\[ \begin{aligned} y_i|\eta_t & \sim \text{Poisson}(E_i\lambda_i)\\ \text{log}(\lambda_i) = \color{#FF6B6B}{\boxed{\eta_i}} & = \color{#FF6B6B}{\boxed{\beta_0 + \beta_1 \mathrm{pm10} + u_i + z_i}} \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + pm10(pm10, model = "linear") 
  ui(region_id, model = "besag", graph = Q) + 
  vi(region_id, model = "iid") 

# define model predictor
eta =  observed ~ Intercept + ui + vi + pm10

# build the observation model
lik = bru_obs(formula = eta, 
              family = "poisson",
              E = expected,
              data = resp_cases)

# fit the model
fit = bru(cmp, lik)

neighbourhood structure

library(spdep)
W.nb <- poly2nb(GGHB.IZ,queen = TRUE)

Respiratory hospitalisation in Glasgow

The Model

\[ \begin{aligned} \color{#FF6B6B}{\boxed{y_i|\eta_t}} & \sim \color{#FF6B6B}{\boxed{\text{Poisson}(E_i\lambda_i)}}\\ \text{log}(\lambda_i) = \eta_i & = \beta_0 + \beta_1 \mathrm{pm10} + u_i + z_i \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + pm10(pm10, model = "linear")
  ui(region_id, model = "besag", graph = Q) + 
  vi(region_id, model = "iid") 

# define model predictor
eta =  observed ~ Intercept + ui + vi + pm10

# build the observation model
lik = bru_obs(formula = eta, 
              family = "poisson",
              E = expected,
              data = resp_cases)

# fit the model
fit = bru(cmp, lik)

neighbourhood structure

library(spdep)
W.nb <- poly2nb(GGHB.IZ,queen = TRUE)

Respiratory hospitalisation in Glasgow

The Model

\[ \begin{aligned} y_i|\eta_t & \sim \text{Poisson}(E_i\lambda_i)\\ \text{log}(\lambda_i) = \eta_i & = \beta_0 + \beta_1 \mathrm{pm10} + u_i + z_i \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + pm10(pm10, model = "linear")
  ui(region_id, model = "besag", graph = Q) + 
  vi(region_id, model = "iid") 

# define model predictor
eta =  observed ~ Intercept + ui + vi + pm10

# build the observation model
lik = bru_obs(formula = eta,  
              family = "poisson",
              E = expected,
              data = resp_cases)

# fit the model
fit = bru(cmp, lik)

neighbourhood structure

library(spdep)
W.nb <- poly2nb(GGHB.IZ,queen = TRUE)

Respiratory hospitalisation in Glasgow

Posterior summaries

	mean	0.025quant	0.975quant
Intercept	−1.52	−1.91	−1.12
pm10	0.09	0.07	0.12
Precision for ui	2.82	2.28	3.45
Precision for vi	21,828.33	1,415.98	85,763.09

Estimated relative risks

pred = predict(fit, resp_cases,
 ~data.frame(
 log_risk = Intercept + ui + vi + pm10,
 risk = exp(Intercept + ui + vi+  pm10)),
      n.samples = 1000)

Improved interpretability: BYM2

BYM2 Model

In the original BYM, the spatially structured component must be scaled so that \(\tau_u\) produces consistent smoothness across different neighborhood structures.

Riebler et al. (2016) proposed a new parametrization of the BYM model based on PC-priors that improves parameter interpretability.

\[ \mathbf{b} = \dfrac{1}{\sqrt{\tau_b}} \left(\sqrt{1-\phi}v^*+\sqrt{\phi}u^*\right). \]

The precision \(\tau_b>0\) controls the marginal variance contribution of the weighted sum \(u^*\) and \(v^*\).
The mixing parameter \(0 \leq \phi \leq 1\) measures the proportion of the marginal variance explained by the structured effect \(u^*\)
- if \(\phi =1\) the model captures only spatially structured variability
- if \(\phi = 0\) it accounts solely for unstructured spatial noise.

BYM2 & PC priors

The BYM2 model allows the specification of Penalized Complexity (PC) priors to control the amount of spatial smoothing and to avoid overfitting.

PC priors shrink the model towards a simpler baseline unless the data provide strong evidence for a more complex structure.
To define the prior for the marginal precision \(\tau_b\) and the mixing parameter \(\phi\), we use the probability statements:

\[ \begin{aligned} P(1/\sqrt{\tau_b} > U) &= \alpha\\ P(\phi < U) &= \alpha \end{aligned} \]

\[ \begin{aligned} P((1/\sqrt{\tau_b}) > 0.5/0.31) &\equiv P(\sigma > 1.61) = 0.01\\ P(\phi < 0.5) &= 2/3 \approx 0.66 \end{aligned} \]

Interpretation:

Prior on \(\tau_b\): low probability of large values
Prior on \(\phi\) assumes that the unstructured random effect accounts for more of the variability than the spatially structured effect.

BYM2 & PC priors

The BYM2 model allows the specification of Penalized Complexity (PC) priors to control the amount of spatial smoothing and avoid overfitting.

PC priors shrink the model towards a simpler baseline unless the data provide strong evidence for a more complex structure.
To define the prior for the marginal precision \(\tau_b\) and the mixing parameter \(\phi\), we use the probability statements:

# define model component
cmp = ~ Intercept(1) +  pm10(pm10, model = "linear") +
  b(region_id, model = "bym2", graph = Q,
               # priors
               hyper = list(theta1 = list("PCprior", c(1.61, 0.01)), # Pr(sd>1) = 0.01 
                            theta2 = list("PCprior", c(0.5, 0.66)))) # Pr(phi<0.5)=0.66

Interpretation:

Prior on \(\tau_b\): low probability of large values
Prior on \(\phi\): assumes that the unstructured random effect accounts for more of the variability than the spatially structured effect.

Take-Home Message

areal data are data on discretised space
spatial relationships are represented through a neighbourhood structure
GMRFs are a computationally efficient representation of spatial structures
these will also be relevant – and very useful – when it comes to approximating spatial structures in continuous space

References

Goicoa, Tomás, Aritz Adin, Marı́a Dolores Ugarte, and James S Hodges. 2018. “In Spatio-Temporal Disease Mapping Models, Identifiability Constraints Affect PQL and INLA Results.” Stochastic Environmental Research and Risk Assessment 32 (3): 749–70.

Lee, Sophie A, Theodoros Economou, and Rachel Lowe. 2022. “A Bayesian Modelling Framework to Quantify Multiple Sources of Spatial Variation for Disease Mapping.” Journal Of The Royal Society Interface 19 (194): 20220440.

Orozco-Acosta, Erick, Aritz Adin, and Marı́a Dolores Ugarte. 2021. “Scalable Bayesian Modelling for Smoothing Disease Risks in Large Spatial Data Sets Using INLA.” Spatial Statistics 41: 100496.

Riebler, Andrea, Sigrunn H. Sørbye, Daniel Simpson, and Håvard Rue. 2016. “An Intuitive Bayesian Spatial Model for Disease Mapping That Accounts for Scaling.” Statistical Methods in Medical Research 25 (4): 1145–65. https://doi.org/10.1177/0962280216660421.

Seaton, Fiona M, Susan G Jarvis, and Peter A Henrys. 2024. “Spatio-Temporal Data Integration for Species Distribution Modelling in r-INLA.” Methods in Ecology and Evolution 15 (7): 1221–32.

Simpson, Daniel, Håvard Rue, Andrea Riebler, Thiago G. Martins, and Sigrunn H. Sørbye. 2017. “Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.” Statistical Science 32 (1): 1–28. https://doi.org/10.1214/16-STS576.