Intro to Spatial Modeling Using INLA

Large Sample Size

============================================================================================================ PREVIOUS: https://rpubs.com/sherloconan/1164312

THIS: https://rpubs.com/sherloconan/1198601

MORE: 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.10 25.06.07 (stable)
• fmesher version 0.1.6 0.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 <- 200

coords <- expand.grid(x=1:n, y=1:n) #n×n lattice
coords$index2 <- coords$index <- 1:nrow(coords) #indices 1–n×n
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))

   

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)

 

scroll to top

   

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 alpha implies 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")

 

scroll to top

   

Results

Besag

 

inla.core.safe: The inla program failed, but will rerun in case better initial values may help. try=1/1

Error: ! vector memory limit of 24.0 Gb reached, see mem.maxVSize()

 

# 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
}

# 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)
model_ICAR$mlik

 

scroll to top

   

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]  2.085677 -3.351189  0.100000  0.000000  0.000000  0.100000
## 
## $theta1$initial
## [1] 2.085677
## 
## 
## $theta2
## $theta2$initial
## [1] -3.351189

# Define the SPDE model
model_SPDE2 <- inla(value ~ 1 + f(index2, model=spde),
                    family="gaussian", data=data,
                    control.compute=list("config"=T, "dic"=T, "waic"=T))
# roughly 2 minutes of run time

summary(model_SPDE2)

## Time used:
##     Pre = 0.914, Running = 120, Post = 5.12, Total = 126 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.006 0.006     -0.005    0.006      0.016 0.006   0
## 
## Random effects:
##   Name     Model
##     index2 SPDE2 model
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  1.00 0.007       0.99     1.00
## Theta1 for index2                       14.53 7.923       3.52    13.33
## Theta2 for index2                       -4.23 2.238      -9.35    -4.01
##                                         0.975quant  mode
## Precision for the Gaussian observations      1.018  1.00
## Theta1 for index2                           33.178  7.91
## Theta2 for index2                           -0.734 -2.72
## 
## Deviance Information Criterion (DIC) ...............: 113371.93
## Deviance Information Criterion (DIC, saturated) ....: 39980.63
## Effective number of parameters .....................: 9.26
## 
## Watanabe-Akaike information criterion (WAIC) ...: 113383.16
## Effective number of parameters .................: 20.60
## 
## Marginal log-Likelihood:  -56708.46 
##  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) -56709.96
## log marginal-likelihood (Gaussian)    -56708.46

 

By default (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

 

scroll to top

   

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_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 = 1.24, Running = 104, Post = 4.71, Total = 110 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.006 0.005     -0.004    0.006      0.016 0.006   0
## 
## Random effects:
##   Name     Model
##     index2 SPDE2 model
## 
## Model hyperparameters:
##                                           mean     sd 0.025quant 0.5quant
## Precision for the Gaussian observations 37.086 14.528     11.169   37.084
## Range for index2                         0.213  0.023      0.171    0.212
## Stdev for index2                         3.828  0.393      3.135    3.801
##                                         0.975quant   mode
## Precision for the Gaussian observations     61.832 33.356
## Range for index2                             0.259  0.211
## Stdev for index2                             4.679  3.732
## 
## Deviance Information Criterion (DIC) ...............: 10973.24
## Deviance Information Criterion (DIC, saturated) ....: 90620.29
## Effective number of parameters .....................: 50618.07
## 
## Watanabe-Akaike information criterion (WAIC) ...: -10121.38
## Effective number of parameters .................: 21701.66
## 
## Marginal log-Likelihood:  -56738.79 
##  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) -56740.28
## log marginal-likelihood (Gaussian)    -56738.79

 

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*}\]

 

scroll to top

   

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.587, Running = 0.695, Post = 0.0689, Total = 1.35 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.006 0.005     -0.004    0.006      0.016 0.006   0
## 
## Model hyperparameters:
##                                         mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 1.00 0.007       0.99     1.00
##                                         0.975quant mode
## Precision for the Gaussian observations       1.02 1.00
## 
## Deviance Information Criterion (DIC) ...............: 113380.95
## Deviance Information Criterion (DIC, saturated) ....: 40004.31
## Effective number of parameters .....................: 1.99
## 
## Watanabe-Akaike information criterion (WAIC) ...: 113380.96
## Effective number of parameters .................: 2.00
## 
## Marginal log-Likelihood:  -56706.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_non$mlik

##                                            [,1]
## log marginal-likelihood (integration) -56706.41
## log marginal-likelihood (Gaussian)    -56706.80

bf10 <- exp(model_SPDE2$mlik - model_non$mlik)
rownames(bf10) <- c("integration", "Gaussian"); colnames(bf10) <- "BF10.est"
bf10

##              BF10.est
## integration 0.0284467
## Gaussian    0.1898565

bf10pc <- exp(model_SPDE_PC2$mlik - model_non$mlik)
rownames(bf10pc) <- c("integration", "Gaussian"); colnames(bf10pc) <- "BF10.est"
bf10pc

##                 BF10.est
## integration 1.947515e-15
## Gaussian    1.274840e-14

 

