Section 8 Appendix 1: Primer on INLA

This section introduces high-resolution mapping based on Bayesian geostatistical techniques to allow for an intuitive understanding of the models—it does not offer an in-depth derivation and explanation of the theory underpinning them. This section is largely based on the books by Zuur et al. (Zuur, Ieno, and Saveliev 2017) and Blangiardo and Cameletti (Blangiardo and Cameletti 2015).

8.1 Area Data versus Point-Referenced Data

There are three types of spatial data (Blangiardo and Cameletti 2015):

• Area or lattice data: the outcome variable to be modeled is a “random aggregate value over an areal unit”(Blangiardo and Cameletti 2015). The difference between area and lattice data is that the latter is aggregated over a regular grid while the former is aggregated over a set of irregular polygons such as administrative boundaries. For instance, income data aggregated at the segmento level are area data.
• Point-referenced (or geostatistical) data: the outcome variable to be modeled is “a random outcome at a specific location” (Blangiardo and Cameletti 2015), such as households’ income measured at the household location georeferenced with latitude and longitude coordinates.
• Spatial point patterns: the outcome variable to be modeled represents the occurrence or not of an event at a given location. While point-referenced data are random outcomes at a specific location, the location of the occurrence or not is itself random in spatial point patterns. For instance, road traffic accidents could be modeled as spatial point patterns.

The EHPM data were collected at the household level but aggregated at the segmento level, so they are area data. Nevertheless, we will also explore a model assuming they are point-referenced data, using the centroids of the segmentos as the location.

8.1.1 Bayesian Approach

In the frequentist approach, the interest lies in estimating the probability of the data, such as on income, given the parameters \(P(data | \beta)\), that is, one assumes the parameters exist. In a Bayesian approach, the interest lies in estimating the probability of the parameters given the data, \(P(\beta | data)\).

The Bayes theorem states that: \[ \begin{aligned} P(A&B)=P(A|B)P(B) \end{aligned} \]

\[ \begin{aligned} P(A|B)=\frac{P(B|A)P(A)}{P(B)} \end{aligned} \] Replacing \(A\) and \(B\) with the quantity of interest gives:

\[ \begin{aligned} P(\beta|data)=\frac{P(data|\beta)P(\beta)}{P(data)} \end{aligned} \] where \(P(\beta|data)\) is the posterior distribution of the parameter \(\beta\), \(P(data|\beta)\) is the likelihood function of the data (e.g., a gamma distribution for income), \(P(\beta)\) is the prior distribution of the \(\beta\) parameter.

This can be expressed as:

\[ \begin{aligned} P(\beta|data)\propto P(data|\beta)P(\beta) \end{aligned} \] which reads as the posterior distribution of \(\beta\) is proportional to \(P(data|\beta)P(\beta)\).

There are various strategies for getting \(P(data|\beta)P(\beta)\), such as a simulation with a Monte Carlo or a Markov Chain Monte Carlo algorithm, or a deterministic approach via Integrated Nested Laplace Approximation (INLA) implemented in the INLA package.

8.2 Model Formulation

This section models the relationship between the outcome variable \(y\), for instance, segmento median income, and a set of covariates \(\mathbf{X}\), with the bold notation denoting that \(\mathbf{X}\) is a matrix with each column corresponding to a covariate.

8.2.1 A Simple Model That Does Not Account for Spatial Dependence

Ignoring, for now, the spatial dimension, the relationship could be modeled as:

\[ \begin{aligned} \mathbf{y}=\beta_{0}+\mathbf{X}\mathbf{\beta}+\mathbf{\epsilon} \end{aligned} \]

\[ \begin{aligned} \mathbf{\epsilon}\sim N(0,\mathbf\Omega) \end{aligned} \]

