如何使用lmer实现多因变量重复测量分析以替代MANOVA？

数据集示例 #

df <- data.frame(
id  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 
       34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 
       19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
       40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 
       29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65), 
collection_point = c(rep(c("Baseline", "Immediate", "3M"), each=28)), 
intervention = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C", 
              "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)), 
scale_A = c(6.5, 7.0, 6.25, 6.0, NA, 7.5, 7.5, 
        8.0, 7.5, 6.75, 7.5, 6.75, 6.75, 6.5, 
        5.75, 6.75, 7.75, 7.5, 7.75, 7.25, 7.75, 
        7.25, 7.25, 5.75, 6.75, NA, 6.75, 7.5, 
        6.75, 7.0, 6.5, 7.0, 7.5, 7.5, 7.5, 
        7.75, 7.25, 7.25, 7.25, 7.5, 6.5, 6.25, 
        6.25, 7.25, 7.5, 6.75, 7.25, 7.25, 7.5, 
        7.25, 7.5, 7.25, NA, 7.0, 7.5, 7.5, 
        6.75, 7.25, 6.5, 7.0, 7.5, 7.5, 7.5, 
        7.75, 7.5, 7.5, 7.5, 7.5, 6.5, 5.75, 
        6.25, 6.75, 7.5, 7.25, 7.25, 7.5, 7.75, 
        7.75, 7.75, 7.5, NA, NA, NA, NA))
scale_B = c(5.0, 6.5, 6.25, 7.0, NA, 5.5, 6.5, 
        6.0, 7.5, 5.75, 6.5, 5.75, 7.75, 6.5, 
        6.75, 7.75, 7.75, 7.5, 7.75, 5.25, 7.75, 
        6.25, 6.25, 6.75, 5.75, NA, 6.75, 6.5, 
        7.75, 6.0, 7.5, 6.0, 7.5, 7.5, 6.5, 
        6.75, 6.25, 6.25, 6.25, 6.5, 6.5, 7.25, 
        7.25, 6.25, 6.5, 7.75, 6.25, 7.25, 6.5, 
        6.25, 6.5, 6.25, NA, 7.0, 6.5, 7.5, 
        7.75, 6.25, 7.5, 6.0, 7.5, 6.5, 6.5, 
        6.75, 6.5, 6.5, 6.5, 7.5, 7.5, 6.75, 
        7.25, 7.75, 6.5, 6.25, 7.25, 6.5, 6.75, 
        6.75, 6.75, 6.5, 5.5, NA, NA, 6.5))
scale_C = c(5.5, 5.0, 7.25, 7.0, 8.0, 5.5, 5.5, 
        8.0, 5.5, 7.75, 5.5, 7.75, 7.75, 7.5, 
        7.75, 7.75, 5.75, 5.5, 5.75, 5.25, 5.75, 
        5.25, 6.25, 7.75, 7.75, NA, 7.75, 5.5, 
        6.75, 6.0, 7.5, 5.0, 5.5, 5.5, 7.5, 
        5.75, 6.25, 5.25, 5.25, 5.5, 7.5, 7.25, 
        7.25, 6.25, 5.5, 7.75, 5.25, 5.25, 7.5, 
        5.25, 6.5, 5.25, 5.0, 5.0, 5.5, 5.5, 
        7.75, 6.25, 7.5, 5.0, 5.5, 5.5, 7.5, 
        5.75, 6.5, 5.5, 5.5, 5.5, 7.5, 7.75, 
        7.25, 7.75, 5.5, 5.25, 5.25, 5.5, 6.75, 
        5.75, 5.75, 5.5, 6.75, NA, 5.75, NA))

字段说明 #

• id：受试者编号
• collection_point：受试者数据采集时间点（重复测量指标）
• intervention：受试者随机分配的干预组（固定效应）
• scale_A/scale_B/scale_C：受试者在各数据采集点填写的不同维度问卷得分（结局指标）

需求背景 #

研究设计为受试者随机分入3个干预组，在3个不同时间点完成A-C共3个同主题不同维度的量表，用于评估干预效果随时间的改善情况。
已实现单结局的lmer拟合：

mixed.lmer.A1<-lmer(scale_A~intervention+(collection_point|id), control =
 lmerControl(check.nobs.vs.nRE = "ignore"), data = df)

尝试运行多因变量lmer代码报错：

mixed.lmer.comb<-lmer(cbind(scale_A, scale_B, scale_C)~intervention+
(collection_point|id), control = lmerControl(check.nobs.vs.nRE = "ignore"), 
data = df)

使用lm运行多因变量分析无法控制collection_point重复测量效应，结果无实际意义，需要实现基于混合效应框架的多因变量拟合分析。

lme4包的lmer()函数原生不支持cbind()格式的多响应变量输入，无法直接运行你写的多因变量代码，可采用以下两种方案实现需求：

方案1：长表重构法（无需额外依赖包，通用性最强） #

将宽格式数据集转换为长格式，新增「量表类型」分类变量，把多结局转化为单结局+分层项的形式，即可用lmer直接拟合，同时控制重复测量效应：

library(tidyr)
library(lme4)
# 转换为长表
df_long <- pivot_longer(df, 
                        cols = starts_with("scale_"),
                        names_to = "scale_type",
                        values_to = "score",
                        values_drop_na = TRUE)
# 拟合多结局混合效应模型
mixed.lmer.comb <- lmer(score ~ intervention * scale_type * collection_point +
                          (1 | id) + (collection_point | id/scale_type),
                        control = lmerControl(check.nobs.vs.nRE = "ignore"),
                        data = df_long)

模型说明：

• 固定效应部分的三向交互项可直接检验不同干预组、不同时间点、不同维度量表的得分差异，完全匹配研究需求
• 随机效应部分控制了受试者个体异质性、同一受试者不同时间点的关联、同一受试者不同量表得分的关联
• 拟合完成后可通过car::Anova(mixed.lmer.comb, type = 3)得到各因子的显著性检验结果，等价于调整了重复测量效应的MANOVA输出

方案2：专用多变量检验方案 #

如果需要输出传统MANOVA格式的多变量检验结果，可使用nlme包的gls()函数搭配car包的Manova()函数实现，该方案对数据缺失值敏感度较高，需提前完成缺失值插补：

library(nlme)
library(car)
# 拟合多变量广义最小二乘模型，控制重复测量相关结构
model_gls <- gls(score ~ intervention * scale_type * collection_point,
                 correlation = corSymm(form = ~ 1 | id/collection_point),
                 data = df_long,
                 na.action = na.omit)
# 输出多变量检验结果
Manova(model_gls, test.statistic = "Pillai")

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