Systematic Review – ページ 4

DTAのメタアナリシス-useRs

useRsのサイトでは、さまざまなメタアナリシスのためのR用のスクリプトを提供しています。

Rのlme4パッケージのglmer()関数を用いる回帰分析でBivariate modelのDiagnostic Test Accuracy (DTA)診断精度研究のメタアナリシスを二項分布Binomial distributionを用いて行うスクリプトを掲載しました。このページの#3-2. Bivariate modelで二項分布を用いる診断精度（Diagnostic Test Accuracy, DTA)研究のメタアナリシスです。また、DTA研究のバイアスリスク評価には、QUADAS-2が用いられることが多いですが、それに準じた評価シートDTA_sheet.xlsxも用意しました。Rのスクリプトを含んでいます。また、madaを用いるスクリプトを含むシートも含めています。

図1. SORC曲線。

図1のようなSROC曲線のプロットを出力します。感度、特異度のForest plotも作成します。コンソールにはさまざまな解析結果と、感度、特異度と95%信頼区間、その下にI2値が出力されます。それとともに、クリップボードに各研究の有病率、感度、特異度と95%信頼区間を格納するので、Excelなどに貼り付けることができます。

続けて2行目を実行すると、感度、特異度の統合値と95%信頼区間、統合値のDiagnostic Odds Ratio （DOR）、陽性尤度比、陰性尤度比とこれらの95%信頼区間およびI二乗値と、RevManへ渡すことによって、SORC曲線を描くパラメータを出力します。それとともにクリップボードにもこれらの値を格納するので、Excelなどに貼り付けることができます。

メタアナリシスの部分のスクリプトは以下の通りです。

ma_res = glmer(formula=cbind(true, n – true ) ~ 0 + sens + spec + (0+sens + spec|Study_ID), data=Y, family=binomial)

文献：
Whiting PF, Rutjes AW, Westwood ME, Mallett S, Deeks JJ, Reitsma JB, Leeflang MM, Sterne JA, Bossuyt PM; QUADAS-2 Group. QUADAS-2: a revised tool for the quality assessment of diagnostic accuracy studies. Ann Intern Med. 2011 Oct 18;155(8):529-36. doi: 10.7326/0003-4819-155-8-201110180-00009. PMID: 22007046.

Cochrane Handbook for Systematic Reviews of Diagnostic Test Accuracy version 2.0, 2022. Link

DTAのメタアナリシス

Diagnostic Test Accuracy (DTA) 診断精度研究のMeta-analysis (MA)メタアナリシスは介入の効果のメタアナリシスとは異なります。手法が異なるだけでなく、結果の解釈や活用法も異なります。

DTA MAの結果、感度・特異度の統合値と95%信頼区間、陽性尤度比、陰性尤度比、診断オッズ比（Diagnostic Odds Ratio, DOR）、Summary Receiver Operating Characteristic (SROC) curve、およびSROCの曲線下面積（Area under the curve, AUC）などの推定値が得られます。

DTA MAの統計学的モデルとして、Reitsma のBivariate model二変量モデル、階層モデルである Rutter & GatsonisのHierarchical Summary Receiver Operating Characteristic (HSROC)モデルの使用が推奨されています。共変量を用いない場合は二変量モデルとHSROCモデルは数学的に同じものです（Arends LR;Harbord RM) 。

DTA MAに関する書籍としては、Biondi-Zoccai G ed. Diagnostic Meta-Analysis: A Useful Tool for Clinical Decision-Making. Springer, Cham, Switzerlandが包括的な内容で、有用と思います。また、Cochrane Handbook for Systematic Reviews of Diagnostic Test Accuracyは2022年度Version 2が発表されており、包括的な内容で、SAS、R用のスクリプトが具体的に解説されており、有用だと思います。

Reitsma JBのオリジナルの論文では、近似正規分布を用いるため、症例数が少なく、感度・特異度が95%程度に高い場合には、誤差が大きくなることをChu Hらが指摘しており、Cochraneはコクランのシステマティックレビューに使用することは推奨しないとしています。また、ゼロイベントのある場合に0.5を加算して補正する方法もずれを生じます。そのため、広く使用されているRのパッケージであるmadaの使用は推奨できないと述べられています。実際に、Chu Hらの論文でシミュレーションの結果を見ると、症例数が25例以下、感度・特異度が95%以上になると二項分布を用いる場合とずれが大きくなるのは確かですが、Chu HのLetterに対するReitsma JBの返答で述べられているように、多くの場合臨床上の問題を生じるほどではないと考えられます。

Cochrane Handbook for Systematic Reviews of Diagnostic Test AccuracyのAppendix 14では同じ対象者でCTとMRIを施行し、診断能を比較した5つの研究のDTA MAの例が記載されており、Rのlme4パッケージのglmer()関数を用いて、GLMM (Generalized Linear Mixed Effects Model)で二項分布による回帰モデルを用いています。Appendix 12では同じ対象者で二つの診断法を実施したのではなく、別の対象者でそれぞれの診断法の感度・特異度を測定した研究をもとに、二つの診断法の診断能を比較するための解析法が記載されています。こちらは、間接的な比較という表現が使われており、同じ対象者で直接比較した研究も含めて解析できる方法が示されています。

