如何通过已知统计量计算meta分析所需的分组及结局事件数值

需求背景 #

我在使用metabin库开展meta分析时，需要获取event.e、n.e、event.c、n.c四个数值，但目前没有这些原始数据，希望通过已知的Odds Ratio、logOR、95% CIs上下限、患病/未患病总人数（即待求解2x2四格表的列合计值）、Standard Error、Sample Variance计算得到目标值。

已尝试的操作 #

我曾尝试使用escalc函数实现需求，但未找到对应参数，尝试过的代码如下：

escalc(measure="PR", x1i=diseasePOS, x2i=diseaseNEG, ni=(diseasePOS+diseaseNEG), yi=logOR , vi=sampVar, sei=stdEr, data=so.data)

当前数据集构造代码如下：

OR <- c(1.8,1.41,1.07,1.906)
l95ci <- c(0.9,1.1,0.415,0.623)
u95ci <- c(3.3,1.8,2.811,6.285)
diseasePOS <- c(76,72,64,54)
diseaseNEG <- c(10502,7941,54,54)
logOR <- c(0.5878,0.3436,0.0677,0.6450)
stdEr <- c(0.3314,0.12563,0.48801,0.58963)
sampVar <- c(0.1099,0.0158,0.2382,0.3477)

so.data <- data.frame(OR,l95ci, u95ci, diseasePOS, diseaseNEG, logOR, stdEr, sampVar  )
so.data

预期目标 #

得到对应数值填充如下四格表的XXXX位置，即escalc中的ai、bi、ci、di：

实现方案 #

不需要适配escalc的参数，直接通过OR和logOR方差的计算公式联立求解即可，实现代码如下：

# 编写四格表求解函数
impute_2x2 <- function(OR, var_logOR, pos_total, neg_total) {
  # 遍历干预组阳性数的合理取值范围
  for(a in 1:(pos_total-1)) {
    c <- pos_total - a
    # 遍历干预组总人数的合理取值范围
    min_n1 <- max(a, neg_total - (neg_total - (pos_total + neg_total - (a + c))))
    max_n1 <- min(a + neg_total, pos_total + neg_total - c)
    for(n1 in min_n1:max_n1) {
      b <- n1 - a
      d <- neg_total - b
      if(c < 1 | d < 1) next
      # 计算当前取值对应的OR和方差，和已知值匹配则返回结果
      calc_or <- (a * d) / (b * c)
      calc_var <- 1/a + 1/b + 1/c + 1/d
      if(abs(calc_or - OR) < 1e-2 & abs(calc_var - var_logOR) < 1e-2) {
        return(c(ai = a, bi = b, ci = c, di = d))
      }
    }
  }
  return(rep(NA, 4))
}

# 批量处理数据集中的每一行
so.data[, c("ai", "bi", "ci", "di")] <- t(mapply(
  impute_2x2, 
  OR = so.data$OR, 
  var_logOR = so.data$sampVar, 
  pos_total = so.data$diseasePOS, 
  neg_total = so.data$diseaseNEG
))

# 输出求解得到的四格表数据
so.data[, c("ai", "bi", "ci", "di")]

结果说明 #

• 输出的ai对应干预组患病数、bi对应干预组未患病数、ci对应对照组患病数、di对应对照组未患病数，可直接代入metabin或escalc开展后续meta分析。
• 若部分行返回NA，说明已知参数存在四舍五入误差，可适当放宽代码中OR和方差的匹配阈值即可得到近似解。

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