Chapter 1 R でネットワークメタ
1.2 インストール
頻度論分析には netmeta パッケージを用いる。
ベイジアン分析には、まず JAGS をインストールする。最新版は 4.3.1。これは R とは関係ない独立したソフトウェアである。
https://mcmc-jags.sourceforge.io/
次に、R から JAGS を使ってベイジアンネットメタを行うためのパッケージとして、gemtc と BUGSnet をインストールする。
install.packages("gemtc")
install.packages("remotes")
library(remotes)
remotes::install_github("audrey-b/BUGSnet@v1.1.0", upgrade = TRUE, dependencies = TRUE)ここで用いていないパッケージも簡単に紹介する。
ベイズ方式として、STAN を用いる場合、Windows であれば RTools、macOS であれば XCode などのプログラミング開発環境が必要になる。
1.3 パッケージ
頻度論分析は netmeta を用いる。
ベイズ分析には、gemtc と BUGSnet を用いる。Hu et al. (2020) は、両者を比較して、「The output given by gemtc is limited. For example, gemtc does not have the option to report the summary effect as either the relative risk or absolute risk. Further, the output is not available in a table format. While BUGSnet is limited to analyzing arm-level data which could be a limitation for veterinary data which is often reported at the contrast level.」と述べている。
本報告書では扱わないパッケージについても簡単に紹介する。頻度論的 NMA のパッケージとしては、pcnetmeta がある。ベイズ的 NMA として、STAN を用いるものに MetaStan と multinma がある。MetaStan は NMA というよりもメタ分析パッケージではあるが、簡易な NMA 機能がある。一方、 multinma は、より高機能であり、CRAN サイトにドキュメントも揃っている (Phillippo et al. 2020)。
https://cran.r-project.org/web/packages/multinma/index.html
ベイズ方式で、BUGS を用いるものに R2OpenBUGS がある。これは BUGS 言語という R とは独立した言語でモデルを記述するものである。LeVasseur et al. (2021) は、ネットワークめた分析にこのパッケージを用い、BUGS 言語を論文の補足資料として提供している。
1.4 アウトカム
ネットワークメタ分析は、複数の治療群を、直接または間接的に比較するメタ分析手法である。ネットワークメタ分析は、大きく分けて頻度論分析とベイズ分析の２手法があり、それぞれ R パッケージが異なる。また、メタ分析する際には効果量として以下のようなものが用いられる。
- 平均差 (mean difference): 2群の差の平均値と標準誤差。体重や HbA1c など、値そのものを共通して使うときに用いる。
- 標準平均差 (standardized mean difference):
- OR/log OR: 頻度論分析の netmeta では、OR/logOR を直接受け付ける。このため、ロジスティック回帰分析で計算された OR/logOR も分析が可能。ベイズ分析パッケージは対応していない。
- 2x2表データ; ２群の発症/非発症の人数データがある場合には、OR または RR が計算できる。OR は、ケース対照や傾向スコアマッチングなど、対照群が恣意的に抽出された際に用いる。そうでなければ RR の方が好ましい。
組み入れている研究が RCT など前後比較の場合は注意が必要。介入前後の平均値では分析できない。(post - pre) の平均値と標準誤差が必要。
1.5 データ
モデルとしては、arm-based model と contrast-based model がある。White IR (2019) によると、両者に優劣はなく、「Models with fixed study intercepts can be recommended and may be implemented with either a CB or an AB model.」としている。
しかし、データの準備方法が異なってくる。また、パッケージによって使用する model が異なる。このため、両者の違いを押さえておきたい。例として、ダイエットのRCT研究で有酸素運動群とレジスタンス運動群があるとしよう。
arm-based model: 平均値と標準偏差（または標準誤差とサンプルサイズ）を群ごとに要約する。ベイス分析の
gemtcとBUGSnet(Béliveau 2019) と頻度論分析のnetmetaはこちらに対応している。contrast-based model: 平均値の差の平均値と標準誤差を要約する。ベイズ分析の
gemtcと頻度論分析のnetmetaはこちらに対応している。必要なのは、２群間の logOR/MD/logHR のいずれか。
contrast-based model において、もし３群（以上）の研究があった場合は、以下のように扱う。３群をそれぞれA群、B群、C群とする。このうち一つの群（C群とする）を参照群とし、A-C と B-C という二つの contrast を作成する。
なお、netmeta の pairwise() で arm-based から contrast-based に変換すすれば contrast-based データとして、BUGSnet などに渡すこともできる。
OR/RR で、元データが「サンプルサイズとイベント数」の場合は、arm-based 形式が良い。まず、netmeta は contrast-based も対応しているものの、イベント数に0があると contrast-based では除かれる。しかし、arm-based で、pairwise() 関数を使って contrast-based に変換すると、+0.5 をして処理するためにペアが除かれることはない。gemtc は contrast-based は全く対応していない。BUGSnet も contrast-based は全く対応していない。
1.5.2 gemtc
gemtc は、arm-based と contrast-based の両方に対応している点ではデータの準備が容易である。
一方で、gemtc はデータフレームの列名を指定する必要がある。
- 平均差の場合: データフレームの列名を [mean, std.err] または [mean, std.dev, sampleSize] とする。
mtc.model()の引数に、likelihood = "normal"を追加する。 - logOR の場合: データフレームの列名を [responders, sampleSize] とする。
mtc.model()の引数に、likelihood = "binom"とlinks = "logit"を追加する。 - logRR の場合: データフレームの列名を [responders, sampleSize] とする。
mtc.model()の引数に、likelihood = "binom"とlinks = "log"を追加する。 - logHR の場合: データフレームの列名を [responders, sampleSize] として、
mtc.model()の引数に、likelihood = "binom"とlinks = "cloglog"を追加する。あるいは、データフレームの列名を [responders, exposure] として、mtc.model()の引数に、likelihood = "binom"とlinks = "log"を追加する。
1.6 メタ分析の流れ
| netmeta | gemtc | BUGSnet | |
|---|---|---|---|
| データ | netmeta() |
なし。データフレームをそのまま使う | data.prep() |
| ネットワークプロット | netgraph() |
plot() |
summary() |
| 要約 | summary() |
summary() |
nettab() |
| ベイズモデル作成 | なし | mtc.model() |
nma.model() |
| モデル実行 | なし | mtc.run() |
nma.run() |
| 収束 | なし | gelman.diag()$mpsrf |
nma.diag() |
| 適合度 | なし | mtc.levplot(mtc.deviance()) |
nma.fit() |
| 一貫性 | decomp.design() |
nma.league()$heatplot |
|
plot(netsplit(), digits=3) |
mtc.nodesplit() |
||
| ランク | netrank() |
rank.probability() |
nma.rank() |
dmetar::sucra(netrank()) |
nma.rank()$sucraplot |
||
| フォレスト | metafor::forest() |
gemtc::forest(relative.effect()) |
nma.forest() |
1.6.1 データ準備
上で詳しく述べたが、データには大きく分けて以下を考慮する必要がある。
- arm-based か contrast-based か
- 入力: 平均値と標準偏差、平均値と標準誤差、観測人数
- 出力: 平均差、標準化平均差、オッズ比、リスク比
gemtc は、データフレームをそのまま分析に用いるが、netmeta と BUGSnet はそれぞれ独自のクラスのオブジェクトに変換して分析する。
1.6.2 Network plot
いずれのパッケージでもネットワークプロットを作成することができる。BUGSnet の Network plot が最も見た目がよさそうである。netmeta は、３群を塗りつぶしで明示したり、3D 表示することもできる。
また、pcnetmeta や viscomp などは独自のネットワークプロット関数を持っている。
1.6.3 要約
いずれのパッケージも、初期段階でペアワイズの要約を表示することができる。
この時にあまりに異質である（\(I^2 > 0.75\) など）ときは、ネットワークメタ分析を行わないという方法論もある。
また、Q 値などから非一貫性を調べることもできる。
1.6.4 モデルの作成と実行
ベイズ的ネットワークメタ分析では、モデルの作成と実行が必要になる。
作成時のパラメータには以下のようなものがある。
likelihoodとlink: アウトカム（SMD、MD、OR、RR）によって決まる。linearModel: 固定効果かランダム効果かを指定する。BUGSnet ではeffect
実行時のパラメータの代表的なものは、以下の通り。
- adapt: 実行した最初の表は一定しない。そのため最初の方のデータは捨てることになる。burnin とも言う。
- iter: 実行回数の総数。adapt 分も含まれるため、adapt より大きい数にする。
- thin: 実行したすべてのシミュレーションを使うのではなく、thin で指定した個数おきのデータを使用する。
まず、Trace 結果をプロットして、adapt が適当かどうかを判断する。最初のほうが一定していない場合、adapt を増やす。
さらに、モデルが収束していることを確認する必要がある。収束していない場合は、iter を増やすなどして収束させる。gelman mpsrf 値が 1 に近いほど収束している。