如何在mgcv::gam的predict.gam()预测中纳入残差方差？

是的，你可以手动将观测水平的残差方差纳入predict.gam()的预测结果中，核心思路是在模型给出的期望值标准误基础上，叠加残差的不确定性，具体步骤如下：

1. 提取模型的残差方差 #

从拟合好的GAM模型对象中，使用sigma(model)可以直接获取残差的标准差（对应你提到的观测水平残差方差的平方根）。

2. 获取期望值的预测结果 #

调用predict.gam()并设置se.fit=TRUE，得到期望值的预测值和其标准误（这部分已经包含了固定系数和分组因子方差参数的不确定性，等价于brms::fitted()）：

pred_exp <- predict(model, newdata = your_new_data, se.fit = TRUE)

3. 计算包含残差方差的总标准误 #

总标准误是期望值标准误的平方与残差标准差的平方之和的平方根：

total_se <- sqrt(pred_exp$se.fit^2 + sigma(model)^2)

4. 构建观测值的预测区间 #

基于总标准误计算预测区间（以95%置信区间为例，使用1.96作为临界值）：

pred_lower <- pred_exp$fit - 1.96 * total_se
pred_upper <- pred_exp$fit + 1.96 * total_se

完整代码示例 #

library(mgcv)
set.seed(123)

# 生成模拟数据并拟合GAM
dat <- gamSim(1, n = 400, dist = "normal", scale = 2)
model <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat)

# 定义新的预测数据
new_dat <- data.frame(x0 = 0.5, x1 = 0.5, x2 = 0.5, x3 = 0.5)

# 获取期望值的预测结果
pred_exp <- predict(model, newdata = new_dat, se.fit = TRUE)

# 提取残差标准差
resid_sd <- sigma(model)

# 计算总标准误和观测值预测区间
total_se <- sqrt(pred_exp$se.fit^2 + resid_sd^2)
obs_pred_lower <- pred_exp$fit - 1.96 * total_se
obs_pred_upper <- pred_exp$fit + 1.96 * total_se

# 输出对比结果
cat("【期望值预测（对应brms::fitted()）】\n")
cat("预测值:", round(pred_exp$fit, 3), "\n")
cat("95%区间:", round(pred_exp$fit - 1.96*pred_exp$se.fit, 3), "~", round(pred_exp$fit + 1.96*pred_exp$se.fit, 3), "\n\n")

cat("【观测值预测（对应brms::predict()）】\n")
cat("预测值:", round(pred_exp$fit, 3), "\n")
cat("95%区间:", round(obs_pred_lower, 3), "~", round(obs_pred_upper, 3), "\n")

注意事项 #

• 如果是广义线性GAM（如泊松、Logistic分布），残差方差并非固定值，而是与预测的期望值相关：
  • 泊松分布：方差 = 期望值，此时总标准误为sqrt(pred_exp$se.fit^2 + pred_exp$fit)
  • 二项分布：方差 = 期望值 × (1 - 期望值)，总标准误为sqrt(pred_exp$se.fit^2 + pred_exp$fit*(1-pred_exp$fit))
• 若模型包含随机效应（如gam(..., random = ~1|group)），predict.gam()的se.fit已经包含了分组因子方差参数的不确定性，无需额外处理。

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