ベイジアンネットワークメタアナリシスのためのRパッケージgemtcで取り扱うデータ形式

ネットワークメタアナリシスでは、3つ以上の介入あるいは治療の効果を複数の研究からまとめて比較し、直接比較だけでなく間接比較のデータも含めて、各ペアを比較する効果指標の統合値と信頼区間(あるいはは確信区間)および効果の大きさの順位を知ることができます。ペア比較メタアナリシスは二つの介入のうち、どちらを選択すべきかという問題に解答を与えてくれますが、ネットワークメタアナリシスは三つ以上の介入のうちどれを選択すべきか、言い換えると最善の介入はどれかという問題に解答を与えてくれる可能性があります。→参考情報

ネットワークメタアナリシス用のソフトウェアとして、頻度論派メタアナリシスのためにはnetmetaというRのパッケージがありますが、gemtcはベイジアンネットワークメタアナリシスのためのRのパッケージで、rjagsというRのパッケージを介してJAGS (Just Another Gibbs Sampler)でMCMC(Malkov Chain Monte Carlo)シミュレーションを実行させます。JAGSとrjagsはあらかじめインストールしておきます。Rで使用する際には、library(rjags);library(gemtc)を最初に実行する必要があるのは、他のRのパッケージと同じです。

MCMCでは事後分布から、必要に応じて、数万から100万程度のサンプリング値を得るので、実行には時間がかかります。事前分布にはほぼ平坦な分布を用いますが、既知の情報に基づいて設定することもできます。

gemtcでは連続変数アウトカムの場合、平均値差Mean Difference (MD)を効果指標として取り扱えますが、最近、標準化平均値差Standardized Mean Difference (SMD)も取り扱えるようになりました。各研究からアームごとの研究ID、治療、平均値、標準偏差、サンプルサイズのデータを抽出します。列名はstudy, treatment, mean, std.dev, sampleSizeとします。同じ研究内では最初の行が参照アームになります。

このようなデータ形式をArm-level dataと呼んでいます。また、あらかじめ計算した平均値差とその標準誤差、あるいは、標準化平均値差とその標準誤差からもネットワークメタアナリシスを実行することもできます。この場合でも、netmetaなどのデータ形式と異なり、1行アームの形式になります。列名はstudy, treatment, diff, std.errとなります。このようなデータ形式をRelative effect dataと呼んでいます。

1つの行に2つのアームを比較した相対的効果指標とその標準誤差が含まれているのですが、参照アームの行も必要で、参照アームのdiffの欄はNAとします。2アームだけしかない場合は、std.errの欄もNAで問題ないのですが、3アーム以上ある場合は、参照アームのstd.errには参照アームの標準誤差の値が必要になります。平均値差であれば、対照アームの標準偏差をサンプルサイズの平方根で割り算した値になりますし、標準化平均値差の場合はサンプルサイズの平方根の逆数になります。参照アームの標準誤差の値がNAだとモデル作成時にエラーになります。

通常は、元データとして、各アームの平均値と標準偏差を用意したほうがいいでしょう。一部の研究で例えば、平均値差と標準誤差、あるいは標準化平均値差と標準誤差のデータしか得られないような場合、他の研究について、これらの相対効果指標と標準誤差を計算して解析するような場合は、Relative effect dataを使うことになるでしょう。

gemtcでは最初に、mtc.network()関数でネットワークオブジェクトを作成しますが、Arm-level dataの場合は、netw=mtc.network(data.ab=data)と記述し、Relative effect dataの場合は、netw=mtc.network(data.re=data)と記述します。

mtc.model()関数は変数netwに対して演算処理を行いますが、likelihoodとlinkの引数に効果指標のタイプに応じた値を設定する必要があります。元データが平均値と標準偏差の場合は、likelihood=”normal”, link=”identity”とすると、統合値は平均値差になります。liklihood=”normal”, link=”smd”とすると統合値は標準化平均値差になります。標準化平均値差はHedge’s gが用いられています。Cohen’s dに対し、サンプルサイズが小さい場合のずれを調整する項が追加されています。

元データが、Relative effect dataで、平均値差と標準誤差の場合、likelihood=”normal”, link=”identity”とすると、統合値は平均値差になります。

元データが、Relative effect dataで、標準化平均値差と標準誤差の場合、likelihood=”normal”, link=”identity”とすると、統合値は標準化平均値差になります。ただし、gemtcのforest()関数で、フォレストプロットを作成すると、ラベルがMean Differenceとなりますので、あとで画像として開いて書き換えるか、forestplotパッケージを使うなどして、別の方法でフォレストプロットを作成する必要があります。

