GAM模型后验模拟绘图异常：抗生素分组预测问题排查

问题背景 #

运行以下GAM模型：

modGI_test <- gam(Met1 ~ Antibiotics + s(Timepoint, k=7) +
                s(Timepoint, by=Antibiotics, k=7, m=1) + 
                s(Timepoint, ID_new, bs="fs", k=7) +
                Age_use + Gender, 
                data=df_full, method="REML")

通过以下代码获取后验模拟信息：

beta <- coef(modGI_test)  
Vb <- vcov(modGI_test) 
Xp <- predict(modGI_test, type = "lpmatrix") 
Beta_sim <- rmvn(n=1000, beta, Vb) 

尝试绘制抗生素分组的时间平滑趋势图，但结果异常（置信区间混乱、曲线无合理时间趋势），绘图代码如下：

##antibiotics yes 
Xp_antibioticsyes <- Xp[, grep("s\(Timepoint\):Antibiotics(Yes)\.\d+", colnames(Xp))]
Beta_antibioticsyes <- Beta_sim[, grep("s\(Timepoint\):Antibiotics(Yes)\.\d+", colnames(Beta_sim))]
Eta_antibioticsyes <- Xp_antibioticsyes %*% t(Beta_antibioticsyes)
Mu_antibioticsyes <- inv_link(modGI_test)(Eta_antibioticsyes)

##antibiotics no 
Xp_antibioticsno <- Xp[, antibiotics_cols_no]
Beta_antibioticsno <- Beta_sim[,  grep("s\(Timepoint\):Antibiotics(No)\.\d+", colnames(Beta_sim))]
Eta_antibioticsno <- Xp_antibioticsno %*% t(Beta_antibioticsno)
Mu_antibioticsno <- inv_link(modGI_test)(Eta_antibioticsno)

global_intercept <- Xp[, "(Intercept)"] %*% t(Beta_sim[, "(Intercept)"])
main_effect <- Xp[, "AntibioticsYes"] %*% t(Beta_sim[, "AntibioticsYes"])

# Add these to the group-specific predictions
Mu_antibioticsyes <- Mu_antibioticsyes + global_intercept + main_effect
Mu_antibioticsno <- Mu_antibioticsno + global_intercept  # No main effect for "No"

##showing only antibiotics yes here
antibiotics_yes <- data.frame(t(apply(Mu_antibioticsyes, 1, function(x) {
  c(mean = mean(x, na.rm = TRUE),
  lowerci = quantile(x, probs = 0.025),
  upperci = quantile(x, probs = 0.975))
})))

antibiotics_yes$Timepoint <- stackoverflow$Timepoint
antibiotics_yes$Group <- rep("Yes", nrow(antibiotics_yes))

 ##doing the same for antibiotics no, rbinding them , then running ggplot with geom_ribbon

补充数据：

structure(list(ID_new = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
  2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 
  4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 
  6L, 6L, 6L), levels = c("1", "3", "4", "5", "6", "7"), class = "factor"), 
  Age_use = c(61L, 61L, 61L, 61L, 61L, 61L, 62L, 51L, 51L, 
  51L, 51L, 51L, 51L, 52L, 30L, 30L, 30L, 30L, 30L, 30L, 31L, 
  50L, 50L, 50L, 50L, 50L, 50L, 51L, 43L, 43L, 43L, 43L, 43L, 
  43L, 44L, 30L, 30L, 30L, 30L, 30L, 30L, 31L), Gender = structure(c(2L, 
  2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
  2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
  2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), levels = c("Female", 
  "Male"), class = "factor"), Antibiotics = structure(c(1L, 
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
  2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("No", 
  "Yes"), class = "factor"), Timepoint = c(1L, 2L, 3L, 4L, 
  5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 
  6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 
  7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L), Met1 = c(0.610312217, -0.892993391, 
  -0.193675438, 0.725636548, 1.028557245, 0.658419815, -0.608733906, 
  0.788168303, 0.031552182, -0.007267893, -0.612162917, -0.016844181, 
  0.773175345, -1.07871784, -0.255036005, 0.282531559, -0.218753849, 
 -0.823197508, 0.415051062, 1.248028647, 0.761003283, -0.413166119, 
  0.384806119, -0.00364492, 0.559821831, 0.091177022, -0.347510957, 
  0.509460406, -0.132033165, 0.479567493, -0.019536566, -0.190650928, 
  0.331425516, 0.691209456, 0.422804123, 0.024092717, -0.294300556, 
  -1.088126749, -1.879112138, 0.277489963, -0.73213862, 1.035703532
  )), row.names = c(NA, -42L), class = "data.frame")

