Intro to Spatial Modeling Using INLA
Zhengxiao Wei, Farouk S. Nathoo
2024-05-02
============================================================================================================ This up-to-date document is available at https://rpubs.com/sherloconan/1164312
Large Sample Size: https://rpubs.com/sherloconan/1198601
Literature Review: https://github.com/zhengxiaoUVic/Spatial
\(\{y(\mathbf{s}_i):\ \mathbf{s}_i\in\mathcal{D}\subseteq\mathbb{R}^d\}\sim\mathcal{GP}(\mu(\mathbf{s}_i),\Sigma(\mathbf{s}_i,\mathbf{s}_j))\) and \(\boldsymbol{\Sigma}={\color{red}{\sigma_s^2}}\!\cdot\!\mathbf{K}+{\color{blue}{\sigma_{\epsilon}^2}}\!\cdot\!\mathbf{I}\), where \(\mathbf{K}\) is a Matérn kernel.
\[\begin{equation*} \tag{1} K_{ij}=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}\delta_{ij}}{l}\right)^\nu\!\cdot K_\nu\left(\frac{\sqrt{2\nu}\delta_{ij}}{l}\right), \end{equation*}\] where \(K_\nu(\cdot)\) is the modified Bessel function of the second kind, \(l\) is the length scale (bandwidth) hyperparameter, and \(\delta_{ij}=\lVert\mathbf{s}_i-\mathbf{s}_j\rVert\).
Gaussian Process (GP) Regression Model
\[\begin{equation*} \tag{2}
\mathbf{y}(\mathbf{s}):=\left(y(\mathbf{s}_1),\dotsb,y(\mathbf{s}_N)\right)^\top\!=\boldsymbol{\mu}(\mathbf{s})+\mathbf{b}(\mathbf{s})+\boldsymbol{\epsilon}(\mathbf{s}),
\end{equation*}\] where \(\mathbf{b}(\mathbf{s})\) is a latent
spatial process capturing spatial dependence. See f() in
the INLA formula.
\[\begin{equation*} \tag{3} \mathbf{b}(\mathbf{s})\sim\boldsymbol{\mathcal{N}}_{\!N}(\mathbf{0},\ {\color{red}{\sigma_s^2}}\!\cdot\!\mathbf{K})\ \perp \!\!\! \perp\ \boldsymbol{\epsilon}(\mathbf{s})\sim\boldsymbol{\mathcal{N}}_{\!N}(\mathbf{0},\ {\color{blue}{\sigma_{\epsilon}^2}}\!\cdot\!\mathbf{I}) \end{equation*}\]
“LGMs (latent Gaussian models): a way to compute, not a way to model” — Håvard Rue
Latent GP \[\begin{equation*} \tag{4} \boldsymbol{\mathcal{N}}_{\!N}\left(\mathbf{y}(\mathbf{s})\mid\boldsymbol{\mu}(\mathbf{s})+\mathbf{b}(\mathbf{s}),\ {\color{blue}{\sigma_{\epsilon}^2}}\!\cdot\!\mathbf{I}\right)\ \times\ \boldsymbol{\mathcal{N}}_{\!N}\left(\mathbf{b}(\mathbf{s})\mid\mathbf{0},\ {\color{red}{\sigma_s^2}}\!\cdot\!\mathbf{K}\right)\ \times\ \pi({\color{red}{\sigma_s^2}},\ l,\ {\color{blue}{\sigma_{\epsilon}^2}}) \end{equation*}\]
Response GP \[\begin{equation*} \tag{5} \boldsymbol{\mathcal{N}}_{\!N}\left(\mathbf{y}(\mathbf{s})\mid\boldsymbol{\mu}(\mathbf{s}),\ {\color{red}{\sigma_s^2}}\!\cdot\!\mathbf{K}+{\color{blue}{\sigma_{\epsilon}^2}}\!\cdot\!\mathbf{I}\right)\ \times\ \pi({\color{red}{\sigma_s^2}},\ l,\ {\color{blue}{\sigma_{\epsilon}^2}}) \end{equation*}\]
The “R-INLA” package constructs a Matérn stochastic partial
differential equation (SPDE) model with the spatial scale \(\kappa=\frac{\sqrt{8\nu}}{\rho}\), where
the spatial range is \(\rho=2 l\quad\)
(this notation should not be confused with the “correlation
coefficient”).
In fact, the corresponding correlation is around
0.14 at the distance \(\rho\).
The joint penalized complexity (PC,
but see ?inla.spde2.pcmatern) hyperprior is \[\begin{equation*} \tag{6}
\pi({\color{red}{\sigma_s}},\
\rho)=\frac{1}{2}d\lambda_\rho\rho^{-1-d/2}\exp\left\{-\lambda_\rho\rho^{-d/2}\right\}\cdot\lambda_{{\color{red}{\sigma_s}}}\exp\left\{-\lambda_{{\color{red}{\sigma_s}}}{\color{red}{\sigma_s}}\right\},
\end{equation*}\] specified to penalize high spatial variances
and large spatial ranges from the base model, thus preventing
over-fitting and over-smoothing.
\(\mathbb{P}({\color{red}{\sigma_s}}>\sigma_0)=p_{{\color{red}{\sigma_s}}}\), e.g., the probability of \({\color{red}{\sigma_s}}\) being greater than 1 is equal to 1% if \(\mathbb{P}({\color{red}{\sigma_s}}>1)=0.01\).
\(\mathbb{P}(\rho<\rho_0)=p_{\rho}\)
Objective: Use the “R-INLA” package to model a Gaussian Markov random field (GMRF) on a square lattice.
- INLA version
24.05.1025.06.07 (stable) - fmesher version
0.1.60.5.0
Synthetic Data
Simulate some spatial data on an \(n\times n\) lattice, assuming a standard normal distribution \(\mathcal{N}(0,1)\) for the response values.
n <- 10
coords <- expand.grid(x=1:n, y=1:n) #10×10 lattice
coords$index2 <- coords$index <- 1:nrow(coords) #indices 1–100
coords[order(coords$x, -coords$y), "index2"] <- 1:nrow(coords) #data arrangement
set.seed(277)
data <- cbind(coords, "value"=rnorm(n*n))
# Plot the simulated data using lattice::levelplot
levelplot(value ~ x*y, data=data, main="Simulated Spatial Data", col.regions=heat.colors(100))
Alternative setup.
# sp::SpatialPixelsDataFrame class
set.seed(277)
spatial_data <- SpatialPixelsDataFrame(points=expand.grid(x=1:n, y=1:n),
data=data.frame("value"=rnorm(n*n),
"index"=1:(n*n),
"index2"=c(t(apply(matrix(1:(n*n),nrow=n),2,rev)))))
# Plot the simulated data using sp::spplot
spplot(spatial_data, "value", xlab="x", ylab="y", main="Simulated Spatial Data")
Data Arrangement
For some models, INLA considers data sorted by column, i.e., a vector with the first column of the grid from top to bottom, followed by the second column and so on (Gómez-Rubio, 2021, ch. 7).
index: lower left corner by rows (internally stored in expand.grid)index2: upper left corner by columns (mapping required by INLA)
plot(0:1,0:1,type='n',bty='n',xlab='',ylab='',xaxt='n',yaxt='n',
main="Location Index",sub="Black (index): lower left by rows\nRed (index2): upper left by columns")
for(i in 1:n) {
for(j in 1:n) {
rect((j-1)/n, 1-i/n, j/n, 1-(i-1)/n)
text((j-0.75)/n, (i-0.25)/n, subset(coords, x==j & y==i)$index)
text((j-0.25)/n, (i-0.75)/n, subset(coords, x==j & y==i)$index2, col="red")
}
}
Read more about aggregated data: https://rpubs.com/sherloconan/1151407.
Modeling
Fit a GMRF model via SPDE and conduct approximate Bayesian inference using the integrated nested Laplace approximation (INLA) method. \[\begin{equation*} \tag{6} f(\mathbf{s})=\sum_{q=1}^{n} w_q\cdot \psi_q(\mathbf{s}),\quad\text{given }(w_1,\dotsb,w_n)^\top\sim\boldsymbol{\mathcal{N}}_{\!n}(\mathbf{0},\ \tilde{\mathbf{\Omega}}^{-1}_w) \end{equation*}\]
Eq. 6 incorporates the piecewise linear basis functions \(\psi_q(\mathbf{s})\) and Gaussian weights \(w_q\) with Markov dependencies determined by a general triangulation of the bounded domain \(\mathcal{D}\) (finite element methods). \(n\) is the number of vertices.
The model formula, e.g., value ~ 1 + covariates + f(),
typically includes the response variable, fixed effects, and random
effects and other latent models specified using the f()
function.
IID: Independent and identically distributed. Let \(\mathbf{K}=\mathbf{I}\) in Eq. 3.
Besag: Used for spatial data that are defined on a lattice. This model is part of the family of intrinsic conditional autoregressive (ICAR) models. Read more.
SPDE: Used for modeling spatial fields that are continuous over an area, particularly useful for unstructured meshes and non-lattice data. Increasing
alphaimplies that the resulting Gaussian random field (GRF) will have more differentiable sample paths, which corresponds to greater smoothness.
We use INLA::inla.mesh.2d to define the adjacency structure
of the lattice, which INLA uses to construct the GMRF. The
max.edge parameter controls the mesh’s granularity (i.e.,
the largest allowed triangle edge length). The unit of measurement for
the edge length in max.edge depends on the units of the
coordinates provided in the loc parameter.
A lower
max.edge=c(inside the boundary triangle, outside the boundary triangle)
indicates a higher resolution.
What is a Mesh?
A mesh is a triangulated representation of a spatial domain \(\mathcal{D}\). It is created by dividing the domain into triangles (or tetrahedrons in 3D). This structure allows the model to approximate the spatial processes over these smaller regions. The mesh helps in defining the precision matrix \(\tilde{\mathbf{\Omega}}\) of the GMRF that represents the spatial random effect. The mesh is passed to the inla.spde2.matern or inla.spde2.pcmatern function to define a spatial model.
Reasons for Triangles Outside the Domain
By including these extra triangles offset, the mesh
ensures more accurate modeling of spatial processes near the boundaries,
improves numerical stability, and enhances the overall quality of the
spatial analysis.
For \(n\times n\) lattice data, a more straightforward approach may involve using an intrinsic GMRF directly. Neighbors are the points directly adjacent on the lattice.
The precision matrix is sparse. \(\tilde{\mathbf{\Omega}}_{ij}=0\) if \(j\notin\{\partial i\cup i\}\), where \(\partial i\) are the indicators of a set of neighbors to location \(\mathbf{s}_i\).
The INLA::inla function is called to fit the model, specifying a Gaussian family and the model components. We also specify control options for the predictor and hyperparameters.
# Define the adjacency structure for the lattice
meshA <- inla.mesh.2d(loc=coords[1:2], max.edge=c(1.5)) #triangulations
meshB <- inla.mesh.2d(loc=coords[1:2], max.edge=c(0.5))
meshC <- inla.mesh.2d(loc=coords[1:2], max.edge=c(1.5, 1.5))
par(mfrow=c(1,3), mar=c(0,0,0,0))
plot(meshA); points(coords, col="red")
plot(meshB); points(coords, col="red")
plot(meshC); points(coords, col="red")
Results
Besag
# Or use the neighborhood structure in an IGMRF model
adj <- matrix(0, n*n, n*n)
for (i in 1:(n*n)) {
xi <- coords$x[i]
yi <- coords$y[i]
neighbors <- with(coords, index[abs(x - xi) + abs(y - yi) == 1]) #rook adjacency
adj[i, neighbors] <- 1
}
adj[1:5,1:5] #neighbors are the points directly adjacent on the lattice## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 1 0 0 0
## [2,] 1 0 1 0 0
## [3,] 0 1 0 1 0
## [4,] 0 0 1 0 1
## [5,] 0 0 0 1 0# Define the ICAR model
model_ICAR <- inla(value ~ 1 + f(index, model="besag", graph=adj,
scale.model=T, constr=T),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
# Summarize the results
summary(model_ICAR)## Time used:
## Pre = 0.773, Running = 0.245, Post = 0.0384, Total = 1.06
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.074 0.106 -0.133 0.074 0.282 0.074 0
##
## Random effects:
## Name Model
## index Besags ICAR model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 9.09e-01 1.28e-01 0.68 9.02e-01
## Precision for index 2.20e+04 2.42e+04 1465.58 1.44e+04
## 0.975quant mode
## Precision for the Gaussian observations 1.18 0.889
## Precision for index 86362.71 3993.208
##
## Deviance Information Criterion (DIC) ...............: 298.20
## Deviance Information Criterion (DIC, saturated) ....: 104.17
## Effective number of parameters .....................: 1.86
##
## Watanabe-Akaike information criterion (WAIC) ...: 298.27
## Effective number of parameters .................: 1.88
##
## Marginal log-Likelihood: -189.53
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_ICAR$mlik## [,1]
## log marginal-likelihood (integration) -189.9294
## log marginal-likelihood (Gaussian) -189.5304
Note: in the unscaled ICAR (the default
unless you set scale.model=T), precision is not equal to
the inverse of a single marginal variance. ICAR is intrinsic/improper
and (after the sum-to-zero constraint) the prior marginal variances
differ by node and depend on the graph.
Matérn SPDE
# Create the spatial structure (SPDE object). See also ?inla.spde.make.index
spde <- inla.spde2.matern(meshC, alpha=2) #smoothness, alpha = nu + d/2
spde$f$hyper.default #default hyperprior## $theta1
## $theta1$prior
## [1] "mvnorm"
##
## $theta1$param
## [1] -1.0104033 -0.2551088 0.1000000 0.0000000 0.0000000 0.1000000
##
## $theta1$initial
## [1] -1.010403
##
##
## $theta2
## $theta2$initial
## [1] -0.2551088# Define the SPDE model
model_SPDE <- inla(value ~ 1 + f(index, model=spde),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_SPDE) # NOT USED!## Time used:
## Pre = 0.686, Running = 0.411, Post = 0.0209, Total = 1.12
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.062 0.115 -0.165 0.061 0.289 0.061 0
##
## Random effects:
## Name Model
## index SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 41502.19 44248.75 3901.42 27993.46
## Theta1 for index -6.52 2.51 -12.29 -6.21
## Theta2 for index 3.27 1.32 1.25 3.10
## 0.975quant mode
## Precision for the Gaussian observations 159174.99 10885.91
## Theta1 for index -2.67 -4.62
## Theta2 for index 6.29 2.27
##
## Deviance Information Criterion (DIC) ...............: -540.49
## Deviance Information Criterion (DIC, saturated) ....: 257.68
## Effective number of parameters .....................: 157.68
##
## Watanabe-Akaike information criterion (WAIC) ...: -579.21
## Effective number of parameters .................: 90.86
##
## Marginal log-Likelihood: -154.71
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_SPDE$mlik # NOT USED!## [,1]
## log marginal-likelihood (integration) -156.3384
## log marginal-likelihood (Gaussian) -154.7138# Does the data arrangement matter? Try index2 (required). Yep.
model_SPDE2 <- inla(value ~ 1 + f(index2, model=spde),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_SPDE2)## Time used:
## Pre = 0.674, Running = 0.341, Post = 0.0235, Total = 1.04
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.052 0.122 -0.19 0.053 0.293 0.053 0
##
## Random effects:
## Name Model
## index2 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 22053.76 2.42e+04 1474.32 14479.41
## Theta1 for index2 -4.12 1.02e+00 -6.38 -4.02
## Theta2 for index2 2.01 5.81e-01 1.06 1.96
## 0.975quant mode
## Precision for the Gaussian observations 86219.24 4021.57
## Theta1 for index2 -2.46 -3.56
## Theta2 for index2 3.31 1.69
##
## Deviance Information Criterion (DIC) ...............: -557.86
## Deviance Information Criterion (DIC, saturated) ....: 224.36
## Effective number of parameters .....................: 124.36
##
## Watanabe-Akaike information criterion (WAIC) ...: -577.99
## Effective number of parameters .................: 81.89
##
## Marginal log-Likelihood: -155.94
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_SPDE2$mlik## [,1]
## log marginal-likelihood (integration) -157.2824
## log marginal-likelihood (Gaussian) -155.9430
By default (double check in progress; Lindgren & Rue, 2015, p. 4),
\(\left(\ln\tau,\ \ln\rho\right)^\top\sim\boldsymbol{\mathcal{N}}_{2}\underset{\hspace{6.5em}\color{gray}{\succ\ 0}}{\left(\begin{pmatrix} \theta_1 \\ \theta_2 \end{pmatrix},\ \begin{pmatrix} \sigma_1^2 & v_{12} \\ v_{12} & \sigma_2^2 \end{pmatrix}\right)}\), \(\qquad\sigma_1^2=\sigma_2^2=1/0.1=10\), and \(v_{12}=0\)
“Theta1 for index” refers to \(\ln\hat{\tau}\), the log (marginal) precision of the spatial random field. Note \(\tau^2=\frac{\Gamma(\nu)}{\Gamma(\alpha)(4\pi)^{d/2}\kappa^{2\nu}}\ {\color{red}{\hat{\sigma}_s^{-2}}}\).
“Theta2 for index” refers to \(\ln\hat{\kappa}\), the log spatial scale (how quickly spatial correlation decays). Note \(\kappa=\frac{2\sqrt{2}}{\rho}\).
The spatial range \(\rho\) describes the distance over which spatial correlation remains significant. Beyond this distance, the correlation between spatial points drops off quickly. A larger \(\rho\) corresponds to slower spatial correlation decay.
Tobler’s First Law of Geography
“Everything is related to
everything else, but near things are more related than distant
things.”
Note: the most principled way to report \({\color{red}{\hat{\sigma}_s^{2}}}\) is to
transform the full posterior of Theta1 and
Theta2 rather than plugging in medians.
spde_res <- inla.spde2.result(model_SPDE2, name="index2", spde=spde, do.transform=T)
spde_res$marginals.variance.nominal # marginal variance
spde_res$marginals.range.nominal # range
Matérn SPDE with PC
# Create the spatial structure (SPDE object with PC prior)
spde_pc <- inla.spde2.pcmatern(meshC,
alpha=2, #smoothness, alpha = nu + d/2
prior.range=c(0.1, 0.01), #P(range < 0.1) = 0.01
prior.sigma=c(1, 0.01) #P(sigma_s > 1) = 0.01
)
# Define the SPDE model
model_SPDE_PC <- inla(value ~ 1 + f(index, model=spde_pc),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_SPDE_PC) # NOT USED!## Time used:
## Pre = 0.681, Running = 0.358, Post = 0.0241, Total = 1.06
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.067 0.124 -0.176 0.067 0.313 0.067 0
##
## Random effects:
## Name Model
## index SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 2.21e+04 2.47e+04 1411.424 1.43e+04
## Range for index 6.32e-01 1.60e-01 0.364 6.16e-01
## Stdev for index 1.60e+00 2.74e-01 1.146 1.57e+00
## 0.975quant mode
## Precision for the Gaussian observations 8.76e+04 3814.819
## Range for index 9.88e-01 0.589
## Stdev for index 2.22e+00 1.502
##
## Deviance Information Criterion (DIC) ...............: -562.88
## Deviance Information Criterion (DIC, saturated) ....: 221.50
## Effective number of parameters .....................: 121.49
##
## Watanabe-Akaike information criterion (WAIC) ...: -581.59
## Effective number of parameters .................: 80.27
##
## Marginal log-Likelihood: -158.47
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_SPDE_PC$mlik # NOT USED!## [,1]
## log marginal-likelihood (integration) -159.7998
## log marginal-likelihood (Gaussian) -158.4753# Does the data arrangement matter? Try index2 (required). Yep.
model_SPDE_PC2 <- inla(value ~ 1 + f(index2, model=spde_pc),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_SPDE_PC2)## Time used:
## Pre = 0.676, Running = 0.349, Post = 0.028, Total = 1.05
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.048 0.132 -0.213 0.048 0.307 0.048 0
##
## Random effects:
## Name Model
## index2 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 2.21e+04 2.48e+04 1417.349 1.43e+04
## Range for index2 7.36e-01 1.90e-01 0.416 7.18e-01
## Stdev for index2 1.46e+00 2.26e-01 1.080 1.43e+00
## 0.975quant mode
## Precision for the Gaussian observations 88050.33 3831.479
## Range for index2 1.16 0.689
## Stdev for index2 1.97 1.372
##
## Deviance Information Criterion (DIC) ...............: -562.55
## Deviance Information Criterion (DIC, saturated) ....: 221.31
## Effective number of parameters .....................: 121.31
##
## Watanabe-Akaike information criterion (WAIC) ...: -580.74
## Effective number of parameters .................: 80.61
##
## Marginal log-Likelihood: -158.29
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_SPDE_PC2$mlik## [,1]
## log marginal-likelihood (integration) -159.6201
## log marginal-likelihood (Gaussian) -158.2880
Recall Eq. 6. See ?inla.spde2.pcmatern. Should not be
confused with the PC
prior for precision. \[\begin{equation*}
\tag{6}
\pi({\color{red}{\sigma_s}},\
\rho)=\frac{1}{2}d\lambda_\rho\rho^{-1-d/2}\exp\left\{-\lambda_\rho\rho^{-d/2}\right\}\cdot\lambda_{{\color{red}{\sigma_s}}}\exp\left\{-\lambda_{{\color{red}{\sigma_s}}}{\color{red}{\sigma_s}}\right\},
\end{equation*}\]
If you want the exact posterior summary (not just transforming point estimates), transform the variance marginal.
spde_res <- inla.spde2.result(model_SPDE_PC2, name="index2", spde=spde_pc, do.transform=T)
spde_res$marginals.variance.nominal
Intercept-Only
# Define the intercept-only model (true model)
model_non <- inla(value ~ 1,
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_non)## Time used:
## Pre = 0.561, Running = 0.19, Post = 0.013, Total = 0.764
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.074 0.106 -0.134 0.074 0.282 0.074 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.907 0.128 0.675 0.901
## 0.975quant mode
## Precision for the Gaussian observations 1.17 0.889
##
## Deviance Information Criterion (DIC) ...............: 298.50
## Deviance Information Criterion (DIC, saturated) ....: 104.34
## Effective number of parameters .....................: 1.99
##
## Watanabe-Akaike information criterion (WAIC) ...: 298.56
## Effective number of parameters .................: 2.00
##
## Marginal log-Likelihood: -159.64
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_non$mlik## [,1]
## log marginal-likelihood (integration) -159.2558
## log marginal-likelihood (Gaussian) -159.6395bf10 <- exp(model_SPDE2$mlik - model_non$mlik)
rownames(bf10) <- c("integration", "Gaussian"); colnames(bf10) <- "BF10.est"
bf10## BF10.est
## integration 7.195635
## Gaussian 40.304056bf10pc <- exp(model_SPDE_PC2$mlik - model_non$mlik)
rownames(bf10pc) <- c("integration", "Gaussian"); colnames(bf10pc) <- "BF10.est"
bf10pc## BF10.est
## integration 0.6947325
## Gaussian 3.8631323
Random Intercepts
Read more about linear mixed-effects models: https://rpubs.com/sherloconan/1184363.
# Define the linear mixed-effects model
model_iid <- inla(value ~ 1 + f(index2, model="iid"
# , hyper=list(prec=list(param=c(0.01, 0.01)))
),
family="gaussian", data=data,
# control.family=list(hyper=list(prec=list(param=c(0.01, 0.01), initial=1))),
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_iid)## Time used:
## Pre = 0.6, Running = 0.257, Post = 0.034, Total = 0.89
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.074 0.106 -0.133 0.074 0.281 0.074 0
##
## Random effects:
## Name Model
## index2 IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 2.24e+04 2.62e+04 1388.791 1.41e+04
## Precision for index2 9.09e-01 1.28e-01 0.677 9.02e-01
## 0.975quant mode
## Precision for the Gaussian observations 91609.62 3724.863
## Precision for index2 1.18 0.893
##
## Deviance Information Criterion (DIC) ...............: -539.04
## Deviance Information Criterion (DIC, saturated) ....: 230.96
## Effective number of parameters .....................: 130.96
##
## Watanabe-Akaike information criterion (WAIC) ...: -558.39
## Effective number of parameters .................: 85.13
##
## Marginal log-Likelihood: -159.69
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_iid$mlik## [,1]
## log marginal-likelihood (integration) -160.1105
## log marginal-likelihood (Gaussian) -159.6902bf10r <- exp(model_SPDE2$mlik - model_iid$mlik)
rownames(bf10r) <- c("integration", "Gaussian"); colnames(bf10r) <- "BF10.est"
bf10r## BF10.est
## integration 16.91352
## Gaussian 42.40262bf10rpc <- exp(model_SPDE_PC2$mlik - model_iid$mlik)
rownames(bf10rpc) <- c("integration", "Gaussian"); colnames(bf10rpc) <- "BF10.est"
bf10rpc## BF10.est
## integration 1.632986
## Gaussian 4.064279
Random Walk
# Random walk of order 2 in two dimensions
model_rw2d <- inla(value ~ 1 + f(index2, model="rw2d", nrow=n, ncol=n),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_rw2d)## Time used:
## Pre = 0.602, Running = 0.25, Post = 0.0336, Total = 0.886
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.074 0.106 -0.133 0.074 0.282 0.074 0
##
## Random effects:
## Name Model
## index2 Random walk 2D
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 9.1e-01 0.13 0.678 9.02e-01
## Precision for index2 2.2e+04 24213.44 1444.672 1.44e+04
## 0.975quant mode
## Precision for the Gaussian observations 1.19 0.889
## Precision for index2 86240.21 3928.356
##
## Deviance Information Criterion (DIC) ...............: 300.10
## Deviance Information Criterion (DIC, saturated) ....: 106.16
## Effective number of parameters .....................: 3.87
##
## Watanabe-Akaike information criterion (WAIC) ...: 300.29
## Effective number of parameters .................: 3.90
##
## Marginal log-Likelihood: -257.16
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_rw2d$mlik## [,1]
## log marginal-likelihood (integration) -257.5552
## log marginal-likelihood (Gaussian) -257.1612
Matérn GRF
Recall Eq. 3.
# Matern2D
model_m2d <- inla(value ~ 1 + f(index2, model="matern2d", nrow=n, ncol=n),
family="gaussian", data=data,
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_m2d)## Time used:
## Pre = 0.604, Running = 0.298, Post = 0.0212, Total = 0.923
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.074 0.106 -0.134 0.074 0.282 0.074 0
##
## Random effects:
## Name Model
## index2 Matern2D model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 9.09e-01 1.28e-01 0.68 9.02e-01
## Precision for index2 2.82e+04 3.39e+04 2318.14 1.79e+04
## Range for index2 1.39e+02 1.66e+02 10.97 8.80e+01
## 0.975quant mode
## Precision for the Gaussian observations 1.18e+00 0.888
## Precision for index2 1.17e+05 6295.873
## Range for index2 5.71e+02 29.876
##
## Deviance Information Criterion (DIC) ...............: 298.56
## Deviance Information Criterion (DIC, saturated) ....: 104.60
## Effective number of parameters .....................: 2.02
##
## Watanabe-Akaike information criterion (WAIC) ...: 298.62
## Effective number of parameters .................: 2.03
##
## Marginal log-Likelihood: -159.80
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')model_m2d$mlik## [,1]
## log marginal-likelihood (integration) -161.1292
## log marginal-likelihood (Gaussian) -159.7985bf01m <- exp(model_non$mlik - model_m2d$mlik)
rownames(bf01m) <- c("integration", "Gaussian"); colnames(bf01m) <- "BF01.est"
bf01m## BF01.est
## integration 6.510369
## Gaussian 1.172373bf01m2 <- exp(model_iid$mlik - model_m2d$mlik)
rownames(bf01m2) <- c("integration", "Gaussian"); colnames(bf01m2) <- "BF01.est"
bf01m2## BF01.est
## integration 2.769751
## Gaussian 1.114350
Discussion
What Went Wrong?
Estimates are sensitive to
the data arrangement: sorting by columns (required) or rows
scaling: e.g., police-reported incidents per 100,000 or 10,000 population; see https://rpubs.com/sherloconan/1187595 “Data \(\times100\)”
the mesh’s edge length:
max.edge=c(inside the boundary triangle, outside the boundary triangle)the model and priors: inla.spde2.matern or inla.spde2.pcmatern
the method: integration or Gaussian for
mliknumerical fluctuations: parallel computations
\(10\times10\) data points. The true data-generating model is the intercept-only model. Columns 3–7 display the posterior median estimates.
| True Value | Besag | SPDE | SPDE + PC | Intercept-Only | Random Intercepts | |
|---|---|---|---|---|---|---|
| Grand Mean | 0 | 0.074 | 0.053 | 0.048 | 0.074 | 0.074 |
| Error Precision | 1 | 0.902 | \(1.45\times10^4\) | \(1.43\times10^4\) | 0.901 | \(1.43\times10^4\) |
| Latent Precision | NA | \(1.44\times10^4\) | \(\exp\{-4.02\}\approx0.018\) | \(1/1.43^2\approx0.49\) | NA | 0.902 |
| Log Marginal Likelihood | -189.929 | -157.282 | -159.620 | -159.256 | -160.120 |
Two methods of computing the log marginal likelihood are
Integration: This approach computes the
mlikby directly integrating the likelihood function over the parameter space. This method is typically more accurate but can be computationally intensive. It’s often used when precision is essential for model comparison.Gaussian: This method approximates
mlikusing a normal approximation. It is generally faster and less computationally demanding than direct integration, but it might be less accurate, particularly in cases where the true likelihood is not well-approximated by a normal distribution (Rue, Martino, & Chopin, 2009, p. 345).
R Script
## Spatial Models using Stochastic Partial Differential Equations
## https://becarioprecario.bitbucket.io/inla-gitbook/ch-spatial.html#spatial-models-using-stochastic-partial-differential-equations
library(INLA)
n <- 10
coords <- expand.grid(x=1:n, y=1:n) #10×10 lattice
coords$index <- NA #indices 1-100
coords[order(coords$x, -coords$y), "index"] <- 1:nrow(coords) #data arrangement
set.seed(277)
data <- cbind(coords, "value"=rnorm(n*n))
# Define a set of basis functions using a two-dimensional mesh
mesh <- inla.mesh.2d(loc=coords[1:2], max.edge=c(1.5, 1.5))
matern <- inla.spde2.matern(mesh, alpha=2) #smoothness, alpha = nu + d/2
matA <- inla.spde.make.A(mesh=mesh, loc=coordinates(coords)) #projector matrix A, mapping the projection of the SPDE to the observed points
sindex <- inla.spde.make.index(name="spatial.field", n.spde=matern$n.spde) #named index vectors for the SPDE model
# Create data structure (no kriging in this case)
data_stack <- inla.stack(data=list(value=data$value),
A=list(matA), #lengths of `A` and `effects` should be equal
effects=list(c(sindex, list(Intercept=1)))) #no covariates in this case
model_SPDE <- inla(value ~ -1 + Intercept + f(spatial.field, model=spde),
family="gaussian", data=inla.stack.data(data_stack, spde=matern),
control.predictor=list(A=inla.stack.A(data_stack), compute=T),
control.compute=list("config"=T, "dic"=T, "waic"=T))
summary(model_SPDE)
model_SPDE$mlik
#Compute statistics in terms or range and variance
spde_est <- inla.spde2.result(inla=model_SPDE, name="spatial.field", spde=matern, do.transf=T)
inla.zmarginal(spde_est$marginals.variance.nominal[[1]]) #(marginal) spatial variance estimate
inla.zmarginal(spde_est$marginals.range.nominal[[1]]) #spatial range estimate
Parting Words
These are basic examples and INLA is capable of much more complex
spatial analyses, including handling non-Gaussian data, incorporating
covariates, and modeling non-stationary processes. See
?inla.spde2.matern: Non-stationary models are supported for
alpha=2 only.
NEXT: https://rpubs.com/sherloconan/1339374
Error in inla.call.builtin() : INLA installation error; no such file
From version 24.05.10, INLA is built for R-4.4.
- R version 4.4.0 (Puppy Cup) has been released on 2024-04-24.
- R version 4.3.3 (Angel Food Cake, 2024-02-29) binary for macOS 11 (Big Sur) and higher and Apple silicon (M1-3) Macs: Download