例えば、netw = mtc.network(data.ab = data)で作成したネットワークオブジェクトの変数名がnetwだとすると、次のステップではmtc.model()関数で、netwに対してモデル作成を行います。

model <- mtc.model(netw, likelihood = lh, link = lk, linearModel=rndfix, n.chain=chain )

変数lh, lkには上記の文字列を格納します。ランダム効果モデルの場合は、rndfix = “random”とします。chainはMCMCを繰り返す回数です、例えば、chain=4とします。

その後、例えば、以下の様にMCMCを実行させます。

burnin = 10000
iteration = 110000
thin = 1
res <- mtc.run(model, n.adapt = burnin, n.iter = iteration, thin = thin)

以上を解説したPDFのスライドファイルをこちらで閲覧できます(4.3MB)。

gemtcはもちろん、リスク比、リスク差、オッズ比、ハザード比も取り扱えます。元データのフォーマットと、likelihood、linkの設定が分からないと、使えないので、今回特に、連続変数アウトカムの場合について、解説しました。実際のRのスクリプトやMCMCの結果の利用法についてはいつか解説したいと思います。

メタアナリシスのためのウェブツール:データ処理をJavaScriptで元データはExcelで

JavaScriptには自然対数関数log()、eの指数関数exp()、平方根sqrt()などの数学の関数を含むMathライブラリがあり、P値を求めたりZ値やt値を求める関数があるjStatという統計解析に必要な関数のライブラリがあります。これらの関数を用いることで、さまざまな統計解析ができます。

また、HTMLのcanvasオブジェクトを用いて、グラフの作成ができます。Forest plot、Funnel plotの描画に使えます。

従って、メタアナリシスをJavaScriptだけで実行することができます。

Excelでメタアナリシスに必要なデータを用意し、ファイルをドラグアンドドロップして、メタアナリシスを実行するウェブツールを作成しました。Forest plot、Prediction Intervalを含むForest plot、Funnel plotを作成します。

Excelシートのデータ形式は、Mindsシートとlabelのカラムを含む2種類に対応します。それらの例を含むExcelファイルはこちらから右クリックしてダウンロードできます(ポップアップメニューから名前を付けてリンクを保存を選択します)。

メタアナリシスは分散逆数法でランダム効果モデルを用い、研究間の分散の計算はRestricted Maximum Likelihood (REML)法を用いています。

Excel sheets to Meta-analysis with Prediction Interval: Web tool by JavaScript
URL: https://stat.zanet.biz/sr/drop_ex_js.htm

図1.メタアナリシスのためのウェブツール:Excel sheets to Meta-analysis with Prediction Interval: Web tool by JavaScript.

Excelファイルを右のDrag and Drop your Excel.xlsx file hereのパネルにドラグアンドドロップして、その上のドロップダウンメニューのシート名の一覧から対象のデータを含むシートを選択し、Do Meta-analysisボタンをクリックします。

図2.対象データを含むシートを選択し、Do Meta-analysisボタンをクリックする。
図3.Prediction Intervalを含むForest plotと含まないForest plotおよびFunnel plotが作成されます。

Do Meta-analysisボタンをクリックして結果が出力された直後の時点では、各研究の効果指標の値と95%信頼区間、統合値の値と95%信頼区間および95%予測区間Prediction Intervalの値がクリップボードに格納されているので、Excelなどに貼り付けて利用できます。

3つの画像のパネルは重なって表示されますが、移動可能です。各画像のパネルを右クリックしてコピーしたり、ファイルとして保存したりできます。また、数値データは上のテキストエリアに出力されるので、Copy to ClipboardボタンをクリックしてExcelなどに貼り付けます。

ドロップダウンメニューから別のシートを選択して、Do Meta-analysisボタンをクリックするとその結果が出力されデータが置き換えられます。

Prediction Intervalの計算はt分布に基づいていますので、tau2が0でも正規分布に基づいて計算している信頼区間と異なる値になります。結果はmetaforあるいはmetaと、同じモデルを用いた場合、同じであることを確認しています。取り扱える効果指標はRR, OR, RD, HR, MD, SMDです。

JavaScriptのコードはdrop_ex_jsとして、GitHubに公開しました。
https://github.com/tmorizane/drop_ex_js.git

