# Fitting Matern spatial random fields in Stan

**URL:** <https://discourse.mc-stan.org/t/fitting-matern-spatial-random-fields-in-stan/25069>\
**Category:** Modeling\
**Tags:** spatial, rstan, specification\
**Created:** [November 1, 2021, 8:34am UTC](https://discourse.mc-stan.org/t/fitting-matern-spatial-random-fields-in-stan/25069 "2021-11-01T08:34:07Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![lionel68](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lionel68/32/1284_2.png) [@lionel68](https://discourse.mc-stan.org/u/lionel68)\
**Post date:** [November 1, 2021, 8:34am UTC](https://discourse.mc-stan.org/t/fitting-matern-spatial-random-fields-in-stan/25069/1 "2021-11-01T08:34:07Z")

</div>

I started [this thread](https://discourse.mc-stan.org/t/spde-models-in-brms/17346) a while ago and given that the topic is moving away from brms, I am starting this new thread.

The idea is to fit spatial random field following a Matern correlation in (R)Stan, based on the results from [this article](https://arxiv.org/pdf/2001.07623.pdf) and using the SPDE approach to estimate the spatial effect followed by INLA we have:

y \sim N(\mu, \sigma)  
\mu = \alpha + A \* \beta  
\beta \sim MVN(0, S^{-1})  
S = \tau^2 \* (\kappa^4 \* C + 2 \* \kappa^2 \* G\_1 + G\_2)

Where \tau and \kappa are the scale and variance of the spatial effect, S is the precision matrix for the spatial effects \beta and A is a matrix projecting from the mesh (where the spatial effect are quantified) to the observations y.

* * *

Writing this model in Stan gives:

```stan
data {
  int<lower=0> N; // number of observations
  int<lower=0> n_knots; // number of interpolation points on the mesh
  matrix[n_knots, n_knots] G[3]; // FEM matrices
  matrix[N, n_knots] A; // projection matrix
  vector[N] y; // the observations
}
parameters {
  real mu; // the intercept 
  real<lower=0> sigma; // residual deviation
  real<lower=0> kappa; // variance of spatial effect
  real<lower=0> tau; // scale of spatial effect
  vector[n_knots] alpha; // std_normal draws for non-centered parametrization
}
transformed parameters {
  matrix[n_knots, n_knots] S; // precision matrix
  vector[n_knots] beta; // spatial coefficient
  S = (tau ^ 2) * (G[1] * kappa ^ 4 + G[2] * 2 * kappa ^ 2 + G[3]);
  beta = cholesky_decompose(S)' \ alpha;
  // implies beta ~ MVN(0, inverse(S))
}
model {
  alpha ~ std_normal();
  mu ~ normal(0, 1);
  sigma ~ lognormal(0, 1);
  kappa ~ lognormal(0, 1);
  tau ~ lognormal(0, 1);
  y ~ normal(mu + A * beta, sigma);
}

```

* * *

After simulating some data in R (see code at the end of the post), the model ended up with divergent transitions and low effective sample sizes. Looking at the pair plots it seems that the scale parameter (\tau) is especially problematic.

 ![spde_stan](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/0b90371a96f3aeefb718b1292b2fdd0916226aa8.png)

Hence my questions:

1. Is the model correctly parametrized? The drawing of the \beta is based on [this post](http://www.statsathome.com/2018/10/19/sampling-from-multivariate-normal-precision-and-covariance-parameterizations/) which uses the inverse of the **upper** cholesky factor of the precision matrix for the non-centered parametrization.
2. How to improve the sampling of \tau?

* * *

Now the R code to generate the data and fit the model:

```
library(rstan)
library(INLA)
library(RandomFields)

# generate 100 x/y coordinates
dat <- data.frame(x = runif(100),
                  y = runif(100))

# make a spatial random field
spat_field <- raster::raster(RFsimulate(RMmatern(nu=1, var = 1, scale = 0.1),
                         x = seq(0, 1, length.out = 100),
                         y = seq(0, 1, length.out = 100),
                         spConform = FALSE))

# draw response values
dat$resp <- rnorm(100, mean = 1 + raster::extract(spat_field, dat), sd = 1)

# the boundary of the spatial field
bnd <- inla.mesh.segment(matrix(c(0, 0, 1, 0, 1, 1, 0, 1), 4, 2, byrow = TRUE))

# a coarse mesh
mesh <- inla.mesh.2d(max.edge = 0.2,
                     offset = 0,
                     boundary = bnd)

# derive the FEM matrices
fem <- inla.mesh.fem(mesh)
# put the matrices in one object
G <- array(c(as(fem$c1, "matrix"),
             as(fem$g1, "matrix"),
             as(fem$g2, "matrix")),
          dim = c(mesh$n, mesh$n, 3))
G <- aperm(G, c(3, 1, 2))

# derive the projection matrix
A <- inla.spde.make.A(mesh, loc = as.matrix(dat[,c("x", "y")]))
A <- as(A, "matrix")

# put all R object into a list
dat_stan <- list(N = nrow(dat),
                 y = dat$resp,
                 n_knots = mesh$n,
                 G = G,
                 A = A)

# fit the model
m <- stan("Programming_stuff/spde_mgcv/spde_stan_new.stan",
          data = dat_stan,
          pars = c("mu", "sigma", "kappa", "tau"))

```

---

<div class="post-metadata">

**Author:** ![nikuehn](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikuehn/32/7595_2.png) [@nikuehn](https://discourse.mc-stan.org/u/nikuehn)\
**Post date:** [November 2, 2021, 12:24pm UTC](https://discourse.mc-stan.org/t/fitting-matern-spatial-random-fields-in-stan/25069/2 "2021-11-02T12:24:14Z")

</div>

I haven’t checked whether the model is correctly parameterized (seems ok), but maybe it helps to parameterize \tau and \kappa in terms of spatial range and standard deviation, as in Equation (4) to (7) of Lindgren and Rue (2015).

Lindgren, F., & Rue, H. (2015). Bayesian Spatial Modelling with R-INLA. _Journal of Statistical Software_ , _63_ (19), 1–25. [Bayesian Spatial Modelling with R-INLA | Journal of Statistical Software](https://doi.org/10.18637/jss.v063.i19)

---

<div class="post-metadata">

**Author:** ![cavieresJJ](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/cavieresjj/32/5730_2.png) [@cavieresJJ](https://discourse.mc-stan.org/u/cavieresJJ)\
**Post date:** [November 2, 2021, 2:54pm UTC](https://discourse.mc-stan.org/t/fitting-matern-spatial-random-fields-in-stan/25069/3 "2021-11-02T14:54:59Z")

</div>

HI!

As say @nikuehn, a best option could be parameterizer the spatial random field using the range (\rho) and the spatial standard deviation (\sigma\_s). You could start trying with a Gamma or Exponential prior for those parameters.

Best wishes