---
title: 因果推論II：モデルを用いた効果推定
tags:  #因果推論  
author: [Tomohiro Shinozaki](https://convert.docswell.com/user/she-knows-a-key)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/LE3NKRDQE5.jpg?width=480
description: 統計数理研究所 医療健康データ科学研究センター 公開講座「疫学・公衆衛生統計」｜12月10日
published: September 08, 22
canonical: https://convert.docswell.com/s/she-knows-a-key/52639Z-2022-09-08-225555
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/LE3NKRDQE5.jpg)

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


# Page. 2

![Page Image](https://bcdn.docswell.com/page/8EDW3YPW7G.jpg)

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


# Page. 3

![Page Image](https://bcdn.docswell.com/page/V7P946NXJ8.jpg)

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


# Page. 4

![Page Image](https://bcdn.docswell.com/page/2JVMX3K3JQ.jpg)

復習：記法
観察変数
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


# Page. 5

![Page Image](https://bcdn.docswell.com/page/5EGGV3QYEL.jpg)

関心のある「平均因果効果」
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 &amp; Robins, 2020
5


# Page. 6

![Page Image](https://bcdn.docswell.com/page/4JQX65K6EP.jpg)

層別できるなら
交絡変数が 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


# Page. 7

![Page Image](https://bcdn.docswell.com/page/K74649PME1.jpg)

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


# Page. 8

![Page Image](https://bcdn.docswell.com/page/LJ1V4K9YEG.jpg)

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


# Page. 9

![Page Image](https://bcdn.docswell.com/page/GJWLXDQ1E2.jpg)

回帰モデル
モデルで 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


# Page. 10

![Page Image](https://bcdn.docswell.com/page/4EZM6GWXJ3.jpg)

回帰モデルの当てはめ
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


# Page. 11

![Page Image](https://bcdn.docswell.com/page/Y76N2YGP7V.jpg)

回帰モデルの予測値
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


# Page. 12

![Page Image](https://bcdn.docswell.com/page/G7552GVP74.jpg)

回帰モデルによる標準化
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


# Page. 13

![Page Image](https://bcdn.docswell.com/page/9J2N4Y3ZER.jpg)

こう考えても同じこと
個人 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


# Page. 14

![Page Image](https://bcdn.docswell.com/page/DEYPMG1NEM.jpg)

個人単位で見た標準化
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 &amp; Keiding, AJE 2011
層別しきれないほどの交絡変数 C にも適用できる
14


# Page. 15

![Page Image](https://bcdn.docswell.com/page/VJN5WG2V78.jpg)

周辺効果、条件付き効果（リスク差）
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


# Page. 16

![Page Image](https://bcdn.docswell.com/page/YE9YXW1VJ3.jpg)

周辺効果、条件付き効果（オッズ比）
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


# Page. 17

![Page Image](https://bcdn.docswell.com/page/GE862N447D.jpg)

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


# Page. 18

![Page Image](https://bcdn.docswell.com/page/LEL92P6VJR.jpg)

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


# Page. 19

![Page Image](https://bcdn.docswell.com/page/4JMZ845NEW.jpg)

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


# Page. 20

![Page Image](https://bcdn.docswell.com/page/PJR65W3479.jpg)

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


# Page. 21

![Page Image](https://bcdn.docswell.com/page/PEXMKP4Z7X.jpg)

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


# Page. 22

![Page Image](https://bcdn.docswell.com/page/3EKG5DZVJD.jpg)

傾向スコアモデル
モデルで 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


# Page. 23

![Page Image](https://bcdn.docswell.com/page/L73NKRXQ75.jpg)

傾向スコアモデルの当てはめ
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


# Page. 24

![Page Image](https://bcdn.docswell.com/page/87DW3Y2WJG.jpg)

傾向スコアモデルの予測値
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


# Page. 25

![Page Image](https://bcdn.docswell.com/page/VJP9461XE8.jpg)

傾向スコアモデルによる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


# Page. 26

![Page Image](https://bcdn.docswell.com/page/2EVMX3Y3EQ.jpg)

傾向スコアモデルによる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


# Page. 27

![Page Image](https://bcdn.docswell.com/page/57GGV39Y7L.jpg)

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


# Page. 28

![Page Image](https://bcdn.docswell.com/page/4EQX65W67P.jpg)

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


# Page. 29

![Page Image](https://bcdn.docswell.com/page/KJ4649KM71.jpg)

正常体温 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


# Page. 30

![Page Image](https://bcdn.docswell.com/page/LE1V4KXY7G.jpg)

手術前変数
Table 1
Female sex — no. (%)
Age (yr)*
BMI (kg/m2)*†
BMI — no. (%)†
Present smoker — no. (%)
Disease
BMI&lt;25
25&lt;BMI&lt;30
30&lt;BMI
Arthritis
Fractures &amp; 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
&lt;.0001
&lt;.0001
&lt;.0001
&lt;.0001
0.52
&lt;.0001
群間で異なる分布
これらの変数によって感染率が異なる可能性が高い
Yamada et al. CID 2020
30


# Page. 31

![Page Image](https://bcdn.docswell.com/page/GEWLXDN172.jpg)

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


# Page. 32

![Page Image](https://bcdn.docswell.com/page/47ZM6G3XE3.jpg)

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


# Page. 33

![Page Image](https://bcdn.docswell.com/page/YJ6N2YXPJV.jpg)

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


# Page. 34

![Page Image](https://bcdn.docswell.com/page/GJ552GXPJ4.jpg)

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


# Page. 35

![Page Image](https://bcdn.docswell.com/page/9E2N4YXZ7R.jpg)

どうやって予測値を得る？
良い方法があります
データセット 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


# Page. 36

![Page Image](https://bcdn.docswell.com/page/D7YPMG8N7M.jpg)

トリック：データを複製
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


# Page. 37

![Page Image](https://bcdn.docswell.com/page/VEN5WGQVJ8.jpg)

回帰モデルの当てはめ
ロジスティックモデル
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


# Page. 38

![Page Image](https://bcdn.docswell.com/page/Y79YXWDVE3.jpg)

� xi ゲット
モデルの予測値 P
回帰モデルの推定には1行目だけ寄与（n = 8841）
2 &amp; 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


# Page. 39

![Page Image](https://bcdn.docswell.com/page/G7862N34JD.jpg)

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


# Page. 40

![Page Image](https://bcdn.docswell.com/page/L7L92PVVER.jpg)

傾向スコアモデル
モデルで 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


# Page. 41

![Page Image](https://bcdn.docswell.com/page/4EMZ84GN7W.jpg)

傾向スコアモデルの当てはめ（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


# Page. 42

![Page Image](https://bcdn.docswell.com/page/PER65WZ4J9.jpg)

傾向スコアモデルの当てはめ（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


# Page. 43

![Page Image](https://bcdn.docswell.com/page/P7XMKP9ZJX.jpg)

傾向スコアモデルの当てはめ（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 &gt; ChiSq
&lt;.0001
0.0007
0.2006
0.1132
0.0003
=α
�0 + α
� 1Male + α
� 2ASA3 + α
� 3DM + α
� 3Age
43


# Page. 44

![Page Image](https://bcdn.docswell.com/page/37KG5DPVED.jpg)

傾向スコアモデルの当てはめ（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


# Page. 45

![Page Image](https://bcdn.docswell.com/page/LJ3NKR4QJ5.jpg)

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


# Page. 46

![Page Image](https://bcdn.docswell.com/page/8JDW3YQWEG.jpg)

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


# Page. 47

![Page Image](https://bcdn.docswell.com/page/VEP946LX78.jpg)

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 &gt; |Z|
&lt;.0001
0.3715
ほかの効果のIPW推定値
95% 信頼区間
効果指標
推定値
リスク差
0.48%
-0.58%
1.54%
リスク比
1.179
0.803
1.731
オッズ比
1.184
0.798
1.759
47


# Page. 48

![Page Image](https://bcdn.docswell.com/page/27VMX3Q37Q.jpg)

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


# Page. 49

![Page Image](https://bcdn.docswell.com/page/5JGGV3WYJL.jpg)

二重ロバスト推定
49


# Page. 50

![Page Image](https://bcdn.docswell.com/page/47QX6536JP.jpg)

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


# Page. 51

![Page Image](https://bcdn.docswell.com/page/KE46491MJ1.jpg)

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


# Page. 52

![Page Image](https://bcdn.docswell.com/page/L71V4KGYJG.jpg)

傾向スコアモデル
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


# Page. 53

![Page Image](https://bcdn.docswell.com/page/G7WLXDK1J2.jpg)

アウトカム回帰モデル
回帰標準化推定量
アウトカム回帰モデル
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


# Page. 54

![Page Image](https://bcdn.docswell.com/page/4JZM6GZX73.jpg)

傾向スコアモデル
+
アウトカム回帰モデル
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


# Page. 55

![Page Image](https://bcdn.docswell.com/page/YE6N2YZPEV.jpg)

傾向スコアモデル
+
アウトカム回帰モデル
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


# Page. 56

![Page Image](https://bcdn.docswell.com/page/GE552GWPE4.jpg)

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


# Page. 57

![Page Image](https://bcdn.docswell.com/page/972N4YQZJR.jpg)

傾向スコアモデル
+
アウトカム回帰モデル
これもぜんぶ二重ロバスト推定量
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 &amp; Robins, Biometrics 2005
Robins et al., Stat Sci 2007
57


# Page. 58

![Page Image](https://bcdn.docswell.com/page/DJYPMGWNJM.jpg)

実例：NHEFSデータ
NHANES I Epidemiologic Follow-Up Study
National Health and Nutrition Examination Survey


www.cdc.gov/nchs/nhanes/nhefs/nhefs.htm
Hernan &amp; 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


# Page. 59

![Page Image](https://bcdn.docswell.com/page/V7N5WG9VE8.jpg)

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


# Page. 60

![Page Image](https://bcdn.docswell.com/page/YJ9YXW2V73.jpg)

主な変数コード
変数名
コーディング
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


# Page. 61

![Page Image](https://bcdn.docswell.com/page/GJ862N54ED.jpg)

解析方針
曝露 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


# Page. 62

![Page Image](https://bcdn.docswell.com/page/LJL92PYV7R.jpg)

IPW推定
ods graphics on;
proc causaltrt data = nhefs_nmv desc;
class exercise active education /desc;
psmodel qsmk(ref=&quot;0&quot;) = 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


# Page. 63

![Page Image](https://bcdn.docswell.com/page/47MZ84NNJW.jpg)

従来のプログラムでは…
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


# Page. 64

![Page Image](https://bcdn.docswell.com/page/P7R65WD4E9.jpg)

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


# Page. 65

![Page Image](https://bcdn.docswell.com/page/PJXMKP1ZEX.jpg)

回帰標準化
ods graphics on;
proc causaltrt data = nhefs_nmv poutcomemod desc;
class exercise active education /desc;
psmodel qsmk(ref=&quot;0&quot;) ;
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


# Page. 66

![Page Image](https://bcdn.docswell.com/page/3JKG5D2V7D.jpg)

回帰標準化
66


# Page. 67

![Page Image](https://bcdn.docswell.com/page/LE3NKRYQE5.jpg)

従来のプログラムでは…
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


# Page. 68

![Page Image](https://bcdn.docswell.com/page/8EDW3Y5W7G.jpg)

AIPW（二重ロバスト）推定量
ods graphics on;
proc causaltrt data = nhefs_nmv ppsmodel poutcomemod desc;
class exercise active education /desc;
psmodel qsmk(ref=&quot;0&quot;) = 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


# Page. 69

![Page Image](https://bcdn.docswell.com/page/V7P946GXJ8.jpg)

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


# Page. 70

![Page Image](https://bcdn.docswell.com/page/2JVMX363JQ.jpg)

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


# Page. 71

![Page Image](https://bcdn.docswell.com/page/5EGGV3NYEL.jpg)

曝露群での効果（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 &amp; 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


# Page. 72

![Page Image](https://bcdn.docswell.com/page/4JQX65Q6EP.jpg)

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 &amp; Matsuyama, Epidemiology 2015
Mao et al.（SMMR 2019）では ATT weight は二重ロバストにできないような
書きぶりだがそんなことはない
CAUSALTRT プロシジャでは非実装
72


# Page. 73

![Page Image](https://bcdn.docswell.com/page/K74649NME1.jpg)

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


# Page. 74

![Page Image](https://bcdn.docswell.com/page/LJ1V4K5YEG.jpg)

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


# Page. 75

![Page Image](https://bcdn.docswell.com/page/GJWLXD51E2.jpg)

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


# Page. 76

![Page Image](https://bcdn.docswell.com/page/4EZM6GNXJ3.jpg)

文献
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 &amp; 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


# Page. 77

![Page Image](https://bcdn.docswell.com/page/Y76N2Y8P7V.jpg)

文献
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


