25.3K Views
September 08, 22
スライド概要
統計数理研究所 医療健康データ科学研究センター 公開講座「疫学・公衆衛生統計」|2021年12月9日
Assoc prof at Tokyo University of Science. PhD in health sciences/MPH at the University of Tokyo. Causal inference in epidemiology/biostatistics.
統計数理研究所 医療健康データ科学研究センター 公開講座「疫学・公衆衛生統計」 因果推論 I 東京理科大学 工学部情報工学科 篠崎 智大 [email protected] 2021年12月9日 15:00~16:30
「因果推論」のアウトライン 因果推論 I (12月9日):因果推論の基礎 効果、交絡、交絡変数、交換可能性、層別解析、回帰、傾向スコア 因果推論 II (12月10日):モデルを用いた効果推定 回帰モデル(条件付き効果、周辺効果、回帰標準化) 傾向スコアモデル(回帰調整、マッチング、逆確率重み付け IPW) 二重ロバスト推定(augmented IPW / IPW回帰標準化) 2
因果効果? ある人(Aさん)が、ある曝露(喫煙)を受けることで、あるイベント(10年後 に死亡)を起こすかどうかを観察した場合 データ(事実) Aさんは喫煙し続けていたら10年後に亡くなった Aさんにとって喫煙は10年後の死亡に効果があった? 3
反事実アウトカムモデル Counterfactual-outcome model 事実データ(factual data) Aさんが2010年1月から喫煙をやめることなく(XA = 1)10年後に死亡(YA = 1) 反事実データ(counterfactual data) Aさんが2010年1月から禁煙を開始したら(XA = 0)10年後に…?(YA = ? if XA = 0) If XA had been 0, what YA would be? 4
潜在アウトカム Potential outcomes Yx:曝露 X = x を受けた場合のアウトカム Aが曝露を受けた場合(XA = 1)のアウトカムYA YAx=1 Aが曝露を受けなかった場合(XA = 0)のアウトカムYA YAx=0 5
潜在アウトカムによる因果効果 Aが曝露を受けた場合(XA = 1)のアウトカムYA YAx=1 Aが曝露を受けなかった場合(XA = 0)のアウトカムYA YAx=0 Aさんに対する因果効果:YAx=1 – YAx=0 データとして得られるのは 曝露ありだった場合(XA = 1):YA = YAx=1 曝露なしだった場合(XA = 0):YA = YAx=0 必ず一方は反事実 6
潜在アウトカムと観察アウトカム Yx=0 Y Causal effect (Yx=1 – Yx=0) 対象者 X Yx=1 A 1 1 1 1 0 B 1 1 0 1 1 C 0 0 0 0 0 D 0 1 0 0 1 E 1 0 1 0 –1 F 0 1 1 1 0 G 0 1 0 0 1 7
観察変数と反事実変数 Yx=1, Yx=0 どちらか一方は実際に観察できない(反事実) 個人 i に対する因果効果 Yix=1 – Yix=0 は調べられない 集団に対する効果 平均因果効果(average causal effect) E(Yx=1) – E(Yx=0) E(Yx=1):集団全員が曝露を受けた場合の期待値 E(Yx=0):集団全員が曝露を受けたなかった場合の期待値 8
平均因果効果 E(Yx=1) – E(Yx=0) 対象者 X Yx=1 Yx=0 Y Causal effect (Yx=1 – Yx=0) A 1 1 1 1 0 B 1 1 0 1 1 C 0 0 0 0 0 D 0 1 0 0 1 E 1 0 1 0 –1 F 0 1 1 1 0 G 0 1 0 0 1 E(Yx=1) – E(Yx=0) = 5/7 – 3/7 = 2/7 でも実際に計算できるのは… 9
Association is not causation Observed Population Unexposed Causation Exposed Association Hernán. JECH 2004 Hernán & Robins, 2020 10
ランダム化 一方の曝露群に予後の悪い人が集まらない Yx (x = 0, 1)と曝露 X が独立 Yx=1, Yx=0 E(Yx=1| X = 1) = E(Yx=1| X = 0) = E(Yx=1) E(Yx=0| X = 1) = E(Yx=0| X = 0) = E(Yx=0) X で決まる2つの集団は交換可能(exchangeable) X = 1 の集団 X = 0 の集団 がおなじ潜在アウトカムの期待値をもつ 11
ランダム化試験では 平均因果効果:E(Yx=1) – E(Yx=0) E(Yx=1) ランダム化されているので、 E(Yx=1) = E(Yx=1 | X = 1) データとして得られるから、 E(Yx=1 | X = 1) = E(Y | X = 1) E(Yx=0) ランダム化されているので、 E(Yx=0) = E(Yx=0 | X = 0) データとして得られるから、 E(Yx=0 | X = 0) = E(Y | X = 0) 平均因果効果 = E(Y | X = 1) – E(Y | X = 0) Association is causation! 12
ランダム化の仮定 ランダム化は実際にできていないが… 曝露 X が Yx とは独立に選ばれているとみなす ただし Yx に関係しない変数にしたがって選ばれていてもよい 交換可能性の仮定とも exchangeability assumption 13
観察研究 ランダム化の仮定は正しくないと分かっている Association is NOT causation 交絡(confounding) 曝露群:集団全員が曝露した場合と異なる E(Y | X = 1) ≠ E(Yx=1) ⇔ E(Yx=1 | X = 1) ≠ E(Yx=1 | X = 0) and/or 非曝露群:集団全員が曝露しなかった場合と異なる E(Y | X = 0) ≠ E(Yx=0) ⇔ E(Yx=0 | X = 1) ≠ E(Yx=0 | X = 0) 14
「条件付き」ランダム化の仮定 リスク因子 C で層別した中で E(Yx=1 | X = 1, C = c) = E(Yx=1| C = c) E(Yx=0 | X = 0, C = c) = E(Yx=0| C = c) 交絡変数(confounders) その層内では「ランダム化」されていると見なせるような変数 交絡変数が全て測定されている仮定ともいう No-unmeasured confounders 交絡変数 C で層別すれば Yx(x = 0, 1)と曝露 X が独立 15
平均因果効果 E(Yx=1) – E(Yx=0) E(Yx=1) = ∑c E(Yx=1|C = c)P(C = c) Cで層別 E(Yx=1|C = c) = E(Yx=1|X = 1, C = c) Cの中では曝露はランダムに決まると仮定 E(Yx=1|X = 1, C = c) = E(Y|X = 1, C = c) データとして観察 ∑c {E(Y|X = 1, C = c)}P(C = c) – ∑c {E(Y|X = 0, C = c)}P(C = c) 層ごとの平均(= 回帰、条件付き期待値)の重み付け平均 層別解析 重みを P(C = c) としたとき、標準化(standardization)という 16
因果効果の識別可能条件 一般的な警句: Association is not causation ランダム化の下では 単純平均の差が因果効果 No-confounders ランダム化と見なせない(交絡がある)状況では 交絡変数が全て測定されているなら、単純平均の差の重み付け平均が因果効果 No-unmeasured confounders 測定されていない交絡変数があると因果効果は求まらない 17
ランダム化すれば交絡がない? 交絡のあるデータ例 : E(Y|X = x) ≠ E(Yx) 対象者 X Yx=1 Yx=0 Y Causal effect (Yx=1 – Yx=0) A 1 1 1 1 0 B 1 1 0 1 1 C 0 0 0 0 0 D 0 1 0 0 1 E 1 0 1 0 –1 F 0 1 1 1 0 G 0 1 0 0 1 どんな X の割り付けをしても E(Y|X = x) = E(Yx) にならなくない…? 18
統計学 meets 因果推論 ランダム化を行っても、個々の有限集団では交絡あり ランダム化している場合 ランダム割り付けの仮想的な繰り返しを標本空間と考えると E(Yx|X = 1, 割り付けパターン) = µ1x これらのパラメータが確率変数 E(Yx|X = 0, 割り付けパターン) = µ0x パラメータが交換可能 exchangeable ベイズ統計と同じ用法 パラメータ (µ1x, µ0x) の確率分布 F(m, m*) が引数に対して対称 個々の標本で µ1x = µ0x を要求するものではない 交絡という用語は実現した割り付けに対して使いたい Greenland & Robins, IJE 1986 Greenland, Robins & Pearl, Stat Sci 1999 19
現実に得られるデータ アウトカム Y=1 Y=0 2 1 1 3 曝露 X=1 X=0 合計 3 4 観察データの背後に反事実アウトカムモデルを仮定 Y = XYx=1 + (1 – X)Yx=0 (☜ causal consistency といいます) 観察データに確率モデルを仮定 例: Y|X = 1 ~ Bin(n, p1), Y|X = 0 ~ Bin(m, p0), 互いに独立 尤度 L(p0, p1) = nCy p1y1 (1 – p1)n – y1 mCy1 p0y0 (1 – p0)m – y0 1 これが正しいためには… 20
統計学 meets 因果推論(つづき) ランダム化していない場合 1. 仮想的母集団(= 理論的な確率分布)からのランダムサンプリング 2. 決定論的(deterministic)でなく確率論的(stochastic)に Yx の値が定まる 3. でもランダムサンプリングなんかしてないし… でも「全員がおなじ確率のメカニズムにしたがう」という仮定が… 標本空間に確率を付与せず、統計的推測をあきらめる Greenland, Epidemiology 1990 反事実モデルと確率モデルは一方が他方を含意する類のものではない 統計的因果推論にはどちらも必要 21
層別解析による平均因果効果推定 22
平均因果効果 E(Yx=1) – E(Yx=0) ランダム化試験では E(Y | X = 1) – E(Y | X = 0) E(Y | X = x):曝露群(x)ごとの平均値で推定 交絡変数が全て測定されている観察研究では ∑c E(Y | X = 1, C = c)P(C = c) – ∑c E(Y | X = 0, C = c)P(C = c) E(Y|X = x, C = c):層(c)×曝露群(x)ごとに推定 P(C = c): 層(c)の大きさ 23
仮想例 層k Sex 1 C X=1 X=0 N Y=1 Age N Y=1 Male Young 1360 577 1640 423 2 Male Old 2800 1709 1200 711 3 Female Young 2660 1707 1140 713 4 Female Old 4020 2978 980 722 10840 6971 4960 2569 Total � E(Y|X = 1) = 64.3% 未調整リスク差 = 64.3% – 51.8% = 12.5% � E(Y|X = 0) = 51.8% 24
Step 1:層ごとの平均値 Age N X=1 Y=1 Male Young 1360 577 1640 423 Male Old 2800 1709 1200 711 層k Sex 1 2 C Risk N X=0 Y = 1 Risk 3 Female Young 2660 1707 1140 713 4 Female 4020 2978 980 722 10840 6971 4960 2569 Total Old 25
Step 1:層ごとの平均値 Age N X=1 Y=1 Male Young 1360 577 0.424 1640 423 0.258 Male Old 2800 1709 0.61 1200 711 0.593 層k Sex 1 2 C Risk N X=0 Y = 1 Risk 3 Female Young 2660 1707 0.642 1140 713 0.625 4 Female 4020 2978 0.741 980 722 0.737 10840 6971 Total Old 4960 2569 � � = 1|X, C) E(Y|X, C) = P(Y 26
Step 2:層ごとの重み C X=1 N Risk X=0 N Risk 層k Sex 1 Male Young 1360 0.424 1640 0.258 2 Male Old 2800 0.61 1200 0.593 Age 3 Female Young 2660 0.642 1140 0.625 4 Female 4020 0.741 980 0.737 Total Old 10840 4960 Nk � P(C) 15800 27
Step 2:層ごとの重み C X=1 N Risk X=0 N Risk Nk � P(C) 層k Sex 1 Male Young 1360 0.424 1640 0.258 3000 2 Male Old 2800 0.61 1200 0.593 4000 0.25 Age 0.19 3 Female Young 2660 0.642 1140 0.625 3800 0.24 4 Female 4020 0.741 980 0.737 5000 0.32 4960 15800 Total Old 10840 � | X = x, C = c) P(C � = c) ∑c E(Y 28
Step 3:重み付け平均 ◆:E(Yx=1) 集団全体が X = 1 だった場合 � | X = 1, C = c)P(C � = c) = 0.424(0.19) + ... + 0.741(0.32) = 62.4% ∑c E(Y ◇:E(Yx=0) 集団全体が X = 0 だった場合 � | X = 0, C = c)P(C � = c) = 0.258(0.19) + ... + 0.737(0.32) = 58.3% ∑c E(Y 調整リスク差 = 62.4% – 58.3% = 4.1% 「標準化リスク差」とも 29
層別解析、Y から見るか? X から見るか? https://filmarks.com/movies/33479 30
因果効果の求め方(1) Y から見る 層別解析 1. 2. 3. 交絡変数 C の組み合わせ c で層に分ける 層ごとにアウトカムの平均 E(Y|X = x, C = c) を計算 層ごとの重み P(C = c) で E(Y|X = x, C = c) を重み付け平均 x = 0, 1 で別々に行うことで集団全体の効果が求められる E(Yx=1) vs. E(Yx=0) 31
因果効果の求め方(2) X から見る 各層の曝露確率 P(X = 1|C = c) を求める 傾向スコア(propensity score:PS) 32
各層の傾向スコア C X=1 N Y=1 X=0 � = 1|C) P(X N Y=1 層k Sex 1 Male Young 1360 577 1640 423 0.453 2 Male Old 2800 1709 1200 711 0.700 3 Female Young 2660 1707 1140 713 0.700 4 Female Old 4020 2978 980 722 0.804 10840 6971 4960 2569 Total Age 傾向スコアがおなじ層をまとめてみよう 33
傾向スコア PS で層別したデータ X=1 Y = 1 Risk � = 1|C) X=0 P(X = PS Y = 1 Risk � P(PS) 層k N 1 1360 577 0.424 1640 423 0.258 0.453 2&3 5460 3416 0.626 2340 1424 0.609 0.700 0.49 4 4020 2978 0.741 980 0.804 0.32 Total 10840 6971 N 722 0.737 4960 2569 0.19 1 34
傾向スコアで「層別」 PS で層別解析 � � = 1, PS)P(PS) = 0.424(0.19) + 0.626(0.49) + 0.741(0.32) = 62.4% ∑PS E(Y|X � � ∑PS E(Y|X = 0, PS)P(PS) = 0.258(0.19) + 0.609(0.49) + 0.737(0.32) = 58.3% C そのもので層別解析した結果とおなじ 35
傾向スコアの推定 C で層別できるなら、層 k の曝露割合 曝露 nk 人、非曝露 mk 人 � = 1| 層 k) = nk/(nk + mk) P(X Sato & Matsuyama, Epidemiology 2003 できないなら、モデル化 (☞次回) 例:ロジスティック回帰 log P(X = 1|C = c) = α0 + α1c1 + α2c2 + α3c3 +… 1 – P(X = 1|C = c) � = 1|C = c) = P(X �0 + α � 1c1 + α � 2c2 + α � 3c3 +…) exp(α �0 + α � 1c1 + α � 2c2 + α � 3c3 +…) 1 + exp(α 36
傾向スコア 交絡変数 C で層別した中での条件付き曝露確率 PS(c) = P(X = 1|C = c) 「対象者の曝露確率」ではない 「交絡変数により定義されるサブグループ」での確率 交絡変数の定義が必要 交絡の定義が必要 交絡変数の組み合わせごとに傾向スコアも異なる 「その人特有の曝露確率」ではない 37
バランシング・スコア Balancing score:b(Z) 何らかの観察変数 Z の実数値関数で、 同じスコア b(Z) の集団の中では曝露 X と Z が独立 = 曝露群間で Z がバランス 最も細かい(finest)バランシング・スコア b(Z): Z の全組み合わせ水準 b(z) = (Z = z となる組み合わせ水準) 最も粗い(coarsest)バランシング・スコア b(z) = P(X = 1 | Z = z) 粗いスコアは細かいスコアの関数で表せる Rosenbaum & Rubin, Biometrika 1983 38
傾向スコアはバランシング・スコア 傾向スコア: PS(C) 交絡変数 C の最も粗いバランシング・スコア PS(c) = P(X = 1 | C = c) 傾向スコア PS(C) が同じ値の集団の中では 曝露 X と交絡変数 C が独立 曝露群間で交絡変数 C のバランスがとれる 39
条件付きランダム化の仮定の下で 条件付きランダム化の仮定 交絡変数 C で層別すれば Yx(x = 0, 1)と曝露 X が独立 傾向スコア PS(C) のバランス性 PS(C) で層別すれば曝露 X と交絡変数 C が独立 PS(C) で層別すれば Yx (x = 0, 1) と曝露 X が独立 Rosenbaum and Rubin, Biometrika 1983 40
傾向スコアのバランス性 C X=1 N Y=1 X=0 � = 1|C) P(X N Y=1 層k Sex 1 Male Young 1360 577 1640 423 0.453 2 Male Old 2800 1709 1200 711 0.700 3 Female Young 2660 1707 1140 713 0.700 4 Female Old 4020 2978 980 722 0.804 10840 6971 4960 2569 Total Age � � P(Male, Old | X = 1) = 2800/5460 = 0.513 = 1200/2340 = P(Male, Old | X = 0) � � P(Female, Young| X = 1) = P(Female, Young| X = 1) (= 0.487) � � P(Male, Young| X = 1) = P(Male, Young| X = 1) (= 0) � � P(Female, Old| X = 1) = P(Female, Old| X = 0) (= 0) 41
傾向スコア、3通りの使い方 42
因果効果の求め方 傾向スコアの推定 C で層別できるなら、層 k の曝露割合 C で層別できないなら、モデルで P(X = 1|C) を近似 3つの方法 層別・回帰調整(stratification/adjustment) マッチング(matching) 重み付け(weighting) 43
傾向スコア利用法(1):層別・調整 adjustment(= conditioning) C で層別できるなら C で層別しても PS(C) = P(X = 1|C)で層別してもおなじ 44
C で層別できないときのPS層別・調整 � � = 1|C = c) モデルから PS(c) = P(X � PS(C) で層別 分位点でサブグループ化 例:0~0.3, 0.3~0.5, …, 0.8~1 各グループ内で比較:Y1k/nk – Y0k/mk k = 1,…,5 各グループの結果を併合 45
傾向スコア層別の注意点 層内でも傾向スコアは連続的に変化 交絡変数のバランスが完全にはとれない residual confounding 代わりに、傾向スコアの共変量調整を行うことも(後述) 回帰モデルで交絡変数の代わりに傾向スコアを入れる 交絡変数でなく傾向スコアで層別した結果を近似 曝露の回帰係数 = 傾向スコアが同じ層内の曝露効果 モデルが正しい必要あり ただし、例外あり 46
傾向スコア利用法(2):マッチング Matching 1. 2. 3. 傾向スコア PS(Ci) をデータ i ごとに求める 曝露群と非曝露群で PS(Ci) の同じ or 近い人をペアにする ペアに選ばれた人だけ解析に利用する 傾向スコアでマッチされた集団 Propensity-matched cohort 各ペアは「同じ交絡変数分布」からのランダムサンプル 個々のペア内で交絡変数が似ている保証はないが… 全体として「交絡がない」とみなして解析 47
C で層別できるとき 傾向スコア PS の値がおなじ層で人数をそろえる 層k 1 2&3 X=1 N Y=1 1360 5460 4 4020 Total 10840 577 3416 ⇒2340/5460 2978 ⇒980/4020 6971 X=0 N Y=1 1640 423 ⇒1360/1640 � = 1|C) P(X 0.453 2340 1424 0.700 980 722 0.804 4960 2569 48
P(X = 1|C) でマッチしたコホート PS matched cohort 層k X=1 N Y=1 X=0 N Y=1 � = 1|C) P(X 1 1360 577 1360 351 2&3 2340 1464 2340 1424 0.700 4 980 726 980 722 0.804 Total 4680 2767 4680 2497 � E(Y|X = 1, PS-matched) = 2767/4680 = 59.1% � E(Y|X = 0, PS-matched) = 2497/4680 = 53.4% 0.453 49
標的集団 効果の定義される集団のこと 層別解析(標準化) � | X = x, C = c)P(C � = c) ∑c E(Y � | X = x, PS = p)P(PS � ∑p E(Y = p) 標的集団:もとの集団全体 E(Yx=1) – E(Yx=0) マッチング � E(Y|X = 1, PS-matched) 標的集団:もとの集団からマッチされた集団 E(Yx=1| PS-matched) – E(Yx=0| PS-matched) 50
交絡変数の分布 P(C) の変化 マッチング前 層k 1 2&3 4 Total PSマッチング後 層k 1 2&3 4 Total X=1 N Y=1 1360 577 5460 3416 4020 2978 10840 6971 X=0 N Y=1 1640 423 2340 1424 980 722 4960 2569 � = 1|C) P(X = PS 0.453 0.700 0.804 X=1 N Y=1 1360 577 2340 1464 980 726 4680 2767 X=0 N Y=1 1360 351 2340 1424 980 722 4680 2497 � = 1|C) P(X = PS 0.453 0.700 0.804 � P(PS) 0.19 0.49 0.32 1 � P(PS) 0.29 0.5 0.21 1 51
PSの分布を変えないマッチング 層k X=1 N Y=1 X=0 N Y=1 1 1360 577 1640 423 2&3 5460 3416 2340 1424 � = 1|C) P(X = PS 0.453 ⇒(980/1640)*(0.19/0.32) 0.700 ⇒(980/2340)*(0.49/0.32) 4 4020 2978 980 Total 10840 6971 4960 722 2569 0.804 � P(PS) 0.19 0.49 0.32 1 52
PSの分布を変えないマッチング 層k X=1 N Y=1 X=0 N Y=1 � = 1|C) P(X = PS 1 1360 577 588 152 2&3 5460 3416 1529 930 0.700 0.49 980 722 0.804 0.32 3097 1804 ⇒588/1360 0.453 � P(PS) 0.19 ⇒1529/5460 4 Total 4020 2978 ⇒980/4020 10840 6971 1 53
PSの分布を変えないマッチング 層k X=1 N Y=1 X=0 N Y=1 0.453 � P(PS) 1 588 249 588 152 2&3 1529 957 1529 930 0.700 0.49 4 980 726 980 722 0.804 0.32 Total 3097 1932 3097 1804 � E(Y|X = 1, PS-matched) = 1932/3097 = 62.4% � E(Y|X = 0, PS-matched) = 1804/3097 = 58.3% � = 1|C) P(X = PS 0.19 1 集団全体を標的とした標準化推定量に一致 54
傾向スコア利用法(3):重み付け Inverse probability weighting(IPW) 個人 i が実際に受けた曝露を受ける確率の逆数で重み付け 曝露群では Wi = 1/PSi 非曝露群では Wi = 1/(1 – PSi) Wi = 1 P(X = xi|C = ci) 重み付けた集団で、そのまま差を計算 Rosenbaum, JASA 1987 Wi で重み付けた集団では 曝露 X と交絡変数 C が独立(statistically exogenous) 条件付きランダム化の仮定と合わせることで Yx(x = 0, 1)と曝露 X が独立(causally exogenous) Hernán, Brumback & Robins, JASA 2001 55
� = 1|C) 1/P(X IPWの計算 Age N X=1 NY=1 Male Young 1360 577 2.21 1640 423 1.83 Male Old 2800 1709 1.43 1200 711 3.33 0.700 Female Young 2660 1707 1.43 1140 713 3.33 0.700 Female 4020 2978 1.24 980 722 5.10 0.804 10840 6971 4960 2569 Sex C Old Total IPW N X=0 NY=1 IPW � = 1|C) P(X 0.453 � = 0|C) 1/P(X 56
IPWで重み付け Sex C X=1 NY=1,IPW IPW NIPW Age NIPW Male Young 3000 1272.8 2.21 3000 Male Old 4000 2441.4 1.43 Female Young 3800 2438.6 Female 5000 3704.0 15800 9856.8 Old Total X=0 � = 1|C) P(X NY=1,IPW IPW 773.8 1.83 0.453 4000 2370.0 3.33 0.700 1.43 3800 2376.6 3.33 0.700 1.24 5000 3683.7 5.10 0.804 15800 9204.1 57
IPWで重み付け Sex C X=1 NY=1,IPW IPW NIPW Age NIPW Male Young 3000 1272.8 2.21 3000 Male Old 4000 2441.4 1.43 Female Young 3800 2438.6 Female 5000 3704.0 15800 9856.8 Old Total X=0 � = 1|C) P(X NY=1,IPW IPW 773.8 1.83 0.453 4000 2370.0 3.33 0.700 1.43 3800 2376.6 3.33 0.700 1.24 5000 3683.7 5.10 0.804 15800 9204.1 � IPW(Y|X = 1) = 9856.8/15800 = 62.4% E � IPW(Y|X = 0) = 9204.1/15800 = 58.3% E 58
統計モデル アウトカム回帰モデル E(Y|X = x, C = c) を近似 傾向スコアモデル P(X = 1|C = c) を近似 59
回帰モデルの係数 例:ロジスティック回帰モデル P(Y = 1|X = x, C = c) = γ0 + γ1x + γ2c log 1 – P(Y = 1|X = x, C = c) exp(γ1) : 「調整」オッズ比 P(Yx=1|C = c) vs. P(Yx=0|C = c) 曝露 X = x を受けているので P(Y = 1|X = x, C = c) = P(Yx=1 = 1|X = x, C = c) C の層内では交絡がないので P(Yx=1 = 1|X = x, C = c) = P(Yx=1 = 1|C = c) C = c という層内での因果効果 = 条件付き効果 60
層ごとの効果 C X=1 N Risk X=0 N Risk Risk Odds ratio Difference 層k Sex 1 Male Young 1360 0.424 1640 0.258 0.166 2.12 2 Male Old 2800 0.61 1200 0.593 0.017 1.08 3 Female Young 2660 0.642 1140 0.625 0.017 1.07 4 Female Old 4020 0.741 980 0.737 0.04 1.02 Total Age 10840 4960 E(Yx=1|C = c) vs. E(Yx=0|C = c) 交絡変数 C = c の層の全員が X = 1 の場合 vs. X = 0 の場合 条件付き効果(⇔ 周辺効果 E(Yx=1) vs. E(Yx=0) ) 61
回帰モデルによる効果推定 あくまで交絡変数 C で層別したい 交絡変数が全て測定されている仮定 できないので、C で層別された層ごとのアウトカム平均を近似 E(Y|X, C) を十分近似できている仮定 回帰係数は条件付き効果 サブグループ効果が等しいと仮定(交互作用項をいれないとき) 周辺効果(平均因果効果)を知るには、もうひと手間必要(☞次回) 集団全体が X = 1 だった場合 vs. X = 0 だった場合 62
傾向スコアモデルによる効果推定 あくまで交絡変数 C で層別したい 交絡変数が全て測定されている仮定 できないので、C で層別された層ごとの曝露確率を近似 P(X = 1|C) を十分近似できている仮定 求まる効果はその後の解析方法しだい 層別・回帰調整(+ 標準化) マッチング IPW 63
傾向スコア PS の層ごとの効果 X=1 N Risk X=0 Risk N Risk Difference 層k P(X=1|C) = PS 1 0.453 1360 0.424 1640 0.258 0.166 2&3 0.700 5460 0.626 2340 0.609 0.017 4 0.804 4020 0.741 980 0.04 Total 10840 0.737 4960 E(Yx=1| PS = p) – E(Yx=0| PS = p) 傾向スコア PS = p の層の全員が X = 1 の場合 vs. X = 0 の場合 64
傾向スコアを回帰モデルで調整 例:ロジスティック回帰モデル P(Y = 1|X = x, PS = p) = γ0 + γ1x + γ2p log 1 – P(Y = 1|X = x, PS = p) exp(γ1) 傾向スコア PS が同じ層内の曝露効果 モデルが正しい必要あり PS の変化に対してリスク変化が一定 PS の値によらず曝露効果一定 ただし、線形モデル、対数線形モデルの場合はモデル誤特定があっても解釈可 Vansteelandt & Daniel, Stat Med 2014 Dukes & Vansteelandt, AJE 2018 65
傾向スコアにより推定される効果 1. 傾向スコアの層別、回帰調整 そのままではPSの値が同じ層での条件付き効果 (回帰)標準化をつかえば周辺効果 2. 傾向スコアマッチング マッチングされた集団における周辺効果 ただし、標的集団の変化に注意 3. 傾向スコアによる重み付け(IPW) 集団全体での周辺効果 他の集団での周辺効果の重みも作れる Sato & Matsuyama (Epi 2003), Li et al. (JASA 2018), Dahabreh et al. (Stat Med 2020) 66
まとめ 観察研究データから因果的命題を導くのはたいへん 因果命題は反事実命題 反事実モデルにおける仮定が必要 仮定さえ置けば、やりようはある 層別解析 回帰 傾向スコア 層別解析をどちらから見るか 仮定に対する感度はどの手法も共有 条件付きランダム化の仮定 交絡変数 C が高次元になるとモデル化の仮定(次回くわしく) 67
Take-home message 交絡調整 = 層別解析 交絡変数でサブグループに分けて(= 層別)比較 交絡調整における2種類の統計モデル アウトカム回帰モデル 傾向スコアモデル E(Y|X, C) P(X = 1|C) モデルをつかう目的 層別できないくらいの交絡変数があるとき 「もし層別できた場合」の結果を近似 アウトカム回帰 傾向スコア 各モデルの推定自体が目的ではない 68
文献 Dahabreh IJ, et al. Extending inferences from a randomized trial to a new target population. Stat Med 2020;39:1999-2014. Dukes O, Vansteelandt S. A note on g-estimation of causal risk ratios. Am J Epidemiol. 2018;187:1079-1084. Greenland S. Randomization, statistics, and causal inference. Epidemiology 1990;1:421-41 Greenland S, Robins JM. Identifiability, exchangeability, and epidemiological confounding. Int J Epidemiol. 1986;15:413419. Greenland S, Robins JM, Pearl J. Confounding and Collapsibility in Causal Inference. Stat Sci. 1999;14:29-46. Hernán MA. A definition of causal effect for epidemiological research. J Epidemiol Community Health. 2004;58:265-271. Hernán MA, Brumback B, Robins JM. Marginal structural models to estimate the joint causal effect of nonrandomized treatments. J Am Stat Assoc 2001;96:440-448 Hernán MA, Robins JM. Causal Inference: What If. Chapman & Hall/CRC, 2020 Li F, et al. Balancing covariates via propensity score weighting. J Am Stat Assoc 2018;113:390-400 Rosenbaum PR. Model-based direct adjustment. J Am Stat Assoc 1987;82:387-394 Rosenbaum PR, Rubin DB. The central role of the propensity score in observational studies for causal effects. Biometrika 1983;70:41-55 Sato T, Matsuyama Y. Marginal structural models as a tool for standardization. Epidemiology 2003;14:680-686 Vansteelandt S, Daniel RM. On regression adjustment for the propensity score. Stat Med. 2014;33:4053-4072. 69
読書案内 Hernán MA, Robins JM. Causal Inference: What If. Chapman & Hall/CRC, 2020. 現時点で最も薦められる因果推論入門書。全編ウェブ公開 Westreich D. Epidemiology by Design: A Causal Approach to the Health Sciences. Oxford Univ Press, 2019. 理論疫学の入門書でもあり、因果推論の入門書でもある。独自の構成だが読みやすく取り扱っている内容のバランスがよい 佐藤, 松山. 交絡という不思議な現象と交絡を取りのぞく解析:標準化と周辺構造モデル. 計量生物学 2011;32:S35–S49. レビュー論文。今回のコース内容の復習にも適している 松山. 18章:因果推論. In: 丹後, 松井編. 新版 医学統計学ハンドブック. 朝倉書店, 2018. 日本語でのレビュー。下記2002年の原稿のアップデート版 佐藤, 松山. 疫学・臨床研究における因果推論. In: 甘利他編. 多変量解析の展開―隠れた構造と因果を推理する. 岩波書店, 2002. 日本語で最初の医学分野での因果推論レビュー。Robins-Greenland的な因果推論 田中. 医学のための因果推論の基礎概念. 計量生物学2019;40:35–62. 割り付けメカニズム(特にランダム化)にもとづく確率モデルを厳格に運用したRubin流の因果推論を、臨床試験の枠組みで解説 黒木学. 構造的因果モデルの基礎. 共立出版, 2017. いわゆるPearl流の因果モデル(構造的因果モデル)の解説書。Rubinの因果モデルなど異なる因果モデル間の考察など類書にない貴重な記載も多い Rothman KJ, Greenland S, Lash TL. Modern Epidemiology, 3rd ed. LWW 2008. 理論疫学のバイブルだが、因果推論の思想が通底しているので疫学外の観察研究データ研究者にも役立つ(4版ではGreenlandが降りました) 田栗, 篠崎.因果推論.In:丹後・松井編. 臨床試験の事典. 朝倉書店,2022予定 反事実因果モデル、交絡、層別解析、回帰モデル、傾向スコア、操作変数法、IPW法について各2~4ページで解説 70
教科書や解説論文もいいですが、クラシックも… 概論 Greenland S. Interpretation and choice of effect measures in epidemiologic analysis. Am J Epidemiol. 1987;125:761-8. Greenland S. Randomization, statistics, and causal inference. Epidemiology. 1990;1:421-429. Greenland S, Brumback B. An overview of relations among causal modelling methods. Int J Epidemiol. 2002;31:1030-1037. Robins JM. Data, design, and background knowledge in etiologic inference. Epidemiology. 2001;12:313-20. Hernán MA, Hernández-Díaz S, Werler MM, Mitchell AA. Causal knowledge as a prerequisite for confounding evaluation: an application to birth defects epidemiology. Am J Epidemiol. 2002;155:176–84. 交絡 Greenland S, Robins JM. Identifiability, exchangeability, and epidemiological confounding. Int J Epidemiol. 1986;15:413-419. Greenland S, Robins JM, Pearl J. Confounding and collapsibility in causal inference. Stat Sci. 1999;14:29-46. Greenland S. Absence of confounding does not correspond to collapsibility of the rate ratio or rate difference. Epidemiology. 1996;7:498-501. 因果DAG Greenland S, Pearl J, Robins JM. Causal diagrams for epidemiologic research. Epidemiology. 1999;10:37-48. Hernán MA, Hernández-Díaz S, Robins JM. A structural approach to selection bias. Epidemiology. 2004;15:615-625. Greenland S. Quantifying biases in causal models: classical confounding vs collider-stratification bias. Epidemiology. 2003;14:300-306. モデル Robins JM, Greenland S. The role of model selection in causal inference from nonexperimental data. Am J Epidemiol. 1986;123:392-402. Greenland S. Summarization, smoothing, and inference in epidemiologic analysis. Scand J Soc Med. 1993;21:227-32. Robins JM, Hernán MA, Brumback B. Marginal structural models and causal inference in epidemiology. Epidemiology. 2000;11:550-560. 71