where \(\mathbf{y}\) is an income vector of size \(n\) for the \(n\) segmentos. With each element \(y_{i}\) of \(i=1, ... ,n\) is the median income level in segmemto \(i\), \(\mathbf{X}\) is the matrix of covariates with \(k\) columns for the \(k\) covariates and \(n\) rows for the \(n\) segmentos, \(\mathbf{\beta}\) is the vector of \(k\) parameters indicating the effect of each covariate on \(y\), and \(\mathbf{\epsilon}\) is the vector of error term of size \(n\), distributed normally with mean \(0\), and variance-covariance defined by the matrix \(\mathbf{\Omega}\).

\[ \begin{aligned} \mathbf\Omega = \left[\begin{array} {rrrr} \sigma^2 & 0 & \dots & 0 \\\ 0 & \sigma^2 & \dots & \vdots \\\ \vdots & &\ddots & 0 \\\ 0 & \dots & 0 & \sigma^2 \end{array}\right]=\sigma^2\mathbf{I} \end{aligned} \]

Each off-diagonal element of \(\mathbf\Omega\) represents the covariance of \(\epsilon_{i}\) and \(\epsilon_{j}\), here set to \(0\). The variance is \(\sigma^2\), which can be more concisely written as \(\sigma^2\mathbf{I}\), where \(\mathbf{I}\) is an identity matrix where the diagonal elements are \(1\) and the off-diagonal elements are \(0\).

The strong assumption of this model is that the term \(\epsilon_{i}\) is independently and identically distributed, often abbreviated as iid. Let us unpack this expression.

Independence. The independence assumption implies that \(\epsilon_{i}\) and \(\epsilon_{j}\) for any two pairs of segmento \(i\) and \(j\) will be independent, even if \(i\) and \(j\) are neighboring segmentos. This is represented by the covariance between \(\epsilon_{i}\) and \(\epsilon_{j}\) for \(i\neq j\) being set to zero in \(\Omega\).

We can expect that it does not hold in the case of average income per segmento. For instance, imagine that a large company is recruiting its workforce in segmentos \(i\) and \(j\) and no data is available in \(\mathbf{X}\) to take this into account. Therefore, \(\epsilon_{i}\) and \(\epsilon_{j}\) will be highly correlated, violating the independence assumption.

Identical distribution. There is no symbol \(i\) in the normal distribution \(N(\mu,\sigma^2)\), meaning that \(\epsilon_{i}\) of all segmentos are assumed to have a mean \(\mu\) and a variance \(\sigma^2\). However, both the mean and the variance might vary across locations. For instance, one can expect that in areas where most of the workforce is employed in agriculture, weather shocks might imply a larger deviation from model predictions than in areas where most of the workforce is employed in the service sector, leading to different variances. Similarly, the model might systematically overpredict or underpredict in some areas, leading to \(\epsilon_{i}\) with \(\mu\) higher or under zero.

The violation of these assumptions leads to incorrect estimation of the precision of the \(\mathbf{\beta}\) parameters.

Furthermore, explicitly modeling the spatial dependence increases the predictive performance of the model; not taking the location into account implies throwing away an important piece of information.

8.2.2 A Spatial Model

To explicitly model the spatial dependence, an additional term is added to the equation.

\[ \begin{aligned} \mathbf{y}=\beta_{0}+\mathbf{X}\mathbf{\beta}+\mathbf{\epsilon}+\mathbf{\upsilon}, \end{aligned} \]

where \(\mathbf{\epsilon}\) is distributed according to the equation presented in the previous subsection and \(\mathbf{\upsilon}\) is a spatially correlated random effect distributed as follows:

\[ \begin{aligned} \mathbf{\upsilon}\sim N(0,\mathbf\Sigma) \end{aligned} \]

\(\Sigma\) is a non-diagonal matrix

