【metafor】メタ解析におけるフォレストプロット解析のやり方【メタ解析】

「メタ解析」と聞いてもどんなことをするかわからない人は多いと思います。この記事では、メタ解析のライブラリであるmetaforを用いてフォレストプロットを作成する方法を学びます。フォレストプロットができると、複数のデータセットを比較することが可能です。ぜひ見てみましょう。

Tested Environment

macOS Monterey (12.4), Quad-Core Intel Core i7, Memory 32GB metafor 4.2-0

公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中

¥3,000 → ¥1,500 今なら50%OFF!!

プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！

技術書ページへ

メタ解析のための統計知識

メタ解析におけるフォレストプロットを学ぶ前に、いくつかの基本的な統計知識を理解しておくことが大切です。これらの用語や概念は、初めて耳にすると少し難しく感じるかもしれません。「そんなものなのだな〜」という程度で収めておくで大丈夫です。

メタ解析のポイント

メタ分析を進めるには、全ての研究から「効果量」を明確に特定しなければなりません。効果量とは、研究の結果の大きさや重要性を示す数値のことを指します。時々、この数値は研究の論文から直接取り出すことができますが、大抵の場合、報告された他のデータからこの数値を計算する必要があります。

この効果量の計算は、メタ分析の最終結果とそのどれだけ解釈できるかに大きな影響を与える可能性があるため、非常に重要です。

比較できる (Comparable):
- 意味: あるテストの結果を他のテストの結果と比較することは、それらが同じ種類の情報を測定している場合にのみ意味があります。
- 例: テストAとテストBはどちらも数学のテストですが、一つは初級レベル、もう一つは上級レベルです。これらのテストの点数を直接比較するのは、それぞれのテストの難易度が異なるため、意味がありません。
計算できる (Computable):
- 意味: 研究の結果から効果量を正確に計算することができなければなりません。これは、メタ分析で研究の結果を合成する際に重要です。
- 例: ある研究では、特定のデータが不足していて効果量を計算することができない場合、その研究をメタ分析に含めることはできません。
信頼できる (Reliable):
- 意味: 計算された効果量は、統計的に信頼性があり、バイアスや誤差が最小限でなければなりません。
- 例: 効果量の計算に使用されるデータや方法に問題がある場合、その効果量は信頼できない可能性があります。例えば、小さなサンプルサイズで計算された効果量は、不正確である可能性が高いです。
解釈できる (Interpretable):
- 意味: 効果量は、研究の目的や質問に対して意味のある形で解釈されなければなりません。
- 例: ある研究では、二つの変数の関連性に関心があります。この場合、効果量として相関係数を使用するのが適切です。しかし、他の種類の効果量を使用すると、結果の解釈が難しくなる可能性があります。

簡単に言うと、研究結果を比較、計算、信頼、および解釈する能力は、メタ分析を行う上での重要な考慮事項です。それぞれの研究とその結果の特性をよく理解することで、正確で信頼性のあるメタ分析の結果を得ることができます。

効果量

効果量(エフェクトサイズ・効果サイズ)とは、実験や調査の結果がどれだけ実質的な意味を持つのかを測定するための統計的なツールの一つです。効果量は、実験や研究の結果がどれだけ大きい（または小さい）かを示す数値です。この値が大きければ大きいほど、研究の結果が実質的に重要であることを示しています。

例えば、二つのグループ間のテストスコアの平均の差（コーエンのd）、二つの変数の関連性（相関係数 r）、またはある疾患に対する二つの異なる治療法の効果の比較（リスク比 RR）などを効果量で表すことができます。

効果量の例

コーエンのd:
- 2つのグループ間の平均の差を、両グループの標準偏差で割ったものです。
- 例：治療前と治療後のテストスコアの差。
相関係数 (r):
- 2つの変数の関連性を示します。
- 例：身長と体重の間の関連性。
リスク比 (RR):
- ある疾患に対する2つの異なる治療法の効果を比較する時に使います。
- 例：ワクチンAとワクチンBの効果の比較。

効果量を知ることで、ある研究が単に統計的に有意であるだけでなく、実際の世界でどれだけ意味があるのかを理解することができます。これにより、研究の結果をより正確に解釈する手がかりが得られます。

ランダム効果モデル

ランダム効果モデルは、異なるグループ間の変動や、階層的、またはグループ化されたデータの変動を考慮するための統計的手法の一つです。例えば、複数の学校で学生のテストスコアを分析する際、各学校の特性（ランダム効果）を考慮する必要があります。

なぜフォレストプロットの際にランダム効果モデルをフィットさせる必要があるのかというと、それは研究間の異質性を考慮するためです。複数の研究やデータセットからのデータを組み合わせる際、それぞれのデータセットは異なる実験条件、患者の特性、測定方法などにより、本質的な異質性を持っています。この異質性は、分析の結果に大きな影響を与える可能性があります。

