カテゴリー: Statistics

Medical statistics principles,
analyes using R, OpenBUGS.

Cochrane Risk of Bias Tool ver.2.0と評価用ウェブツール

ランダム化比較試験のバイアスリスクの評価ツールとしてCochrane Risk of Bias Tool ver.2.0が2019年8月に最終版が完成し、広く使われるようになってきました （Current versionへの Link）。評価ドメインが5つに限定され、概念は同じですが、名称が解説的なものに変わりました。シグナリングクエスチョンに答えてゆくことで、アルゴリズムに従って、評価がLow, Some concerns, Highのいずれかに決まる仕組みになっています。（Version 1.1ではLow, Unclear, High, でした。）アルゴリズムによる自動判定のツールも含んだExcel macroも公開されています。なお、アルゴリズムによる判定と評価者の判定が違う場合は、評価者の判定を優先します。また、クラスターランダム化比較試験、クロスオーバー試験用のツールは別になっています。

シグナリングクエスチョンに答えていけば、判定ができるようになっていますが、アルゴリズムのどこに合致するかを見ないといけませんし、初心者にはシグナリングクエスチョンそのものの意味が分からない場合もあると思います。従来の方法と比べると、判定の手順はかなり労力を要するものになっていると思います。特にアルゴリズムは記憶できる範囲を超える量なので、RoB 2ガイダンスドキュメントを参照しながら判定する必要があり、時間もかかります。

その様な作業を容易にするために、シグナリングクエスチョンの行をクリックすると表示される解説を見ながらシグナリングクエスチョンに答えていくことで、アルゴリズムに従い、自動的に判定するWeb toolを作成しました。（プログラミングはJavaScriptを用いています）。Cochraneのウェブサイトの紹介、RoB 2のガイダンスドキュメントの紹介、RoB 2ウェブツールの紹介、評価シートとの連携、評価シートからR＋metafor, forestplotによるメタアナリシスの実行まで解説した動画（11分13秒）を作成しました。

メタアナリシスをすべきかどうかそれが問題だ

AHRQ Agency for Healthcare Research and Qualityから2018年にQuantitative Synthesis-An Updateというレポートが出されています。（2020年に訂正あり）。執筆者は12名、ピアレビューは15名が担当しています。

メタアナリシスをしようと思った時に、各研究間に異質性がある場合、それらを単にプールして統合値を求めることに当然疑問を抱くと思います。メタアナリシスをすべきかどうかそれが問題だという時に役に立つ情報です。

5章から構成され、以下のテーマが取り上げられています。1. 試験結合の決定、2. エフェクトサイズデータの使用の最適化、3. 研究結合のための統計学的モデルの選択、4. 統計学的異質性の定量化・検定・探索、5. ネットワークメタアナリシス。

今回は、第1章 試験結合の決定の概要を紹介します。メタアナリシスをすべきか決めるためのフローチャートが提示されているので、スライドと解説としてまとめてみました。最後にQuizが2つあります。Link to ”Pooling decision tree”

メタアナリシスをしようと思った時に出てくるさまざまな疑問に答えられる有用な情報になると思います。これに沿って考えることは役に立つと思います。

もう一つ重要なポイントは、この報告書をWorking Groupがどのように作成したかについてです。かれらは、「14カ月にわたって、テレカンファレンスで、①プロジェクトの方向性、スコープ、②仕事の分担・協働、③データの収集・分析について議論し、④ドラフトについて議論し編集を行った。最初の会議では、ドラフトのアウトライン、タイムラインについて議論し、合意形成の方法について同意を得た。」

「大グループは2週おきにテレカンファレンス、章別に小グループに分かれ、それとは別に作業の調整をし、文献レビューの結果について議論し、それぞれの章のドラフトを書いた。後で、全章をひとつにまとめ、2週おきに議論した。」

「文献検索専門家が章ごとの文献検索を行った。それらを全章に渡ってまとめた。16,000の文献を含む、AHRQ SRC Methods Libraryを用いた。追加の文献は、最近のSR、レビュー、エディトリアル、エキスパートのレビューから探した。」

要点は：
・元になる文献のレビューの段階から、チームで議論をしている。
・合意形成の手順についてあらかじめ議論し、それを共有した。
・2週おきの議論を繰り返した。

これほどのエンゲージメントは大変だろうと思いますが、チームワークで作業をする場合これくらいする必要があるだろうと思います。このような手順であれば、論文の解釈についても議論され、参加者の理解が深まり、共有されることになるでしょう。誰か分かっていそうな人が書いた原稿を誰か別の人が査読してＯＫなら終了というのとは全く違うと思います。

