如何为lmer()线性混合效应模型执行Breusch-Pagan异方差检验？

错误原因 #

lmtest包的bptest()函数仅适配传统线性模型（lm生成的S3类对象），而lmer()构建的线性混合效应模型属于S4类的merMod对象，函数内部依赖的$运算符无法直接访问S4对象的组件，因此触发该错误。

替代方法 #

方法1：使用performance包的check_heteroscedasticity()

performance包专门支持混合模型的诊断，内置异方差检验功能：

# 安装并加载包
install.packages("performance")
library(performance)

# 修正模型变量名（原代码中Mg_new应为数据中的Mn_new）
Mg_model <- lmer(Mn_new ~ SEX * S_M + (1|ID), data=data)

# 执行异方差检验
check_heteroscedasticity(Mg_model)

该函数会输出检验统计量、p值及直观的诊断结论。

方法2：使用DHARMa包的testDispersion()

DHARMa包专注于混合模型的残差诊断，通过模拟残差检测异方差：

# 安装并加载包
install.packages("DHARMa")
library(DHARMa)

# 生成模拟残差
sim_res <- simulateResiduals(Mg_model)

# 执行异方差检验
testDispersion(sim_res)

此方法基于模拟残差，更贴合混合模型的统计特性。

方法3：手动构造Breusch-Pagan检验

提取模型条件残差的平方，与原模型固定效应自变量拟合线性模型后，再执行bptest()：

# 提取残差平方和自变量数据
data$resid_sq <- residuals(Mg_model)^2

# 拟合残差平方与固定效应的模型
bp_model <- lm(resid_sq ~ SEX * S_M, data=data)

# 执行Breusch-Pagan检验
bptest(bp_model)

注意：此方法仅针对固定效应部分的异方差，未考虑随机效应影响，适合初步诊断。

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