簡単な例を使って、ランダム効果モデルの動作を実際に見てみましょう。R言語を使用して、グループごとに異なるインターセプトを持つデータを生成し、ランダム効果モデルを使用してこれらのインターセプトを推定します。

この例では、真のインターセプト（x軸）と推定されたインターセプト（y軸）の一致度を視覚的に確認できます。赤い線は、両者が完全に一致する場合の線です。

forestplotの実装方法

まずは全体のコードを紹介します。こちらを実行してください。

# 必要なパッケージのインストール
install.packages("metafor")
library(metafor)

dat.bcg

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg,
              slab=paste(author, year, sep=", ")) # also add study labels
dat

# fit random-effects model
res <- rma(yi, vi, data=dat, test="knha")
res

forest(res, atransf=exp, at=log(c(.05, .25, 1, 4)), xlim=c(-16,6),
       ilab=cbind(tpos, tneg, cpos, cneg), ilab.xpos=c(-9.5,-8,-6,-4.5),
       header="Author(s) and Year", shade="zebra")
text(c(-9.5,-8,-6,-4.5), 15, c("TB+", "TB-", "TB+", "TB-"), font=2)
text(c(-8.75,-5.25),     16, c("Vaccinated", "Control"), font=2)

以下のような出力が出れば成功です！

コードの各部分の詳細について説明します。

install.packages("metafor")
library(metafor)

ここでは、metaforパッケージをインストールし、ライブラリをロードします。

dat.bcg

この行では、dat.bcgデータセットを表示しています。このデータセットは、metaforパッケージに含まれています。

> dat.bcg
   trial               author year tpos  tneg cpos  cneg ablat      alloc
1      1              Aronson 1948    4   119   11   128    44     random
2      2     Ferguson & Simes 1949    6   300   29   274    55     random
3      3      Rosenthal et al 1960    3   228   11   209    42     random
4      4    Hart & Sutherland 1977   62 13536  248 12619    52     random
5      5 Frimodt-Moller et al 1973   33  5036   47  5761    13  alternate
6      6      Stein & Aronson 1953  180  1361  372  1079    44  alternate
7      7     Vandiviere et al 1973    8  2537   10   619    19     random
8      8           TPT Madras 1980  505 87886  499 87892    13     random
9      9     Coetzee & Berjak 1968   29  7470   45  7232    27     random
10    10      Rosenthal et al 1961   17  1699   65  1600    42 systematic
11    11       Comstock et al 1974  186 50448  141 27197    18 systematic
12    12   Comstock & Webster 1969    5  2493    3  2338    33 systematic
13    13       Comstock et al 1976   27 16886   29 17825    33 systematic

dat.bcg　の中身はこの様になっています。

tpos, tneg, cpos, cnegはそれぞれ、治療群と対照群の陽性および陰性の症例数を指します。

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg,
              slab=paste(author, year, sep=", "))

この行では、escalc関数を使用して、各研究の対数リスク比（RR）とその分散を計算します。

ここでdatの中身を見てみると以下のようになっています。yiとviが追加されています。

> dat

   trial               author year tpos  tneg cpos  cneg ablat      alloc      yi     vi 
1      1              Aronson 1948    4   119   11   128    44     random -0.8893 0.3256 
2      2     Ferguson & Simes 1949    6   300   29   274    55     random -1.5854 0.1946 
3      3      Rosenthal et al 1960    3   228   11   209    42     random -1.3481 0.4154 
4      4    Hart & Sutherland 1977   62 13536  248 12619    52     random -1.4416 0.0200 
5      5 Frimodt-Moller et al 1973   33  5036   47  5761    13  alternate -0.2175 0.0512 
6      6      Stein & Aronson 1953  180  1361  372  1079    44  alternate -0.7861 0.0069 
7      7     Vandiviere et al 1973    8  2537   10   619    19     random -1.6209 0.2230 
8      8           TPT Madras 1980  505 87886  499 87892    13     random  0.0120 0.0040 
9      9     Coetzee & Berjak 1968   29  7470   45  7232    27     random -0.4694 0.0564 
10    10      Rosenthal et al 1961   17  1699   65  1600    42 systematic -1.3713 0.0730 
11    11       Comstock et al 1974  186 50448  141 27197    18 systematic -0.3394 0.0124 
12    12   Comstock & Webster 1969    5  2493    3  2338    33 systematic  0.4459 0.5325 
13    13       Comstock et al 1976   27 16886   29 17825    33 systematic -0.0173 0.0714

yi: 効果サイズの推定値です。上のデータの例では、yiは 0.8893 および 1.5854 です。これは対数変換されたリスク比またはオッズ比を表しています。
vi: yiの分散です。上のデータの例では、viは 0.3256 および 0.1946 です。分散の値は、yiの不確実性を表します。viが小さいほど、対応するyiの推定値はより確実です。

res <- rma(yi, vi, data=dat, test="knha")

