Hi all / @paul.buerkner ,

I wonder whether the brms Gaussian Process (GP) support through `gp`

actually supports hierarchical GPs when setting `gr=T`

.

What I mean by hierarchical GP:

Let \vec{M} and \vec{0} be m-dimensional.

So we have

where \mathbf{K}_1 is m\times m dimensional. This is what I would call the population mean GP. Now for each of the N individuals we have that the observation \vec{y}_i (also m-dimensional) of individual i is **i.i.d.** from the following GP:

where \mathbf{K}_2 is another m\times m kernel. Now suppose each individual has observations for the same m-dimensional vector of time points, couldnâ€™t one then do in brms

```
fit <- brms(y ~ gp(time,gr=T) + gp(time,gr=F), data=df)
```

and the term `gp(time,gr=T)`

would correspond to the population GP with kernel \mathbf{K}_1 and `gp(time)`

to the individual GP with kernel \mathbf{K}_2.