因果推論II：モデルを用いた効果推定

統計数理研究所 医療健康データ科学研究センター 公開講座「疫学・公衆衛生統計」 因果推論 II 東京理科大学 工学部情報工学科 篠崎 智大 [email protected] 2021年12月10日 15：00～16：30

「因果推論」のアウトライン 因果推論 I （12月9日）：因果推論の基礎 効果、交絡、交絡変数、交換可能性、層別解析、回帰、傾向スコア 因果推論 II （12月10日）：モデルを用いた効果推定 回帰モデル（条件付き効果、周辺効果、回帰標準化） 傾向スコアモデル（回帰調整、マッチング、逆確率重み付け IPW） 二重ロバスト推定（augmented IPW / IPW回帰標準化） 2

Take-home message（前回） 交絡調整 = 層別解析 交絡変数でサブグループに分けて（= 層別）比較 交絡調整における2種類の統計モデル アウトカム回帰モデル 傾向スコアモデル  E(Y|X, C)  P(X = 1|C) モデルをつかう目的 層別できないくらいの交絡変数があるとき 「もし層別できた場合」の結果を近似   アウトカム回帰 傾向スコア 各モデルの推定自体が目的ではない 3

復習：記法 観察変数 X： 曝露（治療） Y： アウトカム C： 交絡変数 潜在アウトカム変数 Yx=1： X = 1 を受けた場合に観察されるであろうアウトカム Yx=0： X = 0 を受けた場合に観察されるであろうアウトカム X の値に応じて一方だけ観察される E(Y), E(Yx) E(Y|X = x), E(Yx|C = c), E(Y|C = c, X = x) 4

関心のある「平均因果効果」 Observed Population Unexposed Causation E(Y x=0) E(Y x=1) Exposed Association E(Y|X = 0) E(Y|X = 1) Hernán, JECH 2004 Hernán & Robins, 2020 5

層別できるなら 交絡変数が C としてすべて測定されていれば 交絡変数 C で層別して Yx（x = 0, 1）と曝露 X が独立であれば E(Yx) = ∑c E(Y|X = x, C = c)P(C = c) E(Y|X = x, C = c) ：各層（C = c）内の曝露ごと（X = x）の回帰 P(C = c) ：各層（C = c）の集団に占める割合 6

層別解析 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 7

層別解析 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 ∑c E(Y | X = 1, C = c)P(C = c) = 0.424(0.19) + ... + 0.741(0.32) = 62.4% ∑c E(Y | X = 0, C = c)P(C = c) = 0.258(0.19) + ... + 0.737(0.32) = 58.3% 8

