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:

- what we called a global smooth plus random smooths model, or
- 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.