統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

平均値の差のメタアナリシスはどうやるか?

kaiseki daiko banner

いままでオッズ比の統合方法は取り上げてきたが、平均値の差の統合は取り上げてこなかった。

先日ちょっと必要になって、勉強したのでアップしてみる。

いつものメタ・アナリシス入門のS言語からR言語への移植だ。

移植というか、そのままだが。。。

引用書籍

丹後敏郎著 メタ・アナリシス入門

3.2 平均値と標準偏差

新版もある。

統合するのは平均値の差の検定がエンドポイントの解析である試験

データは以下の通り。

mが平均、sが標準偏差、nがサンプルサイズ。

n1 <- c(155,31,75,18,8,57,34,110,60)
m1 <- c(55.0,27.0,64.0,66.0,14.0,19.0,52.0,21.0,30.0)
s1 <- c(47.0,7.0,17.0,20.0,8.0,7.0,45.0,16.0,27.0)
n0 <- c(156,32,71,18,13,52,33,183,52)
m0 <- c(75.0,29.0,119.0,137.0,18.0,18.0,41.0,31.0,23.0)
s0 <- c(64.0,4.0,29.0,48.0,11.0,4.0,34.0,27.0,20.0)

並べてみると以下の通り。

> cbind(n1,m1,s1,n0,m0,s0)
       n1 m1 s1  n0  m0 s0
 [1,] 155 55 47 156  75 64
 [2,]  31 27  7  32  29  4
 [3,]  75 64 17  71 119 29
 [4,]  18 66 20  18 137 48
 [5,]   8 14  8  13  18 11
 [6,]  57 19  7  52  18  4
 [7,]  34 52 45  33  41 34
 [8,] 110 21 16 183  31 27
 [9,]  60 30 27  52  23 20

それぞれの試験の平均値の差と95%信頼区間の算出

# ---- difference in means -----

ad <- m1-m0
s <- ((n1-1)*s1^2+(n0-1)*s0^2)/(n1-1+n0-1)
se <- sqrt((1/n1+1/n0)*s)
tn <- n1+n0
low <- ad-1.96*se
upp <- ad+1.96*se

それぞれの試験の平均値の差と95%信頼区間のグラフ化

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

k <- length(n1)
id <- k:1
plot(ad, id, ylim=c(-4,10), pch=" ", 
     xlim=c(-100,100), yaxt="n",
     ylab="Citation", xlab="Difference in means")
title(main=" Mean difference model ")
symbols(ad, id, squares=sqrt(tn), add=T, inches=0.25)
abline(v=0)
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(50, i, j)
}

f:id:toukeier:20190102184355p:plain

固定効果の計算

固定効果とは、研究の結果が似ていると思えるときに使える方法。

# ----- fixed effects -----

w <- 1/se/se
sw <- sum(w)
adm <- sum(ad*w)/sw
adl <- adm-1.96*sqrt(1/sw)
adu <- adm+1.96*sqrt(1/sw)
q1 <- sum(w*(ad-adm)^2)
df1 <- k-1
pval1 <- 1-pchisq(q1, df1)
q2 <- adm^2*sw
df2 <- 1
pval2 <- 1-pchisq(q2, df2)

変量効果の計算(Dersimonian-Lairdの方法)

変量効果とは、研究同士に無視できない違いがあるときに使う方法。

# ----- random effects -----

tau2 <- (q1-(k-1))/(sw-sum(w*w)/sw)
tau2 <- max(0, tau2)
wx <- 1/(tau2+se*se)
swx <- sum(wx)
admx <- sum(ad*wx)/swx
adxl <- admx-1.96*sqrt(1/swx)
adxu <- admx+1.96*sqrt(1/swx)
qx2 <- admx^2*swx
pvalx2 <- 1-pchisq(qx2, df2)

統合平均値と95%信頼区間を書き入れる

先ほどのグラフに統合平均値と95%信頼区間を書き入れる。

# ----- graph -----

x <- c(adl, adu)
y <- c(-1, -1)
lines(x, y, type="b")
x <- c(adxl, adxu)
y <- c(-2, -2)
lines(x, y, type="b")
abline(v=c(adm, admx), lty=2)
text(60, -1, "Combined: fixed")
text(60, -2, "Combined: random")

f:id:toukeier:20190102184422p:plain

結果を出力する