回帰モデル モデルで E(Y|X = x, C = x) を近似 例：ロジスティック回帰モデル log P(Y = 1|X = x, C = c) = γ0 + γ1x + γ2c1 + γ3c2 +… 1 – P(Y = 1|X = x, C = c) � = 1|X = x, C = c) = � xc = P(Y P exp(�γ0 + �γ1x + �γ2c1 + �γ3c2 +…) 1 + exp(�γ0 + �γ1x + �γ2c1 + �γ3c2 +…) 9

回帰モデルの当てはめ Sex 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 Age 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 X 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 Y 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 N 577 1709 1707 2978 423 711 713 722 783 1091 953 1042 1217 489 427 258 proc logistic desc; model y = x sex age; freq n; output out = pred p = R; run; OR = exp(γ) 95% CI 0.2404 1.272 1.183 1.368 Sex (C1) -0.8674 0.420 0.393 0.449 Age (C2) 0.7425 2.101 1.965 2.247 変数 係数 γ 1（切片） 0.2509 X 10

回帰モデルの予測値 Age N X=1 Risk Male Young 1360 Male Old 層k Sex 1 2 C X=0 Risk � XC P 0.424 � XC P 0.407 1640 0.258 0.351 2800 0.61 0.591 1200 0.593 0.531 N 3 Female Young 2660 0.642 0.620 1140 0.625 0.562 4 Female 4020 0.741 0.774 980 0.737 0.730 Total Old 10840 4960 � xc = P exp(�γ0 + �γ1x + �γ2sex + �γ3age) 1 + exp(�γ0 + �γ1x + �γ2sex + �γ3age) 11

回帰モデルによる標準化 E(Yx) = ∑c P(Y = 1|X = x, C = c)P(C = c) � xcで近似 ただし、P(Y = 1|X = x, C = c) を P � 1c P(C � = c) ∑c P = 0.407(0.19) + ... + 0.774(0.32) = 62.1% � 0c P(C � = c) ∑c P = 0.351(0.19) + ... + 0.730(0.32) = 56.7% 12

こう考えても同じこと 個人 i の交絡変数の値 Ci ごとに P(Y = 1|X = x, C = Ci) を予測 � 1i = 1/{1 + exp(– γ�0 – γ�1 – γ�2C1i – γ�3C2i –…)} x=1：P � 0i = 1/{1 + exp(–�γ0 x=0：P – γ�2C1i – γ�3C2i –…)} P(Y = 1|X = x, C = c) log 1 – P(Y = 1|X = x, C = c) = γ0 + γ1x + γ2c1 + γ3c2 +… � = 1|X = x, C = Ci) = � xi = P(Y P exp(�γ0 + �γ1x + �γ2C1i + �γ3C2i +…) 1 + exp(�γ0 + �γ1x + �γ2C1i + �γ3C2i +…) 13

個人単位で見た標準化 C Nk � xi P 層k Sex 1 Male Young 3000 0.258 2 Male Old 4000 0.593 0.61 3 Female Young 3800 0.625 0.642 4 Female Old 5000 0.737 0.741 Age Total � 1i – P � 0i を全員（n = 15800）で平均 P x=0 x=1 0.424 15800 g-computation とも呼ばれるが古典的な「回帰モデルによる標準化」に外ならない 佐藤, 統計数理 1994; Vansteelandt & Keiding, AJE 2011 層別しきれないほどの交絡変数 C にも適用できる 14

周辺効果、条件付き効果（リスク差） P(Y = 1|X = x, C = c) Risk ： 割合で推定 � XC ： ロジスティック回帰モデルで推定 P C 条件付き効果 X=1 � XC Risk P X=0 � XC Risk P リスク差 � XC Risk P 層k Sex 1 Male Young 0.424 0.407 0.258 0.351 0.166 0.056 2 Male Old 0.61 0.591 0.593 0.531 0.017 0.059 Age 3 Female Young 0.642 0.620 0.625 0.562 0.017 0.058 4 Female 0.741 0.774 0.737 0.730 0.040 0.045 0.624 0.621 0.583 0.567 0.041 0.054 標準化 Old 周辺効果 15

周辺効果、条件付き効果（オッズ比） P(Y = 1|X = x, C = c) Risk ： 割合で推定 � XC ： ロジスティック回帰モデルで推定 P C 条件付き効果 X=1 � XC Risk P X=0 � XC Risk P 層k Sex 1 Male Young 0.424 0.407 0.258 0.351 オッズ比 � XC Risk P 2 Male Old 0.61 0.591 0.593 0.531 1.08 1.27 Age 2.12 1.27 3 Female Young 0.642 0.620 0.625 0.562 1.07 1.27 4 Female 0.741 0.774 0.737 0.730 1.02 1.27 0.624 0.621 0.583 0.567 1.19 1.25 標準化 Old 周辺効果 16

前ページの表より ロジスティック回帰モデル 「条件付きオッズ比が一定」という制約下での近似  交互作用項を入れることで緩めることも可能 オッズ比 ロジスティックモデルの結果   条件付きオッズ比 = 1.27（各層で一定） 周辺オッズ比 = 1.25（！） 併合不能性 non-collapsibility   周辺効果指標が条件付き効果指標の重み付け平均で表せない バイアスではなく、オッズ比（やハザード比）という指標の性質 Greenland, Robins & Pearl, Stat Sci 1999 Shrier & Pang. JCE 2015 17

層別できるなら Part 2 各層の曝露確率 P(X = 1|C = c) を求める 傾向スコア propensity score（PS） 3つの方法（C で層別できたらどれもおなじ） 層別・回帰調整（stratification/adjustment） マッチング（matching） 逆確率重み付け（inverse probability weighting） 18

1/PS IPWの計算 Age N Male Young 1360 577 2.21 1640 423 1.83 0.453 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 Old Total IPW N X=0 NY=1 � = 1|C) P(X = PS IPW X=1 NY=1 C 1/(1 – PS) 19

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 � = 1|C) X=0 P(X = PS 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 20

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 � = 1|C) X=0 P(X = PS 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 21

傾向スコアモデル モデルで P(X = 1|C = c) を近似 例：ロジスティック回帰モデル P(X = 1|C = c) log 1 – P(X = 1|C = c) = α0 + α1c1 + α2c2 + α3c3 +… � c = P(X � = 1|C = c) = PS exp(� α0 + α � 1c 1 + α � 2c 2 + α � 3c3 +…) α0 + α � 1c 1 + α � 2c 2 + α � 3c3 +…) 1 + exp(� 22

傾向スコアモデルの当てはめ Sex 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 Age 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 X 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 Y 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 N 577 1709 1707 2978 423 711 713 722 783 1091 953 1042 1217 489 427 258 proc logistic desc; model x = sex age; freq n; output out = pred p = ps; run; 変数 係数α 1（切片） 0.7331 OR = exp(α) 95% CI Sex -0.7956 0.451 0.421 0.484 Age 0.8012 2.228 2.078 2.389 23

傾向スコアモデルの予測値 Age N X=1 NY=1 N X=0 NY=1 Male Young 1360 577 1640 423 0.484 0.453 Male Old 2800 1709 1200 711 0.677 0.700 Female Young 2660 1707 1140 713 0.700 0.675 Female 4020 2978 980 722 0.804 0.823 10840 6971 4960 2569 Sex C Old Total IPW �c = PS IPW exp(�γ0 + �γ1c) � PS 1 + exp(�γ0 + �γ1c) 24

傾向スコアモデルによるIPW Age N X=1 NY=1 Male Young 1360 577 2.06 1640 423 1.94 0.484 Male Old 2800 1709 1.48 1200 711 3.09 0.677 Female Young 2660 1707 1.48 1140 713 3.08 0.675 Female 4020 2978 1.22 980 722 5.64 0.823 10840 6971 4960 2569 Sex C Old Total �c 1/PS IPW N X=0 NY=1 IPW PS � c) 1/(1 – PS 25

傾向スコアモデルによるIPW Sex C X=1 NY=1,IPW IPW X=0 NY=1,IPW IPW PS Age NIPW Male Young 2807.7 1191.2 2.06 3180.7 820.4 1.94 0.484 Male Old 4137.7 2525.5 1.48 3711.9 2199.3 3.09 0.677 Female Young 3937.9 2527.1 1.48 3512.9 2197.1 3.08 0.675 Female 4886.8 3620.1 1.22 5525.2 4070.6 5.64 0.823 Old Total 15770.0 9863.8 NIPW 15930.7 9287.4 � IPW(Y|X = 1) = 9863.8/15770.0 = 62.5% E � IPW(Y|X = 0) = 9287.4/15930.7 = 58.3% E 26

層別しきれないほどのデータの解析例 27

例：OSSIデータベース 整形外科手術部位感染 Orthopedic surgical site infection 整形外科手術患者 約10,000名 手術前：性・年齢・糖尿病・手術部位、… 手術中：手術時間・出血量・術中体温、… 手術後：予防的抗菌薬・感染・転帰、… 術後30日以内の感染症を比較 Yamada et al. CID 2020 28

正常体温 vs. 低体温 正常体温 低体温 30日以内感染症 あり なし 249 (3.2%) 7584 29 (2.9%) 979 合計 7833 1008 リスク差 0.30%（95%信頼区間 -0.80%, 1.40%） リスク比 1.10 （95%信頼区間 0.76, 1.61） オッズ比 1.11 （95%信頼区間 0.75, 1.64） 術前変数が交絡変数になっているのでは？ 29

手術前変数 Table 1 Female sex — no. (%) Age (yr)* BMI (kg/m2)*† BMI — no. (%)† Present smoker — no. (%) Disease BMI<25 25<BMI<30 30<BMI Arthritis Fractures & dislocations (excluding spine, tumor, and pseudoarthrosis) Spine Tumor (excluding spinal tumors) Others 低体温 正常体温 Total (N=8,841) 5182 (58.6) 66.1 ± 17.1 23.4 ± 4.1 6146 (69.5) 2135 (24.2) 557 (6.3) 1188 (13.4) 2119 (24.0) Hypothermia (N=1,008) 661 (65.6) 68.6 ± 16.6 22.7 ± 3.9 767 (76.2) 195 (19.4) 45 (4.5) 142 (14.1) 313 (31.1) Normothermia (N=7,833) 4521 (57.7) 65.8 ± 17.2 23.5 ± 4.1 5379 (68.7) 1940 (24.8) 512 (6.5) 1046 (13.4) 1806 (23.1) 3419 (38.7) 422 (41.9) 2997 (38.3) 2464 (27.9) 127 (1.4) 712 (8.1) 168 (16.7) 14 (1.4) 91 (9.0) 2296 (29.3) 113 (1.4) 621 (7.9) P Value <.0001 <.0001 <.0001 <.0001 0.52 <.0001 群間で異なる分布 これらの変数によって感染率が異なる可能性が高い Yamada et al. CID 2020 30

たくさんの変数 性、年齢、BMI、糖尿病、手術部位、ASA（米国麻酔科学会）重症度分類、… なかには連続変数も とりあえず4変数を2値化して層別してみよう 24 = 16 の層ができるはず 31

16サブグループに層別 男性 ASA 3以上 糖尿病 年齢 70以上 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 0 1 0 0 0 1 0 1 0 1 1 0 0 1 1 1 1 0 0 0 1 0 0 1 1 0 1 0 1 0 1 1 1 1 0 0 1 1 0 1 1 1 1 0 1 1 1 1 正常体温 感染あり なし 20 1502 83 2136 2 125 9 365 3 47 8 139 1 29 4 48 48 1839 38 808 5 177 11 191 3 43 8 72 2 25 4 38 低体温 感染あり なし 2 198 12 334 1 13 1 66 0 1 3 23 0 2 0 5 3 175 5 87 0 29 1 25 0 0 0 10 0 4 1 7 32

「回帰」の推定 層別解析では層ごとに割合を計算 例： (0, 0, 0, 0) の層での 正常体温でのリスク： 20/(20 + 1502) = 0.013 低体温でのリスク： 2/(2 + 198) = 0.010 しかし… 計算できない層がでてくる  (1, 1, 0, 0) の低体温群など 計算できても1人しかいないサブグループも 層に分けられてはいるけれど、回帰モデルを使って推定を安定させたい Greenland, Stat Med 1991 33

ロジスティック回帰モデル log オッズ = α0 + β×正常体温 + γ1×男性 + γ2 ×ASA高 + γ3×糖尿病 + γ4 ×70歳以上 層×曝露群ごとの回帰（左辺）を数式でまとめて表現（右辺） 推定結果 Wald 95% 信頼限界 カイ二乗値 P-値 -3.4805 238.35 <.0001 -0.2723 0.5125 0.36 0.5486 0.1271 0.1532 0.6514 10.02 0.0015 0.7937 0.1868 0.4276 1.1598 18.05 <.0001 糖尿病 -0.0826 0.1743 -0.4243 0.2591 0.22 0.6355 年齢70歳以上 0.7686 0.136 0.502 1.0352 31.93 <.0001 推定値 標準誤差 切片 -3.9866 0.2582 -4.4928 正常体温 0.1201 0.2002 男性 0.4023 ASA 3以上 34

どうやって予測値を得る？ 良い方法があります データセット OSSI （n = 8841） ID Normothermia HAI Male ASA3 DM Age70 00057 1 0 0 0 0 1 00059 1 0 0 0 1 0 00061 0 0 0 0 0 1 … Normothermia: 正常体温 = 1, 低体温 = 0 HAI （health-associated infection） : 感染あり = 1, なし = 0 35

トリック：データを複製 data OSSI_Dup; set OSSI; do Dup = 1 to 3; if Dup = 2 then do; Normothermia = 1; HAI = .; end; if Dup = 3 then do; Normothermia = 0; HAI = .; end; output; end; run; ID 00057 00057 00057 00059 00059 00059 00061 00061 00061 … Normothermia 1 1 0 1 1 0 0 1 0 HAI 0 . . 0 . . 0 . . Male ASA3 DM Age70 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 Dup 1 2 3 1 2 3 1 2 3 36

回帰モデルの当てはめ ロジスティックモデル proc logistic data = OSSI_Dup desc; model HAI = Normothermia Male ASA3 DM Age70; output out = OSSI_Dup predicted = Px; run; model: ロジスティックモデルの変数を指定  desc (= descending): 値の大きい方の確率をモデル化 output: 新しいデータセットに新しい変数を追加  predicted = Px: 「Px」という変数に予測値を格納 37

� xi ゲット モデルの予測値 P 回帰モデルの推定には1行目だけ寄与（n = 8841） 2 & 3行目はアウトカム HAI を欠測させたのでデフォルトで除外 アウトカムが欠測していても、説明変数に応じて予測値は計算される � 1i 2行目（Normothermia = 1）には P � 0i 3行目（Normothermia = 0）には P ID 00057 00057 00057 00059 00059 00059 00061 00061 00061 … Normothermia 1 1 0 1 1 0 0 1 0 HAI 0 . . 0 . . 0 . . Male ASA3 DM Age70 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 Dup 1 2 3 1 2 3 1 2 3 Px 0.03429 0.03429 0.03053 0.01493 0.01493 0.01326 0.03053 0.03429 0.03053 38

2行目と3行目ごとに平均 proc means data = OSSI_Dup mean; var Px; class Dup Normothermia; run; Dup Normothermia 1 1 2 2 1 3 0 N 7833 1008 8841 8841 Mean 0.0287699 0.0317886 0.0318463 0.0283687 E(Yx=1) vs. E(Yx=0) 信頼区間と p-値の計算 ブートストラップ法 デルタ法 （佐藤, 統計数理 1994 や Greenland, Stat Med 1991） 説明変数が多くなったり連続値を含んでいたりしても使える 39

傾向スコアモデル モデルで P(X = 1|C = c) を近似 例：ロジスティック回帰モデル P(X = 1|C = c) log 1 – P(X = 1|C = c) = α0 + α1c1 + α2c2 + α3c3 +… 共変量の組み合わせ c = (c1, c2, …) ごとに � c = P(X � = 1|C = c) = PS exp(� α0 + α � 1c 1 + α � 2c 2 + α � 3c3 +…) α0 + α � 1c 1 + α � 2c 2 + α � 3c3 +…) 1 + exp(� 共変量の値が個人 i ごとに異なるなら � i = P(X � = 1| C = Ci) = PS exp(� α0 + α � 1C1i + α � 2C2i + α � 3C3i +…) α0 + α � 1C1i + α � 2C2i + α � 3C3i +…) 1 + exp(� 40

傾向スコアモデルの当てはめ（SAS） データセット OSSI （n = 8841） ID 00049 00050 00051 00052 00053 00054 00055 00057 00059 00061 00062 … Normothermia HAI Male ASA3 DM Age 1 0 1 0 0 29 1 0 0 0 0 85 1 1 1 1 0 78 1 0 1 0 0 62 1 0 1 0 0 38 1 0 1 0 0 33 1 0 1 0 1 61 1 0 0 0 0 82 1 0 0 0 1 66 0 0 0 0 0 89 0 1 1 0 0 83 41

傾向スコアモデルの当てはめ（SAS） ロジスティックモデル proc logistic data = OSSI desc; model Normothermia = Male ASA3 DM Age; output out = OSSI predicted = PS; run; model: ロジスティックモデルの変数を指定  desc (= descending): 値の大きい方の確率をモデル化 output: 新しいデータセットに新しい変数を追加  predicted = PS: 「PS」という変数に予測値を格納 42

傾向スコアモデルの当てはめ（SAS） ロジスティックモデル 最尤推定値の分析 パラメータ Intercept Male ASA3 DM Age 自由度 1 1 1 1 1 推定値 2.4984 0.2500 0.1890 -0.1508 -0.0793 P(X = 1|C = c) log 1 – P(X = 1|C = c) 標準誤差 0.1629 0.0740 0.1477 0.0952 0.0221 Wald カイ 2 乗値 235.2193 11.4240 1.6381 2.5096 12.8370 Pr > ChiSq <.0001 0.0007 0.2006 0.1132 0.0003 =α �0 + α � 1Male + α � 2ASA3 + α � 3DM + α � 3Age 43

傾向スコアモデルの当てはめ（SAS） データセット OSSI ID 00049 00050 00051 00052 00053 00054 00055 00057 00059 00061 00062 … �i = PS Normothermia HAI Male ASA3 DM Age 1 0 1 0 0 29 1 0 0 0 0 85 1 1 1 1 0 78 1 0 1 0 0 62 1 0 1 0 0 38 1 0 1 0 0 33 1 0 1 0 1 61 1 0 0 0 0 82 1 0 0 0 1 66 0 0 0 0 0 89 0 1 1 0 0 83 exp(� α0 + α � 1Malei + α � 2ASA3i + α � 3DMi + α � 3Agei) α0 + α � 1Malei + α � 2ASA3i + α � 3DMi + α � 3Agei) 1 + exp(� PS 0.92543 0.86111 0.91044 0.90524 0.92035 0.92321 0.89226 0.86393 0.86109 0.85727 0.88996 44

IPW の計算 data OSSI; set OSSI; IPW = Normothermia/PS + (1 - Normothermia)/(1 - PS); run; ID 00049 00050 00051 00052 00053 00054 00055 00057 00059 00061 00062 … Normothermia HAI Male ASA3 DM Age 1 0 1 0 0 29 1 0 0 0 0 85 1 1 1 1 0 78 1 0 1 0 0 62 1 0 1 0 0 38 1 0 1 0 0 33 1 0 1 0 1 61 1 0 0 0 0 82 1 0 0 0 1 66 0 0 0 0 0 89 0 1 1 0 0 83 PS 0.92543 0.86111 0.91044 0.90524 0.92035 0.92321 0.89226 0.86393 0.86109 0.85727 0.88996 IPW 1.0806 1.1613 1.0984 1.1047 1.0865 1.0832 1.1208 1.1575 1.1613 7.0063 9.0878 45

IPWリスク差 線形回帰モデル ただし、交絡変数は説明変数に入れる必要はない 代わりに重み（IPW）で調整  重みは交絡変数から計算されている ロバスト分散（サンドイッチ分散）を使用  一般化推定方程式（generalized estimating equations）ルーチン proc genmod data = hypothermia ; class v001; model HAI = Normothermia /dist = normal link = identity; weight IPW; 代わりに repeated sub=v001; dist = poisson link= log ☞IPWリスク比 run; dist = bin link= logit ☞IPWオッズ比 46

IPW推定値 パラメータ GEE パラメータ推定値の分析 経験的標準誤差推定 標準誤差 95% 信頼限界 推定値 Intercept 0.0270 Normothermia 0.0048 0.0050 0.0054 0.0172 -0.0058 0.0369 0.0154 Z 5.38 0.89 Pr > |Z| <.0001 0.3715 ほかの効果のIPW推定値 95% 信頼区間 効果指標 推定値 リスク差 0.48% -0.58% 1.54% リスク比 1.179 0.803 1.731 オッズ比 1.184 0.798 1.759 47

解釈上の注意 体温というコントロールの難しい処置 正常体温（X = 1）と低体温（X = 0）という介入は一意的か？    正常体温として具体的に何度に設定するのか どのように体温調整するのか … 因果一致性（causal consistency）が成り立たない  X = x ならば Y = Yx が観察されるというのも「仮定」 どこかに「正常体温にできない」または「低体温にならない」層があると その層で Yx を定義できない 集団全体での E(Yx=1) – E(Yx=0) も定義できない 曝露確率の正値性（positivity）が成り立たない 篠崎, 整形外科 2020 48

二重ロバスト推定 49

平均因果効果の推定方法 モデルの要らない交絡調整 層別解析 層別解析の2つの見方 ： アウトカム回帰と傾向スコア 交絡調整にモデルが必要な状況 アウトカム回帰モデル 傾向スコアモデル モデルを使わざるを得ないなら、どちらをつかうか？ 両方つかう二重ロバスト推定 モデルに必要な仮定は？ 50

AIPW推定量 Augmented IPW IPWにアウトカム回帰を組み込んだ局所セミパラメトリック有効推定量 E(Y | X = x, C = c) のモデル P(X = 1 | C = c) のモデル Tsiatis, 2006 いずれかのモデルが正しければ平均因果効果の一致推定量 二重ロバスト性（double robustness） Bang & Robins, Biometrics 2005 Kang & Shafer, Stat Sci 2007 51

傾向スコアモデル IPW推定量 傾向スコアモデル P(X = 1|C = c) log 1 – P(X = 1|C = c) = α0 + α1c1 + α2c2 + α3c3 +… 個人 i に対して予測値 π�i を得る  π�i = {1 + exp(–� α0 – α � 1C1i – α � 2C2i – α � 3C3i –…)}–1 重み付け平均 µˆ 0,IPW = 1 (1 − X i )Yi i 1 − πˆ i 1 (1 − X i ) i 1 − πˆ i ∑ ∑ µˆ1,IPW = 1 X iYi i πˆ i 1 Xi i πˆ i ∑ ∑ 52

アウトカム回帰モデル 回帰標準化推定量 アウトカム回帰モデル E(Y|X = 1, C = c) 1 – E(Y|X = 1, C = c) = β10 + β11c1 + β12c2 + β13c3 +… E(Y|X = 0, C = c) log 1 – E(Y|X = 0, C = c) = β00 + β01c1 + β02c2 + β03c3 +… log � 0,i, m � 1,i) を得る 個人 i に対して予測値 (m   � 0,i = {1 + exp(–β� 00 – β� 01C1i – β� 02C2i – β� 03C3i –…)}–1 m � 1,i = {1 + exp(–β� 10 – β� 11C1i – β� 12C2i – β� 13C3i –…)}–1 m 平均化 �m � 0) µ �0,reg = E( �m � 1) µ �1,reg = E( 53

傾向スコアモデル + アウトカム回帰モデル AIPW推定量 傾向スコアモデル log  P(X = 1|C = c) 1 – P(X = 1|C = c) = α0 + α1c1 + α2c2 + α3c3 +… π�i = {1 + exp(–� α0 – α � 1C1i – α � 2C2i – α � 3C3i –…)}–1 アウトカム回帰モデル E(Y|X = 1, C = c) 1 – E(Y|X = 1, C = c) = β10 + β11c1 + β12c2 + β13c3 +… E(Y|X = 0, C = c) log 1 – E(Y|X = 0, C = c) = β00 + β01c1 + β02c2 + β03c3 +… log   � 0,i = {1 + exp(–β� 00 – β� 01C1i – β� 02C2i – β� 03C3i –…)}–1 m � 1,i = {1 + exp(–β� 10 – β� 11C1i – β� 12C2i – β� 13C3i –…)}–1 m 54

傾向スコアモデル + アウトカム回帰モデル AIPW推定量 IPW平均 + 回帰標準化 – 回帰のIPW平均  Xi  X Yi + mˆ 1,i − i mˆ 1,i  πˆ i  πˆ i  µˆ1,AIPW = n −1 ∑i  µˆ 0,AIPW = n −1 1 − X i  1− Xi Yi + mˆ 0,i − mˆ 0,i  i  1 − πˆ ˆ 1−πi i   ∑ 二重ロバスト性 （注 ： π と mx は確率変数 C の関数） E(X|C = c) = π(c) と正しく特定  E(Xm1/π) = E[E{Xm1/π|C)}] = E{(m1/π)E(X|C)} = E(m1) E(Y|X = x, C = c) = mx(c) と正しく特定  E(XY/π) = E{E(XY|C)/π} = E[E{X E(Y|X,C) |C}/π] = E[{1⋅m1⋅P(X = 1|C) + 0⋅m0⋅P(X = 1|C)}/π] = E{m1⋅E(X|C)/π} = E(Xm1/π) 55

