メタ分析については、R には多くのパッケージがあり、公式サイト中にも Michael Dewey によるメタ分析のパッケージの紹介ページがあります。
また、Minds 診療ガイドライン作成マニュアル2020 ver. 3.0 には、複数の事例とサンプルコードが解説されています。
R以外のソフトウェアとしては、メタ分析を推奨している Cochrane から RevMan をいうソフトウェアを無償で提供しています。
メタ分析は、参照する文献が用いている統計手法によって違いが出てきます。
実際の論文をもとに、比率 (proportion) のメタ分析を行ってみましょう。ここで、比率 (proportion) とは、ある母集団に対して特定の特徴を持った人の比率です。
なお、rate と proportion の用語は、ここでは深入りしません。
Abate SM, Ahmed Ali S, Mantfardo B, & Basu B (2020). Rate of Intensive Care Unit admission and outcomes among patients with coronavirus: A systematic review and Meta-analysis. PloS one, 15(7), e0235653. https://doi.org/10.1371/journal.pone.0235653
コロナウィルス感染患者の ICU admission のデータです。いわゆる新型コロナ (SARS-COV-2) だけではなく、SARS や MERS のコロナウィルスも含めて、ICU入室率を扱う論文を集めてメタ分析しています。
この論文では、メタ分析を行った際の表データがインターネット上で公開されています。インターネットから、データを直接読み込みましょう。
dfMeta <- readr::read_csv("https://ndownloader.figshare.com/files/23167073")Rows: 38 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): Author, country, corona, Study design
dbl (3): year, event, sample
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.列ごとにデータ型が示されていることが分かります。
ここで、Study design は、列名に空白があるので、バッククォート記号で囲みます。
元データの最後の行は空データで、おそらく元々は統合データを格納するために準備してあったと思われます。この行は不要なので、NA データのある行を削除しました。
R では、データの入っていないことを表す手段が複数あります。NULL や NA です。NA などがあるとエラーが出やすくなります。
dfMeta <- na.omit(dfMeta)
dfMeta$Study_design <- dfMeta$`Study design`メタ分析を行うのは、わずか一行です。なお、ここでは結果を直接表示せず、いったん mpResult に格納しています。
mpResult <- meta::metaprop(data = dfMeta,
event = event,
n = sample,
studlab = paste(Author, year))関数
変数と引数
結果を見てみます。
mpResultNumber of studies: k = 37
Number of observations: o = 32741
Number of events: e = 6220
proportion 95%-CI
Common effect model 0.1900 [0.1858; 0.1943]
Random effects model 0.2834 [0.2057; 0.3765]
Quantifying heterogeneity (with 95%-CIs):
tau^2 = 1.6165; tau = 1.2714; I^2 = 99.0% [98.8%; 99.1%]; H = 9.83 [9.26; 10.44]
Test of heterogeneity:
Q d.f. p-value
Wald 3481.43 36 0
LRT 4446.74 36 0
Details of meta-analysis methods:
- Random intercept logistic regression model
- Maximum-likelihood estimator for tau^2
- Calculation of I^2 based on Q
- Logit transformationこの上までが、R で出力される結果です。
結果においては、異質性 (heterogeneity) 検定の結果が示されています。 Q, \(τ^2\), \(I^2\) などが異質性の指標です。
\(I^2\) の解釈の仕方:
フォレストプロットを作図してみましょう。
meta::forest(mpResult)
あまり美しくないので、オプションを試してみます。
meta::forest(mpResult,
layout = "RevMan5",
comb.fixed = FALSE,
label.right = "Favours control",
col.label.right = "red",
label.left = "Favours experimental",
col.label.left = "green",
prediction = TRUE)
これまで、研究はポジティブな結果が出た時に発表される傾向がありました。これはメタ分析を行う際にはバイアスとなってしまいます。これを出版バイアスと言います。多くの論文では、ファンネルプロットを作成して、視覚的に出版バイアスがあるかどうかを確認するに止まります。
ファンネルプロットを描いてみます。
meta::funnel(mpResult)
メタ分析の論文では、出版バイアスがあるかどうかをファンネルプロットで視覚的に確認したというものが多く見られます。しかしながら、出版バイアスを検証する統計もあります。
metabias(mpResult)Linear regression test of funnel plot asymmetry
Test result: t = 1.40, df = 35, p-value = 0.1698
Bias estimate: 2.8392 (SE = 2.0256)
Details:
- multiplicative residual heterogeneity variance (tau^2 = 94.1829)
- predictor: standard error
- weight: inverse variance
- reference: Egger et al. (1997), BMJこのデータに対し、Egger 検定を行ったところ、
p = 0.1698となり、出版バイアスはあるとは言えません。
R におけるメタ分析は、分析（検定や作図など）がそれぞれ一行でできるほど簡単なものになっています。
なお作図に関しては、出版できるようにするためには大きさなど微調整が必要になります。