从 plink 输出计算优势比置信区间？

机器算法验证 置信区间 遗传学 优势比
2022-03-24 17:34:28

我有plink单倍型分析的输出，但是我没有原始数据。这是基于 Haplotype 的 GLM 关联测试的输出：

SNP1    SNP2    HAPLOTYPE   F   OR  STAT    P
rs1 rs2 22  0.00992 4.23    61.5    4.43E-15
rs1 rs2 12  0.038   1.02    0.217   0.642
rs1 rs2 21  0.00015 5.22E-10    453 1.77E-100
rs1 rs2 11  0.952   0.762   22.9    1.73E-06

以下是 plink 中每一列的解释：

        SNP1    SNP ID of left-most (5') SNP
        SNP2    SNP ID of left-most (3') SNP
   HAPLOTYPE    Haplotype 
           F    Frequency in sample
          OR    Estimated odds ratio
        STAT    Test statistic (T from Wald test)
           P    Asymptotic p-value

问题：根据以上输出，是否可以计算出OR95% 的置信区间？

1个回答

对于置信区间的计算，您需要效果的标准误差，但这些在输出中不可用。但是，可以从 Wald 统计量和优势比估计标准误。

计算如下：

  1. 从优势比取自然对数。这为您提供了逻辑模型的 beta。例如对于表格的第一行： beta=ln(4.23)=1.442
  2. 贝塔的标准误差是通过将贝塔除以沃尔德统计量 (STAT) 的平方根来计算的。然后取结果的绝对值。同样，对于表格的第一行：se=1.442/sqrt(61.5)=0.183。
  3. β 的 95% 置信区间为 beta+/-1.96*se。常数 1.96 来自正态分布。同样，对于第一行数据：1.442-1.96*0.183 ... 1.442+1.96*0.183 = 1.081...1.802。
  4. 最后，您需要将 beta 的置信区间更改为优势比的置信区间。这只需对 beta 的置信区间取幂即可。对于第一行数据：2.71828^1.081 = 2.949 和 2.71828^1.802 = 6.065。

因此，表格第一行的优势比为 4.23，其 95% 置信区间为 2.949-6.065。由于置信区间不包括 1，因此结果具有统计显着性。由于 PLINK 输出的四舍五入，结果可能会出现错误。

这种计算可以在例如 Excel 中实现，但下面也是一个做同样事情的 R 函数（以防万一你也使用 R）。

# The data
or<-structure(list(SNP1 = structure(c(1L, 1L, 1L, 1L), .Label = "rs1", class = "factor"), 
SNP2 = structure(c(1L, 1L, 1L, 1L), .Label = "rs2", class = "factor"), 
HAPLOTYPE = c(22L, 12L, 21L, 11L), F = c(0.00992, 0.038, 
0.00015, 0.952), OR = c(4.23, 1.02, 5.22e-10, 0.762), STAT = c(61.5, 
0.217, 453, 22.9), P = c(4.43e-15, 0.642, 1.77e-100, 1.73e-06
)), .Names = c("SNP1", "SNP2", "HAPLOTYPE", "F", "OR", "STAT", 
"P"), class = "data.frame", row.names = c(NA, -4L))

# The function
orci<-function(or) {
   or$beta<-log(or$OR)
   or$se<-abs(or$beta/sqrt(or$STAT))
   or$lower<-or$beta-1.96*or$se
   or$upper<-or$beta+1.96*or$se
   or$LOWER<-exp(or$lower)
   or$UPPER<-exp(or$upper)
   or$res<-paste(or$OR, " (", round(or$LOWER, 3), "-", round(or$UPPER, 3), ")", sep="")
   return(or)
}

# The calculation
orci(or)

# The result
#SNP1 SNP2 HAPLOTYPE       F       OR    STAT         P         beta         se        lower       upper        LOWER        UPPER                  res
#1  rs1  rs2        22 0.00992 4.23e+00  61.500  4.43e-15   1.44220199 0.18390288   1.08175235   1.8026516 2.949844e+00 6.065710e+00    4.23 (2.95-6.066)
#2  rs1  rs2        12 0.03800 1.02e+00   0.217  6.42e-01   0.01980263 0.04251018  -0.06351733   0.1031226 9.384579e-01 1.108627e+00   1.02 (0.938-1.109)
#3  rs1  rs2        21 0.00015 5.22e-10 453.000 1.77e-100 -21.37335353 1.00420775 -23.34160072 -19.4051063 7.292419e-11 3.736538e-09 0.000000000522 (0-0)
#4  rs1  rs2        11 0.95200 7.62e-01  22.900  1.73e-06  -0.27180872 0.05679965  -0.38313603  -0.1604814 6.817202e-01 8.517337e-01  0.762 (0.682-0.852)
其它你可能感兴趣的问题
上一篇ARIMA 模型中的自由度 下一篇差分设计如何解释时间自相关