I started this thread 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 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:
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.
Hence my questions:
- Is the model correctly parametrized? The drawing of the \beta is based on this post which uses the inverse of the upper cholesky factor of the precision matrix for the non-centered parametrization.
- 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"))