傾向スコアモデル + アウトカム回帰モデル IPW回帰標準化推定量 1. 2. 3. 傾向スコアモデルを推定 � xi, IPW アウトカム回帰モデルを IPW で重み付けて推定 ☞ 予測値 m �m � xi, IPW) ] 全体で重み付けず回帰標準化 ☞ E[ これも二重ロバスト 以下の CAUSALTRT でも実行可 Robins et al., Stat Sci 2007 56

傾向スコアモデル + アウトカム回帰モデル これもぜんぶ二重ロバスト推定量 1. 2. 傾向スコアの推定値 π� i 正準リンク関数 g を用いたアウトカム回帰モデル a. b. c. 3. g{E(Y|X = x, C = c)} = cTβx1 + βx2x/�π – βx3(1 – x)/(1 – π� ) g{E(Y|X = x, C = c)} = cTβx を IPW推定 g{E(Y|X = x, C = c)} = cTβx1 + βx2x�π – βx3(1 – x)(1 – π� )) をIPW推定 �m �m � 1i) – E( � 0i) E( ⇒ β� Bang-Robins ⇒ β� IPW ⇒ β� IPW-PSadj � xi 以外の平均を0にするように β� を選んでいる いずれもAIPW推定量の m なので自動的に二重ロバスト Bang & Robins, Biometrics 2005 Robins et al., Stat Sci 2007 57

実例：NHEFSデータ NHANES I Epidemiologic Follow-Up Study National Health and Nutrition Examination Survey   www.cdc.gov/nchs/nhanes/nhefs/nhefs.htm Hernan & Robins, 2020 ベースラインデータ（1971年） 10年後のフォローアップデータ（1982年） “For the study, medical and behavioral information were collected during an initial physical examination, and again at follow-up interviews approximately one decade later” 58

https://www.icpsr.umich.edu/icpsrweb/ICPSR/index.jsp https://www.hsph.harvard.edu/miguel-hernan/ 59

• https://www.hsph.harvard.edu/miguel-hernan/
• https://www.icpsr.umich.edu/icpsrweb/ICPSR/index.jsp

主な変数コード 変数名 コーディング seqn UNIQUE PERSONAL IDENTIFIER sex 0: MALE 1: FEMALE race 0: WHITE 1: BLACK OR OTHER IN 1971 age AGE IN 1971 education AMOUNT OF EDUCATION BY 1971: 1: 8TH GRADE OR LESS, 2: HS DROPOUT, 3: HS, 4:COLLEGE DROPOUT, 5: COLLEGE OR MORE smokeintensity NUMBER OF CIGARETTES SMOKED PER DAY IN 1971 smokeyrs YEARS OF SMOKING qsmk QUIT SMOKING BETWEEN 1ST QUESTIONNAIRE AND 1982, 1:YES, 0:NO wt82_71 INCREASE IN NUMBER OF CIGARETTES/DAY BETWEEN 1971 and 1982 death DEATH BY 1992, 1:YES, 0:NO exercise IN RECREATION, HOW MUCH EXERCISE? IN 1971, 0:much exercise,1:moderate exercise,2:little or no exercise active IN YOUR USUAL DAY, HOW ACTIVE ARE YOU? IN 1971, 0:very active, 1:moderately active, 2:inactive wt71 WEIGHT IN KILOGRAMS IN 1971 60

解析方針 曝露 X 禁煙 qsmk （1：止めた、0：喫煙） アウトカム Y 10年間の死亡 death （1：死亡、0：生存） 交絡変数 sex、race、age（以上は連続値：2乗項まで調整）、education、smokeintensity、 smokeyrs、exercise、active、wt71 集団全体での効果 E(Yx=1) – E(Yx=0) を推定 IPW アウトカム回帰 二重ロバスト SAS/Proc Causaltrt 61

IPW推定 ods graphics on; proc causaltrt data = nhefs_nmv desc; class exercise active education /desc; psmodel qsmk(ref="0") = sex race age age*age education smokeintensity mokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71 / plots = (psdist weightdist pscovden(effects(age smokeyrs)) ); model death / dist = bin; run; 62

従来のプログラムでは… proc logistic data = nhefs_nmv desc; class exercise active education; model qsmk = sex race age age*age education smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71/ lackfit; output out = nhefs_ps p = ps; run; data nhefs_ps; set nhefs_ps; w = qsmk/ps + (1 - qsmk)/(1 - ps); run; 差を得るため DIST = NORMAL （推定値はIPWRと一致、分布の誤特定下でも分散は妥当） 比が欲しければ DIST = POISSON proc genmod data = nhefs_ps; class seqn; model death = qsmk /dist = normal; weight w; ここでのロバスト分散は傾向スコアモデルの推定を考慮し repeated sub = seqn; ていないため、「近似的」な漸近分散推定量 run; ただし、バイアスは保守的な（過大評価する）方向 63

比が欲しかったら デルタ法 現バージョンの CAUSALTRT ではデフォルトで出力されない 田栗, 2017 “SASによる因果推論：CAUSALTRTプロシジャの紹介” SASではブートストラップ標本を出力可 デフォルトでは変数 _TRTPOM_、_CNTPOM_、_ATE_ に 64

回帰標準化 ods graphics on; proc causaltrt data = nhefs_nmv poutcomemod desc; class exercise active education /desc; psmodel qsmk(ref="0") ; model death = sex race age age*age education smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71 / dist = bin; bootstrap bootci(all) plot seed=1234; run; 65

回帰標準化 66

従来のプログラムでは… data nhefs_nmv; set nhefs_nmv; wt1 = qsmk; wt0 = 1 - qsmk; run; proc genmod data = nhefs_nmv desc; class exercise active education; model death = sex race age age*age education smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71/ dist = bin covb; weight wt1; qsmk = 1 でのみ当てはめ ⇒ E(Yx=1 | C) の予測値を全員に得る output out = pred3 p = pdeath1; run; proc genmod data = pred3 desc; class exercise active education; model death = sex race age age*age education smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71/ dist = bin covb; weight wt0; qsmk = 0 でのみ当てはめ ⇒ E(Yx=0 | C) の予測値を全員に得る output out = pred3 p = pdeath0; run; proc means data = pred3; 分散共分散行列によるデルタ法 var pdeath1 pdeath0; またはブートストラップ法で分散推定 run; 67

AIPW（二重ロバスト）推定量 ods graphics on; proc causaltrt data = nhefs_nmv ppsmodel poutcomemod desc; class exercise active education /desc; psmodel qsmk(ref="0") = sex race age age*age education smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71; model death = sex race age age*age education smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs exercise active wt71 wt71*wt71/ dist = bin; bootstrap bootci(all) plot seed=1234; run; 68

AIPW（二重ロバスト）推定量 69

モデルに含める変数が異なると 「二重ロバスト性」は担保されない 十分な C を条件づけた回帰 十分な C を条件づけた傾向スコア のモデル特定に対する性質 「モデルに含める変数」に対するロバスト性ではない Keil et al. AJE 2018 アドホックな組み合わせ（やらないほうがいいです） 交絡変数の一部を傾向スコアモデルでバランスさせておいて バランスされていない一部を回帰モデルで調整 Shinozaki & Nojima, Epidemiology 2019 70

曝露群での効果（ATT）の推定 average treatment effect on the treated 重み付け（ATT weight, SMR weight, inverse odds weight, …） Xi = 0: Wi = π� i 1 – π� i Xi = 1: Wi = 1 IPW× π� i Sato & Matsuyama, Epidemiology 2003 回帰標準化 n –1 � � 0|X = 1] = n1 ∑ Xi m � 0,i E[m i=1 n –1 � � 1|X = 1] = n1 ∑ Xi m � 1,i E[m i=1 � いずれも E[Y|X = 1] E[Yx=0|X = 1] の推定だけ考えればよい 71

ATTの二重ロバスト推定 � x=0|X = 1] E[Y �m � ATT-weighted[m � ATT-weighted[Y|X = 0] + E[ � 0|X = 1] – E � 0|X = 0] =E Shinozaki & Matsuyama, Epidemiology 2015 Mao et al.（SMMR 2019）では ATT weight は二重ロバストにできないような 書きぶりだがそんなことはない CAUSALTRT プロシジャでは非実装 72

まとめ：交絡調整方法のチョイス 交絡変数 C があるなら、層別解析 交絡変数が C としてすべて測定されている仮定 層別解析ができない場合 層内に十分な人数がいない、個人ごとに C の値が異なる 統計モデル：層別解析の近似手段 回帰モデル 傾向スコアモデル    E(Y|X, C)  P(X = 1|C) 個人ごとに異なる値が得られるが... あくまで交絡変数 C の層ごとの値が欲しい どちらかを使う 両方つかう ⇒ 二重ロバスト推定 73

まとめ：各手法に必要な条件 必要な仮定 1. C で層別できたら X がランダムに決まる 層別解析に共通 2a. E(Y|X, C) を正しく特定 回帰モデル 2b. P(X = 1|C) を正しく特定 傾向スコアモデル 2c. E(Y|X, C) または P(X = 1|C) （または両方）を正しく特定 二重ロバスト推定量 74

Take-home message（再掲） 交絡調整 = 層別解析 交絡変数でサブグループに分けて（= 層別）比較 交絡調整における2種類の統計モデル アウトカム回帰モデル 傾向スコアモデル モデルをつかう目的 層別できないほどの交絡変数があるとき 「もし層別できた場合」の結果を近似   アウトカム回帰 傾向スコア 各モデルの推定自体が目的ではない 75

文献 Bang H, Robins JM. Doubly robust estimation in missing data and causal inference models. Biometrics. 2005;61:962-73. Greenland S. Estimating standardized parameters from generalized linear models. Stat Med. 1991;10:10691074. 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, Robins JM. Causal Inference: What If. Chapman & Hall/CRC, 2020 Kang JDY, Shafer JL, Demystifying double robustness: a comparison of alternative strategies for estimating a population mean from incomplete data. Stat Sci 2007;22:523-539. Keil AP, et al. Resolving an apparent paradox in doubly robust estimators. Am J Epidemiol. 2018;187:891892. Mao H, Li L, Greene T. Propensity score weighting analysis and treatment effect discovery. Stat Methods Med Res. 2019;28:2439-2454. 76

文献 Robins JM, et al. Comment: Performance of double-robust estimators when “inverse probability” weights are highly variable. Stat Sci. 2007;22: 544-559. 佐藤俊哉. 疫学研究における交絡と効果の修飾. 統計数理 1994;42:83–101 篠崎．傾向スコア解析の考え方．整形外科 2020; 71(6):571-576. Shinozaki T, Matsuyama Y. Doubly robust estimation of standardized risk difference and ratio in the exposed population. Epidemiology. 2015;26(6):873-7. Shinozaki Y, Nojima M. Misuse of regression adjustment for additional confounders following insufficient propensity score balancing. Epidemiology. 2019;30:541-548. Shrier I, Pang M. Confounding, effect modification, and the odds ratio: common misinterpretations. J Clin Epidemiol. 2015;68:470-474. Tsiatis AA. Semiparametric Theory and Missing Data. Springer, 2006. Vansteelandt S, Keiding N. Invited commentary: G-computation--lost in translation? Am J Epidemiol. 2011;173:739-42. Yamada et al. Association between normothermia at the end of surgery and postoperative complications following orthopedic surgery. Clin Infect Dis 2020;70:474-482 77