この行では、rma関数を使用して、ランダムエフェクトモデルを適用します。

test="knha": これは、不均一性の検定方法を指定します。"knha"は、Knapp and Hartungの方法を使用して不均一性をテストすることを意味します。

> res

Random-Effects Model (k = 13; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
tau (square root of estimated tau^2 value):      0.5597
I^2 (total heterogeneity / total variability):   92.22%
H^2 (total variability / sampling variability):  12.86

Test for Heterogeneity:
Q(df = 12) = 152.2330, p-val < .0001

Model Results:

estimate      se     tval  df    pval    ci.lb    ci.ub     
 -0.7145  0.1808  -3.9522  12  0.0019  -1.1084  -0.3206  ** 

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

モデルの情報

Random-Effects Model (k = 13; tau^2 estimator: REML): ランダムエフェクトモデルが適用され、13の研究 (k=13) が含まれています。tau^2 の推定にはREML（Restricted Maximum Likelihood）が使用されています。

統計の情報

tau^2: 総不均一性の推定量です。値 (0.3132) が大きいほど、研究間の不均一性が大きいことを意味します。
tau: tau^2の平方根の値です (0.5597)。
I^2: 総不均一性を総変動で割った値で、不均一性の割合を示します。このモデルでは、92.22% です。
H^2: 総変動をサンプリング変動で割った値です (12.86)。

不均一性のテスト

Q(df = 12) = 152.2330, p-val < .0001: Q統計量とその自由度 (df)、p値を表示しています。このp値は非常に低い (< .0001) ため、研究間には有意な不均一性が存在することを示しています。

モデルの結果

estimate: 効果量の推定値です (0.7145)。
se: 推定値の標準誤差です (0.1808)。
tval: t値です (3.9522)。
df: 自由度です (12)。
pval: t検定のp値です (0.0019)。この値が低いため、推定された効果量は統計的に有意です。
ci.lb および ci.ub: 推定値の信頼区間の下限 (1.1084) および上限 (0.3206) です。

最後の行（Signif. codes）は、p値の有意性コードを示しています。このモデルでは、** は p 値が 0.001 以下であることを示しています。

forest(res, atransf=exp, at=log(c(.05, .25, 1, 4)), xlim=c(-16,6),
       ilab=cbind(tpos, tneg, cpos, cneg), ilab.xpos=c(-9.5,-8,-6,-4.5),
       header="Author(s) and Year", shade="zebra")

この行では、forest関数を使用して、フォレストプロットを作成します。

text(c(-9.5,-8,-6,-4.5), 15, c("TB+", "TB-", "TB+", "TB-"), font=2)
text(c(-8.75,-5.25),     16, c("Vaccinated", "Control"), font=2)

最後の2行では、フォレストプロットに追加のテキストを追加しています。

メタ分析においては、個々の研究から「効果サイズ」を計算します。効果サイズは、研究の結果の大きさを表す統計量です。例えば、二つのグループの平均の差や、相関係数、リスク比（RR）などが効果サイズとして使われます。

リスク比の解釈は以下のようになります。

RR = 1: 両方のグループでイベントのリスクは同じです。
RR > 1: 治療群でイベントのリスクが対照群に比べて高いことを示します（治療のネガティブな効果を示すかもしれません）。
RR < 1: 治療群でイベントのリスクが対照群に比べて低いことを示します（治療のポジティブな効果を示すかもしれません）。

フォレストプロットの見方

フォレストプロットを見てみると各プロットの大きさやエラーバーが矢印になっていたりすることに気づくかと思います。それぞれの意味を見ていきましょう。

フォレストプロットの大きさ

各研究のプロットの大きさ（通常は正方形や円で表示される）はその研究の重みを表しています。重みは、メタアナリシスでの各研究の相対的な寄与度を示し、一般的にはその研究のサンプルサイズや分散（または標準誤差）に依存します。

エラーバーが矢印で終わる理由

フォレストプロットにおいて、エラーバーが矢印で終わることがあります。これは、信頼区間がプロットの表示範囲を超えていることを示しています。具体的には、矢印は信頼区間の上限または下限がプロットに表示される領域外にあることを示します。

RE Model（Random Effects Model）のダイヤモンド

ランダム効果モデルに基づいて計算された全体の効果量（または総効果）とその信頼区間を表します。

ダイヤモンドの中央は、全体の効果量の推定値を表します。
ダイヤモンドの幅は、その推定値の信頼区間を表します。

このダイヤモンドは、個々の研究の効果量と比較すると、全体の効果量の視覚的な表現を提供します。フォレストプロットでは、個々の研究の効果量は通常、小さな四角形や丸で表され、それぞれの信頼区間は、それらの横に描かれた水平なエラーバーによって示されます。ダイヤモンドは、これらすべての個別の効果量を統合したものを示すため、通常、プロットの下部に描かれます。