如何调整多变量线性混合模型以获取自变量对各因变量的单独效应？

问题背景 #

你尝试用lmer()构建多变量线性混合模型，将slope_reactivity、peak_reactivity、slope_recovery作为因变量，normalized_CTseverity、gender、age_years、binary_diagnosis作为自变量，study作为随机效应，但原代码仅输出自变量对三个因变量的整体效应，希望在不拆分模型的前提下获取每个因变量的单独效应。

原代码问题说明 #

你当前的写法slope_reactivity + slope_recovery + peak_reactivity ~ ...并不是真正的多变量模型，lmer()会将三个因变量直接相加作为一个复合因变量，因此只能得到整体效应，这是错误的用法。

解决方案 #

以下两种方法都能保留因变量间的关联性，同时输出单个因变量的效应：

方法1：转换为长数据格式建模（推荐，基于基础lme4）

把宽格式数据转成“长格式”，将三个因变量合并为一列分类变量，再通过交互项建模：

1. 转换数据格式：

library(tidyr)
# 把三个因变量转成长数据的一列
df_long <- df_sum_cort %>%
  pivot_longer(
    cols = c(slope_reactivity, slope_recovery, peak_reactivity),
    names_to = "outcome_type",  # 存储因变量名称的列
    values_to = "outcome_value" # 存储因变量数值的列
  )

2. 构建包含交互项的混合模型：

library(lme4)
model_cortisol <- lmer(
  outcome_value ~ outcome_type * (normalized_CTseverity + age_years + gender + binary_diagnosis) + (1 | study),
  data = df_long
)

3. 提取单个因变量的效应：
   用summary()直接查看系数，或者用emmeans包更直观地输出每个自变量对不同因变量的效应：

library(emmeans)
# 查看normalized_CTseverity对每个因变量的效应
emmeans(model_cortisol, ~ normalized_CTseverity | outcome_type)
# 同理替换其他自变量，比如gender
emmeans(model_cortisol, ~ gender | outcome_type)

方法2：使用专门的多变量混合模型包

用mvmer包（基于lme4的扩展）直接构建多变量模型，无需转换数据格式：

library(mvmer)
# 为每个因变量单独指定公式，mvmer会自动考虑残差相关性
model_cortisol <- mvmer(
  list(
    slope_reactivity ~ normalized_CTseverity + age_years + gender + binary_diagnosis + (1 | study),
    slope_recovery ~ normalized_CTseverity + age_years + gender + binary_diagnosis + (1 | study),
    peak_reactivity ~ normalized_CTseverity + age_years + gender + binary_diagnosis + (1 | study)
  ),
  data = df_sum_cort
)
# 查看结果，每个自变量对三个因变量的系数会分别输出
summary(model_cortisol)

注意事项 #

两种方法都能保留因变量间的关联性，避免拆分模型带来的I类错误膨胀。方法1更易上手，依赖的包都是常用工具；方法2是更严谨的多变量模型实现，适合需要直接查看多变量残差相关性的场景。

内容的提问来源于stack exchange，提问作者Roos Eijgenraam