Hyperparameters on group effects covariance matrix

I would like to be able to constrain portions of the covariance matrix of group effects using a hyperparameter. I don’t think that the solution to this problem is the gr function because I’m not trying to represent the covariance structure among the group units, but rather trying to create a custom group-effects covariance matrix sort of like pdClass does for lem4 models.

The specific covariance structure I’m interested in has to do with the time-ordered nature of one of the random effects.

If the data are generated

Sigma <- matrix(c(
  1,      .5,    .25,  .125,  .0625,
  .5,    1,      .5,   .25,   .125,
  .25,    .5,   1,     .5,    .25,
  .125,   .25,   .5,  1,      .5,
  .0625,  .125,  .25,  .5,   1), ncol = 5)

rep_means <- MASS::mvrnorm(1000, mu = rep(0, 5), Sigma = Sigma)

rep_means_contrast <- cbind(rep_means[, 1], apply(rep_means[, 2:5], 2, function(x){
  x - rep_means[, 1]

d <- do.call(rbind, lapply(1:dim(rep_means)[[1]], function(id){
  arow <- rep_means[id, ]
  y1 <- rnorm(length(arow), mean = arow, sd = 1)
  y2 <- rnorm(length(arow), mean = arow, sd = 1)
  return(data.frame(id = id, 
                    rep = as.factor(rep(1:length(arow), 2)), 
                    y = c(y1, y2)))

and if the model is

f <- bf( y ~ 0 + rep + (0 + rep | id) )

I’d like to capture a pooled estimate of the off-diagonals of the group effect of rep using a hyperparameter, as well as provide some prior information that might help speed up sampling. Or in other words, I’d like the model to reflect the possibility that the data generating process is essentially AR(1).

Is there any way to do this in brms? I could imagine altering the model code produced to create a transformed parameter that is the covariance matrix and placing hyperparameters on the cells of that matrix, but it seems like there might be a better way.

1 Like

so I don’t think that there is currently support for arbitrary structures in the covariance matrix (e.g. the way generic3 works in R-INLA). However, AR(1) (and a bunch of other structures) is supported with special syntax, so you might be able to model something very close to this with ar (Set up AR(p) correlation structures — ar • brms)

Maybe something like

y ~ ar(time = rep, gr = id)

would work?

Best of luck with the model!