Bias – DAG Directed Acyclic Graphを用いる分析

バイアスとは

バイアスは「研究結果の系統的な偏り、あるいは、推定の真実からの系統的な偏り」と定義されています。系統的とは？ 偶然による偏りに対して、偶然起きる偏りではないので、系統的な偏りと言います。系統的systematic vs 偶然random という考え方をしているということです。偶然による偏りは、統計学的に説明可能でサンプルサイズが小さいほど大きくなります。バイアスによる偏りは経験的empiricalなデータは限られており、バイアス効果の大きさと方向（過大評価か過小評価か）については評価者が推定せざるをえないことがほとんどです。

実際に得られた結果の効果推定値および信頼区間Confidence intervalから、偶然による偏りだけなのか、バイアスによる偏りなのか、両方が混ざっているのかを見分けるのは難しいです。なお、頻度論派Frequentistの95%信頼区間は同じことを繰り返したら95%の場合はその範囲に本当の値が含まれるということを意味しています。5％の場合は、その範囲外になります。ベイジアンBayesianのアプローチであれば、95%確信区間Credicble intervalは真の値を95%の確率で含む範囲です。

さて、バイアスは多数存在しますが、バイアスの原因となる要素がその研究にあるかどうかを判断することは可能です。それがその研究にあれば、その結果あるいは効果推定値はバイアスのため、偏っている可能性が高くなります。そのようなアプローチが必要になります。

真の値を推定することは、もしバイアス効果の大きさが推定できれば、それによって調整することで、可能になります（Quantitative bias analysis）。例えば、得られた結果のリスク比RR＝0.8で、バイアスの効果がRR＝0.9であれば、真の効果推定値はRR ＝ 0.8/0.9 = 0.89 =exp[ln(0.8)-ln(0.9)] = exp[-0.223-(-0.105)] = exp(-0.118)です。結果RRの分散とバイアス効果のRRの分散を合計するとバイアスで調整した真の効果の推定値の分散はそれらの合計になります。バイアス効果のRRの分散はレビュアが推定値を設定する必要があります（Turner RM 2009）。＊ln 自然対数；exp Exponential （Excelの関数の表示ln(), exp()と同じです, なお、Rでは自然対数はlog()）。

また、バイアスの効果の大きさがどれ位あれば望ましい結果が望ましくない結果に反転するかも知ることができます（Bias adjustment thresholds analysis）（Phillippo DM 2018）。しかしながら、そこまで分析を行うケースは少なく、Cochrane risk of bias toolのように、大きなバイアス（あるいは”実質的なバイアス” material bias）のリスク＝可能性がどれくらいあるかを評価することが一般的です。

バイアスは交絡バイアス、選択バイアス、情報バイアスの３つに分類されるのが一般的です（Lash TL 2021）。

図 バイアスの３つのタイプとDAG。
矢印の始点の方が原因で終点の方が結果を示します。IはIntervention介入またはEと書いてExposure要因暴露、CはCommon cause共通原因＝Confounder交絡因子、OはOutcomeアウトカムを表します。OiはIntermediate outcome中間アウトカム、OeはEndpointエンドポイントです。I*は測定された介入、O*は測定されたアウトカムを表します。Uは測定誤差を引き起こす因子で、多くの場合不明です。
共通原因は交絡バイアス、共通効果は選択バイアス、系統的測定誤差は情報バイアスに対応します。
I←C→O; I→S→Oi←Oeと記述しても共通原因、共通効果であるこをと表現できます。

バイアスの効果を分析する際に、有向非巡回グラフDirected Acyclic Graph (DAG)が用いられることがあります。DAGは因果関係図Causal Diagramと呼ばれることもあります。それぞれの変数がバイアスになりうるかを検討するのに有力なツールです。図には、Luijendijk HJ 2020の論文に基づいて、汎用性のあるDAGの３つのバイアスのタイプについて示してあります。

共通効果→選択バイアス

バイアスの議論の際に、よく引用されるBerkson’s biasは選択バイアスとして知られていますが、共通効果のDAGを使って説明されます(Westreich D 2012)。例えば、クリニック受診患者を対象として糖尿病と認知症の関係を分析した場合、受診の原因が糖尿病の場合もあり認知症の場合もあります。クリニック受診が共通効果になります。クリニックを受診しない患者は選択せず、クリニック受診患者だけを選択して糖尿病と認知症の関係を分析するとバイアスが生じます。このようなバイアスは、前向き研究でも後ろ向き研究でも、観察研究でもランダム化比較試験でも起きる可能性があります。

