R语言加权抽样结果存在偏差，请求排查问题原因

问题背景 #

需求是基于LPP组的nr.genomes概率分布，对CHR组的行进行加权抽样，使抽样结果的nr.genomes分布匹配LPP组，但当前抽样结果存在明显偏差。现有代码流程如下：

现有代码 #

# 数据框列名定义
colnames(table) <- c('PF', 'verticality', 'nr.genomes', 'mean.branch', 'tot.branch','nr.species', 'replicon')

# 从LPP组计算nr.genomes的概率分布
freq <- table(table$nr.genomes[(table$replicon=='LPP')])
prob <- as.data.frame(prop.table(freq))
colnames(prob) <- c('nr.genomes', 'prob')

# 合并概率到原数据框并筛选CHR组
table.probs = merge(x = table, y = prob, by = "nr.genomes")
table.chr <- table.probs[(table.probs$replicon=='CHR'),]

# 多轮抽样
df.final <- data.frame() # 补充初始化
for(i in 1:100) {
  sampled_rows <- sample(1:NROW(table.chr), size = 1097, prob = table.chr$prob)
  sampled_df <- table.chr[sampled_rows, ]
  run <- rep(as.character(i), 1097)
  df <- cbind(run, sampled_df)
  df.final <- rbind(df.final, df)
}

偏差原因分析 #

1. 概率权重逻辑错误
   你从LPP组计算的是每个nr.genomes值的全局占比，但直接将这个概率赋值给CHR组中所有对应nr.genomes的行，会导致抽样时的实际权重是「该nr.genomes的概率 × 该类别在CHR组中的行数」。举个例子：

   • LPP组中nr.genomes=5的概率是0.3，nr.genomes=10的概率是0.3
   • CHR组中nr.genomes=5有200行，nr.genomes=10只有20行
   • 抽样时，nr.genomes=5的总权重是200×0.3=60，nr.genomes=10的总权重是20×0.3=6
     最终抽样结果中nr.genomes=5的占比会远高于0.3，完全偏离LPP组的目标分布。

2. 抽样对象逻辑混淆
   你的目标是让抽样结果的nr.genomes分布匹配LPP组，但当前代码是对行按nr.genomes对应的概率抽样，而不是先对**nr.genomes类别**按目标分布分配抽样数量，再从对应类别中抽取行。

修正方案 #

正确的思路是先按LPP组的分布确定每个nr.genomes类别需要抽取的行数，再从CHR组的对应类别中随机抽取对应数量的行：

# 1. 从LPP组获取目标分布
lpp_data <- table[table$replicon == 'LPP', ]
target_dist <- prop.table(table(lpp_data$nr.genomes))
target_dist_df <- as.data.frame(target_dist)
colnames(target_dist_df) <- c('nr.genomes', 'prob')

# 2. 准备CHR组数据
chr_data <- table[table$replicon == 'CHR', ]

# 3. 初始化最终结果框
df.final <- data.frame()

# 4. 多轮抽样
total_sample <- 1097
n_runs <- 100

for(i in 1:n_runs) {
  # 计算每个nr.genomes需要抽取的数量（四舍五入处理，确保总数符合要求）
  target_counts <- round(total_sample * target_dist_df$prob)
  # 调整总数，避免因四舍五入导致偏差
  if(sum(target_counts) != total_sample) {
    diff <- total_sample - sum(target_counts)
    # 找到概率最大的类别调整数量
    adjust_idx <- which.max(target_dist_df$prob)
    target_counts[adjust_idx] <- target_counts[adjust_idx] + diff
  }
  
  sampled_list <- list()
  # 对每个nr.genomes类别抽取对应数量的行
  for(j in 1:nrow(target_dist_df)) {
    current_nr <- target_dist_df$nr.genomes[j]
    current_count <- target_counts[j]
    # 筛选CHR组中该nr.genomes的行
    subset_rows <- chr_data[chr_data$nr.genomes == current_nr, ]
    # 随机抽取（如果行数不足则全部抽取）
    if(nrow(subset_rows) >= current_count) {
      sampled_subset <- subset_rows[sample(nrow(subset_rows), current_count), ]
    } else {
      sampled_subset <- subset_rows
      warning(paste("CHR组中nr.genomes=", current_nr, "的行数不足，已全部抽取"))
    }
    sampled_list[[j]] <- sampled_subset
  }
  
  # 合并单次抽样结果
  sampled_df <- do.call(rbind, sampled_list)
  sampled_df$run <- as.character(i)
  # 合并到最终结果
  df.final <- rbind(df.final, sampled_df)
}

额外说明 #

• 如果CHR组中某些nr.genomes类别不存在，需要提前从目标分布中移除该类别并重新归一化概率
• 四舍五入的调整逻辑可以优化为：用floor计算基础数量后，将剩余数量随机分配到多个类别，避免只调整最大概率的类别

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