scroll to top

   

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.666, Running = 1.96, Post = 1.5, Total = 4.12 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.006 0.005     -0.004    0.006      0.016 0.006   0
## 
## Random effects:
##   Name     Model
##     index2 IID model
## 
## Model hyperparameters:
##                                         mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 2.22 0.629       1.24     2.13
## Precision for index2                    1.95 0.492       1.16     1.90
##                                         0.975quant mode
## Precision for the Gaussian observations       3.69 1.97
## Precision for index2                          3.08 1.79
## 
## Deviance Information Criterion (DIC) ...............: 103365.20
## Deviance Information Criterion (DIC, saturated) ....: 60553.97
## Effective number of parameters .....................: 20550.42
## 
## Watanabe-Akaike information criterion (WAIC) ...: 103432.72
## Effective number of parameters .................: 15818.14
## 
## Marginal log-Likelihood:  -56715.03 
##  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) -56715.55
## log marginal-likelihood (Gaussian)    -56715.03

bf10r <- exp(model_SPDE2$mlik - model_iid$mlik)
rownames(bf10r) <- c("integration", "Gaussian"); colnames(bf10r) <- "BF10.est"
bf10r

##             BF10.est
## integration 266.1449
## Gaussian    715.9357

bf10rpc <- exp(model_SPDE_PC2$mlik - model_iid$mlik)
rownames(bf10rpc) <- c("integration", "Gaussian"); colnames(bf10rpc) <- "BF10.est"
bf10rpc

##                 BF10.est
## integration 1.822078e-11
## Gaussian    4.807331e-11

 

scroll to top

   

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.622, Running = 23.7, Post = 2.61, Total = 27 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.006 0.005     -0.004    0.006      0.016 0.006   0
## 
## Random effects:
##   Name     Model
##     index2 Random walk 2D
## 
## Model hyperparameters:
##                                              mean       sd 0.025quant  0.5quant
## Precision for the Gaussian observations      1.00 7.00e-03       0.99      1.00
## Precision for index2                    117860.96 4.36e+04   52714.17 111089.14
##                                         0.975quant     mode
## Precision for the Gaussian observations   1.02e+00     1.00
## Precision for index2                      2.22e+05 98606.76
## 
## Deviance Information Criterion (DIC) ...............: 113397.27
## Deviance Information Criterion (DIC, saturated) ....: 40024.15
## Effective number of parameters .....................: 22.89
## 
## Watanabe-Akaike information criterion (WAIC) ...: 113397.24
## Effective number of parameters .................: 22.84
## 
## Marginal log-Likelihood:  -102892.26 
##  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) -102892.7
## log marginal-likelihood (Gaussian)    -102892.3

 

scroll to top

   

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.623, Running = 26.9, Post = 1.71, Total = 29.2 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.005 0.005     -0.005    0.005      0.016 0.005   0
## 
## Random effects:
##   Name     Model
##     index2 Matern2D model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for the Gaussian observations 1.00e+00 7.00e-03      0.989     1.00
## Precision for index2                    1.79e+05 1.35e+07      0.000     7.67
## Range for index2                        7.00e+08 3.99e+10     12.755  3357.80
##                                         0.975quant  mode
## Precision for the Gaussian observations   1.02e+00  1.00
## Precision for index2                      6.18e+05  0.00
## Range for index2                          1.00e+08 13.90
## 
## Deviance Information Criterion (DIC) ...............: 113380.79
## Deviance Information Criterion (DIC, saturated) ....: 40037.19
## Effective number of parameters .....................: 4.06
## 
## Watanabe-Akaike information criterion (WAIC) ...: 113380.86
## Effective number of parameters .................: 4.12
## 
## Marginal log-Likelihood:  -56705.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_m2d$mlik

##                                            [,1]
## log marginal-likelihood (integration) -56707.37
## log marginal-likelihood (Gaussian)    -56705.29

bf01m <- exp(model_non$mlik - model_m2d$mlik)
rownames(bf01m) <- c("integration", "Gaussian"); colnames(bf01m) <- "BF01.est"
bf01m

##             BF01.est
## integration 2.616026
## Gaussian    0.220960

bf01m2 <- exp(model_iid$mlik - model_m2d$mlik)
rownames(bf01m2) <- c("integration", "Gaussian"); colnames(bf01m2) <- "BF01.est"
bf01m2

##                 BF01.est
## integration 2.796120e-04
## Gaussian    5.859562e-05

 

scroll to top

   

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 mlik

• numerical fluctuations: parallel computations

 

\(200\times200\) data points. The true data-generating model is the intercept-only model. Columns 3–6 display the posterior median estimates.

	True Value	SPDE	SPDE + PC	Intercept-Only	Random Intercepts
Grand Mean	0	0.006	0.006	0.006	0.006
Error Precision	1	1	37.100	1	1.59
Latent Precision	NA	\(\exp\{14.14\}\approx1.38\times10^{6}\)	\(1/3.801^2\approx0.069\)	NA	2.72
Log Marginal Likelihood		-56709.97	-56740.28	-56706.41	-56715.66

   

Two methods of computing the log marginal likelihood are

• Integration: This approach computes the mlik by 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 mlik using 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).

 

scroll to top

   

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 <- 200

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

 

scroll to top

   

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