错误原因分析 #

1. 正则表达式匹配错误：mgcv中by=Antibiotics生成的列名格式是s(Timepoint):AntibioticsYes（无括号包裹Yes），原代码里的Antibiotics(Yes)是错误的，且括号未正确转义，导致提取不到正确的平滑项列。
2. 遗漏全局时间平滑项：模型中的s(Timepoint)是所有样本的基础时间趋势，原代码未将其加入到分组预测中，导致趋势缺失核心部分。
3. 固定效应缺失：未加入Age_use和Gender的固定效应贡献，预测值不完整。
4. 个体随机效应处理混乱：若要绘制群体层面的分组趋势，应排除个体特异性的s(Timepoint, ID_new, bs="fs")项；若要包含，需正确提取，但原代码未处理这部分，导致预测值混入个体变异或缺失。
5. 数据引用错误：stackoverflow$Timepoint是无效的变量名，应使用原数据df_full$Timepoint。

修正后的代码 #

1. 定义各效应的列名匹配规则 #

# 全局时间平滑项
global_time_cols <- grep("^s\\(Timepoint\\)\\.\\d+$", colnames(Xp))
# 抗生素Yes组的时间平滑项（by项）
abx_yes_time_cols <- grep("^s\\(Timepoint\\):AntibioticsYes\\.\\d+$", colnames(Xp))
# 固定效应列：截距、抗生素主效应、年龄、性别
fixed_cols <- c(grep("^\\(Intercept\\)$", colnames(Xp)),
                grep("^AntibioticsYes$", colnames(Xp)),
                grep("^Age_use$", colnames(Xp)),
                grep("^GenderMale$", colnames(Xp)))
# 个体随机效应（如果绘制群体图则排除，个体图则保留）
indiv_time_cols <- grep("^s\\(Timepoint\\,ID_new\\)\\.\\d+$", colnames(Xp))

2. 计算群体层面的分组预测（排除个体随机效应） #

# 构建预测矩阵：固定效应 + 全局时间 + 分组时间
# --- 抗生素Yes组 ---
Xp_yes <- Xp[, c(fixed_cols, global_time_cols, abx_yes_time_cols)]
Beta_yes <- Beta_sim[, c(fixed_cols, global_time_cols, abx_yes_time_cols)]
Eta_yes <- Xp_yes %*% t(Beta_yes)
Mu_yes <- inv_link(modGI_test)(Eta_yes)

# --- 抗生素No组 ---
# No组没有by项，所以是固定效应 + 全局时间
Xp_no <- Xp[, c(fixed_cols[!grepl("AntibioticsYes", colnames(Xp)[fixed_cols])], 
                global_time_cols)]
Beta_no <- Beta_sim[, colnames(Xp_no)]
Eta_no <- Xp_no %*% t(Beta_no)
Mu_no <- inv_link(modGI_test)(Eta_no)

3. 计算均值与置信区间并整理数据 #

# 处理Yes组
yes_stats <- data.frame(t(apply(Mu_yes, 1, function(x) {
  c(mean = mean(x, na.rm = TRUE),
    lowerci = quantile(x, probs = 0.025),
    upperci = quantile(x, probs = 0.975))
})))
yes_stats$Timepoint <- df_full$Timepoint
yes_stats$Group <- "Yes"

# 处理No组
no_stats <- data.frame(t(apply(Mu_no, 1, function(x) {
  c(mean = mean(x, na.rm = TRUE),
    lowerci = quantile(x, probs = 0.025),
    upperci = quantile(x, probs = 0.975))
})))
no_stats$Timepoint <- df_full$Timepoint
no_stats$Group <- "No"

# 合并数据
plot_data <- rbind(yes_stats, no_stats)
# 按时间点和分组排序
plot_data <- plot_data[order(plot_data$Timepoint, plot_data$Group), ]

4. 绘制ggplot图 #

library(ggplot2)
ggplot(plot_data, aes(x = Timepoint, y = mean, color = Group)) +
  geom_ribbon(aes(ymin = lowerci, ymax = upperci, fill = Group), 
              alpha = 0.2, color = NA) +
  geom_line(linewidth = 1) +
  labs(x = "时间点", y = "Met1预测值", color = "抗生素使用", fill = "抗生素使用") +
  theme_minimal()

额外说明 #

• 如果需要绘制个体层面的时间趋势，只需在预测矩阵中加入indiv_time_cols即可。
• 对于全局s(Timepoint)的单独绘图，只需提取global_time_cols对应的列，结合截距和其他固定效应计算预测值即可。

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