问询GAM模型中两组张量积曲面差异区域的统计检验方法

核心思路：聚类置换检验是合理选择 #

你提到的未校正置信区间确实容易因多重比较导致假阳性，而聚类置换检验正是针对这类二维/空间曲面差异检验的最优方法之一——它考虑了网格点之间的空间相关性，通过置换分组标签生成零分布，能有效控制整体Ⅰ类错误率，比单纯的p值校正更适合定位连续的差异区域。

具体实施步骤（R环境） #

1. 提取原始差异曲面 #

首先用mgcv的predict()生成两组的拟合曲面，计算差值：

# 构造a和b的网格数据
newdata <- expand.grid(
  a = seq(min(dat$a), max(dat$a), length.out = 100),
  b = seq(min(dat$b), max(dat$b), length.out = 100),
  factor = unique(dat$factor),
  # 协变量设为均值或典型值
  covariate = mean(dat$covariate),
  id = NA  # 设为NA会预测边际效应（积分掉随机效应）
)
# 假设你已拟合好模型gam_obj
# gam_obj <- gam(y ~ te(a,b,by=factor) + covariate + s(id,bs="re"), data=dat, method="REML")
# 预测两组的拟合值及标准误
pred <- predict(gam_obj, newdata = newdata, type = "response", se.fit = TRUE)
newdata$fit <- pred$fit
newdata$se <- pred$se.fit

# 拆分两组并计算差值
pred_group1 <- subset(newdata, factor == levels(dat$factor)[1])
pred_group2 <- subset(newdata, factor == levels(dat$factor)[2])
diff_surface <- data.frame(
  a = pred_group1$a,
  b = pred_group1$b,
  diff = pred_group2$fit - pred_group1$fit,
  diff_se = sqrt(pred_group1$se^2 + pred_group2$se^2)
)
# 计算每个网格点的t统计量
diff_surface$t_val <- diff_surface$diff / diff_surface$diff_se

2. 设计置换方案（适配随机效应） #

因为模型包含受试者ID随机效应，置换必须保持受试者内的分组结构：

• 如果factor是组间变量（每个受试者只属于一组）：按受试者为单位置换分组标签，用permute::shuffleSet生成块置换索引：

  library(permute)
  set.seed(123)
  # nperm为置换次数，blocks按受试者ID分组
  perm_indices <- shuffleSet(n = length(unique(dat$id)), nperm = 1000, replace = FALSE)
  

• 如果factor是组内变量（每个受试者有两组观测）：在每个受试者内部置换观测的factor标签。

3. 执行聚类置换检验 #

循环生成置换后的模型，计算检验统计量（这里用最大t值，也可选用聚类总t值）：

library(mgcv)
# 初始化存储置换后的最大t值
perm_max_t <- numeric(1000)

for (i in 1:1000) {
  # 生成置换后的分组标签（组间情况）
  perm_id_groups <- unique(dat$id)[perm_indices[i,]]
  perm_dat <- dat
  perm_dat$factor <- factor(ifelse(perm_dat$id %in% perm_id_groups, 
                                   levels(dat$factor)[1], levels(dat$factor)[2]))
  
  # 重新拟合模型
  perm_gam <- gam(y ~ te(a,b,by=factor) + covariate + s(id,bs="re"), 
                  data=perm_dat, method="REML")
  
  # 预测置换后的差异曲面并计算t值
  perm_pred <- predict(perm_gam, newdata = newdata, se.fit = TRUE)
  newdata$perm_fit <- perm_pred$fit
  perm_pred1 <- subset(newdata, factor == levels(dat$factor)[1])
  perm_pred2 <- subset(newdata, factor == levels(dat$factor)[2])
  perm_diff_t <- (perm_pred2$perm_fit - perm_pred1$perm_fit) / 
    sqrt(perm_pred1$se^2 + perm_pred2$se^2)
  
  # 记录最大t值
  perm_max_t[i] <- max(abs(perm_diff_t))
}

4. 确定显著差异区域 #

从置换分布中取临界值，标记原始t值超过临界值的区域：

# 双侧检验取97.5%分位数作为临界值
crit_val <- quantile(perm_max_t, 0.975)
# 标记显著区域
diff_surface$significant <- abs(diff_surface$t_val) > crit_val
# 可视化差异曲面+显著区域
library(ggplot2)
ggplot(diff_surface, aes(a, b, fill = diff)) +
  geom_tile() +
  geom_tile(data = subset(diff_surface, significant), fill = "red", alpha = 0.3) +
  scale_fill_viridis_c() +
  labs(title = "差异曲面（红色为显著差异区域）")

替代方法：参数化Bootstrap（gratia包） #

如果觉得置换检验计算量太大，可以用gratia包的difference()函数，它支持基于模型的参数化Bootstrap生成校正置信区间，自动处理多重比较：

library(gratia)
# 计算两组曲面的差异，带bootstrap置信区间
diff_obj <- difference(gam_obj, 
                       smooth = "te(a,b)", 
                       by = "factor",
                       n = 100,  # 网格点数
                       n_boot = 1000)
# 可视化差异曲面+显著区域（置信区间不包含0的区域）
draw(diff_obj, ci_level = 0.95, show_significant = TRUE)

关键注意事项 #

• 置换次数：至少1000次，复杂模型可通过foreach+doParallel实现并行计算。
• 聚类定义：若选用聚类总t值作为检验统计量，需明确连续区域规则（如4邻域/8邻域），避免孤立点被误判。
• 随机效应处理：绝对不能打乱所有观测的factor标签，必须保持受试者内的结构，否则会破坏模型的随机效应假设。

内容的提问来源于stack exchange，提问作者Josué García Arch