For this example, I am using the data "appendix_example1_wide.SUPP.FINAL.csv" posted here.
In the paper, the authors use the to MCMCglmm package fit a multivariate multilevel model. Of particular interest, they use the following code to estimate covariances between random effects (e.g., random intercepts from each outcome):
prior1 <- list(R = list(V = .2*diag(2), nu = 3), G = list(G1 = list(V = diag(c(.2,.2,.1,.1)), nu = 5))) mult.mcmc <- MCMCglmm( cbind(y1, y2) ~ -1 + trait + trait:tx + trait:time + trait:time:tx, data = data1, random = ~ us(-1 + trait + trait:time):pid, rcov = ~ us(trait):units, family = c("gaussian","gaussian"), nitt = 25000, burnin = 5000, prior = prior1 ) summary(mult.mcmc)
Is it possible to specify the same model using the brms package? Here's my attempt (that did not even closely reproduce the estimates from the paper):
mod1 <- brm(cbind(y1,y2) ~ 1 + time*tx + (1 + time|pid),data=data)