共通効果で条件付けされるすなわちconditioned on （分類される）ある層だけを対象として選択したり、共通効果で調整するとバイアスが生じます。すなわち、図に示す変数Sに基づく層のひとつを分析する、あるいは変数SでIとOeの関係を調整した分析を行うとバイアスが生じます。このようなバイアスの結果は過大評価になる場合も、過小評価になる場合もあり、例えば、上記の例だと糖尿病は認知症のリスクを高めることはないという結果が得られる可能性があります。

SがOiの影響をうけて変動する結果生じるバイアスとして、脱落によるバイアス、症例減少バイアスAttrition biasがあります。その場合、介入群の方で脱落がより多いというような、Iからの影響も受けます。

コンシールメントがないためにランダム化が歪んで起きるバイアスも共通効果で説明できるはずです。皆さんも考えてみて下さい。

共通効果はCollider 合流因子とも呼ばれ、共通効果によるバイアスは選択バイアスに相当します。選択バイアスをCollider biasと呼ぶこともあります。

共通原因→交絡バイアス

共通原因は交絡因子に相当し、介入とアウトカムの両方に関係がある因子です。因果関係はI←C→Oの方向です。共通原因＝交絡因子の影響によって介入が変動し、アウトカムも変動する場合に、介入のアウトカムに対する効果を単純に分析すると交絡バイアスが生じます。交絡バイアスは分析の時点で層別分析や多変量回帰分析などで調整することが可能ですが調整の程度はさまざまです。

例えば、盲検化がされていないため“別の治療を受ける”という“治療企図からの乖離”が起きると、対象者の介入の内容を変え、治療の効果＝アウトカムを変えてしまうので、交絡因子になります。共通原因により生じるバイアスは交絡バイアスに相当します。

系統的測定誤差→情報バイアス

系統的測定誤差あるいは誤分類はさまざまな原因で起きます。図に示すUは、未知の因子ですが、それが影響して測定誤差が生じた介入がI*、アウトカムがO*です。測定誤差の原因がIにある場合も（赤い矢印）、Oにある場合も（緑の矢印）あります。

例えば、盲検化がされていないために、アウトカム測定者が患者が受けている治療を知ることができるため、介入に有利な測定結果を出してしまい、過大評価の結果が得られた場合、系統的測定誤差の赤の矢印の因果関係が作動したための情報バイアスの一例になります。このようなバイアスは介入が新しい治療法でアウトカム測定者がより高い効果を期待しているような場合に起き得ます。このようなバイアスは検出バイアスに相当します。

矢印の始点の方が原因で終点の方が結果を示します。IはIntervention介入またはEと書いてExposure要因暴露、CはCommon cause共通原因＝Confounder交絡因子、OはOutcomeアウトカムを表します。OiはIntermediate outcome中間アウトカム、OeはEndpointエンドポイントです。I*は測定された介入、O*は測定されたアウトカムを表します。Uは測定誤差を引き起こす因子で、多くの場合不明です。

DAGの各変数

これら、I, C, O, Oi, Oe, I*, O*, Uは変数を表し、例えば、I=1は介入あり、I=0は介入なし、のような値が設定され、例えばO＝１治癒、O＝０非治癒とすると、条件付き確率の式を用いてP(O=1|I=1)と記述すると、Iが介入ありの場合のOが治癒となる確率を表し、介入群の治癒確率を表すことになります。

図中ボックスで囲んであるのが実際の分析対象となる変数です。Sは選択を表す変数ですが、Sのボックスの意味は、介入と中間アウトカムに基づく対象者の除外が研究デザインあるいは分析のしかたで起きることを示しています。Sの値が０なら対照、1なら介入という設定や、0は非脱落、1は脱落、０は報告する、1は報告しないというような値を設定できます。

左下の共通原因のDAGの青い点線はI→C→Oという本来ないはずの因果関係がBack door pathバックドア経路として開かれるということを示しています。

系統的測定誤差のDAGは必ずしも矢印のすべてが同時に起きる事象ということではなく、どの矢印が有効かは、分析対象の研究によって異なってきます。

観察研究の場合は、IをE要因曝露に置き換えます。

