Hyperparameters on group effects covariance matrix

Interfaces brms
,
1

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

set.seed(12)
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
2

Hi,
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!

3

Thanks! I actually have multiple observations in each repetition (experimental conditions) and a categorical variable that codes the condition, so just using ar will not work. However, I am going to look at a two-step approach where I sample from

bf(y ~ 0 + rep + rep:cond_code + (0 + rep || id) + (0 + rep:cond_code || id))

and model the random effects using ar like

bf(RE | se(RE_sd) ~ ar(time = rep_time, gr = id))

where RE is each id's posterior mean for rep:cond_code and RE_sd is each id's posterior standard deviation. This two-step process seems a bit hackish though.

1 Like