在R中使用prop.test()获取多分组的比例与置信区间

需求说明 #

• 按县统计多种药物的阳性率，真实数据包含更多药物与区县
• NA值代表未送检，不计入样本总数
• 必须使用prop.test(x, y, conf.level = 0.95)计算每个县的阳性比例及95%置信区间

模拟数据 #

county <- c("Erie", "Orange", "Erie", "Orange", "Erie", "Orange", "Erie", "Orange", "Erie", "Orange", "Erie", "Orange")
drug1 <- c("Positive", "Negative", "Negative", "Positive", NA, "Negative", "Positive", "Negative", "Positive", "Negative", "Positive", "Negative")
drug2 <- c("Positive", NA, "Negative", "Negative", "Negative", "Positive", "Positive", NA, "Negative", "Negative", "Negative", "Positive")
data <- data.frame(county, drug1, drug2)

问题 #

不能直接用以下简单的group_by()+summarise()方式计算比率，需要改用prop.test()实现：

data |>
  group_by(county) |>
  summarise(drug1 = sum(drug1=="Positive", na.rm = TRUE)/
              sum(!is.na(drug1)),
            drug2 = sum(drug2=="Positive", na.rm = TRUE)/
              sum(!is.na(drug2))
            ) |>
  ungroup()

解决方案 #

可以通过数据格式整理+分组调用prop.test实现，以下是两种可行方法：

方法1：宽转长后分组处理 #

先将宽格式数据转为长格式，再按县和药物分组，调用prop.test提取比例与置信区间：

library(dplyr)
library(tidyr)

# 转为长格式并剔除未送检的NA值
data_long <- data |>
  pivot_longer(cols = starts_with("drug"), names_to = "drug", values_to = "result") |>
  filter(!is.na(result))

# 分组计算阳性率及置信区间
result_df <- data_long |>
  group_by(county, drug) |>
  summarise(
    pos_count = sum(result == "Positive"),
    total_count = n(),
    prop = prop.test(pos_count, total_count, conf.level = 0.95)$estimate,
    conf_low = prop.test(pos_count, total_count, conf.level = 0.95)$conf.int[1],
    conf_high = prop.test(pos_count, total_count, conf.level = 0.95)$conf.int[2],
    .groups = "drop"
  )

print(result_df)

方法2：保留宽格式批量处理 #

如果希望结果保持原数据的宽格式，可使用across()批量处理每个药物列：

library(dplyr)

result_wide <- data |>
  group_by(county) |>
  summarise(
    across(
      starts_with("drug"),
      ~{
        pos <- sum(.x == "Positive", na.rm = TRUE)
        total <- sum(!is.na(.x))
        test_res <- prop.test(pos, total, conf.level = 0.95)
        tibble(
          prop = test_res$estimate,
          conf_low = test_res$conf.int[1],
          conf_high = test_res$conf.int[2]
        )
      },
      .names = "{.col}_{.value}"
    ),
    .groups = "drop"
  )

print(result_wide)

结果说明 #

两种方法都会输出：

• 县名（county）
• 药物标识（或对应列名）
• 阳性率（prop）
• 95%置信区间下限（conf_low）与上限（conf_high）

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