Matern fields over a mesh: spde()

spde(fem, gr = node) in a formula adds a gaussian Matern field through the SPDE representation of Lindgren, Rue and Lindstrom (2011): on a fixed mesh the field's precision is Q = tau^2 (kappa^4 M0 + 2 kappa^2 M1 + M2), with the mesh's three finite-element matrices as fixed data. It is a predictor term, not a bar term, and it is available once frmtmb.spatial is loaded: library(frmtmb.spatial). Without the package, frmtmb::frm() refuses the term and names this package. spde() is frmtmb's own spelling; brms has no such term.

Value

spde() is a formula term, not a function: frmtmb::frm() reads it at parse time, and the value it contributes is a random-effect block of the model. This page documents the term and returns nothing.

Arguments of the term

fem

The mesh's finite-element matrices, as a list holding c0, g1 and g2 (fmesher::fm_fem()'s names) or M0, M1 and M2 (INLA's). They must be square and of one size, one row per mesh node. Put the list in data2 = list(fem = fem). Mesh construction stays outside the package. alpha = 2, the order the three-matrix form covers, is the only one.

gr

The node of each observation, by ROW NUMBER: whole numbers in 1..nrow(M0), as an integer or numeric column, a character column of those numbers, or a factor whose levels are those numbers. The matrices carry no dimnames, so nothing else can reconcile a set of labels with them, and any other labeling is refused. A node no row visits keeps its column and its prior.

What the fit reports

frmtmb::confint_varcorr() reports the marginal standard deviation sd(spde) and the range range(spde), under the planar alpha = 2 identities range = sqrt(8) / kappa and sigma^2 = 1 / (4 pi kappa^2 tau^2). Every mesh node is a level of the block, so newdata names a node by its row number too. A general projector matrix, which spreads an observation over the vertices of the triangle it falls in, is not supported yet.

See also

vignette("spatial", package = "frmtmb.spatial"), and frmtmb::frmtmb_register_covstruct() for the seam the term is registered through.

Examples

# the alpha = 2 finite-element triple of a regular 1-D mesh
nn <- 10
h <- 0.5
C0 <- diag(rep(h, nn))
C0[1, 1] <- C0[nn, nn] <- h / 2
G <- matrix(0, nn, nn)
for (i in seq_len(nn - 1)) {
  G[i, i] <- G[i, i] + 1 / h
  G[i + 1, i + 1] <- G[i + 1, i + 1] + 1 / h
  G[i, i + 1] <- G[i + 1, i] <- -1 / h
}
fem <- list(c0 = C0, g1 = G, g2 = G %*% solve(C0) %*% G)
set.seed(5)
d <- data.frame(node = rep(seq_len(nn), each = 5), x = rnorm(5 * nn))
d$y <- 0.3 + 0.4 * d$x + sin(d$node / 2) + rnorm(nrow(d), 0, 0.3)
fit <- frm(bf(y ~ x + spde(fem, gr = node)), family = gaussian(),
           data = d, data2 = list(fem = fem))
confint_varcorr(fit)
#>                  block        term  type  estimate       lwr        upr
#> 1 spde(fem, gr = node)    sd(spde)    sd 0.3054464 0.1517612  0.6147655
#> 2 spde(fem, gr = node) range(spde) range 4.3155470 1.5616748 11.9256237