統計ER

R, EZR, SPSS, KH Coder を使ったデータ分析方法を紹介するブログ。ニッチな内容が多め

相関係数のメタアナリシスはどうやるか?

相関係数を統合したい場合はどうやるか?

相関係数を変換したり、重みがちょっと変わっているだけで、オッズ比を統合する方法とほとんど同じ。

統計ソフトRでやってみた。

引用書籍

ネタ本はこちら。4.2 相関係数から。

新版はこちら。

個々の研究の相関係数と95%信頼区間の準備

使うデータはこちら。r が相関係数。n がサンプルサイズ。z がFisherの分散安定化変換後の値。Hyperbolic arc-tangent 関数 atanh() を使うと簡単に変換できる。z に変換した後の標準誤差は  \sqrt{\frac{1}{n-3}} になる。これが絶妙。zl, zu は z の95%信頼区間下限と上限。rl, ru は zl, zuを逆変換した(元の世界に戻した)値。

r <- c(0.307,-0.01,0.300,0.119,0.194,0.248)
n <- c(107,1524,154,6165,4138,1559)
z <- atanh(r)

se <- sqrt(1/(n-3))

Za <- qnorm(1-0.05/2)
zl <- z-Za*se
zu <- z+Za*se

rl <- (exp(2*zl)-1)/(exp(2*zl)+1)
ru <- (exp(2*zu)-1)/(exp(2*zu)+1)

Fisherの分散安定化変換の逆変換はtanh()で行ける

Hyperbolic arc-tangentの逆は、hyperbolic tangentだった。なので、rl, ru の計算式はtanh()で置き換えられる。以降はtanh()を使う。

> (exp(2*zl)-1)/(exp(2*zl)+1)
[1]  0.12439197 -0.06018300  0.14890459  0.09431717  0.16450217  0.20082642
>  tanh(zl)
[1]  0.12439197 -0.06018300  0.14890459  0.09431717  0.16450217  0.20082642

個々の研究の相関係数と95%信頼区間をグラフにする

個々の研究をグラフにするには以下の通り。

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

k <- length(n)
id <- k:1
plot(r, id, ylim=c(-4, k+1), pch=" ", xlim=c(-1,1), yaxt="n", 
     ylab="Citation", xlab="Correlation coefficients")
title(" Correlation coefficient model ")
symbols(r, id, squares=sqrt(n-3), add=T, inches=0.25)
abline(v=0)
for (i in 1:k){
 j <- k-i+1
 x <- c(rl[i], ru[i])
 y <- c(j, j)
 lines(x, y, type="l")
 text(-0.5, i, j)
}

グラフはこんな風に出力される。

f:id:toukeier:20190104121705p:plain

固定効果で相関係数を統合する

 w_i = n_i -3 で重みづけ平均を取ると統合相関係数の推定値(Fisherの分散安定化変換後)。tanh()で戻すと求める点推定値になる。

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

w <- n-3
sw <- sum(w)
varz <- sum(w*z)/sw
varr <- tanh(varz)
varrl <- tanh(varz-Za*sqrt(1/sw))
varru <- tanh(varz+Za*sqrt(1/sw))
q1 <- sum(w*(z-varz)^2)
df1 <- k-1
pval1 <- 1-pchisq(q1, df1)
q2 <- varz^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)
varzd <- sum(wx*z)/swx
varrd <- tanh(varzd)
varrdl <- tanh(varzd-Za*sqrt(1/swx))
varrdu <- tanh(varzd+Za*sqrt(1/swx))
qx2 <- varzd^2*swx
pvalx2 <- 1-pchisq(qx2, df2)

固定効果と変量効果の統合値をグラフに書き入れる

固定効果(グラフ中はCombined: fixed)と変量効果(グラフ中はCombined: random)をグラフに書き入れる。

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

x <- c(varrl, varru)
y <- c(-1, -1)
lines(x, y, type="b")
x <- c(varrdl, varrdu)
y <- c(-2, -2)
lines(x, y, type="b")
abline(v=c(varr, varrd), lty=2)
text(-0.5, -1, "Combined: fixed")
text(-0.5, -2, "Combined: random")

グラフはこんな感じになる。

f:id:toukeier:20190104122328p:plain

制限付き最尤推定量を用いた方法(REML)で相関係数を統合する

Restricted Maximum Likelihood estimator(REML)を適用した方法。以下のスクリプトで繰り返し収束計算を行い解を求める。 \hat{\tau}^2 の初期値にはDerSimonian-Lairdの方法で求めた値を用いる。

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

intau <- tau2
tau <- intau
#
nrep <- 10
newt <- 1:nrep
for (i in 1:nrep){
 wb <- 1/(tau+se*se)
 zmb <- sum(wb*z)/sum(wb)
 qf <- k/(k-1)*(z-zmb)^2-se*se
 dkx <- (-1*sum(z*wb*wb)+zmb*sum(wb*wb))/sum(wb)
 qf2 <- -2*k/(k-1)*(z-zmb)*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)
zRM <- sum(z*wg)/swg
rRM <- tanh(zRM)
rRMl <- tanh(zRM-Za*sqrt(1/swg))
rRMu <- tanh(zRM+Za*sqrt(1/swg))
qx2RM <- zRM^2*swg
pvalx2RM <- 1-pchisq(qx2RM, df2)
tau2 <- tau

REMLの結果をグラフに書き入れる

# ----- REML graph -----

x <- c(rRMl, rRMu)
y <- c(-3, -3)
lines(x, y, type="b")
abline(v=rRM, lty=2)
text(-0.5, -3, "Combined: REML")

グラフはこちら。95%信頼区間は若干広くなるが、統計学的有意ではあった。

f:id:toukeier:20190104122855p:plain

結果の数値を並べて比べてみる。

点推定値は固定効果→DerSimonian-Lairdの方法→REMLとだんだんに大きくなり、95%信頼区間もだんだん広くなっていった。最終的に、統計学的有意ではあった( p_{REML} = 0.0003

> ## Fixed effects
> list(round(c(rFE=varr, LL=varrl, UL=varru, Q1=q1, df1=df1, p1=pval1, Q2=q2, df2=df2, p2=pval2),4))
[[1]]
     rFE       LL       UL       Q1      df1       p1       Q2      df2       p2 
  0.1463   0.1298   0.1627  76.8331   5.0000   0.0000 295.8090   1.0000   0.0000 

> 
> ## Random effects
> list(round(c(rDL=varrd, LL=varrdl, UL=varrdu, Q2DL=qx2, df2=df2, p2DL=pvalx2, "tau^2"=intau),4))
[[1]]
    rDL      LL      UL    Q2DL     df2    p2DL   tau^2 
 0.1741  0.0966  0.2496 19.0392  1.0000  0.0000  0.0078 

> 
> ## REML
> list(round(c(rRM=rRM, LL=rRMl, UL=rRMu, Q2RM=qx2RM, df2=df2, p2RM=pvalx2RM, "tau^2"=tau2),4))
[[1]]
    rRM      LL      UL    Q2RM     df2    p2RM   tau^2 
 0.1800  0.0831  0.2735 13.0604  1.0000  0.0003  0.0130 

まとめ

相関係数の統合を統計ソフトRでやってみた。相関係数にFisherの分散安定化変換をほどこして、 n-3 を重みにして重みづけ平均を取るのが中核だ。あとはオッズ比の統合(漸近分散法)と同じ方法であった。相関係数が関心事の研究を統合するときにはこの方法が使える。

漸近分散法によるオッズ比の統合は以下の記事をどうぞ。

toukeier.hatenablog.com