著者が作成したのは、
drop_ex_js.htm
drop_ex_js.js
meta-an2-js.js
meta-an2-pi.js
です。
以下のJavaScriptのライブラリを用いています。jStat https://github.com/jstat/jstat.git
Math.js https://mathjs.org/
xlsx.full.min.js https://www.cdnpkg.com/xlsx/file/xlsx.full.min.js/?id=78603
jquery.min.js, jquery-ui.min.js https://jquery.com/download/

Rのmeta, metaforパッケージでメタアナリシスを実行する:Prediction Intervalも出力

Rのsource()関数を使うと、ウェブサーバーに置いたRスクリプトのファイルを直接読み込んでローカルで実行させることができます。例えば、以下のスクリプトは1行で書かれていますが、これをRのエディターに書き込んで実行すると、クリップボードのデータを変数exdatに格納し、source()関数に書かれているURLからmeta_lo_meta.Rというファイルに書かれているスクリプトをダウンロードしてRで実行します。このファイルには、変数exdatに対する処理が記述されており、metaパッケージを利用して、メタアナリシスを実行し、Forest plot、Funnel plotを作成するRスクリプトが記述されていますが、それらをRで実行します。

分散逆数法によるランダム効果モデルのメタアナリシスで、研究間の分散の計算にはRestricted Maximum Likelihood(REML)法を用います。

exdat = read.table("clipboard", sep = "\t", header = TRUE);source("https://stat.zanet.biz/useRs/scripts/meta_lo_meta.R");

メタアナリシスに必要なデータをExcelで用意し、Excelシートからデータを読み込んで、metaあるいはmetaforで処理できる形式にreformatし、クリップボードにコピーする操作をJavaScriptで実行することができます。なお、Macの場合は、exdat = read.table(“clipborad”, sep=”\t”, header = TRUE)の部分を以下の様に書き換えてください。つまり、”clipboard”をpipe(“pbpaste”)に書き換えてください。

exdat=read.delim(pipe("pbpaste"),sep="\t",header=TRUE);source("https://stat.zanet.biz/useRs/scripts/meta_lo_meta.R");

Excelファイルをドラグアンドドロップすることで、各シートのデータを読み込み、ドロップダウンメニューのシート名の一覧からシートを選択し、Reformat Dataでmetaあるいはmetaforで処理できる形式にreformatするウェブページに、上記のRスクリプトを選択するパネルを設定したウェブページを作成しました。Reformat Data from Excel and Do Meta-analysis URL: https://stat.zanet.biz/sr/drop_ex_cb.htm

図1.Reformat Data from Excel and Do Meta-analysis in R.

メタアナリシスに必要なデータを用意するExcelシートにはMindsの評価シートをそのまま使うこともできますし、labelのカラムを含む一定の形式のシートを使うこともできます。それらの例を含むExcelファイルはこちらから右クリックしてダウンロードできます(ポップアップメニューから名前を付けてリンクを保存を選択します)。Excelで開いて、データ部分を書き換えて使用してください。効果指標のタイプはExcelシート中で指定します。RR, OR, RD, HR, MD, SMDを取り扱えます。

使い方ですが、最初にRを起動しておきます。そして上記のURLから、Reformat Data from Excel and Do Meta-analysisのページを開き、metaforを使うか、metaを使うかを選択します。

metaforを選択すると以下のスクリプトが上段のテキストエリアに表示されます。

exdat = read.table("clipboard", sep = "\t", header = TRUE);source("https://stat.zanet.biz/useRs/scripts/forest_metafor_bs2.R");source("https://stat.zanet.biz/useRs/scripts/add_pi_forest.R");

右上のCopy & Paste in R Editorのボタンをクリックし、Rに戻り、ファイルメニューから新しいスクリプトを選択して、開かれたRエディタのフィールドに貼り付けます。

図2.Rのエディタにスクリプトを貼り付け(Ctr+V)た状態。

ウェブページに戻り、Excelファイルを右下のDrag and Drop your Excel.xlsx file hereのパネルにドラグアンドドロップします。その上のドロップダウンメニューにシート名の一覧が見られるようになりますので、メタアナリシスを実行したいシート名を選択してください。

図3.Excelファイルをドラグアンドドロップした後、ドロップダウンメニューのシート名の一覧。

続いて、Reformat Dataボタンをクリックすると、下のテキストエリアにリフォーマットしたデータが表示され、同時にデータはクリップボードにコピーされます。

図4.Excelシートのデータのリフォーマット。

そのまま、Rに戻り、先ほどの1行のスクリプトにカーソルが置かれた状態で、実行ボタンをクリックします。もし、途中でコピー操作を別の対象で行ったような場合は、Copy to Clipboardボタンをクリックして同じ操作をしてください。