感度のロジットと偽陽性率のロジットを従属変数とし、各研究の感度のロジットと偽陽性率ロジットで回帰モデルを作っています。ロジットはオッズの自然対数、つまりlogit(se)=ln[se/(1-se)], logit(sp)=ln[sp/(1-sp)], logit(fpr)=-logit(sp)です。

具体的には図1に示すようなデータに対して、回帰分析を行い、ロジットのExponentialでオッズに変換し、さらにオッズ/（１＋オッズ）で割合に変換して感度・特異度の値を求めています。図1の例は同一症例でCTとMRIを施行しており、これら二つの検査法の感度・特異度の統計学的な比較が行われています。

通常、元になるデータは各研究IDとTP, FP, FN, TN （True Positive, False Positive, False Negative, True Negative)の人数のデータと検査法のデータです。それらから、図1の形式のデータフレームを作成し、以下のスクリプトで解析を行います。

###Comparison of sensitivity and specificity between two tests done in the same subjects###
###Random-effects meta-analysis with bivariate model using binomial distribution###
library(lme4)
library(lmtest)
###Y is a data frame as shown in Fig. 1###
(B = glmer(formula=cbind(true, n – true) ~ 0 + seCT + seMRI + spCT + spMRI + (0+sens + spec|Study_ID), data=Y, family=binomial))
(C = glmer(formula=cbind(true, n – true) ~ 0 + sens + spCT + spMRI + (0+sens + spec|Study_ID), data=Y, family=binomial))
###Is there a statistically significant difference in sensitivity between CT and MRI?
lrtest(B,C)

###Is there a statistically significant difference in specificity between CT and MRI?
lrtest(B,D)

実際には各研究IDとTP, FP, FN, TN のデータから図1の形式のデータフレームを作成する部分のスクリプトも必要ですが、ここではlmer()関数の回帰モデルと図1のデータラベルの関係を考えて、回帰分析とメタアナリシスの関係を考えるきっかけになればと思います。

なお、ここでは触れませんでしたが、感度・特異度の統合値と95%信頼区間の計算はCTとMRIについて別々に回帰分析を行った結果から算出します。

文献：
Chu H, Cole SR: Bivariate meta-analysis of sensitivity and specificity with sparse data: a generalized linear mixed model approach. J Clin Epidemiol 2006;59:1331-2 author reply 1332-3. doi: 10.1016/j.jclinepi.2006.06.011 PMID: 17098577

Reitsma JB, Glas AS, Rutjes AW, Scholten RJ, Bossuyt PM, Zwinderman AH: Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. J Clin Epidemiol 2005;58:982-90. doi: 10.1016/j.jclinepi.2005.02.022 PMID: 16168343

Rutter CM, Gatsonis CA: A hierarchical regression approach to meta-analysis of diagnostic test accuracy evaluations. Stat Med 2001;20:2865-84. PMID: 11568945

Harbord RM, Deeks JJ, Egger M, Whiting P, Sterne JA: A unification of models for meta-analysis of diagnostic accuracy studies. Biostatistics 2007;8:239-51. doi: 10.1093/biostatistics/kxl004 PMID: 16698768

Arends LR, Hamza TH, van Houwelingen JC, Heijenbrok-Kal MH, Hunink MG, Stijnen T: Bivariate random effects meta-analysis of ROC curves. Med Decis Making 2008;28:621-38. doi: 10.1177/0272989X08319957 PMID: 18591542

Rのsource()関数の応用

統計解析プラットフォームとして広く使われており、プログラミング言語であるRにはsource(URL)という関数があります。Rのスクリプトをウェブサーバーにアップロードしておき、そのファイルのURLをsource()関数の引数として設定し、Rでそれを実行させると、ウェブサーバーからスクリプトが読み込まれて実行されます。

解析対象のデータはクリップボードにコピーしておいて、それをread.delim()関数で変数に読み込んで、その変数に対して処理を行うスクリプトを用意しておくことで、各自用意したデータを解析するシステムが作れます。解析したいデータはCSVファイルで用意して、それを読み込ませて解析する方法も可能です。

source()関数で読み込むRのスクリプトはウェブサーバーに置いてありますが、それを読み込んで実行させるスクリプトは各自のPCにインストールしたRで実行します。それを読み込んで実行させるスクリプトは、各自で保存しておくこともできますし、別のウェブページから提供することも可能です。useRsもそのようなウェブページのひとつですが、別の方法でも同じことが可能です。

例えば、メタアナリシスのさまざまなモデルに対応したスクリプトを用意して、ウェブページで選択して、実行させるというようなことができます。

システマティックレビューに必要なExcelシートを集め、メタアナリシスのためのRのスクリプトを含めたBookを作りました。2023_excel_book_for_sr_v.0.96.xlsxから自分の必要なシートをシートのタブを右クリックして、My SR Book.xlsxにコピーを追加し、目次のシートに各シートの名前を入力してそれらへのリンクを設定し、利用することができます。

