如何用R循环或数组批量计算多变量的风险比（Risk Ratio）

问题背景 #

现有野餐疫情后的患者隐私保护虚构数据集，包含7名患者的患病状态（sick）及5种食物的食用情况（apple、banana、strawberry、yogurt、sandwich），所有分类变量取值均为yes/no。手动逐个计算单食物风险比效率极低，需实现批量处理所有食物变量，输出规整统一的风险比结果。

数据集示例：

# A tibble: 7 × 7
  patient  sick  apple banana strawberry yogurt sandwich
  <chr>    <chr> <chr> <chr>  <chr>      <chr>  <chr>   
1 patient1 yes   no    no     yes        no     no      
2 patient2 yes   yes   no     no         yes    no      
3 patient3 no    yes   yes    yes        yes    yes     
4 patient4 yes   no    no     yes        yes    no      
5 patient5 no    no    no     no         yes    no      
6 patient6 no    no    no     no         no     no      
7 patient7 yes   no    no     no         yes    yes   

此前手动计算apple风险比的流程为分组统计频数、构建2×2矩阵、调用epitools::riskratio.wald()，但该方式无法适配大量食物的批量处理需求。

批量处理解决方案 #

借助tidyverse生态中的purrr迭代工具，配合自定义函数实现批量计算，自动整理每个食物的风险比、置信区间及P值为统一表格。

步骤1：加载依赖包 #

library(tidyverse)
library(epitools)

步骤2：定义批量计算函数 #

自定义函数接收食物变量名，完成频数统计、矩阵构建、风险比计算，并提取关键结果：

calculate_risk_ratio <- function(food_var) {
  # 分组统计患病与食用情况的频数
  freq_table <- picnic %>%
    group_by(sick, .data[[food_var]]) %>%
    summarise(n = n(), .groups = "drop") %>%
    pivot_wider(names_from = sick, values_from = n, values_fill = 0)
  
  # 构建符合epitools要求的2×2矩阵（行：未食用/食用；列：未患病/患病）
  epi_matrix <- matrix(c(freq_table$no[freq_table[[food_var]] == "no"],
                         freq_table$yes[freq_table[[food_var]] == "no"],
                         freq_table$no[freq_table[[food_var]] == "yes"],
                         freq_table$yes[freq_table[[food_var]] == "yes"]),
                       nrow = 2, byrow = TRUE,
                       dimnames = list(Predictor = c("Not Exposed", "Exposed"),
                                       Outcome = c("Not Sick", "Sick")))
  
  # 计算风险比
  rr_result <- riskratio.wald(epi_matrix)
  
  # 提取并整理关键结果
  tibble(
    food = food_var,
    risk_ratio = rr_result$measure["Exposed", "estimate"],
    ci_lower = rr_result$measure["Exposed", "lower"],
    ci_upper = rr_result$measure["Exposed", "upper"],
    p_value_chi = rr_result$p.value["Exposed", "chi.square"],
    p_value_fisher = rr_result$p.value["Exposed", "fisher.exact"]
  )
}

步骤3：批量处理所有食物变量 #

指定目标食物变量列表，用map_dfr()迭代计算并合并结果：

# 定义需要处理的食物变量
food_vars <- c("apple", "banana", "strawberry", "yogurt", "sandwich")

# 批量计算并合并为统一表格
rr_summary <- map_dfr(food_vars, calculate_risk_ratio)

# 查看结果
rr_summary

示例输出 #

运行后将得到规整的结果表格：

# A tibble: 5 × 6
  food       risk_ratio ci_lower ci_upper p_value_chi p_value_fisher
  <chr>           <dbl>    <dbl>    <dbl>       <dbl>          <dbl>
1 apple           0.833    0.175     3.96       0.809              1
2 banana          0.5      0.0556    4.50       0.536              1
3 strawberry      2.67     0.713     9.97       0.147              0.2
4 yogurt          1.2      0.387     3.71       0.737              1
5 sandwich        1.5      0.333     6.75       0.589              1

说明 #

• 矩阵构建严格遵循riskratio.wald()的输入要求，确保风险比计算逻辑正确
• map_dfr()自动将单食物结果合并为tibble，方便后续分析或导出
• 结果包含风险比、95%置信区间及卡方检验、Fisher精确检验的P值，满足流行病学分析需求

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