バイアスは数多く存在し、観察研究はバイアスの影響を受けやすいですが、ランダム化比較試験もバイアスに無縁ではなく、多くのバイアスの影響を受ける可能性があります。

ランダム化比較試験のバイアス評価について、特にCochrane risk of bias tool ver. 2.0 RoB 2)を中心にスライドと解説の資料を作成しました。RoB 2を用いたランダム化比較試験のエビデンス評価の作業をする際に参考にしてください → Link

文献：
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 PubMed

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 PubMed

Lash TL, VanderWeele TJ, Haneuse S, Rothman KJ: Modern Epidemiology (FORTH EDITION). 2021,  Wolters Kluwer, PA, USA. Amazon

Luijendijk HJ, Page MJ, Burger H, Koolman X. Assessing risk of bias: a proposal for a unified framework for observational studies and randomized trials. BMC Med Res Methodol. 2020 Sep 23;20(1):237. doi: 10.1186/s12874-020-01115-7. PubMed

Westreich D: Berkson’s bias, selection bias, and missing data. Epidemiology 2012;23:159-64. doi: 10.1097/EDE.0b013e31823b6296 PMID: 22081062 PubMed

Hernán MA, Monge S: Selection bias due to conditioning on a collider. BMJ 2023;381:1135. doi: 10.1097/EDE.0000000000000031 PMID: 37286200 PubMed

バイアスに関する論文はたくさんありますが、役立ちそうな文献を少しあげておきます：
Hernán MA, Monge S: Selection bias due to conditioning on a collider. BMJ 2023;381:1135. doi: 10.1136/bmj.p1135 PMID: 37286200 PubMed

Sjölander A: Selection Bias with Outcome-dependent Sampling. Epidemiology 2023;34:186-191. doi: 10.1097/EDE.0000000000001567 PMID: 36722800 PubMed

Lu H, Cole SR, Howe CJ, Westreich D: Toward a Clearer Definition of Selection Bias When Estimating Causal Effects. Epidemiology 2022;33:699-706. doi: 10.1097/EDE.0000000000001516 PMID: 35700187 PubMed

Multi-nomial logistic regression

応答変数（従属変数）が二値変数、例えばある疾患・そうでないで、いくつかの説明変数からそれを推定するモデルを作成したい場合、ロジスティック回帰分析を用いることができます。説明変数は連続変数でも名義変数でも使えます。

応答変数が3つ以上ある場合はMulti-nomial logistic regression analysis多項ロジスティック回帰分析を用いることができます。例えば、診断したい病態が3つあり、各症例で3種類の検査のデータがあるような場合、このデータを解析して、各検査の値に対する係数を算出し、別の症例や検証用の症例のデータに対して、その症例に対する疾患確率を計算することができます。診断したい疾患が3つであれば、それぞれに対する疾患確率の値を3つ得ることができます。実臨床では、ある疾患・そうでないの2つを鑑別することよりも、3つ以上の複数の疾患のどれかを診断するような状況の方がはるかに多いと思います。前者は感度・特異度のパラダイムで処理できますが、後者の場合、感度・特異度の考え方では対処が難しくなります。

想定される疾患・病態が3つ以上あり、検査や診断法が2つ以上あるような場合、Multi-nomial logistic regressionに基づくモデルは役に立つはずです。

Rのパッケージでnnetを用いると、その解析ができます。解析に必要なデータは図1のようなものになります。

図1．Multi-nomial logistic regression用のデータの例。

このようなデータをExcelなどで用意しておきます。なお、ここに示すのは架空のデータです。

nnetはBrian Ripley, William VenablesによるRのパッケージです。 Link

以下Rでスクリプトを実行していきます。変数名は図1のものを用いていますが、実際に図1のデータを解析しているわけではありません。スクリプトと手順を示すのが目的です。それぞれ自分のデータで必要に応じて変数名など変更して実行してください。

#Install packages needed.

packneed=c(“nnet”,”tcltk2″);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.”)}

今回取り上げたパッケージ以外に、NeuralNetTools、EffectStars2というパッケージも有用です。今回Multi-nomial logistic regression analysisをnnetパッケージのmultinom()関数で行っていますが、nnetはそれだけでなく、Neural networkを用いる解析もできます。

#Read in packages.

library(nnet)
library(tcltk2)

#Read in data.

train.dat=read.delim(“clipboard”,sep=”\t”,header=TRUE)

