I am struggling a bit in R. I want to make a multivariate linear mixed model with the variables slope_reactivity, peak_reactivity, and slope_recovery as dependent variables and normalized_CTseverity, gender, age_years and binary_diagnosis and dependent variables. Study is the random effect.

I got as far as this code:

model_cortisol <- lmer(slope_reactivity + slope_recovery + peak_reactivity ~ normalized_CTseverity + age_years + gender + binary_diagnosis + (1 | study), data = df_sum_cort)

However this only gives me the effect of normalized_CTseverity etc on slope_reactivity, slope_recovery and peak_reactivity together as a whole. I also want the effect of the independent variables on the dependent variables seperate from eachother.

How do I do this? Do I change something in the code of my model or can I do some sort of post hoc test? I don’t want to do seperate models as the slope_reactivity, slope_recovery and peak_reactivity are strongly linked.