メタアナリシス―オッズ比をマンテル・ヘンツェルの方法で統合する方法

オッズ比をマンテル・ヘンツェル（Mantel-Haenszel）の方法で統合する方法の解説。

マンテル・ヘンツェル法は、２ｘ２の分割表を統合する方法で、層別解析の方法。

メタアナリシスに応用する方法。

＞＞もう統計で悩むのを終わりにしませんか？

↑1万人以上の医療従事者が購読中

オッズ比のメタアナリシス―マンテル・ヘンツェル法の準備
オッズ比メタアナリシスのグラフ表示
オッズ比のメタアナリシスをmetaforパッケージで行う
まとめ
関連記事
参考書籍

オッズ比のメタアナリシス―マンテル・ヘンツェル法の準備

データは以下の通り。

a <- c(3,7,5,102,28,4,98,60,25,138,64,45,9,57,25,65,17)
n1 <- c(38,114,69,1533,355,59,945,632,278,1916,873,263,291,858,154,1195,298)
c <- c(3,14,11,127,27,6,152,48,37,188,52,47,16,45,31,62,34)
n0 <- c(39,116,93,1520,365,52,939,471,282,1921,583,266,293,883,147,1200,309)

dat <- data.frame(a,n1,c,n0)

まず分析の準備。

下図と同じように変数名を作成。

オッズ比は、 $\frac{ad}{bc}$ で計算。

ai <- dat$a
bi <- dat$n1 - dat$a
ci <- dat$c
di <- dat$n0 - dat$c
tn <- dat$n1 + dat$n0
n1 <- dat$n1
n0 <- dat$n0
m1 <- ai + ci
m0 <- bi + di

各研究のオッズ比、標準誤差、対数オッズ比、95%信頼区間下限・上限、各研究の重みを計算する。

or <- ai*di/bi/ci
se <- sqrt(tn/bi/ci)
lgor <- log(or)
low <- exp(lgor-1.96*se)
upp <- exp(lgor+1.96*se)
w <- 1/se/se

各研究のオッズ比、95%信頼区間を並べてみる。

> round(cbind(ORi=or, LLi=low, ULi=upp),4)
         ORi    LLi    ULi
 [1,] 1.0286 0.1920 5.5103
 [2,] 0.4766 0.2211 1.0274
 [3,] 0.5824 0.2274 1.4912
 [4,] 0.7818 0.6064 1.0079
 [5,] 1.0719 0.6125 1.8760
 [6,] 0.5576 0.1789 1.7377
 [7,] 0.5991 0.4726 0.7594
 [8,] 0.9244 0.6241 1.3692
 [9,] 0.6543 0.4051 1.0568
[10,] 0.7155 0.5799 0.8826
[11,] 0.8078 0.5610 1.1633
[12,] 0.9618 0.6162 1.5015
[13,] 0.5525 0.2730 1.1184
[14,] 1.3252 0.8614 2.0387
[15,] 0.7252 0.4236 1.2416
[16,] 1.0558 0.7349 1.5169
[17,] 0.4893 0.2986 0.8020

統合オッズ比を推定し、統合オッズ比の95%信頼区間を計算する。

#----- Mantel-Haenszel Odds Ratio -----
mhor <- sum(ai*di/tn)/sum(bi*ci/tn)

#----- Variance of Mantel-Haenszel Odds Ratio -----
P <- (ai+di)/tn
Q <- (bi+ci)/tn
R <- ai*di/tn
S <- bi*ci/tn

mhv <- (sum(P*R)/sum(R)^2 + sum(P*S+Q*R)/sum(R)/sum(S) + sum(Q*S)/sum(S)^2)/2

#----- 95% confidence interval -----
mhorl <- exp(log(mhor)-1.96*sqrt(mhv))
mhoru <- exp(log(mhor)+1.96*sqrt(mhv))

list(round(c(ORmh=mhor, LL=mhorl, UL=mhoru),4))

統合オッズ比と95%信頼区間は以下の通り。

> list(round(c(ORmh=mhor, LL=mhorl, UL=mhoru),4))
1
  ORmh     LL     UL 
0.7816 0.7058 0.8655

均質性の検定（重みは漸近分散法と同じ計算方法）、有意性の検定（Petoの方法と同じ）を行う。

k <- length(ai)

w1 <- 1/(1/ai+1/bi+1/ci+1/di)

q1 <- sum(w1*(lgor-log(mhor))^2)
df1 <- k-1
pval1 <- 1-pchisq(q1, df1)

