R中基于分层子组的Bootstrap抽样代码正确性验证请求

问题背景 #

我有一个覆盖波罗的海南部的数据集，包含以下采样特征：

• 部分站点每年采样一次，仅采集一个抓取样本
• 部分站点仅采样一次且仅一个抓取样本
• 部分站点连续多年采样，每次采集多个抓取样本

数据集核心字段：

• Stations：站点的唯一名称
• Stationname：对应采样时间的站点名
• Stationname_Hol：每个抓取样本对应的站点名
• bqi：目标分析变量

需求目标：

1. 随机选择一个站点，再从该站点的采样时间中随机选一个Stationname，接着从该Stationname对应的抓取样本中随机选一个Stationname_Hol，记录对应的bqi值
2. 每次抽取69个（与站点总数相同）上述方式生成的bqi值，重复该过程10000次
3. 对每次抽取的69个bqi值计算第20百分位数，最终求这10000个百分位数的中位数

数据集概况 #

• 270条Stationname_Hol（抓取样本）记录
• 96条Stationname（采样时间）记录
• 69个唯一的Stations（站点）
• 对应bqi数值

现有代码 #

Ref_EIG4a <- subset(Data_EIG4a, Within_RefArea=="Referenz", select = c(stationsnamen_hol, stationsnamen, Stations, bqi))
    
for (i in 1:10000) { 
  for (j in 1:69) {
    station <- sample(unique(Ref_EIG4a$Stations), 1)
    stationname <- sample(unique(Ref_EIG4a$stationsnamen[Ref_EIG4a$Stations == station]), 1)
    hol <- sample(unique(Ref_EIG4a$stationsnamen_hol[Ref_EIG4a$stationsnamen == stationname]), 1)
    if (j==1) {
      picked_bqi <- (Ref_EIG4a$bqi[Ref_EIG4a$stationsnamen_hol == hol])
    }
    else {
      picked_bqi <- c(picked_bqi, Ref_EIG4a$bqi[Ref_EIG4a$stationsnamen_hol == hol])
    }
  }
  if (i==1) {
    Q20_BQI <- quantile(picked_bqi, probs = 0.2)
  }
  else {
    Q20_BQI <- c(Q20_BQI, quantile(picked_bqi, probs = 0.2))
  }
}

MD_Q20_Boot_Ref_EIG4a <- median(Q20_BQI)  #Median berechnen 

代码正确性分析 #

核心逻辑匹配度 #

现有代码的逻辑是：

• 外层循环执行10000次抽样过程
• 内层循环每次从所有站点中随机选1个（有放回），再依次选择对应采样时间、抓取样本，最终拼接出69个bqi值
• 对每次的69个值计算20%分位数，最后取中位数

如果你的需求是允许每次抽取的69个bqi值对应重复站点（即每次抽样是有放回地从站点池选69次），那么代码的核心逻辑是正确的。但如果你的需求是每次抽取的69个bqi值必须覆盖所有69个站点（每个站点恰好选一次），那么代码存在错误——内层循环的站点选择是有放回的，会导致部分站点被重复选中，部分站点未被选中。这种情况下，应该遍历所有站点，再对每个站点执行后续的采样时间、抓取样本选择。

潜在问题 #

1. 效率低下：
   • 每次循环中用c()拼接picked_bqi和Q20_BQI向量，会频繁复制内存，10000次循环会显著拖慢运行速度，建议预先分配向量空间。
2. 抽样的明确性：
   • sample()函数默认是无放回抽样，但当仅抽取1个元素时，有放回和无放回效果一致；如果后续调整抽样数量，需明确指定replace参数。
3. 变量命名一致性：
   • 数据集字段名存在大小写不一致（如Stations vs stationsnamen），虽不影响运行，但降低了代码可读性。

优化后的代码 #

场景1：允许重复站点（与原代码逻辑一致，优化效率） #

# 预处理数据
Ref_EIG4a <- subset(Data_EIG4a, Within_RefArea == "Referenz", 
                    select = c(stationsnamen_hol, stationsnamen, Stations, bqi))

# 预先分配结果向量，提升运行效率
Q20_BQI <- numeric(10000)

for (i in 1:10000) {
  picked_bqi <- numeric(69)  # 预先分配当前轮次的bqi存储向量
  for (j in 1:69) {
    # 随机选择站点
    station <- sample(unique(Ref_EIG4a$Stations), 1)
    # 从该站点的采样时间中随机选一个
    stationname <- sample(unique(Ref_EIG4a$stationsnamen[Ref_EIG4a$Stations == station]), 1)
    # 从该采样时间的抓取样本中随机选一个
    hol <- sample(unique(Ref_EIG4a$stationsnamen_hol[Ref_EIG4a$stationsnamen == stationname]), 1)
    # 提取对应bqi值
    picked_bqi[j] <- Ref_EIG4a$bqi[Ref_EIG4a$stationsnamen_hol == hol]
  }
  # 计算当前轮次的20%分位数并存储
  Q20_BQI[i] <- quantile(picked_bqi, probs = 0.2)
}

# 计算最终中位数
MD_Q20_Boot_Ref_EIG4a <- median(Q20_BQI)

场景2：必须覆盖所有69个站点（每个站点选一次） #

# 预处理数据
Ref_EIG4a <- subset(Data_EIG4a, Within_RefArea == "Referenz", 
                    select = c(stationsnamen_hol, stationsnamen, Stations, bqi))

# 获取所有唯一站点列表
all_stations <- unique(Ref_EIG4a$Stations)
Q20_BQI <- numeric(10000)

for (i in 1:10000) {
  picked_bqi <- numeric(69)
  # 遍历所有站点（若需随机顺序，可替换为for (station in sample(all_stations))）
  for (j in seq_along(all_stations)) {
    station <- all_stations[j]
    stationname <- sample(unique(Ref_EIG4a$stationsnamen[Ref_EIG4a$Stations == station]), 1)
    hol <- sample(unique(Ref_EIG4a$stationsnamen_hol[Ref_EIG4a$stationsnamen == stationname]), 1)
    picked_bqi[j] <- Ref_EIG4a$bqi[Ref_EIG4a$stationsnamen_hol == hol]
  }
  Q20_BQI[i] <- quantile(picked_bqi, probs = 0.2)
}

MD_Q20_Boot_Ref_EIG4a <- median(Q20_BQI)

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