Bootstrap重采样处理分组嵌套数据拟合tidy回归模型报错问询

问题原因 #

第一次报错原因 #

你将nest()后的分组汇总表直接传入bootstraps()，此时重采样的对象是分组行而非每个分组内的原始观测，得到的splits对象中仅包含分组字段和存储子数据的data列，不存在conc字段，因此抛出变量未找到错误。

第二次调整后无有效模型的原因 #

map(boots, "splits", ~lm(...))的调用逻辑错误：boots列每个元素是单组的bootstrap结果表，提取splits后得到的是该组对应的1000个rsplit对象的列表，你需要对每个rsplit对象单独迭代拟合模型，现有写法仅迭代到bootstrap结果表层级，未深入到单个rsplit对象。

解决代码 #

library(tidyverse)
library(tidymodels)

# 示例数据
dat <- 
  structure(list(site = c("mb", "mb", "mb", "mb", "mb", "mb", "mb", 
  "mb", "sp", "sp", "sp", "sp", "sp", "sp", "sp", "sp"), year = c(2015, 
  2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
  2015, 2015, 2015, 2015), yday = c(15, 15, 35, 35, 48, 48, 69, 
  69, 15, 15, 37, 37, 49, 49, 69, 69), samp_depth_cat2 = structure(c(1L, 
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Mid-2", 
  "Bottom"), class = "factor"), analyte = c("NO3", "NO3", "NO3", 
  "NO3", "NO3", "NO3", "NO3", "NO3", "NH4", "NH4", "NH4", "NH4", 
  "NH4", "NH4", "NH4", "NH4"), conc = c(44.8171069465267, 44.7775358035268, 
  33.3678662097523, 33.0710828871279, 25.8427604055115, 26.9309658742058, 
  23.7585524380667, 17.5240386949382, 8.35832733633183, 9.29280745341615, 
  10.0797380595417, 10.2322058970515, 13.7930668951239, 15.6226805882773, 
  25.3003042764332, 16.8723637466981)), row.names = c(NA, -16L), class = c("tbl_df", 
  "tbl", "data.frame"))

set.seed(27)

# 修正后的运行代码
lm_boot <- dat %>% 
  # 按分组维度嵌套
  group_by(site, year, samp_depth_cat2, analyte) %>% 
  nest() %>% 
  # 对每个分组的独立数据集单独做bootstrap重采样
  mutate(boots = map(data, ~bootstraps(.x, times = 1000, apparent = TRUE))) %>% 
  # 展开bootstrap结果，每个重采样样本占单独一行
  unnest(boots) %>% 
  # 对每个重采样的split对象拟合模型，需用analysis()提取重采样后的数据集
  mutate(
    model = map(splits, ~lm(conc ~ yday, data = analysis(.x))),
    coef_info = map(model, tidy)
  )

# 提取所有重采样的回归系数
boot_coefs <- lm_boot %>% 
  unnest(coef_info)

# 按分组计算回归系数的百分位置信区间
percentile_intervals <- int_pctl(
  .data = lm_boot,
  statistics = coef_info,
  .by = c(site, year, samp_depth_cat2, analyte)
)

# 查看结果
percentile_intervals

核心调整点 #

• 分组嵌套后仅对每个分组内部的原始观测做bootstrap，保证组间独立性，符合分组分析的逻辑
• 用rsample内置的analysis()函数从split对象中提取重采样后的分析数据集，避免直接传入split对象本身导致的变量识别错误
• 调用int_pctl时通过.by参数指定分组维度，直接输出每个分组对应的截距、斜率的置信区间结果

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