求可计算线性混合效应偏残差的R函数及termplot报错解决方法

你遇到的报错是因为termplot在处理lmer生成的merMod类模型时，内部调用的predict.merMod函数不支持type="terms"参数，导致无法生成所需的拟合项来计算偏残差。下面提供三种实用的解决方法：

方法一：手动计算偏残差，自定义绘图 #

偏残差的核心逻辑是原始残差 + 该自变量对应的固定效应拟合值，我们可以手动拆分模型结果来计算，再用基础绘图或ggplot2实现可视化：

library(lme4)
library(ggplot2)
library(gridExtra)

# 拟合混合效应模型
lmer.cars <- lmer(mpg ~ disp + drat + qsec + (1|cyl), data=mtcars)

# 提取模型的固定效应设计矩阵（去除截距项）
X <- model.matrix(lmer.cars)[, -1]
# 提取对应自变量的固定效应系数
fixef_coefs <- fixef(lmer.cars)[-1]
# 计算每个自变量的固定效应拟合值
term_fits <- X %*% fixef_coefs

# 提取模型残差
model_resids <- residuals(lmer.cars)

# 计算每个自变量的偏残差
partial_resids <- as.data.frame(model_resids + term_fits)
colnames(partial_resids) <- colnames(X)
# 保留原数据的自变量值
partial_resids <- cbind(partial_resids, mtcars[, colnames(X)])

# 用ggplot2批量生成偏残差图
plot_list <- lapply(colnames(X), function(var) {
  ggplot(partial_resids, aes(x = .data[[var]], y = .data[[var]])) +
    geom_point(color = "purple") +
    geom_smooth(method = "loess", se = FALSE, color = "black") +
    labs(x = var, y = "偏残差") +
    theme_bw()
})

# 2x2排列图形
grid.arrange(grobs = plot_list, ncol = 2)

方法二：用sjPlot包快速生成 #

sjPlot包对混合效应模型的可视化支持很友好，plot_model函数可以直接生成偏残差图，无需手动计算：

library(sjPlot)
library(lme4)

lmer.cars <- lmer(mpg ~ disp + drat + qsec + (1|cyl), data=mtcars)

# 绘制偏残差图，展示原始数据点
plot_model(lmer.cars, type = "resid", show.data = TRUE, colors = "purple") +
  theme_bw()

方法三：用effects包适配混合模型 #

你之前尝试的effects包其实支持merMod模型，只需通过Effect函数提取单变量效应后，结合残差绘制偏残差图：

library(effects)
library(lme4)

lmer.cars <- lmer(mpg ~ disp + drat + qsec + (1|cyl), data=mtcars)

# 设置2x2绘图布局
par(mfrow = c(2,2))

# 遍历每个自变量绘制偏残差图
for (var in c("disp", "drat", "qsec")) {
  # 提取该变量的效应对象
  eff_obj <- Effect(var, lmer.cars)
  # 计算偏残差：残差 + 该变量的拟合效应（去除截距）
  partial_resid <- residuals(lmer.cars) + (eff_obj$fitted - fixef(lmer.cars)[1])
  # 绘制效应曲线，再叠加偏残差点
  plot(eff_obj, main = var, ylab = "偏残差")
  points(mtcars[[var]], partial_resid, col = "purple")
}

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