観察研究で治療とアウトカムの関連を推定する時、測れた共変量を全て調整しても、喫煙量や全身状態のように診療記録へ十分な精度で残らない要因が治療選択とアウトカムの両方に効いていれば、推定値は真の値から離れたままになります。考察へ「未測定交絡の影響は否定できない」と書いて終える論文は多いのですが、その一文はバイアスの向きも大きさも読者へ伝えないため、定量的バイアス分析と陰性対照は、その空白に数値を入れる際に用いられます。
未測定交絡の想定を1組だけ置く単純バイアス分析、想定に幅を持たせる確率的バイアス分析、必要な交絡の強さを1つの数へまとめるE値、そして治療効果がないはずのアウトカムを使う陰性対照は、同じ未測定交絡を扱いながら返す量が異なり、その判断の一端を本記事が担えると幸いです。
なお、$A$は治療の有無、$Y$は主要アウトカム、$U$は測っていない交絡因子で、どれも1か0を取る変数です。観察されたリスク比を$\mathrm{RR}^{\mathrm{obs}}$、$U$の影響を除いた後のリスク比を$\mathrm{RR}^{\mathrm{adj}}$、$U$がアウトカムに与えるリスク比を$\mathrm{RR}_{UY}$において、調整によって推定量そのものが変わる話は生存時間解析における因果推論で扱っていますのでご確認ください。
観察された1.804が含んでいるもの
乱数で生成した20,000例のデータセットを使って話を進めます。治療群9,999例、対照群10,001例で、未測定交絡因子$U$は治療群の50%、対照群の20%が持ち、$U$を持つ人はアウトカムが3倍起きやすく、治療そのものの効果はリスク比1.30と置いてあります。
この設定で観察された主要アウトカムは、治療群13.3%（1,329例）、対照群7.4%（737例）でした。リスク比は1.804、95%信頼区間は1.656〜1.965となり、$U$が見えていない解析者の手元には、真の1.30よりはるかに強い有害方向の関連だけが残ります。
もっとも、信頼区間の幅が狭い点には注意が必要です。標本誤差は20,000例で十分に小さくなっていますが、交絡による偏りは症例数を増やしても縮まないため、区間の狭さは結論の確からしさを保証しません。
単純バイアス分析：想定を1組置いて補正する
定量的バイアス分析の最も簡単な形は、バイアスパラメータを1組だけ決めて観察値を割り戻す方法です。Greenland（1996）が整理した補正式では、$U$の保有割合と$U$のアウトカムへの影響から、群間のリスクが何倍ズレるかを表すバイアス因子を作ります。
$$\mathrm{BF} = \frac{\mathrm{RR}_{UY}\Pr(U=1 \mid A=1) + \Pr(U=0 \mid A=1)}{\mathrm{RR}_{UY}\Pr(U=1 \mid A=0) + \Pr(U=0 \mid A=0)} \tag{1}$$
式(1)は、$U$の分布の違いだけで治療群と対照群のリスクが何倍ズレるかを表しています。分子は治療群、分母は対照群について、$U$を持つ人と持たない人のリスクを保有割合で加重した値です。
$$\mathrm{RR}^{\mathrm{adj}} = \frac{\mathrm{RR}^{\mathrm{obs}}}{\mathrm{BF}} \tag{2}$$
式(2)は、観察されたリスク比のうち$U$の分布差で説明できる分をバイアス因子として取り除く操作です。ここで補正されるのは偏りだけで、標本誤差は手つかずのまま残ります。
実際に$\mathrm{RR}_{UY}=3.0$、保有割合0.50対0.20を入れるとバイアス因子は1.429となり、補正リスク比は1.804を1.429で割った1.263になりました。設計で置いた真の値1.30に近い値です。
この補正が正しく効いていることは、別の計算経路でも確かめられます。$U$が見えている前提で$U$の周辺分布へ標準化したリスク比は1.304で、補正式から得た1.263との差は0.042にとどまりました。
もっとも、実務で$\mathrm{RR}_{UY}$と保有割合を正確に知っている状況はまずありません。単純バイアス分析の価値は正解を当てることではなく、どの想定を置けば結論が変わるのかを読者へ開示する点にあります。
想定に幅を持たせる確率的バイアス分析
そこで、バイアスパラメータを一点で決め打ちせず、分布から繰り返し引いて補正を反復する方法が使われます。Lash et al.（2014）が確率的バイアス分析と呼んでいるのがこの手順で、バイアスパラメータの不確かさと標本誤差を同じ枠内で扱えます。
この手順に沿って、今回は$\mathrm{RR}_{UY}$を2.0〜4.0、治療群の保有割合を0.40〜0.60、対照群を0.12〜0.28の一様分布から20,000回引き、あわせて観察リスク比の標本誤差も対数スケールの正規乱数で加えました。補正後のリスク比の中央値は1.283、シミュレーション区間は0.990〜1.562で、引いた結果の96.9%が1を超えました。
なお、この区間を信頼区間のように読むのは避けた方がよいでしょう。0.990という下限は「データが示した下限」ではなく、置いた一様分布の端に引きずられた値であり、範囲を0.5だけ広げれば下限もそのぶん動きます。
したがって、確率的バイアス分析を報告する時は、区間の数値より先に、どの範囲をどの分布で置いたかとその根拠を書く必要があります。Lash et al.（2014）が繰り返し求めているのもこの点で、バイアスパラメータの出所を示さない感度分析は、精密に見えるだけで検証できません。
E値は想定の代わりに下限を返す
想定を置くことをためらう場面では、逆向きの問い方があります。VanderWeele & Ding（2017）のE値は、観察された関連を説明し切るために未測定交絡因子が治療とアウトカムの双方と持つ必要のある最小の関連の強さを返します。
$$E = \mathrm{RR}^{\mathrm{obs}} + \sqrt{\mathrm{RR}^{\mathrm{obs}}\left(\mathrm{RR}^{\mathrm{obs}}-1\right)} \tag{3}$$
式(3)は、観察されたリスク比だけからE値が決まることを表しています。バイアスパラメータを一つも置かずに計算できる代わりに、返ってくるのは補正後の効果量ではなく、必要な交絡の強さの下限です。
実際に、今回の1.804に対するE値は3.008、信頼区間下限1.656に対しては2.697でした。未測定交絡因子が治療の有無とアウトカムの両方に3倍以上の関連を持っていなければ、1.804という関連は交絡だけでは説明し切れない、という読み方になります。
ここで読み違えやすいのは、E値が大きいほど結論が頑健だと単純に判断してしまうことです。E値は補正後の推定値を返さないため、3.008という数が臨床的に大きいかどうかは、その領域で実際に観察される交絡の強さと照らして初めて決まり、この区別はE値とtipping pointを並べた生存時間解析における感度分析でも同じ形で出てきます。
陰性対照アウトカムが残りの交絡を見せる
ここまでは想定を置く論理でしたが、データそのものに交絡の痕跡を探す方法もあります。Lipsitch et al.（2010）は、治療の影響を受けないと考えられるアウトカムや曝露を陰性対照と呼び、そこで関連が出たら交絡か選択の偏りが残っている証拠になると論じました。
この考え方に沿えば、治療効果がゼロのはずのアウトカムでリスク比が1から離れた時、そのズレは交絡が作ったものになります。今回のデータでは、治療の影響を受けない設計にした陰性対照アウトカムのリスクが治療群14.1%（1,412例）、対照群10.5%（1,048例）で、リスク比は1.348、95%信頼区間1.250〜1.452でした。
しかも、真の効果がゼロと分かっている量で1.348が出ている以上、主要アウトカムの1.804をそのまま治療の効果として読むことはできません。
さらに、このズレは想定の検算にも使えます。陰性対照で観察された1.348を式(1)に当てはめて$\mathrm{RR}_{UY}$について解くと2.508となり、設計で陰性対照アウトカムに置いた2.50とほぼ一致しました。
陰性対照を補正へ使う時の条件
陰性対照の使い方には、交絡を検出する使い方と、ズレの大きさを推定値の補正へ使う使い方があります。Shi et al.（2020）が整理しているように、後者は近接因子を用いた枠組みとして発展していますが、成り立つために必要な条件は検出より格段に厳しくなります。
この補正へ進むには、陰性対照が治療から因果的な影響を一切受けないことに加えて、主要アウトカムと同じ交絡構造を共有していることまで仮定する必要があります。診療記録から後付けで選んだ陰性対照でこの仮定を主張するのは、実務上かなり難しくなります。
そのため、補正まで踏み込む使い方を主解析に据えることは適切ではありません。陰性対照アウトカムの選定が事前に登録されていて、かつ交絡構造の仮定を臨床的な根拠とともに書ける場合を除き、検出のツールとして使い、補正は感度分析の一部に留める方が査読でも説明が通ります。
逆に、検出としての使い方には無理な仮定がほとんど要りません。Arnold & Ercumen（2016）は無作為化試験でも陰性対照アウトカムが役立つ場面を挙げていて、割付が保たれていれば1付近に収まるはずの値がズレていれば、運用上の問題を早い段階で拾えます。
論文でどこまで書くか
定量的バイアス分析は、考察の一文を数値へ置き換えるために書きます。Lash et al.（2014）は、報告の際にバイアスパラメータの値とその出所、置いた分布の形、補正前後の推定値を並べて示すよう求めています。
この求めに沿って今回の例を並べると、観察されたリスク比1.804、想定を1組置いた補正値1.263、想定に幅を持たせた時の中央値1.283と区間0.990〜1.562、E値3.008、陰性対照アウトカムのリスク比1.348が並びます。数値の並びだけでも、結論が想定にどれだけ依存しているかが読み手に伝わります。
この並びをどう作るかについて、私なら、まず陰性対照アウトカムを見ます。想定を置く前にこのデータで交絡がどちらへどれだけ効いているかの手がかりが得られ、そこで得た値は確率的バイアス分析の範囲設定の根拠として論文に書けます。
なお、傾向スコアで測れた共変量を揃えた後も、この記事の数値は同じ形で必要になります。傾向スコアマッチングも周辺構造モデルも、測れた変数の情報しか使えません。
Rで確かめる
ここまでの数値は、baseRとstatsだけで再現できます。乱数シードを42に固定し、未測定交絡因子を持つ20,000例と、治療の効果を入れない陰性対照アウトカムを生成しました。
set.seed(42)
n_total <- 20000
treatment <- rbinom(n_total, size = 1, prob = 0.5)
# 未測定交絡因子Uは治療群に50%、対照群に20%
prev_u <- ifelse(treatment == 1, 0.50, 0.20)
confounder_u <- rbinom(n_total, size = 1, prob = prev_u)
risk_outcome <- 0.05 * 1.30^treatment * 3.00^confounder_u
outcome <- rbinom(n_total, size = 1, prob = risk_outcome)
# 陰性対照アウトカムのリスクにtreatmentは入れない
risk_nc <- 0.08 * 2.50^confounder_u
negative_control <- rbinom(n_total, size = 1, prob = risk_nc)次に、補正式とE値の定義をそのまま関数にして、観察されたリスク比へ当てはめました。
obs_rr <- mean(outcome[treatment == 1]) / mean(outcome[treatment == 0])
se_log <- sqrt(1 / sum(outcome[treatment == 1]) - 1 / sum(treatment == 1) +
1 / sum(outcome[treatment == 0]) - 1 / sum(treatment == 0))
bias_factor <- function(rr_uy, p_u1, p_u0) {
(rr_uy * p_u1 + (1 - p_u1)) / (rr_uy * p_u0 + (1 - p_u0))
}
bf_correct <- bias_factor(3.00, 0.50, 0.20)
e_value <- function(rr) rr + sqrt(rr * (rr - 1))
round(c(obs_rr, bf_correct, obs_rr / bf_correct,
e_value(obs_rr), e_value(exp(log(obs_rr) - 1.96 * se_log))), 3)#> [1] 1.804 1.429 1.263 3.008 2.697この出力の読み方として押さえておきたいのは、補正リスク比1.263とE値3.008が違う単位の量だという点です。前者はリスク比のスケールに乗った補正後の効果量、後者は交絡因子との関連の強さであり、両者を並べて大小を比べる意味はありません。
続いて、バイアスパラメータを一様分布から引き、補正を20,000回繰り返しました。標本誤差も同じ枠で扱うため、観察リスク比の対数へ正規乱数を加えています。
n_draws <- 20000
draw_rr_u <- runif(n_draws, min = 2.0, max = 4.0)
draw_p_u1 <- runif(n_draws, min = 0.40, max = 0.60)
draw_p_u0 <- runif(n_draws, min = 0.12, max = 0.28)
draw_obs <- exp(rnorm(n_draws, mean = log(obs_rr), sd = se_log))
draw_adj <- draw_obs / bias_factor(draw_rr_u, draw_p_u1, draw_p_u0)
hist(draw_adj, breaks = 60, col = "grey85", border = "white",
main = "Probabilistic bias analysis of an unmeasured confounder",
xlab = "Bias-adjusted risk ratio", ylab = "Number of draws")
round(c(median(draw_adj), quantile(draw_adj, c(0.025, 0.975)),
mean(draw_adj > 1)), 3)#> [1] 1.283 0.990 1.562 0.969
また、図1から読み取れるのは、補正後の値が1点に定まらないことです。裾がどこまで伸びるかは引いたパラメータの範囲で決まるため、範囲を変えた時に結論がどう動くかも併せて見ておきます。
一方、図2で確かめられるのは、交絡の強さだけでは補正値が決まらないことです。$U$の保有割合が治療群と対照群で近いほど曲線は寝て、$\mathrm{RR}_{UY}$をいくら大きくしても観察値はほとんど動きません。
自分の研究で何を判断するか
感度分析の数値を読む時に最初に決めるのは、その研究の結論が定性的なものか定量的なものかです。治療に害があるかどうかだけを扱う研究なら、補正後も96.9%の引きが1を超えたという情報は判断材料になりますが、効果量を報告したい研究では、1.263と1.562のどちらを書くかという問題が残ります。
この選択に正解はありませんが、決め方はあります。事前に定めた解析計画で主たる補正の想定を1組決めておき、確率的バイアス分析はその想定の周りを見る補助として位置づければ、結果を見てから都合のよい数値を選んだ、という疑いは残りません。
その判断を実際のデータで詰める段になると、未測定交絡の想定範囲を文献からどう持ってくるか、陰性対照として何を選べるかが具体的な検討事項になります。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
- Lash TL, Fox MP, MacLehose RF, Maldonado G, McCandless LC, Greenland S. Good practices for quantitative bias analysis. International Journal of Epidemiology. 2014; 43(6): 1969-1985.
- Greenland S. Basic Methods for Sensitivity Analysis of Biases. International Journal of Epidemiology. 1996; 25(6): 1107-1116.
- Lipsitch M, Tchetgen Tchetgen E, Cohen T. Negative controls: a tool for detecting confounding and bias in observational studies. Epidemiology. 2010; 21(3): 383-388.
- Arnold BF, Ercumen A. Negative Control Outcomes: A Tool to Detect Bias in Randomized Trials. JAMA. 2016; 316(24): 2597-2598.
- VanderWeele TJ, Ding P. Sensitivity Analysis in Observational Research: Introducing the E-Value. Annals of Internal Medicine. 2017; 167(4): 268-274.
- Shi X, Miao W, Tchetgen Tchetgen E. A Selective Review of Negative Control Methods in Epidemiology. Current Epidemiology Reports. 2020; 7(4): 190-202.
