如何在R中对mgcv包拟合的非线性回归执行Bootstrap？

当然可以对mgcv生成的广义可加模型（GAM）执行Bootstrap操作，针对你遇到的问题，以下是具体的实现方案：

1. 基础统计量的Bootstrap（以伪R²为例） #

你之前用boot包时的核心问题是统计量函数里误拟合了线性模型（lm），应该替换为gam。以下是修正后的代码：

# 加载必需的包
library(mgcv)
library(boot)

# 自定义统计量函数：计算GAM的伪R²（解释方差比例）
rsq_gam <- function(data, indices) {
  # 按Bootstrap抽样的索引提取数据
  d <- data[indices, ]
  # 拟合与原模型结构一致的GAM
  fit <- gam(y ~ s(x1, bs = "tp") + s(x2, bs = "tp") + s(x3, bs = "tp"), data = d)
  # 计算伪R²：1 - (模型偏差 / 空模型偏差)
  return(1 - deviance(fit) / deviance(gam(y ~ 1, data = d)))
}

# 设置随机种子保证结果可重复
set.seed(123)
# 执行Bootstrap，R为抽样次数
boot_results <- boot(data = your_dataset, statistic = rsq_gam, R = 1000)

# 查看Bootstrap结果
print(boot_results)
# 绘制统计量的分布直方图
plot(boot_results)
# 提取95% BCA置信区间（偏态数据下更稳健）
boot.ci(boot_results, type = "bca")

2. 平滑项的Bootstrap p值生成 #

如果你目标是给每个平滑项生成稳健的p值（应对残差偏态），有两种高效的方法：

方法1：手动用boot包实现Case Bootstrap #

通过重复抽样观测值，生成平滑项统计量的分布，再与原模型统计量对比计算p值：

# 自定义函数：提取每个平滑项的F统计量
smooth_f_stats <- function(data, indices) {
  d <- data[indices, ]
  fit <- gam(y ~ s(x1, bs = "tp") + s(x2, bs = "tp") + s(x3, bs = "tp"), data = d)
  # 从模型摘要中提取平滑项的F值
  return(summary(fit)$s.table[, "F"])
}

# 执行Bootstrap
set.seed(456)
boot_smooth <- boot(data = your_dataset, statistic = smooth_f_stats, R = 2000)

# 提取原模型的平滑项F值
orig_f_values <- summary(nonlinear)$s.table[, "F"]

# 计算每个平滑项的Bootstrap p值：统计量大于原统计量的比例
bootstrap_pvals <- sapply(1:length(orig_f_values), function(i) {
  mean(boot_smooth$t[, i] >= orig_f_values[i])
})

# 给p值命名，对应每个平滑项
names(bootstrap_pvals) <- rownames(summary(nonlinear)$s.table)
print(bootstrap_pvals)

方法2：用mgcv内置的置换检验 #

mgcv自带permutation.test.gam()函数，专门针对GAM平滑项做置换检验，无需手动写Bootstrap逻辑，更适合残差偏态的场景：

# 对已拟合的nonlinear模型执行置换检验
perm_test_results <- permutation.test.gam(nonlinear, nsim = 2000)
print(perm_test_results)

这个函数会直接输出每个平滑项的置换p值，完全规避了正态假设的限制，计算效率也更高。

3. 关于boot.pval包的问题 #

boot_summary()依赖car包的Boot函数，但car::Boot并未实现针对mgcv GAM模型的残差Bootstrap支持，因此无法直接使用，只能通过上述手动实现或mgcv内置工具解决。

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