R语言分层Bootstrap回归自动化咨询：多分组批量处理方案

我已经能熟练实现基础的Bootstrap回归了，下面是我的示例代码：

library(boot)
sampledata = data.frame(y = sample(0:20, r = T, size = 5000), x1 = runif(5000), x2 = runif(5000), group1 = sample(1:3, r = T, size = 5000), group2 = sample(0:1, r = T, size = 5000))
# 用于获取模型系数的函数
rsq <- function(formula, data, indices) {
  d <- data[indices,]
  fit <- lm(formula, data=d)
  return(coef(fit))
}
# 250次重复抽样的Bootstrap
results <- boot(data=sampledata, statistic=rsq, R=250, formula=y~x1+x2)
df <- data.frame(results$t)
names(df) <- names(results$t0)

但现在我需要针对数据的子分组自动化执行这个流程，目前手动拆分6个分组分别运行Bootstrap再合并结果的方式太繁琐了，而且实际工作里分组类别更多。我注意到boot包有strata参数，但一直没法正常使用它。另外我还有个疑问：到底应该先按分组拆分数据再做Bootstrap抽样，还是先抽样再拆分？想请教下怎么写代码实现所有分组的Bootstrap估计自动化，处理这种分层模型的需求。

手动处理的核心代码示例 #

sampledata1 = subset(sampledata, group1 == 1 & group2 == 0)
# 其余5个分组的拆分代码类似
results <- boot(data=sampledata1, statistic=rsq, R=250, formula=y~x1+x2)
df1 <- data.frame(results$t)
names(df1) <- names(results$t0)
df1$group1 = 1
df1$group2 = 0
# 其余分组的Bootstrap运行代码类似
DATA = rbind(df1, df2, df3, df4, df5, df6)

我自己试着写了个函数，但好像没法正常工作：

rsq <- function(formula, data, indices) {
  dataSUB <- data %>% filter(GROUP1 == group1, GROUP2 == group2)
  d <- dataSUB[indices,]
  fit <- lm(formula, data=d)
  return(coef(fit))
}

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