Hi,
Let’s say I have 10 correlated measurements of some quantity in 300 participants, and I want to identify characteristics associated with high values. I could run 10 independent models, or use a multivariate model. But I thought of another trick (not sure if it’s common): put the data in “long” format, and model the quantity of any measurement, with random intercepts and random slopes by type of measurement.
One advantage is that I only use one model, with the random intercept becoming the mean value of each type of measurement. Another advantage is that, with the random slope, I am able to estimate an “average” association of each of the covariates with all 10 measurements, while at the same time being able to look at specific associations. I find this average effect very useful to communicate complex results.
I made some simulations and it seems to work exactly as intended:
# settings
n_patients = 300
n_groups = 10
set.seed(123)
mu = runif(n=n_groups,min=1,max=20) # mean value per group
sigma = runif(n=n_groups,min=0.5,max=2) # std dev per group
p_female = .5 # proportion of females
beta_female = rnorm(n=n_groups,mean=2,sd=2) # beta of females per group
rho = randcorr(p = n_groups)
Sigma = diag(sigma) %*% rho %*% diag(sigma)
# simulate
set.seed(123)
betas = tibble(group=paste0("V",1:n_groups),beta_female=beta_female)
data_long =
MASS::mvrnorm(n_patients, mu, Sigma) %>%
as_tibble() %>%
mutate(female=rbinom(n_patients,1,p_female)) %>%
pivot_longer(1:n_groups,names_to = "group") %>%
left_join(betas,by="group") %>%
mutate(value=value+female*beta_female)
# model
modB4 = brm(value ~ -1 + female + (1 + female|group), data=data_long)
summary(modB4)
ranef(modB4)
I guess I could do something similar with a multivariate normal model, but the thing is that my real application is a bit more complex with data imputation and a hurdle-gamma model, and I don’t think it would work as well.
The only thing that I could not find is how to extract the variance-covariance matrix for the random intercept (so basically Sigma
or rho
in my code). I found a similar question here, but no answer either. Is there a way to do that? I would like to avoid recoding everything in stan.
Thanks in advance!