図5.Rでmetafor, forestplotパッケージでメタアナリシスを実行した結果。

Rでmetafor, forestplotパッケージでメタアナリシスを実行した結果が出力されます。予測区間Prediction Intervalを含まないForest plot、それを含むForest plot、およびFunnel plotが作成され、クリップボードには各研究の効果指標と95%信頼区間、95%予測区間の値が格納されています。そのまま、元のExcelシートにクリップボードのデータを貼り付けることができます。

各プロットはコピーしてExcel、PowerPoint, Wordなどに貼り付けて利用することができます。また、Rのファイルメニューからさまざまな形式の画像ファイルとして保存することもできます。

さて、Reformat Data from Excel and Do Meta-analysis のウェブページで、metaを選択した場合は、Layout:メニューからForest plotの形式を選択できます。

図6.metaを選択した場合は、Layout:メニューからForest plotの形式を選択できます。

例えば、RevMan5を選択した場合の例です:

図7.metaでRevMan5のLayoutの例。

予測区間を含むForest plotとそれを含まないForest plot、Funnel plotが作成され、クリップボードには各研究の効果指標と95%信頼区間、95%予測区間の値が格納されています。なお、metaでは研究数が10件以上ないと、Beggの検定、Eggerの検定は実行されません。

metaでLayoutをJAMAにした別のデータの例です:

図8.metaでJAMAのLayoutの例。

クリップボードを介して、スクリプトやデータを取り込むので、途中で別の対象に対してコピー操作をすると入れ替わることに注意してください。

metaでLayoutをRevMan5にする場合:

exdat = read.table("clipboard", sep = "\t", header = TRUE);source("https://stat.zanet.biz/useRs/scripts/meta_lo_revman.R");

metaでLayoutをJAMAにする場合:

exdat = read.table("clipboard", sep = "\t", header = TRUE);source("https://stat.zanet.biz/useRs/scripts/meta_lo_JAMA.R");

いずれのスクリプトも、Macの場合は、”clipboard”をpipe(“pbpaste”)に書き換えてください。

Reformat Data from Excel and Do Meta-analysis URL: https://stat.zanet.biz/sr/drop_ex_cb.htmを使う利点は、複数のシートを含むExcelファイルを1回読み込ませると、各シートのデータに対してメタアナリシスをセルの範囲を選択することなく、連続して、実行できる点です。益と害の複数のアウトカムに対するメタアナリシスを行う場合、それぞれのアウトカムに対してひとつのシートにデータが入力されていると思いますが、それらに対するメタアナリシスを続けて実行できます。

ランダム効果モデルによるメタアナリシスの予測区間Prediction Interval

メタアナリシスでは複数の研究の効果推定値の統合値とその95%信頼区間が計算されます。ランダム効果モデルでは研究間のばらつきの指標としてτ2値が計算され、予測区間Prediction Intervalの計算は統合値の分散にτ2(タウ二乗)を加算した値を元に計算されます。そのため予測区間は信頼区間よりも幅広くなり、τ2が大きいほど、またI2が大きいほどその差が大きくなります。ただし、予測区間に意味があるのはランダム効果モデルの場合です。

また信頼区間は正規分布を元に計算するのが一般的で、95 %信頼区間であればZ値1.96を統合値の標準誤差に掛け算して±して計算します。一方、予測区間の場合も、正規分布を用いることが可能ですが、正規分布ではなくt分布に基づいて計算することが一般的で、研究数-1の自由度のt分布におけるp = 0.975(累積確率が0.975)のt値を統合値の標準誤差に掛け算して±して計算します。Rのメタアナリシス用のパッケージであるmetaでのmetabin等の関数のデフォルトの設定は、t分布を用いる様になっています。

信頼区間に正規分布、予測区間にt分布を用いると、τ2値が0の場合でも、これら二つの区間は違う値になり、t分布の方が幅が広くなります。もし、両方とも正規分布あるいはt分布を用いれば、τ2値が0の場合は、同じ値になるはずです。また、τ2値の計算にDersimonian-Laird法を用いる場合と、REML(Restricted Maximum Likelihood)法を用いる場合でも異なる値になる可能性があります。(現在は、REML法の使用が推奨されています。)

Rのメタアナリシス用のパッケージであるmetaforでは予測区間の計算はrma関数の結果に対してpredict関数で処理することで行われますが、デフォルトでは正規分布に基づいて計算されるようになっています。metaforの場合は、rma関数で引数test=”t”に設定すると、信頼区間がt分布に基づいて計算され、その結果をpredict関数で処理するとt分布に基づく予測区間が得られます。信頼区間をデフォルトの設定で正規分布に基づいて計算した場合、予測区間をt分布に基づいて計算するには自分で計算する必要があります。

