MacでR, JAGS, rjagsを使うNetwork meta-analysisをやってみる

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のパッケージであるrjagsgemtcをインストールして、動かして見ました。

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アームで同時に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を結ぶ線が太く、そこに多くの研究がある事が示されています。2と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などによって偏りが生じている可能性をチェックする必要があると言えます。

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

New PubMed 2020 April 1

現在2020年4月1日の時点でのNew PubMedについての情報です。

検索結果の画面で、Display Optionsを見ると、Sort byはBest matchがデフォルトになっています。Display optionsのボタンは画面右上のSearchの下にあります。(現在Display optionsの表示が以下の歯車のボタンに変更されています。)

これをMost recentに変えることができます。今までのデフォルトはMost recentでした。Most recentにした場合、すぐ右に昇順か降順かを指定する階段のようなアイコンのボタンがありますので、目的に応じて設定します。

また、1ページに表示する文献数はデフォルトで10になっていますが、Per pageからこれも変更することができます。

検索結果をファイルとしてダウンロードして保存する場合、Saveボタンをクリックします。Formatには当初RIS形式がありましたが、現時点ではなくなっています。FormatのPubMedが従来のMEDLINE形式に相当します。ファイルとしてダウンロードすると、ファイル名の拡張子が.nbibとなっていますが、通常のテキストファイルで、拡張子.txtと同じです。

また、CSV形式でダウンロードすると、Excelなどで開くことができます。

一行にひとつの文献の情報が含まれています。PMID, Title, Authors, Citation, First Author, Journal/Book, Publication Year, Create Date, PMCID, NIHMS ID, DOIのラベルがついています

この中で、DOIはDigital Object Identifierの略ですが、スラッシュの前がDOI財団が付与した番号、後ろがジャーナルがつけた番号です。DOIは文献に1対1で対応しており、重複はありません。DOIをコピーしてGoogleあるいはBingで検索すると、検索結果の一覧の上位にその文献へのリンクが表示されますので、直接ジャーナルのその論文を開くことができます。また、ブラウザーのURLのフィールドにhttps://doi.org/に続けて、DOIを書いて、Enter keyを押すと、直接ジャーナルのその文献を開くことができます。

Send toからもCitation mangerを選択して、Create fileでもPubMed形式のファイルとしてダウンロードできるようになっています。

各文献の情報の下に表示されているCiteをクリックすると引用のための形式で小さなフィールドに表示されます。CopyボタンをクリックするとPCのクリップボードにコピーされるので、執筆中の原稿にそのまま貼り付けることができます。形式も、4種類用意されています。

New PubMedに関する以前の投稿はこちらです。

SARS-CoV-2に関する文献検索

PubMedを以下の検索式で検索してみます。Searchの部分をクリックするとPubMedが開かれて検索結果が表示されます。新規タブで新規ウインドウに表示されます。

検索式1 :フリー・フルテキストでヒトを対象にした英語と日本語の2019-2020年の文献に限定しています。重要な論文でもフリー・フルテキストでない論文は出てきません。 -> Search 1
検索式2:それも含めた検索はこちら。->Search 2

検索式1: (sars-cov-2 OR covid-19 OR “new coronavirus” OR “novel coronavirus”) AND 2019:2020[dp] AND humans[mh] AND (english[la] OR japanese[la]) AND “loattrfree full text”[sb]
または
検索式2: (sars-cov-2 OR covid-19 OR “new coronavirus” OR “novel coronavirus”) AND 2019:2020[dp] AND humans[mh] AND (english[la] OR japanese[la])

humans[mh]を除いた検索式は以下の通りになります。ヒトを対象にした研究に限定しないので、引き出される文献数が多くなります。Searchの部分をクリックしてください。
検索式3:フリー・フルテキストで英語と日本語の2019-2020年の文献です。->Search 3
検索式4:フリー・フルテキスト以外の論文も含める検索はこちら。->Search 4
検索式5:出版年月日とMeSH用語による制限のない、引き出される文献数が最も多い検索式はこちら。 ->Search 5

検索式3: (sars-cov-2 OR covid-19 OR “new coronavirus” OR “novel coronavirus”) AND 2019:2020[dp] AND (english[la] OR japanese[la]) AND “loattrfree full text”[sb]
または
検索式4: (sars-cov-2 OR covid-19 OR “new coronavirus” OR “novel coronavirus”) AND 2019:2020[dp] AND (english[la] OR japanese[la])
または
検索式5:sars-cov-2 OR covid-19 OR “new coronavirus” OR “novel coronavirus”

Searchの部分をクリックした場合、その検索式をPubMed側にGETメソッドで送り、結果が表示されるので、その時点の最新の検索結果が表示されます。

さて、検索式1の結果を見ると、中国の武漢の医師を中心としたグループからGRADEアプローチに準じた”Rapid advice guideline”が発表されています。

Jin YH, Cai L, Cheng ZS, Cheng H, Deng T, Fan YP, Fang C, Huang D, Huang LQ, Huang Q, Han Y, Hu B, Hu F, Li BH, Li YR, Liang K, Lin LK, Luo LS, Ma J, Ma LL, Peng ZY, Pan YB, Pan ZY, Ren XQ, Sun HM, Wang Y, Wang YY, Weng H, Wei CJ, Wu DF, Xia J, Xiong Y, Xu HB, Yao XM, Yuan YF, Ye TS, Zhang XC, Zhang YW, Zhang YG, Zhang HM, Zhao Y, Zhao MJ, Zi H, Zeng XT, Wang YY, Wang XH; , for the Zhongnan Hospital of Wuhan University Novel Coronavirus Management and Research Team, Evidence-Based Medicine Chapter of China International Exchange and Promotive Association for Medical and Health Care (CPAM). A rapid advice guideline for the diagnosis and treatment of 2019 novel coronavirus (2019-nCoV) infected pneumonia (standard version). Mil Med Res. 2020 Feb 6;7(1):4. doi: 10.1186/s40779-020-0233-6. PMID: 32029004; PMCID: PMC7003341.

スピードがすごいですね。診療ガイドライン作成方法論についてはだいぶ前から学習を進めていた、あるいは方法論の専門家が育成されていたということなのでしょう。

PubMedを開くとCenters for Disease Control and Prevention (CDC)のCoronavirus (COVID-19)に関するページへのリンクとNational Institutes of Health (NIH)のCoronavirus (COVID-19)の研究に関するページへのリンクが表示されます。

バイアスの効果とバイアス調整

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

バイアスのモデル化については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