観察研究でCox回帰を当てはめ、年齢と重症度と併存疾患を調整したうえでハザード比0.68が出たとします。表に並べた交絡因子は入れ終えたつもりでも、治療の選択に関わった生活機能や紹介経路が測定されていなければ、この値だけで因果効果を確定できません。感度分析は、見えていない交絡や欠測や打ち切りの仮定をどこまで動かすと主解析の読み方が変わるのかを、主張ではなく数字として残す作業です。
欠測データで扱った完全ケース解析と多重代入も、感度分析としてよく用いられるペアです。ここでは未測定交絡に絞り、E値とtipping pointがそれぞれ何を答えているのかを、ランダムに生成したデータセットで確かめます。
主解析の値が動く理由
同じ患者データでも、解析対象を完全ケースに絞るか欠測値を多重代入するかでハザード比は変わります。また、治療の選択に関わった因子が測定されていなければ、交絡因子を調整したCox回帰でも推定値に偏りが残ります。主解析と別の値を並べることが感度分析なのではなく、どの仮定を動かし、その仮定が研究のどの部分に関わるのかを先に決めることが出発点です。
解析の順序も結果の読まれ方を左右します。主解析の結果を見てから感度分析を探し始めると、都合のよい設定を選びやすくなるので、解析計画の段階で完全ケース、代入の方法、打ち切りの扱い、未測定交絡の範囲まで書いておきます。
ここで混同しやすい点があります。主解析と感度分析の推定値が近いことは、その設定の範囲では結論が大きく動かなかったという情報にすぎず、未測定交絡が存在しないことの証拠にはなりません。値が離れた場合も感度分析が失敗したのではなく、結論がどの仮定に依存しているのかが見えたことになります。
Cox回帰は、打ち切りが独立していること、比例ハザード性、調整した共変量で交絡を扱いきれていることを、同じモデルの中で暗黙に使っています。情報的打ち切りで問題になるのは追跡から外れる仕組みであり、未測定交絡で問題になるのは治療の選択とイベントの両方に関わる因子です。どちらもハザード比を動かしますが、動かしている仮定は異なります。そのため、打ち切りの重み付けを行った結果を未測定交絡への感度分析と呼ぶことはできません。
E値が置く未測定交絡の強さ
VanderWeeleとDingが提案したE値は、測定されていない因子が治療とアウトカムの双方にどれだけ強く関連すれば、観測された関連をリスク比1まで動かせるかを表す数値です。保護的に見えるハザード比では逆数を使い、イベントがまれだという近似のもとで次の形になります。
$$E(\mathrm{HR}) = \mathrm{HR}^{-1} + \sqrt{\mathrm{HR}^{-1}\left(\mathrm{HR}^{-1}-1\right)}, \qquad 0 < \mathrm{HR} < 1$$
$\mathrm{HR}$は観測されたハザード比、$E(\mathrm{HR})$はその関連をリスク比1まで説明するために必要な関連の強さです。治療とアウトカムの両方に対して、これだけの強さで関連する因子が隠れていれば、観測された値は交絡だけで作れることになります。
点推定値のE値だけを書くと、必要な交絡の強さを大きく見せることになります。図1で分かるとおり、上限0.936に対応するE値は1.34で、点推定値の半分ほどです。論文には点推定値と、因果効果なしに最も近い信頼区間限界の両方のE値を並べます。限界が既に1をまたいでいる場合は、E値を1と書くのではなく、推定の不確実性が効果なしを含んでいると書きます。
E値は未測定交絡の有無を検定する値ではありません。また、ハザード比をそのままリスク比として扱う近似には、追跡終了時のイベントがまれだという条件が付きます。イベントが多い研究、比例ハザード性が崩れている研究、時間依存交絡を扱う研究では、E値だけを結論の根拠にすることは適切ではありません。
偏りの向きを取り違えない
ここからは、治療とイベントの両方に関わる二値の因子を乱数で作り、主解析からは意図的に外したデータセットで確かめます。1500人のうち157人にイベントが起き、イベント割合は10.5%です。生成時に置いた治療の効果はハザード比0.70で、この因子はイベントのハザードを2倍にし、治療を受ける確率も上げるようにしています。
主解析のハザード比は0.682（95%信頼区間0.497〜0.936）、この因子を加えた参照モデルでは0.567（0.408〜0.788）でした。因子の分布を見ると、治療群では56.0%、対照群では28.7%が保有しており、イベントへのハザード比は1.94です。危険因子が治療群に多いので、主解析のハザード比は1へ引き寄せられていました。
ここが読み違えやすいところです。E値2.29が答えているのは、0.682を1まで動かすにはどれだけ強い交絡が必要かという問いで、そのためには危険因子が対照群に多い必要があります。ところが今回のデータで実際に働いていた交絡は逆向きで、効果を消す方向ではなく弱める方向に効いていました。E値の数字が大きいことは、研究にある交絡が弱いことを保証しません。向きの違う交絡を測っているためです。
補正式で参照モデルを再現する
向きと強さを自分で置いて計算するなら、LinとPsatyとKronmalの補正式が使えます。未測定因子が二値で、治療群での保有割合を$p_t$、対照群での保有割合を$p_c$、イベントへのハザード比を$\gamma$とすると、次のようになります。
$$\mathrm{HR}_{\text{adj}} = \mathrm{HR}_{\text{obs}} \times \frac{p_c(\gamma – 1) + 1}{p_t(\gamma – 1) + 1}$$
$\mathrm{HR}_{\text{obs}}$は主解析で得られたハザード比、$\mathrm{HR}_{\text{adj}}$はその因子を調整できた場合のハザード比です。分子と分母が入れ替わる形なので、危険因子が治療群に多ければ補正後の値は小さくなり、対照群に多ければ大きくなります。
今回のデータに$p_t = 0.560$、$p_c = 0.287$、$\gamma = 1.94$を入れると、0.682は0.567になりました。参照モデルの推定値と小数第3位まで一致します。真の値が分かる場面でこの一致を確かめておくと、実データで同じ式を使う時に、式が何を計算しているのかを説明できます。
tipping pointは範囲で読む
E値が必要な強さを1つの数に要約するのに対し、tipping point解析は保有割合の差とイベントへの関連を別々に動かし、どの組み合わせでハザード比1に届くかを面として扱います。補正式を$\mathrm{HR}_{\text{adj}} = 1$について解けば、境界の$p_c$は直接求まります。
数値で並べると、イベントへのハザード比が2.0の因子なら対照群に55.9ポイント多く、2.5なら40.4ポイント、3.0なら32.6ポイント、3.5なら28.0ポイント多い必要がありました。1.5の場合は102.5ポイントとなり、保有割合が100%を超えるので到達しません。つまり、弱い因子がどれだけ偏っていても0.682は説明できず、イベントのハザードを2倍以上にする因子が、対照群に半分近く多く偏っている状態を考えて初めて結論が変わります。
この範囲は臨床的にあり得る値と照らして読みます。既に測定してある重症度指標が治療の選択に及ぼす関連より大きな仮定を置くのであれば、その仮定が成り立つ場面を本文で説明します。私なら、E値を先に出して既知の交絡因子の関連と比べ、そのあとで研究で測れなかった因子を1つ挙げ、図3の縦軸と横軸がその因子で何を意味するのかを書きます。設定値を0.1刻みで広く並べるだけでは、読者はどの行を重く見ればよいのか分かりません。
Rで確認する
データを作り、未測定因子を外した主解析と、含めた参照モデルを当てはめます。
set.seed(42)
library(survival)
n <- 1500
u <- rbinom(n, size = 1, prob = 0.45)
treatment <- rbinom(n, size = 1, prob = plogis(-0.5 + 1.3 * u))
event_time <- rexp(n, rate = 0.00008 * exp(log(0.70) * treatment + log(2) * u))
censor_time <- rexp(n, rate = 0.0008)
followup <- pmin(event_time, censor_time)
event <- as.integer(event_time <= censor_time)
main_fit <- coxph(Surv(followup, event) ~ treatment)
reference_fit <- coxph(Surv(followup, event) ~ treatment + u)#> events main reference
#> 157.000 0.682 0.567E値と補正式を関数にします。補正式は保有割合と関連の強さを引数に取ります。
e_value <- function(hr) { rr <- if (hr < 1) 1 / hr else hr; rr + sqrt(rr * (rr - 1)) }
adjust <- function(hr, p_treated, p_control, gamma) {
hr * (p_control * (gamma - 1) + 1) / (p_treated * (gamma - 1) + 1)
}
tipping <- function(hr, p_treated, gamma) {
((1 / hr) * (p_treated * (gamma - 1) + 1) - 1) / (gamma - 1) - p_treated
}
sapply(c(1.5, 2.0, 2.5, 3.0, 3.5), function(g) tipping(0.682, 0.20, g))#> 1.025 0.559 0.404 0.326 0.280
#> E値: 点推定 2.293、信頼区間上限 1.340
#> 補正式 adjust(0.682, 0.560, 0.287, 1.942) = 0.567tippingの返り値は保有割合の差なので、1.025のように1を超えた値は到達しないことを表します。読み違えやすいのはE値と補正式の関係で、E値は治療との関連とアウトカムとの関連を同じ強さに揃えた場合の最小値です。片方だけが強い組み合わせは図3の曲線上に無数にあり、E値はそのうちの1点にすぎません。
論文にどう書くか
STROBEの説明文書は、対象集団の選び方、変数定義、欠測データ、交絡因子、測定誤差などを感度分析で扱えるとしたうえで、行った感度分析を記載するよう求めています。主解析、動かした仮定、設定した範囲、各設定でのハザード比と95%信頼区間を並べ、事後的に追加した分析があればその理由も分けて書きます。
未測定交絡ではE値とtipping pointのどちらか一方だけで済ませず、研究で想定できる因子に応じて選びます。単一の要約値が必要ならE値、交絡因子の頻度とイベントへの関連を具体的に議論するならtipping pointが読みやすい方法です。傾向スコアマッチングで共変量のバランスを確認したあとも、測定されていない因子は別に扱います。
結果が1をまたぐ設定だけを抜き出して書くより、あらかじめ決めた範囲を表か図で全て載せる方が読み手に伝わります。そして、偏りの向きも明記します。向きを書かずに数字だけ並べると、今回のように帰無方向へ働いていた交絡と、効果を作り出す方向の交絡が、同じ「未測定交絡あり」として読まれてしまいます。
ここまで読んで、自分のデータでどの仮定を動かせばよいのか迷われた方もいるかもしれません。実際には、測れなかった因子の候補、その因子の想定される頻度、主解析の推定精度まで含めて決めることになります。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
次に確かめること
今回の補正式は未測定因子が1つで二値、そして比例ハザード性が成り立つ場面を前提にしています。因子が連続量だったり複数あったりする場合は、同じ式をそのまま当てはめられません。また、何を因果効果として推定したいのかが決まっていなければ、感度分析で動かす仮定も定まらないため、因果推論で扱った推定対象の設定を先に確かめておくと迷いが減ります。
参考文献
1. VanderWeele TJ, Ding P. Sensitivity analysis in observational research: introducing the E-value. Annals of Internal Medicine. 2017;167(4):268-274.
2. Lin DY, Psaty BM, Kronmal RA. Assessing the sensitivity of survival analyses to unmeasured confounders. Biometrics. 1998;54(3):948-963.
3. Ding P, VanderWeele TJ. Sensitivity analysis without assumptions. Epidemiology. 2016;27(3):368-377.
4. von Elm E, Altman DG, Egger M, Pocock SJ, Gotzsche PC, Vandenbroucke JP. The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. PLoS Medicine. 2007;4(10):e296.
5. VanderWeele TJ, Mathur MB. Commentary: developing best-practice guidelines for the reporting of E-values. International Journal of Epidemiology. 2020;49(5):1495-1497.
