カテゴリー: Statistics

Rはフリーの統計解析プラットフォームおよびプログラミング言語ですが、今は相当普及していて、アカデミアでも、使っている人が大勢いると思います。The Comprehensive R Archive NetworkではWindows版のみならずMac版、Linux版が提供されています。普段Windows10のノートPCを使うことが多いのですが、周りにはMacを使う人も多く、私自身も以前はMacを主に使っていました。今も時々iMacを使うことがあります。しばらくMacでRを使うことはなかったのですが、今日はMacでRの最新バージョンをインストールして、さらにXQuartzをインストールし、Bayesian推定に使われるJAGS (Just Another Gibbs Sampler)をインストールし、Rのパッケージであるrjagsとgemtcをインストールして、動かして見ました。

Network meta-analysisネットワークメタアナリシスにはすでに紹介したことのある、OpenBUGSが使われる事が多いのですが、ネットワークメタアナリシス用のRのパッケージであるgemtcとpcnetmetaはRからrjagsを介してJAGSを動かしてMarkov Chain Monte Carlo (MCMC) simulationを実行するようになっています。gemtcはコントラストベースcontrast-basedのネットワークメタアナリシス、pcnetmetaはアームベースarm-basedのネットワークメタアナリシスができるRのパッケージです。

rjagsとJAGSを使う事で、WindowsだけでなくMacでも同じようにデータ解析ができるというのが最大の利点だろうと思います。OpenBUGSはもともとWinBUGSと呼ばれていたもので、Windowsでしか動きませんが、JAGSは同じコードがMacでもWindowsでも動くので、Macの利用者のことを考えてこのような選択になったのだろうと思います。Windowsであれば、RからOpenBUGSを動かしてデータのやり取りをするにはBRugsというパッケージがあります。

Mac用のRはR-4.0.0.pkgでバージョンは4.0.0になっています(2020.5.22時点)。CRANのメインページから左サイドバーのMirrorsを開き、Japanのサイトのいずれかを選びます。Download R for (Mac) OS Xのページを開き、R-4.0.0.pkgをダウンドードして他のMac用ソフトと同じようにインストールします。続いて、同じページにあるXQuartzのリンクを開き、XQuartz-2.7.22.dmgをダウンロードしてインストールします。

JAGSは上記のリンクを開いて、Downloadsの項のfiles pageのリンクを開き、JAGSのリンクを開き、4.xを開き、MacOSを開いて、JAGS-4.4.0.dmgをダウンロードしてインストールします。JAGSはRとは独立したソフトウェアです。

これで後はRを起動して、Rのパッケージでであるrjagsとgemtcをインストールします。個別に、自分でRの”パッケージとデータ”メニュー から”パッケージインストーラ”を選んで、そこからインストールすることもできますし、以下のスクリプトをRのファイルメニューから新規文書で開かれるエディタに貼り付けて実行（command+return)することでもインストールする事ができます。その行にカーソルを置いてcommand+returnで実行されます。（Windowsだとctrl+rまたは実行ボタン）。packneed=c( )のパッケージ名を変更する事で他のパッケージのインストールにも使えます。

packneed=c(“gemtc”,”rjags”);current=installed.packages();addpack=setdiff(packneed,rownames(current));url=”https://cran.ism.ac.jp/”;if(length(addpack)>0){install.packages(addpack,repos=url)};if(length(addpack)==0){print(“Already installed.”)}

ここではuseRsで私が紹介しているgemtcを動かすスクリプトを実行した結果を示します。”#4. GeMTCを用いるベイジアンネットワークメタアナリシス”のScriptをクリックして出てくるバーのD：の右の四角をクリックするとクリップボードに内容がコピーされるので、Excelを開いて貼り付け、データを確認して下さい。サンプルデータの形式でデータを用意すれば、（上が二値変数、下が連続変数のデータサンプルです）。二値変数アウトカムでも連続変数アウトカムでもどちらも解析可能です。効果指標の値が小さいほうが望ましい場合は、上の#4-1のScript、値が大きい方が望ましい場合は、下の#4-2のScriptを使います。W:がWindows用、M:がMac用のスクリプトになっています。