\[ \begin{aligned} \mathbf\Sigma = \left[\begin{array} {rrrr} \sigma_{u}^2 &\phi_{1,2}^2 & \dots & \phi_{1,n}^2 \\\ \vdots & \sigma_{u}^2 & \dots & \vdots \\\ \vdots & &\ddots & \phi_{n-1,n-1}^2 \\\ \dots & \dots & \dots & \sigma_{u}^2 \end{array}\right] \end{aligned} \]

with \(\sigma_{u}^2\) as the variance of the \(\upsilon\) and \(\phi_{i,j}=corr(\upsilon_{i},\upsilon_{j})\neq0\) is the correlation between the spatial random effect at locations \(i\) and \(j\).

The challenge now is to estimate the \(\phi_{i,j}\) parameters. In the present case, the matrix \(\mathbf\Sigma\) has a dimension of \(1664*1664\), that is, 2,768,896 parameters to estimate.

8.2.3 Estimating the Model with Areal Data

When working with areal data, quantifying proximity with a Euclidean distance is not adequate. First, one would need to figure out which point of the area should be considered when measuring the distance. One option could be to take some measure of the center of each segmento; however, larger segmentos would then appear more remote than smaller ones. The preferred approach with areal data is hence slightly different: proximity is defined in terms of which segmentos share a border.

For simplicity’s sake, let us assume for the moment that the link function is a simple identity function, such that \(g(\mu(s_{i}))=\mu(s_{i})\). Therefore, \(\mu(s_{i})\) can be expressed as

\[ \begin{aligned} \mu(s_{i})=\eta_{i}=\mathbf{X(s_{i})}\mathbf{\beta}+\upsilon_{i}. \end{aligned} \]

A conditional autoregressive model (CAR) correlation function is used to model the spatial random effect \(\upsilon_{i}\) (when dealing with point-referenced data, a Matern correlation and its SDPE expression are used instead, see the next subsection).

The spatial dependency is assumed to be Markovian in nature, that is, the spatial correlation can be summarized by the spatial correlation between direct neighbors.

The distribution of each \(\upsilon_{i}\) conditional on all other \(\upsilon\) is expressed as

\[ \begin{aligned} \upsilon_{i}|\upsilon_{-i} \sim N(\sum_{j\neq i}^{N}c_{i,j}\upsilon_{j}, \sigma_{i}^{2}), \end{aligned} \]

where \({-i}\) means all except \(i\). In the simplest model, called the intrinsic CAR, \(c_{i,j}\) is \(1\) if \(j\) is a neighbor of \(i\) and \(0\) otherwise. The conditional mean of each \(\upsilon_{i}\) is hence an average of its direct neighbor.

The joint distribution of all \(\upsilon_{i}\) is also Gaussian and written as

\[ \begin{aligned} \mathbf{\upsilon} \sim N(0, (\mathbf{I-C})^{-1}\mathbf{M}), \end{aligned} \]

where \(\mathbf{I}\) is a matrix of \(1\) and \(0\), \(\mathbf{C}\) contains the \(c_{i,j}\), and \(\mathbf{M}\) is the covariance matrix. Various CAR models specify various function for \(\mathbf{C}\) and \(\mathbf{M}\).

The homogeneous CAR model set is

\[ \begin{aligned} \mathbf{C}=\phi\mathbf{A} \space \space \space \space \space \mathbf{M}=\mathbf{I}\sigma^{2}_{CAR}, \end{aligned} \]

where \(\mathbf{A}\) is \(1\) if the areas are neighboring and \(0\) otherwise, and \(\phi\) needs to be estimated: the larger it is, the more the random effects are correlated. The intrinsic CAR models set \(\phi\) to 1.

One of the issues with CAR models is that they lead to overfitting, that is, models perform well on the training data set but poorly on out-of-sample predictions on the test set for validation purposes.

