【学習動機】
メタアナリシスの本や記事を読んでも, 実践例は全て過去のメタアナリシスの例をやるだけで, 人工データで手法が上手くいっていることの確認をしていることがほとんどない. 自分の頭の整理も兼ねて, これから気になったことを少しずつ試したり, 調べたりしていく. コードはRで書く. ここでは, パッケージの使い方の紹介などは基本的にはしないつもりだ.
今回は, 離散データのケース. 2*2の分割表のリスク比(効果量:Effect Size)を推定して統合する実験をする.
丹後『新版 メタ・アナリシス入門 ─エビデンスの統合をめざす統計手法─』のp75あたりを主に参考にしている. これを読みながら目を通すと分かりやすいかもしれない. 他にもメタアナリシスの本(和書はほとんどないので, だいたい洋書)を読んでいるところで, そのうち, どの本がどのような感じか書くつもりだ.
【学習内容】
パッケージの読み込みと乱数シードの決定.
library(tidyverse)
set.seed(55555)2*2分割表データを固定効果モデルで生成(リスク比の推定)
まずは, 真のリスク比を設定し, 固定効果モデルを真のモデルとしてデータを生成して, 固定効果モデルでパラメータの推定をし, 均質性の検定を行う.
真のリスク比を0.8, 統合する研究数を10, 各研究のTreat群とControl群の人数を30,
Control群のイベント発生率は0.3, Treat群は0.3*0.8=0.24とする.
イメージとしては, 治療を行うとイベントが起きにくくなるということ.
ES_RR <- 0.8
Nstudy <- 10
nT <- nC <- 30
# Treat群とControl群のイベントの起こる確率
pC <- rep(0.3, Nstudy)
pT <- pC * ES_RR各研究について, イベント数データを生成する.
# Treatment群のイベント数
XT <- rbinom(Nstudy, nT, prob = pT)
# Comparisonのイベント数
XC <- rbinom(Nstudy, nC, prob = pC)
df <-
data.frame(
nT = nT,
nC = nC,
XT = XT,
XC = XC
)
df## nT nC XT XC
## 1 30 30 11 11
## 2 30 30 8 10
## 3 30 30 6 10
## 4 30 30 11 12
## 5 30 30 4 12
## 6 30 30 6 9
## 7 30 30 3 11
## 8 30 30 7 8
## 9 30 30 4 8
## 10 30 30 8 8リスク比の推定値, 対数リスク比推定値, 標準誤差, 信頼区間, 固定効果モデルのときのウェイトを計算.
df <-
df %>% mutate(
RR_hat = (XT*nC)/(XC*nT),
LRR_hat = log(RR_hat),
SE = sqrt((nT-XT)/(XT*nT) + (nC-XC)/(XC*nC)),
CI_L = LRR_hat - 1.96*SE,
CI_U = LRR_hat + 1.96*SE,
w = 1/(SE^2)
)
df %>% select(-nT, -nC)## XT XC RR_hat LRR_hat SE CI_L CI_U w
## 1 11 11 1.0000000 0.00000000 0.3393398 -0.6651061 0.66510605 8.684211
## 2 8 10 0.8000000 -0.22314355 0.3979112 -1.0030495 0.55676243 6.315789
## 3 6 10 0.6000000 -0.51082562 0.4472136 -1.3873643 0.36571302 5.000000
## 4 11 12 0.9166667 -0.08701138 0.3279874 -0.7298667 0.55584400 9.295775
## 5 4 12 0.3333333 -1.09861229 0.5163978 -2.1107519 -0.08647264 3.750000
## 6 6 9 0.6666667 -0.40546511 0.4594683 -1.3060230 0.49509274 4.736842
## 7 3 11 0.2727273 -1.29928298 0.5979764 -2.4713167 -0.12724927 2.796610
## 8 7 8 0.8750000 -0.13353139 0.4485426 -1.0126749 0.74561213 4.970414
## 9 4 8 0.5000000 -0.69314718 0.5552777 -1.7814915 0.39519713 3.243243
## 10 8 8 1.0000000 0.00000000 0.4281744 -0.8392219 0.83922186 5.454545統合リスク比(固定効果モデルのあてはめ)と信頼区間を計算(pはpoolの意)
RRp <- exp(sum(df$w * df$LRR_hat)/sum(df$w))
RRp## [1] 0.7261641RRp_CI_L <- exp(log(RRp) - 1.96*sqrt(1/sum(df$w)))
RRp_CI_L## [1] 0.5564975RRp_CI_U <- exp(log(RRp) + 1.96*sqrt(1/sum(df$w)))
RRp_CI_U## [1] 0.947559異質性について調べるため, 均質性の検定を行う. 統計量Q1を計算.
Q1 <- sum(df$w * (df$LRR_hat - log(RRp))^2)
Q1## [1] 7.808073Q1は自由度が研究数-1のカイ二乗分布に従う.
1-pchisq(Q1, df=Nstudy-1)## [1] 0.5535978リスク比について有意性の検定を行う.
Q2 <- (log(RRp)^2*sum(df$w))
Q2## [1] 5.5542171-pchisq(Q2, df=1)## [1] 0.01843621いちおう、このケースで, 変量効果モデルのあてはめをDerSimonian-Lairdの方法で行い, 変量効果の分散の推定値を確認する. 0になっているはず.
tau2_hat <- max(0,(Q1-(Nstudy-1))/(sum(df$w)-(sum(df$w^2)/sum(df$w))))
tau2_hat## [1] 02*2分割表データを変量効果モデルで生成(リスク比の推定)
次に,真のリスク比の従う正規分布を設定し, 変量効果モデルを真のモデルとしてデータを生成して, 固定効果モデルでパラメータの推定をし, 均質性の検定を行った後, 変量効果モデルをあてはめる.
真のリスク比の分布 N(0.8, 0.3^2)で, 各研究の真のリスク比はその分布に従って生成されているとする.
ES_RR_mu <- 0.8
ES_RR_tau <- 0.3
ES_RR_tau2 <- ES_RR_tau^2
Nstudy <- 10
nT <- nC <- 30
ES_RR <- rnorm(Nstudy, mean = ES_RR_mu, sd = ES_RR_tau)
pC <- 0.3
pT <- pC * ES_RR各研究について, イベント数データを生成する.
# Treatment群のイベント数
XT <- rbinom(Nstudy, nT, prob = pT)
# Comparisonのイベント数
XC <- rbinom(Nstudy, nC, prob = pC)
df <-
data.frame(
nT = nT,
nC = nC,
XT = XT,
XC = XC
)
df## nT nC XT XC
## 1 30 30 3 11
## 2 30 30 7 10
## 3 30 30 6 8
## 4 30 30 7 11
## 5 30 30 9 7
## 6 30 30 3 7
## 7 30 30 6 8
## 8 30 30 12 8
## 9 30 30 3 15
## 10 30 30 5 8リスク比の推定値, 対数リスク比推定値, 標準誤差, 信頼区間, 固定効果モデルのときのウェイトを計算.
df <-
df %>% mutate(
RR_hat = (XT*nC)/(XC*nT),
LRR_hat = log(RR_hat),
SE = sqrt((nT-XT)/(XT*nT) + (nC-XC)/(XC*nC)),
CI_L = LRR_hat - 1.96*SE,
CI_U = LRR_hat + 1.96*SE,
w = 1/(SE^2)
)
df %>% select(-nT, -nC)## XT XC RR_hat LRR_hat SE CI_L CI_U w
## 1 3 11 0.2727273 -1.2992830 0.5979764 -2.4713167 -0.1272493 2.796610
## 2 7 10 0.7000000 -0.3566749 0.4197505 -1.1793859 0.4660360 5.675676
## 3 6 8 0.7500000 -0.2876821 0.4743416 -1.2173917 0.6420276 4.444444
## 4 7 11 0.6363636 -0.4519851 0.4087781 -1.2531903 0.3492200 5.984456
## 5 9 7 1.2857143 0.2513144 0.4327835 -0.5969413 1.0995702 5.338983
## 6 3 7 0.4285714 -0.8472979 0.6399405 -2.1015812 0.4069855 2.441860
## 7 6 8 0.7500000 -0.2876821 0.4743416 -1.2173917 0.6420276 4.444444
## 8 12 8 1.5000000 0.4054651 0.3763863 -0.3322521 1.1431823 7.058824
## 9 3 15 0.2000000 -1.6094379 0.5773503 -2.7410444 -0.4778314 3.000000
## 10 5 8 0.6250000 -0.4700036 0.5082650 -1.4662031 0.5261958 3.870968統合リスク比(固定効果モデルのあてはめ)と信頼区間を計算(pはpoolの意)
RRp <- exp(sum(df$w * df$LRR_hat)/sum(df$w))
RRp## [1] 0.7099767RRp_CI_L <- exp(log(RRp) - 1.96*sqrt(1/sum(df$w)))
RRp_CI_L## [1] 0.5301898RRp_CI_U <- exp(log(RRp) + 1.96*sqrt(1/sum(df$w)))
RRp_CI_U## [1] 0.9507291変量効果の分散(研究間のバラツキ, ES_RR_tau2)を推定するため, 統計量Q1を計算.
Q1 <- sum(df$w * (df$LRR_hat - log(RRp))^2)
Q1## [1] 13.99194検定の結果も見ておく.
1-pchisq(Q1, df=Nstudy-1) ## [1] 0.1226122変量効果の分散を推定.
tau2_hat <- max(0,(Q1-(Nstudy-1))/(sum(df$w)-(sum(df$w^2)/sum(df$w))))
tau2_hat## [1] 0.1245095変量効果の分散の真値は0.09であり, 近いところを推定しているだろうか.
異質性の尺度I2
I2 <- (Q1 - (Nstudy - 1))/Q1*100
I2## [1] 35.67726DerSimonian-Lairdの方法で使うウェイトを計算
df <-
df %>% mutate(
w_RE = 1/(SE^2 + tau2_hat)
)統合リスク比(変量効果モデルのあてはめ)
RRp <- exp(sum(df$w_RE * df$LRR_hat)/sum(df$w_RE))
RRp## [1] 0.6749969RRp_CI_L <- exp(log(RRp) - 1.96*sqrt(1/sum(df$w_RE)))
RRp_CI_L## [1] 0.4666199RRp_CI_U <- exp(log(RRp) + 1.96*sqrt(1/sum(df$w_RE)))
RRp_CI_U## [1] 0.9764281有意性の検定
Q2 <- (log(RRp)^2*sum(df$w_RE))
Q2## [1] 4.3540611-pchisq(Q2, df=1)## [1] 0.03692081まとめ
人工データを作って実験してみた. 今回の例は一応, 上手くいっているが, 乱数シードを変えると, 変量効果モデルの分散が0と推定されてしまったり, 異質性がないと判断されてしまうことが多々あった.
変量効果と他の要素のバランスが重要なのかもしれない. 分散要素の比重だとは思うのだけれど, どのうようになるとどうのような結果が起きるかがまだよく分かっていない.
モデルを理解したと思ってコードを書いてはいるのだが, 正しく理解できているかまだ自信がない(メタアナリシスの数理モデルが明確に記述されている本があまり見当たらないため). 間違いに気付いたら書き換える予定だ.