正規分布の場合の計算式は以下のとおりですが、t分布の場合は、z1-α/2をt1-α/2に置き換えて計算します。t分布の自由度は、研究数をkとした場合k – 1です。95%予測区間の場合は、α=0.05です。

metaforの場合、例えばrma関数で処理した結果を変数resに格納したとすると、t分布に基づく予測区間の計算は以下のスクリプトで計算できます。resから統合値b、その標準誤差se、τ2値 tau2、研究数kを得て、qt関数で得たt値 tvalを合わせて95%予測区間を計算します。この例は、リスク比の場合なので、b、se、τ2値は対数スケール上で計算されているので、変数convlogを1に設定し、上限値uppert、下限値lowertからExponentialを計算し、predupt、predlowtを計算しています。

#Prediction Interval based on t-distribution.
mu_hat = res$b
se_mu = res$se
tau2 = res$tau2
k = res$k
tval = qt(0.975, df = k - 1)
lowert = mu_hat - tval * sqrt(se_mu^2 + tau2)
uppert = mu_hat + tval * sqrt(se_mu^2 + tau2)
predlwt=lowert
predupt=uppert
if(convlog==1){
predlwt=exp(lowert)
predupt=exp(uppert)
}

すでに述べたように、τ2値が0の場合は、95%信頼区間も95%予測区間も同じ正規分布、あるいはt分布を用いて計算すれば同じ値になるはずです。しかし、例えばmetaの場合であれば、デフォルトでは信頼区間は正規分布で予測分布はt分布なので、同じ値にはなりません。

Cochrane Handbook for Systematic Reviews of Interventionsによるとメタアナリシスの結果を提示する際に、統合値の信頼区間(Confidence Interval)に加えて予測区間(Prediction Interval)を含めることは、必須とまでは明記されていませんが、特にランダム効果モデルを用いたメタアナリシスにおいては、それを含めることが強く推奨されており、適切でない解釈を避けるために重要であるとされています。ただし、ソースは、研究数が適切(例:5つ以上)で、ファンネルプロットに明確な非対称性がない場合に、予測区間の使用を推奨しています。そして、研究数が少ない場合、予測区間は見せかけ上、広く見えたり狭く見えたりする可能性があり、問題となることがあるとされています。(10.10.4.2 Interpretation of random-effects meta-analyses)

また、metaやmetaforでは、頻度論派Frequentistの方法でメタアナリシスを実行していますので、信頼区間は真の効果推定値の分布を表しているわけではありません。同じことを繰り返したら、そのたびに異なる信頼区間の値が得られますが、95%の場合は、真の値がその間に含まれているという意味です。例えば、リスク比が1.0以上になる確率を知りたい場合、頻度論派の信頼区間に基づいて推測するのではなく、ベイジアンメタアナリシスの確信区間に基づいて推測する必要があります。前者の場合、信頼区間と確信区間がほぼ同じということを前提にしていることになります。

文献:

Higgins JP, Thompson SG, Spiegelhalter DJ: A re-evaluation of random-effects meta-analysis. J R Stat Soc Ser A Stat Soc 2009;172:137-159. doi: 10.1111/j.1467-985X.2008.00552.x PMID: 19381330
URL: https://pubmed.ncbi.nlm.nih.gov/19381330/

Cochrane Handbook for Systematic Reviews of Interventions
URL: https://www.cochrane.org/authors/handbooks-and-manuals/handbook

metafor https://cran.r-project.org/web/packages/metafor/index.html
Viechtbauer W (2010). “Conducting meta-analyses in R with the metafor package.” Journal of Statistical Software, 36(3), 1–48. doi:10.18637/jss.v036.i03.

meta https://cran.r-project.org/web/packages/meta/index.html
Balduzzi S, Rücker G, Schwarzer G (2019). “How to perform a meta-analysis with R: a practical tutorial.” Evidence-Based Mental Health, 153–160.

ランダム効果モデルのメタアナリシスの結果を臨床に適用する場合、Prediction Intervalも含め結果の解釈には注意が必要です。YouTube チャンネル IZ statに「メタアナリシスの臨床応用:ランダム効果モデルの解釈-ダイアログ- 」の動画を公開しています。

IZ stat URL:https://www.youtube.com/channel/UCqzzJbfQwKDvValptCAM4gg