Multivariate model correlations between responses and categorical variables

Modeling
1

I’m working on a gaussian multivariate model with four response variables and one categorical variable as a explanatory variable. My interest here are the correlations between the response variables per each level of the Treatment variable. For example:
cor1: Y1-Trt1 X Y2-Trt1,
cor2: Y1-Trt2 X Y2-Trt2

I use a grouping structure of gr(site:plot, by = “Treatment”, id = “V”, cor = TRUE). So a nested random effect correlated between the response variables with a separate variance-covariance matrix for each level (three levels) of Treatment variable.

f.1 <- bf(y1 ~ 1 + Trt + gr(site:plot, by = "Trt", id = "V", cor = TRUE))
f.2 <- bf(y2 ~ 1 + Trt + gr(site:plot, by = "Trt", id = "V", cor = TRUE))
f.3 <- bf(y3 ~ 1 + Trt + gr(site:plot, by = "Trt", id = "V", cor = TRUE))
f.4 <- bf(y4 ~ 1 + Trt + gr(site:plot, by = "Trt", id = "V", cor = TRUE))

What I get from this definition is:

cor: y1_Intercept:"Trt"Trt,y2_Intercept:"Trt"Trt)
cor: y1_Intercept:"Trt"Trt,y3_Intercept:"Trt"Trt)
cor: y2_Intercept:"Trt"Trt,y3_Intercept:"Trt"Trt)
cor: y1_Intercept:"Trt"Trt,y4_Intercept:"Trt"Trt)
cor: y2_Intercept:"Trt"Trt,y4_Intercept:"Trt"Trt)
cor: y3_Intercept:"Trt"Trt,y4_Intercept:"Trt"Trt)

The estimates for each of these correlations correspond to grouping

gr(site:plot, id = "V", cor = TRUE)).

So this is probably not the way to get these correlations i’m looking for…? What would be the proper way to do this?

2

There’s a large community of brms users here, but in case you don’t get a response soon enough, you could try explaining the data/variables and what kind of relationship you are trying to establish between them – you already did that a bit in the first sentence, but it’s very brief and a bit cryptic.

If you can, you could also describe the model you are using in more general terms (e.g. through mathematical notation). That could help other people (like me) who may not know what bf or gr mean.

3

by should not be character vector, I think you simply want

f.1 <- bf(y1 ~ 1 + Trt + gr(site:plot, by = Trt, id = "V", cor = TRUE))

I think you are effectively fitting “separate” VCVs for the single level "Trt" instead of the multiple levels of Trt. See the example in ?gr.