A solution is to smooth the pattern of the spatial random effects. With a lower \(\sigma^{2}_{CAR}\), the \(\upsilon\) will vary less from neighbor to neighbor, leading to less overfit. To push the \(\sigma^{2}_{CAR}\) lower, it is possible to work on the priors. The standard way of doing this is by using penalized complexity priors whereby one specifies the probability that \(\sigma^{2}_{CAR}\) will be larger than a given number. To further limit the risk of overfitting, decompose the spatial random effects into one part that is spatially correlated and another that is pure noise. This is called the modified Besag-York-Mollie (Simpson and al. 2017), which we will be using in the next section.

8.2.4 Estimating the Model with Point-Referenced Data

We now assume that the segmentos’ median incomes are point-referenced data.

We reconsider \(\mathbf{y}\) as a random variable distributed at \(n\) locations \(s_{1}, ...,s_{n}\). The income data can be considered as a sample of \(y(s_{1}), ..., y(s_{n})\) from the random process \(\mathbf{y(s)}\).

Assuming that the \(y(s_{i})\) are normally distributed, we have

\[ \begin{aligned} y(s_{i})\sim N(\mu(s_{i}),\sigma^2), \end{aligned} \]

where \(\mu(s_{i})\) is expressed as a function of a structured additive predictor \(\eta_{i}\), such that \(g(\mu(s_{i}))=\eta_{i}\). For the sake of simplicity, let us assume for the moment that the link function is a simple identity function such that \(g(\mu(s_{i}))=\mu(s_{i})\). \(\mu(s_{i})\) can be expressed as

\[ \begin{aligned} \mu(s_{i})=\eta_{i}=\mathbf{X(s_{i})}\mathbf{\beta}+\upsilon_{i}+\epsilon_{i}. \end{aligned} \]

Assuming that the \(\upsilon_{i}\) is normally distributed, it makes it a Gaussian field (\(GF\)):

\[ \begin{aligned} \upsilon_{i} \sim GF(\mathbf{0},\mathbf{\Sigma}). \end{aligned} \]

Assuming that the covariance is Markovian, meaning that only the spatial correlation can be summarized by the spatial correlation between direct neighbors, then the \(\upsilon_{i}\) is distributed according to a Gaussian Markov random field (GMRF).

\[ \begin{aligned} \upsilon_{i} \sim GRMF(\mathbf{0},\mathbf{\Sigma}). \end{aligned} \]

As \(\mathbf{\Sigma}\) can be large (a matrix with more than 2 million elements in the case of the segmentos), a deterministic structure is imposed on it to speed the estimation.

When working with point-referenced data, the variance-covariance \(\mathbf{\Sigma}\) is expressed in terms of the Matern correlation function.

\[ \begin{aligned} \mathbf{\Sigma}=\sigma_{\upsilon}^2 cor_{Matern}(\upsilon(s_{i}),\upsilon(s_{j})). \end{aligned} \]

The advantage of the Matern correlation function is that only a few parameters need to be estimated.

To further simplify the computation, the Matern correlation can be re-expressed with an SPDE. Once the SPDE equation is solved, the SPDE parameters can be used to solve the Matern correlation parameters.

To solve the SPDE, a spatial mesh is created on the point-referenced data. For each node of the mesh, we get a value \(w_{k}\) via the finite element approach. The \(w_{k}\) also form a GMRF. Once we have the \(w_{k}\), we can calculate the \(\upsilon_{i}\) as a weighted sum of \(w_{k}\) time \(a_{i,k}\)s, where \(a_{i,k}\) is the distance of \(s{i}\) to node \(k\). This allows us to obtain the posterior distribution of the \(\upsilon_{i}\)s.

References

Blangiardo, M., and M. Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-Inla. John Wiley & Sons.

Simpson, D., and H. Rueand A. Riebler et al. 2017. “Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.” Statistical Science 32 (1): 1–28.

Zuur, A. F., E. N. Ieno, and A. A. Saveliev. 2017. Beginner’s Guide to Spatial, Temporal, and Spatial-Temporal Ecological Data Analysis with R-Inla. Vol. 1. Highland Statistics Ltd.