R语言mgcv包GAM模型参数项绘图原始数据提取方法问询

mgcv自带的plot.gam()函数返回的列表仅包含平滑项的绘图数据，分类变量这类参数项的对应数据可通过以下两种方式获取：

方法1：使用termplot()直接提取所有项数据 #

termplot()是R内置的回归项绘图工具，支持gam类模型，关闭绘图开关即可返回所有项（含参数项、平滑项）的完整绘图数据：

# 提取所有项的绘图数据，不直接绘图
term_data <- termplot(mod, se = TRUE, plot = FALSE, data = data)

# 提取分类参数项pred2的绘图数据
param_pred2 <- term_data$pred2
print(param_pred2)

输出内容包含分类水平对应的x轴坐标、拟合值、上下95%置信区间边界，和mgcv默认绘图的参数项结果完全匹配。

方法2：手动计算（自定义自由度更高） #

如果需要更灵活的控制（比如调整置信区间水平、控制其他协变量的取值），可以直接用模型预测的方式生成参数项的绘图数据，和默认绘图逻辑完全一致：

# 构造预测用的新数据集：参数项取所有水平，其余协变量取均值/参考水平
newdata <- expand.grid(
  pred2 = levels(data$pred2),
  pred1 = mean(data$pred1),
  pred3 = mean(data$pred3)
)

# 计算拟合值和标准误
pred_res <- predict(mod, newdata = newdata, se.fit = TRUE, type = "terms")

# 提取参数项pred2的拟合值、置信区间
param_pred2_manual <- data.frame(
  pred2_level = newdata$pred2,
  fit = pred_res$fit[, "pred2"],
  se = pred_res$se.fit[, "pred2"],
  lower_ci = pred_res$fit[, "pred2"] - 1.96 * pred_res$se.fit[, "pred2"],
  upper_ci = pred_res$fit[, "pred2"] + 1.96 * pred_res$se.fit[, "pred2"]
)
print(param_pred2_manual)

# 如果需要提取参数项对应的偏残差（和plot(mod, residuals=T)里的散点对应）
partial_resid <- residuals(mod, type = "partial")[, "pred2"]

该方法生成的数据可以直接用于ggplot2等工具自定义绘图，可根据需求调整置信区间的显著性水平、其余协变量的固定取值。

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