R语言基于bootstrap方法比较两个分箱数据集均值的相关咨询

1 异常值处理方法 #

你当前场景下优先推荐IQR分位数截断法，鲁棒性不受极端值影响，具体逻辑：

• 对每个分箱单独计算第一四分位数Q1、第三四分位数Q3，得到四分位距IQR=Q3-Q1
• 判定超出[Q1 - 1.5*IQR, Q3 + 1.5*IQR]的数值为异常值，可选择直接删除或替换为对应边界值
• 若不想损失样本量，也可直接用分箱内的中位数替代均值做统计，中位数本身对异常值不敏感

2 完整实现代码 #

2.1 依赖包加载与异常值处理函数定义 #

# 加载依赖包，初学者直接安装tidyverse即可：install.packages("tidyverse")
library(tidyverse)

# 自定义IQR法异常值清洗函数，输入数值向量，返回清洗后的向量
clean_outlier <- function(x) {
  q1 <- quantile(x, 0.25, na.rm = T)
  q3 <- quantile(x, 0.75, na.rm = T)
  iqr <- q3 - q1
  lower <- q1 - 1.5*iqr
  upper <- q3 + 1.5*iqr
  x[x < lower | x > upper] <- NA
  # 这里返回删除异常值后的向量，要替换的话把上面NA改成lower或upper即可
  return(na.omit(x))
}

2.2 Bootstrap抽样与均值计算 #

# 先代入你提供的示例数据集
dataset1=structure(list(genenames = c("data1", "data2", "data3", "data4", "data5", "data6"), 
                        bin1 = c(0,20,9,0,2,0), 
                        bin2 = c(5,20,8,30,10,0), 
                        bin3 = c(0,0,1,1,3,0),
                        bin4 =c(6, 20, 10, 5, 0, 1),
                        bin5 =c(10,15,30,10,9, 4)), 
                   class = "data.frame", row.names = c(NA, -6L))

dataset2=structure(list(genenames = c("data10", "data11", "data12", "data13", "data14", "data15"), 
                        bin1 = c(0,30,0,0,20,0), 
                        bin2 = c(0,0,8,10,20,0), 
                        bin3 = c(0,10,19,15,3,10),
                        bin4 =c(30, 0, 0, 25, 0, 20),
                        bin5 =c(0,5,0,20,30, 29)), 
                   class = "data.frame", row.names = c(NA, -6L))

# 参数设置，你实际使用时把n_sample改成100，n_iter改成1000即可
set.seed(123) # 固定随机种子保证结果可复现
n_sample <- 3 # 示例数据只有6个样本，所以抽样量设为3，实际使用时改100
n_iter <- 100 # 重复抽样次数，实际使用建议≥1000

# 定义bootstrap计算函数
boot_mean <- function(data, n_sample, n_iter) {
  # 先去掉genenames列，对每个分箱做异常值清洗
  clean_data <- data[,-1] %>% mutate(across(everything(), clean_outlier))
  # 重复抽样计算均值
  res <- replicate(n_iter, {
    # 对每个分箱有放回抽样n_sample个，计算均值
    clean_data %>% 
      summarise(across(everything(), ~mean(sample(.x, n_sample, replace = T), na.rm = T)))
  })
  # 整理结果：返回每个分箱的bootstrap均值、95%置信区间上下限
  res %>% 
    t() %>% 
    as.data.frame() %>% 
    pivot_longer(everything(), names_to = "bin", values_to = "mean") %>% 
    group_by(bin) %>% 
    summarise(
      boot_mean = mean(mean),
      ci_low = quantile(mean, 0.025),
      ci_high = quantile(mean, 0.975)
    )
}

# 分别计算两个数据集的bootstrap结果
boot1 <- boot_mean(dataset1, n_sample, n_iter) %>% mutate(dataset = "数据集1")
boot2 <- boot_mean(dataset2, n_sample, n_iter) %>% mutate(dataset = "数据集2")
boot_all <- bind_rows(boot1, boot2)

2.3 分箱均值对比可视化 #

ggplot(boot_all, aes(x = bin, y = boot_mean, color = dataset, group = dataset)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  # 可选添加置信区间带
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high, fill = dataset), alpha = 0.2, color = NA) +
  labs(x = "分箱", y = "Bootstrap校正后均值", color = "数据集", fill = "数据集") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

3 Bootstrap使用注意事项 #

• 必须设置set.seed()固定随机种子，否则每次运行结果不一致，无法复现
• 抽样必须使用有放回抽样（即sample函数的replace = T），这是bootstrap方法的标准要求，无放回抽样属于子采样，会低估结果的方差
• 重复抽样次数建议设置在1000~5000区间，次数太少会导致均值和置信区间估计不稳定，次数太多会不必要的增加运算耗时
• 异常值清洗必须在抽样之前完成，不要先抽样再做异常值处理，否则异常值会进入抽样池干扰估计结果
• 若两个数据集同个分箱的样本量差异超过20%，建议抽样量设置为较小样本量的50%，不要超过最小样本量的80%，避免出现抽样偏差

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