観察データで治療効果を推定するとき、重症な患者ほどその治療を受けやすいという構造があります。単純に治療群と非治療群の生存を比較すると、実際には有効な治療であっても、有害に見えることさえあります。
このような交絡への対処として、回帰調整や傾向スコアを用いた方法が使われます。どちらも交絡を調整するための手法ですが、得られるハザード比は必ずしも一致しません。これは、そもそも、推定している効果の量が異なるためであり、どちらか一方が間違っているからではありません。
2000例の架空データを使って、この違いを確認します。真の因果効果は、条件付きハザード比0.70となるように設定しています。
3つの数字が出る
coxph(Surv(tm, st) ~ trt, data = d) # 調整なし
coxph(Surv(tm, st) ~ trt + sev + age, data = d) # 回帰調整
ps <- predict(glm(trt ~ sev + age, family = binomial, data = d), type = "response")
w <- ifelse(d$trt == 1, mean(d$trt)/ps, (1-mean(d$trt))/(1-ps))
coxph(Surv(tm, st) ~ trt, data = d, weights = w, robust = TRUE)
調整なしで1.66となるのは、交絡の影響です。重症な患者に治療が集中しているため、治療群の予後が悪く見えています。ここまでは想定どおりです。
問題は、回帰調整の0.71とIPTWの0.76がなぜ異なるのかという点です。単なる推定誤差として片づけることもできますが、両者が異なる理由には、統計的にはっきりした意味があります。
条件付き効果と周辺効果
回帰調整で得られるのは、条件付きハザード比です。「重症度と年齢が同じ患者同士を比較したとき、治療によってハザードがどの程度変わるか」という効果を表します。Cox比例ハザードモデルで扱ってきた調整済みハザード比が、この考え方に対応します。
一方、IPTWで得られるのは、周辺的な効果を表すハザード比です。「この集団全体に治療を行った場合と、集団全体に治療を行わなかった場合で、ハザードがどう変わるか」という効果です。個々の患者の重症度を直接揃えるのではなく、重み付けによって治療群と非治療群の共変量分布を揃え、集団全体として比較できる状態を作ります。
ハザード比には非崩壊性という性質があるため、条件付きハザード比と周辺ハザード比は、交絡を完全に調整できたとしても一致するとは限りません。同じデータを無作為化した状況で計算すると、条件付き効果は0.694、周辺効果は0.762となりました。図1の2本の破線がこれに対応します。予後因子による異質性が大きいほど、条件付きハザード比と周辺ハザード比の差も大きくなります。
したがって、論文で報告されたハザード比を比較するときには、まず「何を推定しているハザード比なのか」を確認する必要があります。回帰調整を使った論文とIPTWを使った論文で数字が違っていても、それだけで結果が矛盾しているとは言えません。
傾向スコアで何が起きているか
傾向スコアは、観測された共変量から治療を受ける確率を推定したものです。
$$e(\mathbf{x}) = P(A = 1 \mid \mathbf{X} = \mathbf{x}) \tag{1}$$
この傾向スコアを使って重みを作ると、もともと治療を受けにくかったのに治療を受けた患者や、治療を受けやすかったのに受けなかった患者を重く数えることになります。結果として、治療群と非治療群の共変量分布が揃った仮想的な集団を作ることができます。
傾向スコアを使った解析では、共変量のバランスと傾向スコアの重なりを確認することが重要です。バランスが取れていなければ、測定した交絡因子が十分に調整されていません。重なりがなければ、そもそも両群を比較できる患者が存在しないことになります。
重なりの確認は、正値性の仮定を確認することにもつながります。傾向スコアが0や1に近い患者が多いと、重みが極端に大きくなり、一部の患者が解析結果を大きく左右する可能性があります。
極端な重みが残る場合には、上位1%などで重みを切り詰める方法や、傾向スコアが両群で重なっている範囲に解析対象を限定する方法があります。ただし、後者では解析対象となる集団そのものが変わります。したがって、「誰についての効果を推定しているのか」も変わることに注意が必要です。
標準誤差にも注意が必要です。重み付けした解析では通常の分散推定では不十分になることがあるため、ロバスト分散を使います。Rのcoxphではrobust=TRUEを指定します。
マッチングという別の使い方
傾向スコアは、重み付け以外にも利用できます。たとえば、傾向スコアが近い患者同士を組にするマッチング、傾向スコアの範囲で層に分ける層別化、傾向スコアを共変量として回帰モデルに入れる方法があります。
マッチングは直感的で、解析の考え方を説明しやすい方法です。ただし、適切なマッチ相手が見つからなかった患者は解析から外れるため、対象集団が変わります。治療群の患者を対象にマッチングを行えば、「治療を受けた患者における効果」を推定する形になり、集団全体を対象としたIPTWの推定量とは異なります。したがって、マッチングによって除外された患者の数と特徴を報告する必要があります。
マッチング後の解析では、マッチした患者同士が独立ではない点にも注意が必要です。クラスターを持つデータと同様に、ロバスト分散や層別Cox回帰を使って、この依存性を考慮します。この処理が省略されている論文もあります。
調整しても消えないもの
ここまでの方法は、観測された交絡因子によって交絡を十分に調整できるという仮定に立っています。しかし、未測定交絡が存在しないことをデータから検証することはできません。傾向スコアによって測定した共変量のバランスがどれだけ改善しても、それはあくまでモデルに入れた変数についての話です。
図2の右で標準化差が0.01まで低下しているのは、モデルに入れた重症度と年齢についてです。カルテに記録されていない全身状態や患者の希望が治療選択に影響しているのであれば、それらによる交絡は残ります。
したがって、傾向スコアを用いた解析を無作為化試験の代わりとして扱うのは適切ではありません。原稿を確認するときには、未測定交絡についてどのように考察しているかを確認する必要があります。E値のような感度分析を用いて、観測された関連を説明するためにどの程度の未測定交絡が必要なのかを示すと、議論を具体化できます。
さらに、時間依存の交絡が存在する場合には、ここまでの方法だけでは対応できません。治療が途中で変更され、その治療選択が過去の検査値に依存し、その検査値自体が過去の治療の影響を受けているような構造です。この場合、検査値を通常の共変量として調整すると過調整になる可能性があり、調整しなければ交絡が残ります。時変共変量で触れた問題であり、周辺構造モデルなどの方法が必要になります。
未測定交絡をE値で評価する
未測定交絡が存在するかどうかは検証できませんが、どの程度の未測定交絡があれば観察された結果を説明できるのかを評価することはできます。E値は、そのための指標の一つです。
観測されたハザード比が0.76だったとき、これを1まで押し戻すために必要な未測定交絡の強さは、リスク比のスケールで約1.95になります。これは、未測定の交絡因子が治療の有無と転帰の両方に、それぞれ少なくとも約1.95倍の関連を持ち、かつその関連が測定済みの共変量では説明されない程度に残っている必要がある、という意味です。
この数字を測定済みの共変量による効果と並べると、未測定交絡の影響をより具体的に考察できます。この例では、重症度のハザード比が2.5前後でした。重症度と同程度に強い未測定の予後因子が残っていれば結論が変わる可能性がありますが、そのような強い因子がカルテに記録されずに残っているかどうかは、臨床的な知識も踏まえて判断する必要があります。
E値はEValueパッケージなどを使って計算できます。信頼区間の境界についてもE値を計算して報告すると、推定値だけでなく統計的不確実性も含めて評価できます。
正値性が崩れているとき
傾向スコアが0や1に近い患者では、その患者に対応する比較対象が存在しない可能性があります。たとえば、90歳で重度の腎障害がある患者が、その治療を受けたことが1例もなければ、その患者と同じような集団における治療効果をデータから直接推定することはできません。
このような状態を正値性の違反と呼びます。特に、単なる偶然ではなく、臨床的な禁忌などによって特定の患者には一方の治療しか行われない場合には、構造的な正値性違反が生じます。重みを切り詰めても、この問題そのものが解決するわけではありません。
対処の一つは、両群に十分な患者が存在する傾向スコアの範囲に解析対象を限定することです。ただし、その場合には対象集団が変わるため、結論をどの集団に適用できるのかを明記する必要があります。
実際の解析では、まず傾向スコアの分布を治療群と非治療群に分けて確認します。図2の左のように十分な重なりがあれば解析を進めやすくなりますが、分布が大きく偏っている場合には、重みの計算以前に、このデータで因果効果を推定できるのかを検討する必要があります。
報告の仕方
まず、推定した効果の種類を明示します。条件付きハザード比なのか、周辺ハザード比なのかを区別し、英語で記載する場合もconditionalとmarginalを使い分けます。図1のように両者は異なる値になり得るため、推定量の定義を示さずに数字だけを並べても、読者は適切に解釈できません。
次に、傾向スコアモデルに投入した変数を全て記載し、それらを選択した理由も説明します。治療と転帰の双方に関係する交絡因子を中心に設定し、治療開始後に生じる変数を安易に調整変数へ入れないことが重要です。
共変量のバランスも必ず評価します。標準化差の表や図を使って、重み付け前後を比較します。バランスの評価にp値を使うことは一般には勧められません。p値は症例数に強く依存するため、大規模なデータでは臨床的には小さな差でも統計的に有意になるからです。
重みの分布も報告します。最大値や上位パーセンタイル、重みを切り詰めた場合にはその閾値を示します。図3のような図があれば、極端な重みを持つ少数例が推定結果を支配していないかを読者が確認できます。
そして、未測定交絡についての議論を記載します。どのような因子が測定されていないのか、それが治療選択と転帰のどちらにどの方向の影響を与え得るのかを考察します。可能であれば感度分析も併せて示します。
次に確かめること
調整済みハザード比の解釈についてはCox比例ハザードモデル、ハザード比と非崩壊性の関係についてはハザード関数で扱っています。ハザード比以外の効果指標を使うことで、この問題を避けられる場合もあります。たとえば、RMSTの差は崩壊性を持つため、条件付き効果と周辺効果の違いを考えるうえで別の選択肢になります。
時間依存の交絡を扱う周辺構造モデルや、観察データを仮想的な臨床試験として設計し直すTarget Trial Emulationについては、それぞれ別の記事で扱います。
どの推定量を主解析に置くか、傾向スコアにどの変数を入れるか、未測定交絡をどこまで考慮するかは、研究の目的とデータの構造を見て決める必要があります。Dr.データサイエンスでは、こうした因果推論の解析設計についてのご相談を承っています。
参考文献
Hernan MA, Robins JM. Causal Inference: What If. Chapman and Hall/CRC; 2020.
Austin PC. An introduction to propensity score methods for reducing the effects of confounding in observational studies. Multivariate Behavioral Research. 2011;46(3):399-424.
Hernan MA. The hazards of hazard ratios. Epidemiology. 2010;21(1):13-15.
Robins JM, Hernan MA, Brumback B. Marginal structural models and causal inference in epidemiology. Epidemiology. 2000;11(5):550-560.
VanderWeele TJ, Ding P. Sensitivity analysis in observational research: introducing the E-value. Annals of Internal Medicine. 2017;167(4):268-274.
