如何采用精确法正确计算疫苗效力的置信区间

方法说明 #

队列研究中疫苗效力计算的精确法，核心是基于人年发病率服从泊松分布的假设，先计算发病率比（RR）的精确置信区间，再通过公式VE = (1 - RR)*100%转换得到VE的置信区间，和单组率计算使用的Clopper-Pearson法逻辑同源，但需要额外纳入两组人年的权重参数。之前使用的简化代码省略了人年权重步骤，仅在两组人年完全相等时成立，无法适配不等规模队列场景。
关于置信水平参数：双侧α=5%的95%置信区间，等价于两侧各分配2.5%尾概率的区间，和单侧α=2.5%的置信界直接对应，参数设置会直接影响区间宽度，计算时可根据需求灵活调整。

可复现R实现代码 #

以下代码完全适配两组人年不等的场景，支持自定义单侧/双侧置信水平，无需依赖额外第三方包即可运行：

# 定义精确率比计算函数
calc_ve_exact <- function(vax_cases, unvax_cases, vax_py, unvax_py, conf.level = 0.95, side = "two.sided") {
  # 计算点估计
  rate_vax <- vax_cases / vax_py
  rate_unvax <- unvax_cases / unvax_py
  rr_point <- rate_vax / rate_unvax
  ve_point <- (1 - rr_point) * 100
  
  # 确定尾概率
  if(side == "two.sided") {
    alpha <- 1 - conf.level
    alpha_lo <- alpha/2
    alpha_hi <- alpha/2
  } else if(side == "lower") { # 输出VE单侧下界
    alpha_lo <- 1 - conf.level
    alpha_hi <- 0
  } else if(side == "upper") { # 输出VE单侧上界
    alpha_lo <- 0
    alpha_hi <- 1 - conf.level
  }
  
  total_cases <- vax_cases + unvax_cases
  # 二分法求解RR下界（对应VE上界）
  rr_lo <- uniroot(function(rr) {
    p <- (vax_py * rr) / (vax_py * rr + unvax_py)
    pbinom(vax_cases, total_cases, p) - alpha_hi
  }, interval = c(1e-8, 100))$root
  # 二分法求解RR上界（对应VE下界）
  rr_hi <- uniroot(function(rr) {
    pbinom(vax_cases - 1, total_cases, p = (vax_py * rr)/(vax_py * rr + unvax_py), lower.tail = F) - alpha_lo
  }, interval = c(1e-8, 100))$root
  
  # 转换为VE区间
  ve_lo <- (1 - rr_hi) * 100
  ve_hi <- (1 - rr_lo) * 100
  
  return(list(
    VE = round(ve_point, 1),
    CI_lower = round(ve_lo, 1),
    CI_upper = round(ve_hi, 1)
  ))
}

# 代入验证数据集计算
result <- calc_ve_exact(
  vax_cases = 176,
  unvax_cases = 221,
  vax_py = 11793,
  unvax_py = 5809,
  conf.level = 0.95,
  side = "two.sided"
)
print(result)

验证结果 #

代入提供的示例数据集运行上述代码，输出结果为：

• 疫苗效力点估计：60.8%
• 双侧95%置信区间：52.0% – 68.0%
  和参考结果完全匹配。

参数调整说明 #

• 若需计算单侧α=2.5%的VE下界，将参数设置为conf.level = 0.975, side = "lower"即可
• 若需计算单侧α=2.5%的VE上界，将参数设置为conf.level = 0.975, side = "upper"即可
• 其他置信水平需求直接修改conf.level参数即可适配

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