この例は連続変数のアウトカムで、下のサンプルデータを用いています。Node-splittingモデルのプロットも出力するので、直接比較と間接比較の効果推定値のインコヒレンスincoherenceの評価がやりやすくなると思います。

自分でデータを用意するときは、Excelで同じ形式で用意し、データ範囲を選択し（この例ではセルA1からM11）、コピー操作をしてから、Rに戻り以下のいずれかのスクリプトを実行します（エディターのスクリプト上にカーソルを置いて、command+enterを押す）。上のスクリプトは値が小さい方が望ましい効果指標の場合、下のスクリプトは値が大きい方が望ましい効果指標の場合です。これらのスクリプトはuseRsでScriptをクリックして出てくるバーのM:の右の四角をクリックするとクリップボードにコピーされます。（Windowsの場合はW:の右の四角をクリックします）。

exdat=read.delim(pipe(“pbpaste”),sep=”\t”,header=TRUE);source(“http://zanet.biz/med/stat/ma/nma_gemtc_simpleL.R”)

exdat=read.delim(pipe(“pbpaste”),sep=”\t”,header=TRUE);source(“http://zanet.biz/med/stat/ma/nma_gemtc_simpleH.R”)

これらのスクリプトでは、私が利用しているウェブサーバーに置いてあるスクリプトのファイルをsource( )関数で読み込んで実行することでクリップボードにコピーされたデータに対して、データ処理、解析が行われるようになっています。スクリプトのファイルはhttpプロトコールでインターネットを介してダウンロードされ、そのままRで実行されます。Rは自分のPCに置かれているスクリプトファイルだけでなく、インターネットを介し別のサーバーに置かれているスクリプトを実行できる仕組みがあり、またデータもそのような仕組みで取り扱え、遠隔での共同作業がやりやすいようになっている点が素晴らしいと思います。

最後に、二値変数のアウトカムの場合のデータサンプルとそのネットワークグラフを示します。

解析データは1行に一つの硏究のデータが含まれています。データの内t[,1], t[,2], t[,3]には治療番号のデータが入力されています。治療番号1が参照治療comparatorになります。必要に応じて、カラムや行を増やして使う事ができますna[]は治療アーム数です。

上の例では、2行目と3行目の研究は3アームで同時に３つの治療が比較されたランダム化比較試験です。それ以外の研究は2アームなのでt[,3]はNAになっています。行3,と行22の2つの研究は治療2と3を直接比較した研究です。治療2と3を治療1を介して間接比較できる研究は治療1と3を直接比較した行2, 4〜10, 13〜15, 18〜21の15研究と治療1と治療2を直接比較した行11,12,17の3研究があります。これら治療1と2の直接比較と治療2と3の直接比較から得られる効果推定値を差を求めることで、治療2と3の間接比較の効果推定値が得られることになります。そして、これら直接比較の効果推定値と間接比較の効果推定値の一致の程度を評価する事が必要になるので、ネットワークメタアナリシスに基づくエビデンスの確実性の評価は作業量が多くなり大変です。

ネットワークグラフで閉じられたループがあると、直接比較の効果推定値だけでなく間接比較の効果推定値が得られる事がわかります。この例では、2と3を結ぶ線より、3と1を結ぶ線が太く、そこに多くの研究がある事が示されています。２と1を結ぶ線もあるので、2と3を比較する間接比較のデータが得られ、2と3を比較するネットワーク効果推定値には間接比較のデータがかなり影響するであろうと考えられます。実際、Node-splittingモデルのプロットを見ると、3 vs 2のnetwork estimateはindirectのestimateの影響で、オッズ比が0.93 (0.22-3.9)から1.4 (0.63-3.3)と大きな値になっています。

なお、この解析結果は#4-1のデータを#4-2のスクリプトで解析したものです。つまり、オッズ比の値が大きい方が望ましい結果です。Rankogramを見ると、治療4が最も治療効果が高い確率が最も高い事が言えます。しかし、上のNode-splittingモデルの4 vs 3の比較を見ると、間接比較の点推定値がかなり大きく、直接比較との間の検定ではP=0.0819ですが、バイアス、非直接性、非移行性intransitivityなどによって偏りが生じている可能性をチェックする必要があると言えます。

このようにネットワークメタアナリシスに基づくエビデンスの確実性の評価はかなりの作業量が必要となる事がわかります。ネットワークがもっと大きく、複雑になるとペア比較の数も多くなり、ネットワークの構造が複雑になり、間接比較と直接比較の効果推定値を比較する事も大変な作業になります。

バイアスの効果はリスク比あるいはオッズ比などで表すことができます。バイアスの効果を定量的に推定できるのであれば、実際に得られた効果推定値をそれで調整して真の値により近づけることができます。

バイアスのモデル化についてはTurner RM 2009らのJournal of the Royal Statistical Society Series Aに発表された論文が精細な内容で参考になります。個別研究のバイアスドメイン・項目を評価しバイス調整をしたうえで、複数の研究のメタアナリシスの統合値を得る手法です。

バイアスの効果とバイアス調整について図にしてみました。

バイアスの効果には大きさと方向（過大評価、過小評価）と不確実性の3つの要素があります。

バイアスの効果を定量的解析する場合には、例えば、リスク比あるいはオッズ比を効果指標にした場合、自然対数に変換すると、正規分布に従うとみなせ、研究で得られた効果推定値は真の値とバイアス効果の加算された値になります。バイアス調整をする場合は、推定されるバイアス効果の大きさをリスク比で表し、その自然対数を引き算することで、真の値が得られはずです。分散は両者の合計になります。多変量正規分布の原理です。

メタアナリシスで得られた統合値の95%信頼区間の上限値（あるいは下限値）がNull effectあるいは臨床的閾値を超えるほどの大きさのバイアスの効果は簡単に計算できます。一般的に行われているバイアスリスクRisk of biasの評価結果でバイアスリスクが高く、バイアスの大きさがこの値よりも大きな効果を持っていると判断できる場合は、エビデンスの確実性が低くなると考えることができます。もし、そのようなバイアス効果の大きさがあり得ないほど大きいのであれば、ある程度のバイアスリスクがあっても、エビデンスの確実性は高いと考えることができるでしょう。

文献：
Turner RM, Spiegelhalter DJ, Smith GC, Thompson SG: Bias modelling in evidence synthesis. J R Stat Soc Ser A Stat Soc 2009;172:21-47. PMID: 19381328

Turner RM, Lloyd-Jones M, Anumba DO, Smith GC, Spiegelhalter DJ, Squires H, Stevens JW, Sweeting MJ, Urbaniak SJ, Webster R, Thompson SG: Routine antenatal anti-D prophylaxis in women who are Rh(D) negative: meta-analyses adjusted for differences in study design and quality. PLoS One 2012;7:e30711. PMID: 22319580

Darvishian M, Gefenaite G, Turner RM, Pechlivanoglou P, Van der Hoek W, Van den Heuvel ER, Hak E: After adjusting for bias in meta-analysis seasonal influenza vaccine remains effective in community-dwelling elderly. J Clin Epidemiol 2014;67:734-44. PMID: 24768004

Dias S, Sutton AJ, Welton NJ, Ades AE: Evidence synthesis for decision making 3: heterogeneity–subgroups, meta-regression, bias, and bias-adjustment. Med Decis Making 2013;33:618-40. PMID: 23804507

Thompson S, Ekelund U, Jebb S, Lindroos AK, Mander A, Sharp S, Turner R, Wilks D: A proposed method of bias adjustment for meta-analyses of published observational studies. Int J Epidemiol 2011;40:765-77. PMID: 21186183

Wilks DC, Mander AP, Jebb SA, Thompson SG, Sharp SJ, Turner RM, Lindroos AK: Dietary energy density and adiposity: employing bias adjustments in a meta-analysis of prospective studies. BMC Public Health 2011;11:48. PMID: 21255448

Wilks DC, Sharp SJ, Ekelund U, Thompson SG, Mander AP, Turner RM, Jebb SA, Lindroos AK: Objectively measured physical activity and fat mass in children: a bias-adjusted meta-analysis of prospective studies. PLoS One 2011;6:e17205. PMID: 21383837

Schnell-Inderst P, Iglesias CP, Arvandi M, Ciani O, Matteucci Gothe R, Peters J, Blom AW, Taylor RS, Siebert U: A bias-adjusted evidence synthesis of RCT and observational data: the case of total hip replacement. Health Econ 2017;26 Suppl 1:46-69. PMID: 28139089

Doi SA, Barendregt JJ, Onitilo AA: Methods for the bias adjustment of meta-analyses of published observational studies. J Eval Clin Pract 2013;19:653-7. PMID: 22845171

McCarron CE, Pullenayegum EM, Thabane L, Goeree R, Tarride JE: Bayesian hierarchical models combining different study types and adjusting for covariate imbalances: a simulation study to assess model performance. PLoS One 2011;6:e25635. PMID: 22016772

McCarron CE, Pullenayegum EM, Thabane L, Goeree R, Tarride JE: The importance of adjusting for potential confounders in Bayesian hierarchical models synthesising evidence from randomised and non-randomised studies: an application comparing treatments for abdominal aortic aneurysms. BMC Med Res Methodol 2010;10:64. PMID: 20618973

Phillippo DM, Dias S, Welton NJ, Caldwell DM, Taske N, Ades AE: Threshold Analysis as an Alternative to GRADE for Assessing Confidence in Guideline Recommendations Based on Network Meta-analyses. Ann Intern Med 2019;170:538-546. PMID: 30909295

Caldwell DM, Ades AE, Dias S, Watkins S, Li T, Taske N, Naidoo B, Welton NJ: A threshold analysis assessed the credibility of conclusions from network meta-analysis. J Clin Epidemiol 2016;80:68-76. PMID: 27430731

Phillippo DM, Dias S, Ades AE, Didelez V, Welton NJ: Sensitivity of treatment recommendations to bias in network meta-analysis. J R Stat Soc Ser A Stat Soc 2018;181:843-867. PMID: 30449954

model {
for(i in 1:tS) {
p[i]<-phi(mu[t[i]]+ vi[s[i], t[i]]) # model
r[i]~dbin(p[i], totaln[i]) # binomial likelihood
}
for(j in 1:sN){
vi[j, 1:tN]~dmnorm(mn[1:tN], T[1:tN,1:tN]) # multivariate normal distribution
}
invT[1:tN, 1:tN]<-inverse(T[,])
for (j in 1:tN){
mu[j]~dnorm(0, 0.001)
sigma[j]<-sqrt(invT[j, j])
probt[j]<-phi(mu[j]/sqrt(1+invT[j, j]))
#population-averaged treatment specific event rate
}
T[1:tN,1:tN]~dwish(R[1:tN, 1:tN], tN) # Wishart prior
for (k in 1:tN) {
rk[k]<- tN+1-rank(probt[],k) # ranking
best[k]<-equals(rk[k],1) # prob {treatment k is best}
}
for (j in 1:tN){ # calculation of RR, RD and OR
for (k in (j+1):tN){
RR[j, k] <- probt[k]/probt[j]
RD[j, k] <- probt[k]-probt[j]
OR[j, k] <- probt[k]/(1-probt[k])/probt[j]*(1-probt[j])
}
}
}
#DATA
list(tN=4,sN=24,tS=50,mn=c(0, 0, 0, 0),R=structure(.Data=c(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1),.Dim=c(4,4)))
s[] t[] r[] totaln[]
1 1 9 140
1 3 23 140
1 4 10 138
2 2 11 78
2 3 12 85
2 4 29 170
3 1 75 731
3 3 363 714
4 1 2 106
4 3 9 205
5 1 58 549
5 3 237 1561
6 1 0 33
6 3 9 48
7 1 3 100
7 3 31 98
8 1 1 31
8 3 26 95
9 1 6 39
9 3 17 77
10 1 79 702
10 2 77 694
11 1 18 671
11 2 21 535
12 1 64 642
12 3 107 761
13 1 5 62
13 3 8 90
14 1 20 234
14 3 34 237
15 1 0 20
15 4 9 20
16 1 8 116
16 2 19 146
17 1 95 1107
17 3 134 1031
18 1 15 187
18 3 35 504
19 1 78 584
19 3 73 675
20 1 69 1177
20 3 54 888
21 2 20 49
21 3 16 43
22 2 7 66
22 4 32 127
23 3 12 76
23 4 20 74
24 3 9 55
24 4 3 26
END