如何为lme4线性混合效应模型计算稳健标准误？

步骤1：安装并加载所需包 #

install.packages(c("lme4", "merDeriv", "clubSandwich"))
library(lme4)
library(merDeriv)
library(clubSandwich)

步骤2：确认原模型拟合 #

先确保你已经正确拟合了线性混合效应模型：

# 你的原模型代码
TD_PE_analysis1 <- lmer(PE ~ learn_prof * trial * condition + (1|subject) + (1|image), data=dat)

步骤3：提取模型得分矩阵（merDeriv） #

merDeriv的get_scores()函数会提取模型的得分信息，这是计算稳健标准误的核心输入：

model_scores <- get_scores(TD_PE_analysis1)

步骤4：计算稳健方差-协方差矩阵（clubSandwich） #

用clubSandwich的vcovCR()函数生成稳健VC矩阵。针对你的研究设计，按受试者（subject）聚类是最合理的选择（同一受试者的60轮观测存在非独立性）：

# 基于受试者聚类的CR2型稳健标准误
vcov_robust <- vcovCR(TD_PE_analysis1, cluster = dat$subject, type = "CR2")

• type = "CR2"：小样本下表现最优的稳健标准误类型，推荐使用；若需其他类型，可替换为"CR1"或"CR0"。
• 若需考虑image的非独立性，可尝试cluster = interaction(dat$subject, dat$image)，但优先推荐受试者层面聚类。

步骤5：输出带稳健标准误的模型结果 #

用coef_test()函数输出系数的统计检验结果，包含稳健标准误、z/t值和p值：

# 用z检验输出结果
coef_test(TD_PE_analysis1, vcov = vcov_robust, test = "z")

# 若偏好t检验（基于聚类自由度），可改为：
# coef_test(TD_PE_analysis1, vcov = vcov_robust, test = "t")

关键说明 #

• 你通过Fligner检验发现的方差不齐问题，稳健标准误能有效缓解异方差对统计推断的干扰，这个解决方案适配你的问题。
• 随机效应部分仍沿用原模型的设定，稳健标准误仅修正固定效应的推断结果。

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