I am attempting to use the R package “brms” to evaluate a bivariate mixed-model regression model (specifically, looking at continuous outcomes for both left and right eyes over time). In applying your model to my dataset, I had several questions:
I’m unsure how the correlation/dimensional dependency between the two dimensions are handled?
Can the covariance structure for both within each dimension and between the two dimensions be modified or specified by the user? For example, following a traditional repeated measures linear mixed model, we can specify a covariance structure such as compound symmetry or auto-regressive. However, when we extend this to a bivariate outcome, this structure essentially comes a cube (if I’m understanding correctly). What variance/covariance structure does your modeling assume and can the user modify this cube?
Lastly, when I model my data with your package, the model summary provides different fixed effect coefficient values for each dimension (left eye and the right eye for my data), such as the impact of baseline age is different for right eye outcomes than left eye outcomes). Is there a way to use this package to obtain global coefficient values for these fixed effects, rather than dimension-specific values? I have included my code below, if that is helpful, where DeltaSEQOD and DeltaSEQOS are continuous measures, sex, ethnicity, and race are all categorical variables, and age and the other two baseline measures are continuous. We have 5 follow-up visits so currently treating visit as continuous with a random intercept and slope for visits.
Hi @morgan.
I am not an expert on this but have a little experience so I’ll try my best to answer your questions:
I think when you run your code you likely see a warning message along the lines of: “Setting ‘rescor’ to TRUE by default for this model” - yes? To explain this a bit - there is a brmsformula option called set_rescor (rescor → residual correlations) that you did not specify in your code, but is being added by default. So in effect your formula would be: formula = formulaOD + formulaOS + set_rescor(TRUE) If you look at the helpfile ( ?set_rescor() ) you see the following explanation:
also should you want you can turn it off by using + set_rescor(FALSE). (Note: it is good practice to expliciilty include + set_rescor() since there is a warning notice that the default behaviour will change to FALSE in future.)
I don’t fully understand your question here in that I don’t know what you mean by a cube structure. However this might be helpful → you can modify the residual correlation structure to include correlation between given levels of a grouping variable across outcomes using the notation: (Visit | p | PtID)) in both your outcomes. “p” is just a label - you can use any name you like. The effect is the brms will model group level effects as correlated where the same ( | p | ) formula structure appears. You can find some explanation of this here: Estimating Multivariate Models with brms and more in the help file help("brmsformula") if you scroll down to “Group-level terms”.
I don’t know if you can obtain “global coefficient values” rather than dimension specific ones with brms. However, it is worth asking the question: do you really want these? For example if your outcome variables DeltaSEQOD and DeltaSEQOS happen to have really different scales how would you interpret a global coefficient across both? Or even if they are on the same scale, perhaps the interdependence of the covariates differ between outcomes - if you use global coefficients you will be averaging over such subtleties.
Finally if you wish to understand how the correlation structures are modelled in detail I suggest to take a look at the generated stancode from your models. For example you could generate code as follows and check the differences between different model specifications:
library(brms)
library(palmerpenguins)
# fit some GP models
bform1 <- bf(mvbind(bill_length_mm, flipper_length_mm) ~ sex +
(1|species))
bform2 <- bf(mvbind(bill_length_mm, flipper_length_mm) ~ sex +
(1|p|species))
# view the Stan code generated for different models
stancode(bform1 + set_rescor(TRUE), data = penguins)
stancode(bform1 + set_rescor(FALSE), data = penguins)
stancode(bform2 + set_rescor(TRUE), data = penguins)
# Run the models if you want
fit1 <- brm(bform1, data = penguins, chains = 2, cores = 2)
fit1a <- brm(bform1 + set_rescor(FALSE), data = penguins, chains = 2, cores = 2)
fit2 <- brm(bform2 + set_rescor(TRUE), data = penguins, chains = 2, cores = 2)