outは、固定効果の統合平均値と95%信頼区間、outq1は均質性の検定(検定統計量、自由度、P値)、outq2は有意性の検定(検定統計量、自由度、P値)、outRは変量効果の統合平均値と95%信頼区間、out2Rは変量効果の有意性の検定、tau2は {\tau}^2の推定値。

> # ----- output variables -----
> out
[1] -3.4939 -5.0265 -1.9612
> outq1
[1] 241.059   8.000   0.000
> outq2
[1] 19.963897  1.000000  0.000008
> outR
[1] -14.0972 -24.4456  -3.7488
> outq2R
[1] 7.129061 1.000000 0.007584
> tau2
[1] 218.7216

list()を使うと、少し見やすくできる。

> ## fixed effects
> list(round(c(adFE=adm, LL=adl, UL=adu, Q1=q1, df1=df1,p1=pval1, Q2=q2, df2=df2, p2=pval2),4))
[[1]]
    adFE       LL       UL       Q1      df1       p1       Q2      df2       p2 
 -3.4939  -5.0265  -1.9612 241.0590   8.0000   0.0000  19.9639   1.0000   0.0000 

> 
> ## random effects
> list(round(c(adDL=admx, LL=adxl, UL=adxu, Q2DL=qx2, df2=df2, p2DL=pvalx2, "tau^2"=tau2),4))
[[1]]
    adDL       LL       UL     Q2DL      df2     p2DL    tau^2 
-14.0972 -24.4456  -3.7488   7.1291   1.0000   0.0076 218.7216 

より柔軟な変量モデル―制限付き最尤推定量(REML)

変量モデルのDersimonian-Lairdの方法は簡便な方法だが、より柔軟で理論的に自然な方法が制限付き最尤推定量(REML, restricted maximum likelihood estimator)だ。現代はPCを使って反復計算が簡単にできるようになったので、REMLの適用が望ましい。Newton-Raphson法を用いて繰り返し収束計算で {\tau}^2の推定値を求めている。

# ----- REML method -----

intau <- tau2
tau <- intau
#
nrep <- 10
newt <- 1:nrep
for (i in 1:nrep){
 wb <- 1/(tau+se*se)
 admb <- sum(ad*wb)/sum(wb)
 qf <- k/(k-1)*(ad-admb)^2-se*se
 dkx <- (-1*sum(ad*wb*wb)+admb*sum(wb*wb))/sum(wb)
 qf2 <- -2*k/(k-1)*(ad-admb)*dkx
 h <- sum(wb*wb*(qf-tau))
 dh <- sum(-2*wb*wb*wb*(qf-tau)+wb*wb*(qf2-1))
 newt[i] <- tau - h/dh
 rel <- abs((newt[i]-tau)/tau)
 tau <- newt[i]
}
wg <- 1/(tau+se*se)
swg <- sum(wg)
adRM <- sum(ad*wg)/swg
adRMl <- adRM-1.96*sqrt(1/swg)
adRMu <- adRM+1.96*sqrt(1/swg)
qx2RM <- adRM^2*swg
pvalx2RM <- 1-pchisq(qx2RM, df2)
x <- c(adRMl, adRMu)
y <- c(-3, -3)
lines(x, y, type="b")
abline(v=adRM, lty=2)
text(60, -3, "Combined: REML")
tau2 <- tau
outRM <- round(c(adRM, adRMl, adRMu), 4)
outq2RM <- round(c(qx2RM, df2, pvalx2RM), 6)

結果は、以下の通り。Dersimonian-Lairdの方法よりもさらに95%信頼区間が広がっていて、統計学的有意でなくなっている。

> outRM
[1] -15.1213 -32.6671   2.4244
> outq2RM
[1] 2.853310 1.000000 0.091186
> tau2
[1] 684.829

> list(round(c(adRM=adRM, LL=adRMl, UL=adRMu, Q2RM=qx2RM, df2=df2, p2RM=pvalx2RM, "tau^2"=tau2),4))
[[1]]
    adRM       LL       UL     Q2RM      df2     p2RM    tau^2 
-15.1213 -32.6671   2.4244   2.8533   1.0000   0.0912 684.8290 

REMLの結果をグラフに足すと以下のようになる。 f:id:toukeier:20190102185519p:plain

まとめ

試験のエンドポイントが平均値の差の検定である試験を統合する方法を書籍に沿って実施してみた。

図を見ると明らかに試験同士が大きく異なっていて、無視できない差があると理解できる。

試験同士に無視できない差がある場合は、変量効果モデル(REML)を使うことが望ましい。