如何在GAM中分析两因子与连续变量的三向交互及组间差异

研究背景 #

我们的研究要检验碳源添加（factor1：amended/unamended）、温度（factor2：high/low）两个因子及其交互对土壤CO₂通量随时间变化的影响，实验采用完全析因设计，36个土壤样本在9天内重复测量。目前已拟合两个GAM模型，但都无法直接判断四个处理组的平滑曲线是否存在两两显著差异，以下是具体解决办法：

针对模型1（分组平滑项GAM）的事后检验 #

模型1用了有序因子f1f2_o的分组平滑项，要做两两比较，按以下步骤操作：

1. 提取各组拟合曲线与置信区间
   先构造包含所有处理组和连续day序列的新数据集，再用predict()获取每个组在不同day的拟合值和标准误，进而计算置信区间：

   # 构造用于预测的新数据集
   new_data <- expand.grid(
     day = seq(min(data$day), max(data$day), length.out = 100),
     f1f2_o = unique(data$f1f2_o),
     ID = sample(unique(data$ID), 1) # 随便选一个ID即可，随机效应不影响组间比较
   )
   # 获取预测值和标准误
   pred_result <- predict(m1, newdata = new_data, se.fit = TRUE)
   new_data$fit_value <- pred_result$fit
   new_data$se <- pred_result$se.fit
   new_data$ci_lower <- new_data$fit_value - 1.96 * new_data$se
   new_data$ci_upper <- new_data$fit_value + 1.96 * new_data$se
   

2. 直接做平滑项的两两比较
   使用mgcv包的compare_smooths()函数（需mgcv版本≥1.8-38），直接对分组平滑项进行两两显著性检验：

   library(mgcv)
   # 指定要比较的分组平滑项
   smooth_comp <- compare_smooths(m1, smooth = "s(day):f1f2_o")
   # 查看两两比较的显著性结果
   print(smooth_comp)
   

   这个函数会输出每两组平滑曲线差异的统计量和p值，能直接看出哪几组的曲线形状存在显著差异。

针对模型2（sz样条交互项GAM）的事后检验 #

模型2用了sz样条构建三向交互，要拆解两两处理组的差异，有两种可行思路：

1. 重构分组变量后拟合新模型
   先在数据集中创建包含四个处理组的组合变量，再拟合一个类似模型1的分组平滑GAM，最后用compare_smooths()做比较：

   # 创建组合分组变量
   data$treatment_group <- interaction(data$factor1, data$factor2)
   # 重新拟合分组平滑GAM
   m3 <- gam(flux ~ treatment_group + s(day, k = 7) + s(day, by = treatment_group, k = 7) + s(ID, bs = "re"),
             data = data, method = "REML")
   # 对新模型的分组平滑项做两两比较
   comp_m3 <- compare_smooths(m3, smooth = "s(day):treatment_group")
   print(comp_m3)
   

2. 构建对比矩阵做精准检验
   若只想检验特定两组的差异，可通过构建对比矩阵实现。比如检验amended_high和unamended_low的曲线差异：

   # 查看模型参数命名，定位对应组的平滑项参数位置
   coef_names <- names(coef(m2))
   idx_ah <- grep("amended.high", coef_names)
   idx_ul <- grep("unamended.low", coef_names)
   # 构建对比向量：让两组参数相减，检验差值是否为0
   contrast_vec <- rep(0, length(coef(m2)))
   contrast_vec[idx_ah] <- 1
   contrast_vec[idx_ul] <- -1
   # 执行显著性检验
   diff_test <- linearHypothesis(m2, contrast_vec)
   print(diff_test)
   

   这种方法需要熟悉模型参数的命名规则，适合针对性检验特定组的差异。

可视化辅助验证 #

不管用哪种统计方法，都建议把各组的平滑曲线和置信区间画出来，直观观察差异：

library(ggplot2)
ggplot(new_data, aes(x = day, y = fit_value, color = f1f2_o)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper, fill = f1f2_o), alpha = 0.2) +
  labs(x = "天数", y = "CO₂通量", color = "处理组", fill = "处理组") +
  theme_bw()

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