Hieararchical 2D+1D space-time modelling using brms

Interfaces brms
1

I just posted a similar thread asking about approaching this topic using raw Stan & Gaussian Processes, but I wanted to break out this brms-specific query as a separate topic.

[background recap]
I have space-time data (x,y,time) from the neuroimaging field where the data is collected hierarchically such that there are multiple “trials” of data from each of multiple individuals, and I’d like to achieve inference on the 2D+1D smoothed hypersurface averaged-across-trials-and-individuals while permitting systematic deviations from this surface in each individual (and trials just being random noise deviates from the individual’s surface), ideally with shrinkage applied to the magnitude of the individual deviations.

[brms-specific query]
I see that there’s already been discussion of non-hierarchical 2D+1D space-time models, using splines in brms. I’m less familiar with brms than stan, so I’m curious if others might be able to weigh in on whether it would be straightforward to add hierarchy to the smoothing in the manner I describe above.

(tagging @ucfagls and @MilaniC, who seem to have pertinent expertise)

1 Like
2

You can do this with various combinations of splines. How you do it will depend on whether you think the 2d surface is isotropic or not in terms of the wiggliness. Either way you’ll need to use tensor products of smooths, as unless I am mistaken there isn’t an efficient way to mingle/interact a spline built with mgcv’s spline constructors (which brms uses) and brms more efficient random effects.

Colleagues and I described how this would work in a recent paper focused on mgcv, but the model construction applies equally to brms. See doi.org/10.7717/peerj.6876 and especially the hypothetical bird movement examples (except note that you have to use t2() for tensor products in brms because if the way it sets up the mixed model form of the GAM.)

Briefly, if you want a more random effect approach (“random surfaces” varying in time) you could fit models in two modes:

  1. what we called a global smooth plus random smooths model, or
  2. just the random smooths model

These models would assume the same general degree of wiggliness for each random surface as they would involve one set of smoothness parameters for all the smooths. If you want to allow for different degrees of wiggliness between individuals then you need to use the by variable alternatives.

Anyway, the first model would involved something like

y ~ t2(x, y, time) +
    t2(x, y, time, ind, c('cr', 'cr', 'cr', 're'), full = TRUE)

where the first smooth is the average time-varying surface and the second smooth is the random time-varying surfaces for individuals.

If you just want the random surfaces and don’t need the global smooth (i.e. you don’t plan on predicting for another individual) then the model simplifies to:

y ~ t2(x, y, time, ind, c('cr', 'cr', 'cr', 're'), full = TRUE)

If you think space is isotropic then you can use tensor interactions of a 2d thin plate or duchon spline, say, with a 1d time spline, for example:

t2(x, y, time, bs = c('ds', 'cr'), d = c(2, 1)) +
t2(x, y, time, ind, bs = c('ds', 'cr', 're'), d = c(2, 1, 1))

or (if you don’t want the global smooth) just

t2(x, y, time, ind, bs = c('ds', 'cr', 're'), d = c(2, 1, 1))

where the spatial term is a 2d Duchon spline as indicated by the argument d.

Throughout, you’ll want to be setting k to reasonable values for all the smooths and each marginal smooth of a tensor product.

4 Likes
3

This is awesome, thanks!

1 Like
4

Oh, @ucfagls any reason you suggest the “cr” basis? Way back when I used to be more familiar with mgcv I had the feeling that “ts” was better for generalization thanks to the shrinkage penalty, but even then I wasn’t confident in my expertise in this area.

5

They’re a lot more efficient when setting up the basis for the model matrix than thin plate splines and if you have a enough data for splines comprising three covariates you’ll likely notice the difference just in terms of generating the model matrix.

The default in these smooths is "cr" so when I needed to specify one basis to "re", I had to use something for the others and wanted to keep with the defaults.

A better option is to use select = TRUE in mgcv but I don’t know if that applies at all to brms models. The shrinkage splines are a good second alternative, but they do imply that the model should shrink the wiggly parts more than the linear parts (strictly the functions in the penalty null space) because of how the extra penalty is constructed.

Wouldn’t the priors on the terms have a similar effect to shrink the null space functions too?

2 Likes
6

Yes, good point.