リスク比、オッズ比、リスク差、ハザード比、平均値差、標準化平均値差 (RR, OR, RD, HR, MD, SMD)のメタアナリシスをRのmetaforパッケージで、診断精度研究のメタアナリシスをRのmadaパッケージを利用して実行します。Forest plotの作成にはRのforestplotパッケージを用いる場合もあります。Ｒをインストールし、必要なパッケージをＲでインストールしてから使います。

正規分布に従う複数の変数に重みの値を掛け算した値の総和の分散

Linear combination of n variables with normal distributions with weights. 重みを係数coefficientsと言い換えてもいいです。

すぐにはピンとこないかもしれませんが、例えば、日本人の夫婦の身長の合計、つまり二人の身長の合計の平均値と分布を知りたいとします。それぞれ夫と妻の身長の分布が正規分布に従っているとします。日本人の夫の身長の平均値と分散、分散は標準偏差の二乗です（個々の値と平均値の差の二乗値の平均値が[標本]分散です）、が分かっていて、妻の身長の平均値と分散が分かっているとします。夫と妻のペアはランダムな組み合わせだとすると、（実際には背の高い妻は背の高い夫がいるというようなある程度の相関があるかもしれませんが、まずは妻と夫の身長の間にはそのような相関が無い、つまり共分散が0と仮定しておきます）、夫の身長と妻の身長の合計値の分布はどうなるでしょうか？平均値はそれぞれの平均値の和になり、分散はそれぞれの分散の和になります。

日本人の妻の集団からランダムに一人抜き出し、日本人の夫の集団からランダムに一人抜き出し、かれらの身長を測って、合計するということを繰り返した場合、その合計値の分布はどうなるかという風に考えてみて下さい。

このような計算がどいう時に使われているかというと、突然話が飛びますが、例えば分散逆数法Inverse-variance methodによるメタアナリシスの際に統合値の分散を求める際に使われています。計算法を図1に示します。今度は、妻と夫の身長という二つの変数ではなく、研究の数分の変数を扱います。

一つの研究の効果推定値（リスク比、オッズ比、ハザード比の自然対数、連続変数であれば平均値）に対して、その分散の値の逆数を重みとして掛け算して、その総和を重みの総和で割り算すると統合値が得られます。

その統合値の分散を計算するにはどうするか？分散の逆数の総和の逆数をΣの内側に移動させることができるので、各研究の効果推定値に掛け算される値は、その総和で各研究の重みの値、つまり、各研究の効果推定値の分散の逆数を割り算した値になります。これが、係数として各研究の効果推定値に掛け算されているとみなせます。そして、統合値の値の分散はこの係数の二乗値を各研究の効果推定値の分散に掛け算した値の総和になります。

ただし、各研究の効果推定値は独立していて、相関が無いことが前提です。また、それぞれの効果推定値は正規分布に従うことが前提です。つまり、リスク比、オッズ比、ハザード比の自然対数、連続変数であればその値が、正規分布に従うことを前提としています。

さて、正規分布に従う複数の変数に重みの値を掛け算した値の総和の分散のより一般化した計算法を図2に示します。

図2．正規分布に従う複数の変数に重みの値（係数a）を掛け算した値の総和の分散の計算法。

各変数のペアで相関がない場合で、この図で示すCov、つまり、共分散の値が0の場合が、上で述べた計算です。図2の2つ目と3つ目の式でCov（X_i,X_j）=0となるので、各研究の効果推定値に係数の二乗を掛け算した値の総和が統合値の分散になります。また、係数a_iが１で共分散が０の場合は、それぞれの変数の分散を合計すればX_iの合計の分散になります。最初に述べた、妻と夫の身長の和の分布の分散の計算の場合はこれに相当します。

メタアナリシスの分散逆数法の統合値の分散は各研究の分散の逆数の総和の逆数になりますが、実は、ここで二つの図で示した計算式から証明することができます。図３、図4にそれを示します。

図4．分散逆数法メタアナリシスの統合値の分散の計算。各研究の効果推定値の分散の逆数の総和の逆数が統合値の分散になることを示す。

また、これらの計算には、分散共分散行列を計算に用いることもでき、変数間に相関がある場合にも対応できます。その際には、行列計算の知識が必要になります。

今回解説した、係数を掛け算した正規分布に従う変数の平均値の分散の計算は価値観で重みづけした効果推定値の総和、すなわち正味の益（net benefit, benefit-harm balance）の分散の計算でも、バイアス効果で調整した統合値の分散の計算でも用いられます。

一般化した言い方をすると、”正規分布に従う複数の変数に重みの値を掛け算した値の総和の分散”ということになります。重みづけ平均値とあわせて理解しておく必要があると思います。

そして得られた分散の値の平方根に1.96を掛け算してブラスマイナスすると95%信頼区間が得られます。さらに、例えば正味の益が0以上の確率やある閾値以上の確率を計算することもできます。