如何基于R数据集生成按termx分组的每万人患病率及95%Wilson置信区间表

DevOps 小满

2026-9-25

实现步骤

你可以借助R的tidyverse生态搭配置信区间计算包完成需求，完整操作代码如下：

1. 安装并加载依赖包

未安装过所需包先执行安装：

install.packages(c("tidyverse", "epitools", "knitr"))

加载包：

library(tidyverse)
library(epitools)
library(knitr)

2. 调整孕周分组的因子顺序（可选）

如果原数据中termx的因子水平顺序和你需要的孕周从小到大顺序不一致，先运行这一步调整列顺序：

b.pcalc <- b.pcalc %>%
  mutate(termx = factor(termx, levels = c(
    "<28 weeks", "28 to <32 weeks", "32 to <37 weeks", 
    "37-38 weeks", "39-40 weeks", "41-<42 weeks", "42+ weeks"
  )))

3. 数据重塑+统计计算

# 宽表转长表，统一处理所有出生缺陷变量
b.pcalc_long <- b.pcalc %>%
  pivot_longer(
    cols = -termx,
    names_to = "Birth defects",
    values_to = "status"
  ) %>%
  mutate(status = as.numeric(as.character(status)))

# 分组计算患病率和95% Wilson置信区间
stat_result <- b.pcalc_long %>%
  group_by(termx, `Birth defects`) %>%
  summarise(
    n_case = sum(status),
    n_total = n(),
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(
    # 计算每万人患病率
    prev_10k = n_case / n_total * 10000,
    # 计算Wilson置信区间并转成每万人单位
    ci = list(binomial.wilson(n_case, n_total, conf.level = 0.95)),
    ci_low = ci[1, "lower"] * 10000,
    ci_high = ci[1, "upper"] * 10000,
    # 拼接成要求的展示格式，保留1位小数
    show_text = sprintf("%.1f (%.1f,%.1f)", prev_10k, ci_low, ci_high)
  ) %>%
  # 转回宽表匹配输出格式
  pivot_wider(
    id_cols = `Birth defects`,
    names_from = termx,
    values_from = show_text
  )

如果你不想使用epitools包，也可以用binom包的binom.confint函数计算Wilson置信区间，只需要把计算ci的部分替换为：
ci = list(binom.confint(n_case, n_total, conf.level = 0.95, method = "wilson"))
ci_low = ci$lower * 10000
ci_high = ci$upper * 10000