R中使用tidymodels为分组数据框列表回归结果实现Bootstrap抽样

实现逻辑 #

在每个year+group拆分后的子组内做有放回自助抽样，不需要修改原有分组拆分逻辑，只需要对每个子组新增bootstrap抽样流程即可，以下是两种实现方案：

方案1：基于现有代码的最小改动 #

和你原有写法逻辑完全兼容，改动量最小：

library(dplyr)

# 原始数据生成逻辑不变
year <- rep(2014:2018, length.out=10000)
group <- sample(c(0,1,2,3,4,5,6), replace=TRUE, size=10000)
value <- sample(10000, replace=T)
female <- sample(c(0,1), replace=TRUE, size=10000)
smoker <- sample(c(0,1), replace=TRUE, size=10000)
dta <- data.frame(year=year, group=group, value=value, female=female, smoker=smoker)

# 设定bootstrap抽样次数，可根据算力和精度需要调整
R <- 100

# 1. 按year+group拆分子组，逻辑和原有代码完全一致
table_list <- dta %>%
  group_by(year, group) %>%
  group_split()

# 2. 对每个子组做R次bootstrap抽样、拟合模型、生成预测
boot_result <- lapply(table_list, function(sub_df) {
  # 对当前子组生成R个bootstrap样本
  lapply(1:R, function(i) {
    # 有放回抽样，样本量和原个子组大小一致
    boot_idx <- sample(nrow(sub_df), size = nrow(sub_df), replace = TRUE)
    boot_df <- sub_df[boot_idx, ]
    # 拟合probit模型
    mod <- glm(smoker ~ female, data = boot_df, family = binomial(link = "probit"))
    # 生成预测值：如果需要对bootstrap样本本身预测，把newdata改为boot_df即可
    pred <- predict.glm(mod, newdata = sub_df, type = "response")
    return(pred)
  }) %>% do.call(cbind, .) # 把R次预测按列合并，每列对应一次bootstrap结果
})

# 最终boot_result是和原table_list等长的列表，每个元素是子组所有观测的R次bootstrap预测值矩阵
# 可以直接按需计算自助法统计量，比如均值、95%置信区间：
boot_pred_mean <- lapply(boot_result, rowMeans)
boot_pred_ci <- lapply(boot_result, function(x) t(apply(x, 1, quantile, probs = c(0.025, 0.975), na.rm = TRUE)))

方案2：tidymodels原生实现 #

符合tidymodels工具链规范，方便后续扩展模型验证、参数调优等功能：

library(dplyr)
library(tidyr)
library(rsample)
library(purrr)

# 原始数据不变
year <- rep(2014:2018, length.out=10000)
group <- sample(c(0,1,2,3,4,5,6), replace=TRUE, size=10000)
value <- sample(10000, replace=T)
female <- sample(c(0,1), replace=TRUE, size=10000)
smoker <- sample(c(0,1), replace=TRUE, size=10000)
dta <- data.frame(year=year, group=group, value=value, female=female, smoker=smoker)

R <- 100 # 自助抽样次数

# 嵌套分组+分组内bootstrap全流程
boot_res <- dta %>%
  nest(data = -c(year, group)) %>% # 按year+group嵌套生成子数据集
  mutate(
    # 对每个子组生成R个bootstrap样本
    boots = map(data, ~ bootstraps(.x, times = R, apparent = FALSE)),
    # 对每个bootstrap样本拟合模型
    models = map(boots, ~ map(.x$splits, function(split) {
      glm(smoker ~ female, data = analysis(split), family = binomial(link = "probit"))
    })),
    # 生成对原样本的预测
    pred = map2(models, data, ~ map_dfc(.x, function(mod) predict(mod, newdata = .y, type = "response")))
  )

# 提取所有预测结果：boot_res$pred就是每个子组的R次bootstrap预测值矩阵

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