位置尺度广义可加模型RMSE估计的不确定性Bootstrap方法问询

问题背景 #

拥有按基因型（Genotype，因子型）、单株（Individual，因子型）分组的植物株高（Height，数值型，单位cm）随日序（Doy，数值型）观测数据（每单株每日1条记录），已用位置尺度广义可加模型（location-scale GAM）计算RMSE，但无法通过Bootstrap估计RMSE的不确定性。此前尝试的循环抽样代码运行缓慢且每次返回相同RMSE，抽样逻辑存在问题。

问题分析 #

• 抽样逻辑错误：原代码按Doy和Genotype分组抽样，破坏了单株的时间序列完整性，且未实现“按基因型抽取固定数量单株”的需求。
• 固定种子导致结果一致：CVgam中固定seed=29，每次交叉验证的分组完全相同，因此RMSE无变化。
• 未初始化结果向量：RMSE未预先定义长度，可能导致性能问题。
• 模型重复拟合效率低：循环内重复拟合复杂GAM模型，未做优化。

修正方案 #

1. 正确的Bootstrap抽样逻辑 #

按Genotype分组，每个基因型内抽取固定数量（如n=4）的单株，保留这些单株的所有观测数据，确保单株时间序列完整。

2. 动态设置随机种子 #

每次循环使用不同随机种子，保证抽样和交叉验证的随机性。

3. 预初始化结果向量 #

提前定义RMSE向量，提升运行效率。

4. 可选：并行加速循环 #

利用foreach和doParallel包并行运行循环，减少总耗时。

修正后的代码 #

基础版（单线程） #

# 初始化RMSE结果向量
n_boot <- 100  # Bootstrap次数，可根据需求调整
RMSE <- numeric(n_boot)

# 定义每个基因型抽取的单株数量
n_ind_per_geno <- 4

set.seed(123)  # 全局种子保证可重复性
for (i in 1:n_boot) {
  # 按基因型分组，抽取固定数量单株
  sampled_individuals <- data %>%
    group_by(Genotype) %>%
    slice_sample(n = n_ind_per_geno, replace = TRUE) %>%  # replace=TRUE允许同一单株多次被抽中
    pull(Individual) %>%
    unique()
  
  # 提取抽样得到的单株的所有观测数据
  datax <- data %>% filter(Individual %in% sampled_individuals)
  
  # 拟合位置尺度GAM模型
  model <- gam(list(Height ~ s(Doy, bs = 'ps', by = Genotype) +
                      s(Doy, Individual, bs = "re") +
                      Genotype, 
                    ~ s(Doy, bs = 'ps', by = Genotype) +
                      s(Doy, Individual, bs = "re") +
                      Genotype),
               family = gaulss(),
               method = "REML",
               data = datax)
  
  # 提取位置项公式并做交叉验证
  form <- formula.gam(model)
  CV <- CVgam(form[[1]], datax, nfold = 10, debug.level = 0, 
              method = "GCV.Cp", printit = FALSE, cvparts = NULL, gamma = 1)
  
  # 计算并存储RMSE
  RMSE[i] <- sqrt(CV$cvscale[1])
}

# 查看RMSE的分布
summary(RMSE)
hist(RMSE)

并行加速版 #

library(foreach)
library(doParallel)

n_boot <- 100
n_ind_per_geno <- 4
RMSE <- numeric(n_boot)

# 设置并行核心数（根据电脑配置调整，一般为核心数-1）
cl <- makeCluster(4)
registerDoParallel(cl)

set.seed(123)
RMSE <- foreach(i = 1:n_boot, .combine = 'c', .packages = 'mgcv') %dopar% {
  # 抽样逻辑同基础版
  sampled_individuals <- data %>%
    group_by(Genotype) %>%
    slice_sample(n = n_ind_per_geno, replace = TRUE) %>%
    pull(Individual) %>%
    unique()
  
  datax <- data %>% filter(Individual %in% sampled_individuals)
  
  # 拟合模型
  model <- gam(list(Height ~ s(Doy, bs = 'ps', by = Genotype) +
                      s(Doy, Individual, bs = "re") +
                      Genotype, 
                    ~ s(Doy, bs = 'ps', by = Genotype) +
                      s(Doy, Individual, bs = "re") +
                      Genotype),
               family = gaulss(),
               method = "REML",
               data = datax)
  
  form <- formula.gam(model)
  CV <- CVgam(form[[1]], datax, nfold = 10, debug.level = 0, 
              method = "GCV.Cp", printit = FALSE, cvparts = NULL, gamma = 1)
  
  sqrt(CV$cvscale[1])
}

stopCluster(cl)

# 分析结果
summary(RMSE)
hist(RMSE)

关键说明 #

• 抽样逻辑：先按基因型抽取单株，再提取这些单株的所有观测，保证了单株时间序列的完整性，符合分层数据的Bootstrap要求。
• 随机种子：移除CVgam中的固定seed，全局设置一个种子保证整体可重复性，同时每次循环的抽样和交叉验证都是随机的。
• 并行加速：foreach结合doParallel可大幅减少循环耗时，尤其当Bootstrap次数较多时效果明显。

内容的提问来源于stack exchange，提问作者Bertold Mariën