o <- ai
e <- (ai+bi)*(ai+ci)/tn
v <- ((ai+bi)*(ci+di)/tn)*((ai+ci)*(bi+di)/tn)/(tn-1)

q2 <- (abs(sum(o-e))-0.5)^2/sum(v)
df2 <- 1
pval2 <- 1-pchisq(q2, df2)

list(round(c(Q1=q1, df1=df1,P1=pval1, Q2=q2, df2=df2, P2=pval2),4))

結果はこちら。

> list(round(c(Q1=q1, df1=df1,P1=pval1, Q2=q2, df2=df2, P2=pval2),4))
1
     Q1     df1      P1      Q2     df2      P2 
21.4811 16.0000  0.1607 22.2667  1.0000  0.0000

マンテル・ヘンツェルでの均質性の検定は、漸近分散法の重みを使った上述の方法よりもBreslow-Day検定のほうが適切。

DescToolsパッケージの BreslowDayTest() で実施できる。

２ｘ２の分割表の準備が手間がかかる。

array()で準備している。

all <- c(0,0,0,0)
for (i in 1:17){
  new <- c(ai[i], ci[i], bi[i], di[i])
  all <- c(all, new)
  }
out <- all[5:72]

out.bd <- array(out,
 dim=c(2,2,17),
 dimnames=list(exposure=c("Beta brokade", "Control"),
                 event=   c("Yes","No"),
                 site=    c(paste("site", 1:17, sep=""))
               )
)

library(DescTools)
BreslowDayTest(out.bd)
BreslowDayTest(out.bd, correct=T)

Breslow-Day検定の結果はこちら。

連続性補正の有無では結果変わらず。

> BreslowDayTest(out.bd)

        Breslow-Day test on Homogeneity of Odds Ratios

data:  out.bd
X-squared = 21.714, df = 16, p-value = 0.1527

> BreslowDayTest(out.bd, correct=T)

        Breslow-Day Test on Homogeneity of Odds Ratios (with Tarone correction)

data:  out.bd
X-squared = 21.714, df = 16, p-value = 0.1527

オッズ比メタアナリシスのグラフ表示

個々の研究のオッズ比95%信頼区間のグラフに、マンテル・ヘンツェルの方法で推定した、統合オッズ比と95%信頼区間を乗せて図示する。

# ------------- individual graph ----------------

id <- k:1

plot(exp(lgor), id, ylim=c(-2,20),
log="x", xlim=c(0.1,10), yaxt="n", pch="",
ylab="Citation", xlab="Odds ratio")

title(main=" Mantel-Haenszel method ")

symbols(exp(lgor), id, squares=sqrt(tn),
add=TRUE, inches=0.25)

for (i in 1:k){
j <- k-i+1
x <- c(low[i], upp[i])
y <- c(j, j)
lines(x, y, type="l")
text(0.1, i, j)
}

# -------------- Combined graph --------------

mhorx <- c(mhorl, mhoru)
mhory <- c(-1, -1)
lines(mhorx, mhory, type="o", lty=1, lwd=2)
abline(v=c(mhor), lty=2)
abline(v=1)
text(0.3, -1, "Combined")

＞＞もう統計で悩むのを終わりにしませんか？

↑1万人以上の医療従事者が購読中

オッズ比のメタアナリシスをmetaforパッケージで行う

metaforパッケージのrma.mh()を使った方法。

library(metafor)

rma.mh.res <- rma.mh(ai=a, bi=n1-a, ci=c, di=n0-c, data=dat)

rma.mh.res

結果は、estimate, ci.lb, ci.ub。
均質性の検定は、Test for Heterogeneity。
有意性の検定は、Cochran-Mantel-Haenszel Test。
Breslow-Day検定は、Tarone's Test for Heterogeneity。

それぞれチェック。

> rma.mh.res

Fixed-Effects Model (k = 17)

Test for Heterogeneity: 
Q(df = 16) = 21.4811, p-val = 0.1607

Model Results (log scale):

estimate      se     zval    pval    ci.lb    ci.ub
 -0.2464  0.0520  -4.7370  <.0001  -0.3484  -0.1445

Model Results (OR scale):

estimate   ci.lb   ci.ub
  0.7816  0.7058  0.8655

Cochran-Mantel-Haenszel Test:    CMH = 22.2667, df = 1,  p-val < 0.0001
Tarone's Test for Heterogeneity: X^2 = 21.7143, df = 16, p-val = 0.1527