#Excelで必要なデータ範囲をコピーしてからRに戻りこのスクリプトを実行します。クリップボードを介して、データが変数train.datに格納されます。Macの場合は、train.dat=read.delim(pipe(“pbpaste”),sep=”\t”,header=TRUE) です。

#Set a reference as a factor variable.

train.dat$Diagnosis = as.factor(train.dat$Diagnosis)
train.dat$Diagnosis = relevel(train.dat$Diagnosis, ref = “D1”)

ベースとなるDiagnosisを設定します。その際に、数値ではなくfactorとして取り扱うので、as.factor()関数を用いています。ref引数で設定するDiagnosisは通常は対照となるような診断・病態です。

#Multi-nomial logistic regression analysis.

multinom_model = multinom(Diagnosis ~ Test1+Test2+Test3, data = train.dat)

診断 ~ 説明変数1 + 説明変数2 + 説明変数3の様にformulaを記述します。

#Show the results: Intercepts, coefficients for explanatory variables.

summary(multinom_model)

#P values for explanatory variables.

z=summary(multinom_model)$coefficients/summary(multinom_model)$standard.errors
P_value=(1-pnorm(abs(z),0,1))*2
P_value

summary()関数ではP値は計算されないので、このようなスクリプトで計算しています。

#Odds ratios for explanatory variables (exponentials of the coefficients).

exp(coef(multinom_model))

#Prediction the diagnosis with explanatory variables for each case with the training data.

(predDiagnosis=predict(multinom_model, newdata=train.dat, type=”probs”))

ここでは、トレーニングに用いたデータに対して、Predictionを行っています。適合度が高いので、正診率は高くなります。

#Prepare the result data to be saved.

results=cbind(train.dat,predDiagnosis)

#Save the result data as a text file (tab-separated).

filnam=tclvalue(tkgetSaveFile(initialfile=”.txt”,filetypes=”{{Text Files} {.txt}} {{All files} *}”));if(filnam!=””){write.table(results,filnam,sep=”\t”,row.names=FALSE,col.names=TRUE)}

ファイル名を設定して任意のフォルダーにテキストファイルとして、元のデータとPredictionの結果を合わせてタブ区切りのテーブルの形式で保存します。

#Calcuating disease probabilities for individual explanatory variable data.
Set data for the case:
an example

Test1=3
Test2=14
Test3=5
case=cbind(Test1,Test2,Test3)
row.names(case)[1]=”case”

任意の検査結果の場合の、診断の推定を行います。

#Output Diagnosis class.

(CasePredicted=predict(multinom_model, newdata = case, type=”class”))

3つの診断の内、疾患確率が一番高いものが出力されます。

#Output probabilities.

(CasePredicted=predict(multinom_model, newdata = case, type=”probs”))

3つの診断のそれぞれの疾患確率が出力されます。

診断のPredictionは上記のpredict()関数を用いてできますが、自分でスクリプトを書いて計算したい場合は以下の様にします。

Multi-nomial logistic regressionの結果、Reference diagnosis以外のDiagnosisに対する切片、係数の値が得られます。これらの値を用いて、各症例の疾患確率を計算します。切片、係数の値は、coef(multinom_model)でデータフレームとして得られます。そこから、必要な切片、係数の値を読み出して、その症例のそれぞれの検査結果の値にそれぞれの対応する係数を掛け算して、総和を求め、さらに切片の値を加算します。得られた値のExponentialを求めます。それを各Diagnosisについて計算し、それらの総和を求め、さらに1を加算します。その値で、各DiagnosisのExponentialの値を割り算すると、各Diagnosisに対する疾患確率が得られます。加算する1はReference diagnosisの切片と係数を0に設定するためで、0のExponentialが1になるためです。

検査が3つ、すなわち、説明変数が3つの場合で、その症例の検査結果がx1, x2, x3でああれば、
診断iについて、yi = ai + b1i*x1 + b2i*x2 + b3i*x3
P(診断i) = exp(yi)/(1 + Σexp(yi))

通常のLogistic regressionの場合は、y = a + b1*x1 + b2*x2 + b3*x3で、
P(診断)=1/(1 + exp(-y)) ですから計算法が違います。

文献：
Ciaburro G, Venkateswaran B: Neural Networks with R: Smart models using CNN, RNN, deep learning, and artificial intelligence principles. 2017, Packt Publishing, Birmingham, UK.

Xu J: Modern Applied Regression. 2023, CRC press, Taylor & Francis Group.

Multi-nomial logistic regression with R. R-bloggers.