因果推論のための分析手法は様々ありますが、回帰モデルを使った主なアプローチのRでの実装方法とその推定結果の比較をします。
モチベーション的な部分は以下をご参照ください。
シミュレーションデータを使って、各手法がどのような（主にモデリングに関する）仮定に基づいているのか、それが結果の違いにどのように影響しているのかをみていきます。
なお、Rマークダウンで書いたものをコピペしただけなのであまり見た目は美しくないですが、予め ご了承ください。
そのうち書籍化予定ですので、そのときにはもっと見やすく＆詳しく説明します・・・
library(tidyverse)
library(broom)
library(gtsummary)
library(MatchIt)
library(geepack)
library(ltmle)シナリオ
「禁煙がその後の体重に与える影響を知りたい」というモチベーションのもと、もともと喫煙者の集団から「最近禁煙をはじめた人」「喫煙を継続している人」でその後の体重を比較した。交絡因子は年齢・性別・飲酒の有無のみとする。
擬似データの概要:
変数
-
Y: アウトカム（連続変数）
- フォローアップ終了時の体重
-
A: 曝露（二値変数）
- ベースラインでの禁煙
-
C1: 交絡因子（連続変数）
- 年齢 (中心化済み; 平均 = 0, SD = 5)
-
C2 and C3: 交絡因子（二値変数）
- 男性 (C2=1; サンプルの40%) and 飲酒 (C3=1; 30%)
真のアウトカムモデル:
\[ E[Y|A,C1,C2,C3]=\beta_0+\beta_1A+\beta_2C_1+\beta_3C_1^2+\beta_4C2+\beta_5C3+\beta_6A*C_2 \]
真の曝露モデル:
\[ logit(Pr[A=1|C_1,C_2,C_3]) = \gamma_0 + \gamma_1C_1 + \gamma_2C_1^2 +\gamma_3C_2 +\gamma_4C_3 +\gamma_5C_2*C_3 \]
set.seed(0)
n.obs = 10000 #サンプルサイズの設定
#---- 真のアウトカムモデルのパラメーター ----
b0 = 60
b1 = 5
b2 = -0.3
b3 = -0.1
b4 = 8
b5 = 3
b6 = 2
#---- 真の曝露オッズモデルのパラメーター ----
g0 = log(0.20/(1-0.20))
g1 = log(1.01)
g2 = log(1.005)
g3 = log(0.6)
g4 = log(0.5)
g5 = log(0.8)
#アウトカム期待値を計算する関数
## 上で指定したモデルパラメーターを使用
mean_out <- function(C1, C2, C3, exposure){
b0 + b1*exposure + b2*C1 + b3*I(C1^2) + b4*C2 + b5*C3 + b6*exposure*C2
}
#曝露確率を計算する関数
## 上で指定したモデルパラメーターを使用
prob_exp <- function(C1, C2, C3){
exp(g0 + g1*C1 + g2*I(C1^2) + g3*C2 + g4*C3 + g5*C2*C3)/(1 + exp(g0 + g1*C1 + g2*I(C1^2) + g3*C2 + g4*C3 + g5*C2*C3))
}
#データをシミュレート
df.sim <- tibble("ID" = seq(from = 1, to = n.obs, by = 1),
"C1" = rnorm(n = n.obs, mean = 0, sd = 5),
"C2" = rbinom(n = n.obs, size = 1, p = 0.4),
"C3" = rbinom(n = n.obs, size = 1, p = 0.3),
"Pexposure" = prob_exp(C1, C2, C3),
"Exposure" = rbinom(n = n.obs, size = 1,
prob = Pexposure),
"Outcome" = as.numeric(mean_out(C1,C2,C3, Exposure))) #
#データの分布を確認
df.sim %>%
mutate_at(c("C2","C3","Exposure"),factor) %>%
select(Outcome,Exposure,C1,C2,C3) %>%
tbl_summary(
by = Exposure,
statistic = list(c(Outcome, C1) ~ "{mean} ({sd})")
) %>%
add_p() %>%
add_overall() %>%
modify_spanning_header(c("stat_1","stat_2") ~ "**Exposure**")| Characteristic | Overall, N = 10,0001 | Exposure | p-value2 | |
|---|---|---|---|---|
| 0, N = 8,3551 | 1, N = 1,6451 | |||
| Outcome | 62.5 (5.8) | 61.9 (5.4) | 65.3 (6.6) | <0.001 |
| C1 | 0.1 (5.0) | 0.0 (4.9) | 0.4 (5.5) | 0.005 |
| C2 | <0.001 | |||
| 0 | 6,042 (60%) | 4,885 (58%) | 1,157 (70%) | |
| 1 | 3,958 (40%) | 3,470 (42%) | 488 (30%) | |
| C3 | <0.001 | |||
| 0 | 6,992 (70%) | 5,635 (67%) | 1,357 (82%) | |
| 1 | 3,008 (30%) | 2,720 (33%) | 288 (18%) | |
|
1 Mean (SD); n (%) 2 Wilcoxon rank sum test; Pearson's Chi-squared test |
||||
曝露因子の有無でC1~C3の分布が結構違う。交絡が生じてそうですね。
真の治療効果:
-
女性における禁煙が体重に与える効果(Conditional Effect): \(\beta_1 = 5\) (kg)
-
男性における禁煙が体重に与える効果(Conditional Effect): \(\beta_1+\beta_6 = 5 + 2 = 7\) (kg)
-
（男女合わせた）母集団全体における効果(Marginal Effect): \(5*0.6+7*0.4=5.8\) (kg)
- 女性60%, 男性40%の集団
重回帰分析
アウトカムに対する重回帰モデルを用いてConditional Effectを推定する
正しく設定されたモデル
df.sim %>%
lm(Outcome ~ Exposure*C2 + C1 + I(C1^2) + C3, data = .) %>%
tidy(conf.int = TRUE)## # A tibble: 7 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 60.0 1.22e-14 4.92e15 0 60.0 60.0
## 2 Exposure 5.00 2.32e-14 2.16e14 0 5.00 5.00
## 3 C2 8. 1.56e-14 5.12e14 0 8. 8.
## 4 C1 -0.3 1.41e-15 -2.12e14 0 -0.3 -0.300
## 5 I(C1^2) -0.1 2.06e-16 -4.86e14 0 -0.1 -0.100
## 6 C3 3. 1.55e-14 1.94e14 0 3.00 3.
## 7 Exposure:C2 2.00 4.11e-14 4.87e13 0 2.00 2.Exposureの推定値=5, 掛け算項の推定値=2となり、真のモデルと一致。
調整なしモデル
df.sim %>%
lm(Outcome ~ Exposure,
data = .) %>%
tidy(conf.int = TRUE)## # A tibble: 2 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 61.9 0.0616 1006. 0. 61.8 62.1
## 2 Exposure 3.32 0.152 21.9 6.62e-104 3.03 3.62Exposureに対する推定値3.32となり、真の効果(5, 7, 5.8)から大きくずれている。C1~C3を調整しない交絡バイアス。
調整あり&モデルの誤設定1
df.sim %>%
lm(Outcome ~ Exposure + C1 + C2 + C3,
data = .) %>%
tidy(conf.int = TRUE)## # A tibble: 5 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 57.6 0.0538 1071. 0. 57.5 57.7
## 2 Exposure 4.88 0.0939 52.0 0. 4.70 5.07
## 3 C1 -0.307 0.00691 -44.5 0. -0.321 -0.294
## 4 C2 8.16 0.0706 116. 0. 8.02 8.30
## 5 C3 2.96 0.0755 39.1 1.02e-311 2.81 3.10交絡バイアスはないが、アウトカムモデルが誤設定
-
C1のFunctional formが線形
- 本当は二次関数形
-
Exposure*C2がない
- 男女で効果が違うのにConstant effectを仮定
Exposureの推定値=4.882となり、真の効果(5, 7, 5.8)からずれている。モデルの誤設定によるバイアス。
調整あり＆モデルの誤設定2
df.sim %>%
lm(Outcome ~ Exposure + C1 + I(C1^2) + C2 + C3,
data = .) %>%
tidy(conf.int = TRUE)## # A tibble: 6 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 59.9 0.00583 10264. 0 59.9 59.9
## 2 Exposure 5.63 0.00939 599. 0 5.61 5.65
## 3 C1 -0.300 0.000689 -435. 0 -0.301 -0.299
## 4 I(C1^2) -0.100 0.000100 -998. 0 -0.100 -0.0999
## 5 C2 8.29 0.00704 1178. 0 8.28 8.30
## 6 C3 3.00 0.00753 398. 0 2.98 3.01交絡バイアスはないが、アウトカムモデルが誤設定
- Exposure*C2がない
Exposureの推定値=5.63となり、真の効果(5, 7, 5.8)からずれている。モデルの誤設定によるバイアス。
標準化(G-formula/G-computation)
正しく設定されたモデル
# オリジナルデータのコピーを作成
df.sim.a1 <- df.sim %>%
mutate(Outcome = NA,
Exposure = 1) #全員にExposure = 1を割付
df.sim.a0 <- df.sim %>%
mutate(Outcome = NA,
Exposure = 0) #全員にExposure = 0を割付
df.sim.combined <-
bind_rows(df.sim.a1,df.sim.a0)
# オリジナルデータにアウトカムモデルをフィット
## 重回帰分析と同じ
gcomp.fit <- df.sim %>%
lm(Outcome ~ Exposure*C2 + C1 + I(C1^2) + C3, data = .)
# コピーデータを使ってアウトカムの予測
df.sim.combined$pred <- predict(gcomp.fit, newdata = df.sim.combined)
# ATE推定値: A=1の行の予測値平均とA=0の行の予測値平均の差
df.sim.combined %>%
group_by(Exposure) %>%
summarise(
mean.Y = mean(pred)
) %>%
pivot_wider(
names_from = Exposure,
names_glue = "mean.Y.{Exposure}",
values_from = mean.Y
) %>%
mutate(
ATE = mean.Y.1-mean.Y.0
)## # A tibble: 1 x 3
## mean.Y.0 mean.Y.1 ATE
## <dbl> <dbl> <dbl>
## 1 61.6 67.4 5.79ATEの推定値=5.79なので真のMarginal Effect = 5.8とほぼ一致。
モデルの誤設定
-
C1のFunctional formが線形
- 本当は二次関数形
-
Exposure*C2がない
- 男女で効果が違うのにConstant effectを仮定
gcomp.fit <- df.sim %>%
lm(Outcome ~ Exposure + C1 + C2 + C3, data = .)
df.sim.combined$pred <- predict(gcomp.fit, newdata = df.sim.combined)
df.sim.combined %>%
group_by(Exposure) %>%
summarise(
mean.Y = mean(pred)
) %>%
pivot_wider(
names_from = Exposure,
names_glue = "mean.Y.{Exposure}",
values_from = mean.Y
) %>%
mutate(
ATE = mean.Y.1-mean.Y.0
)## # A tibble: 1 x 3
## mean.Y.0 mean.Y.1 ATE
## <dbl> <dbl> <dbl>
## 1 61.7 66.6 4.88ATEの推定値=4.88となり真のMarginal Effect = 5.8と不一致。モデルの誤設定によるバイアス。
そもそもConstant effectを仮定しているので、標準化をしても同じモデルを使った重回帰分析と推定結果は同じ。
傾向スコア(PS)
# 曝露モデルのフィット
PS.fit.correct <-
df.sim %>%
glm(Exposure ~ C1 + I(C1^2) + C2*C3, family = "binomial", data=.)
PS.fit.misspecified <-
df.sim %>%
glm(Exposure ~ C1 + C2 + C3, family = "binomial", data=.)
# 傾向スコアの計算
df.sim$PS.correct <-
predict(PS.fit.correct, type = "response")
df.sim$PS.misspecified <-
predict(PS.fit.misspecified, type = "response")層化
傾向スコアの十分位点を使って10の層にサンプルを分割。層ごとにアウトカムの群間比較(A=1vs.A=0)。
df.sim %>%
mutate(PS.correct.strata = gtools::quantcut(PS.correct, 10)) %>% #10分位点
group_by(Exposure,PS.correct.strata) %>% #曝露値・PS層ごとにグループ
summarise(mean.Y = mean(Outcome)) %>% #各グループごとのアウトカム平均を計算
pivot_wider( # wideデータに変換
names_from = Exposure,
names_glue = "mean.Y.{Exposure}",
values_from = mean.Y
) %>%
mutate(
ATE.strata = mean.Y.1-mean.Y.0, # 層ごとのATEを計算
ATE = mean(ATE.strata) # 層ごとのATEの平均
)## # A tibble: 10 x 5
## PS.correct.strata mean.Y.0 mean.Y.1 ATE.strata ATE
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 [0.0554,0.0729] 69.7 76.8 7.05 5.64
## 2 (0.0729,0.106] 62.9 68.3 5.38 5.64
## 3 (0.106,0.129] 62.9 68.2 5.32 5.64
## 4 (0.129,0.135] 67.0 74.3 7.26 5.64
## 5 (0.135,0.155] 65.0 71.2 6.28 5.64
## 6 (0.155,0.198] 60.3 66.3 6.04 5.64
## 7 (0.198,0.202] 59.9 65.0 5.07 5.64
## 8 (0.202,0.212] 59.1 64.1 4.99 5.64
## 9 (0.212,0.235] 57.4 62.4 5.02 5.64
## 10 (0.235,0.608] 51.9 55.9 4.02 5.64- 層ごとのATE推定値が大きく異なる。傾向スコア（層）による効果修飾があったと考えるのがリーズナブル
- PS層のなかで調整しきれない交絡(residual confounding)がある可能性
回帰モデルでの調整
アウトカムに対する重回帰分析でC1~C3を調整する代わりに傾向スコアを調整。
df.sim %>%
lm(Outcome ~ Exposure + PS.correct, data=.) %>%
tidy()## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 73.7 0.0809 912. 0
## 2 Exposure 5.63 0.0819 68.8 0
## 3 PS.correct -74.0 0.464 -159. 0Exposureの推定値=5.63で真のMarginal Effect = 5.8と微妙に違う。モデルの誤設定によるバイアス。
モデルの仮定：
-
PSのFunctional formが線形
-
Exposure*PSがない
- 「PSによる効果修飾なし」の仮定
g-formulaとの合わせ技
g-formulaと合わせることでPSによる効果修飾をモデルしつつ、Marginal effectを推定可能。
df.sim.a1 <- df.sim %>%
mutate(Outcome = NA,
Exposure = 1) #全員にExposure = 1を割付
df.sim.a0 <- df.sim %>%
mutate(Outcome = NA,
Exposure = 0) #全員にExposure = 0を割付
df.sim.combined <-
bind_rows(df.sim.a1,df.sim.a0)
gcomp.PS.fit <- df.sim %>%
lm(Outcome ~ Exposure*PS.correct, data = .) # PSによる交互作用をモデル
df.sim.combined$pred <- predict(gcomp.PS.fit, newdata = df.sim.combined)
df.sim.combined %>%
group_by(Exposure) %>%
summarise(
mean.Y = mean(pred)
) %>%
pivot_wider(
names_from = Exposure,
names_glue = "mean.Y.{Exposure}",
values_from = mean.Y
) %>%
mutate(
ATE = mean.Y.1-mean.Y.0
)## # A tibble: 1 x 3
## mean.Y.0 mean.Y.1 ATE
## <dbl> <dbl> <dbl>
## 1 61.6 67.5 5.90ATE推定値 = 5.9で真のMarginal Effect = 5.8に近づいた。一致しないのはPSのFunctional formの仮定によるものだと考えられる。
マッチング
fit.match <- matchit(Exposure ~ C1 + I(C1^2) + C2*C3, method = "nearest", data=df.sim)
df.matched <- match.data(fit.match)
plot(fit.match, type = "hist")
df.matched %>%
mutate_at(c("C2","C3","Exposure"),factor) %>%
select(Outcome,Exposure,C1,C2,C3) %>%
tbl_summary(
by = Exposure,
statistic = list(c(Outcome, C1) ~ "{mean} ({sd})")
) %>%
add_p() %>%
add_overall() %>%
modify_spanning_header(c("stat_1","stat_2") ~ "**Exposure**")| Characteristic | Overall, N = 3,2901 | Exposure | p-value2 | |
|---|---|---|---|---|
| 0, N = 1,6451 | 1, N = 1,6451 | |||
| Outcome | 62 (7) | 60 (6) | 65 (7) | <0.001 |
| C1 | 0.4 (5.6) | 0.3 (5.6) | 0.4 (5.5) | 0.7 |
| C2 | 0.7 | |||
| 0 | 2,305 (70%) | 1,148 (70%) | 1,157 (70%) | |
| 1 | 985 (30%) | 497 (30%) | 488 (30%) | |
| C3 | 0.7 | |||
| 0 | 2,721 (83%) | 1,364 (83%) | 1,357 (82%) | |
| 1 | 569 (17%) | 281 (17%) | 288 (18%) | |
|
1 Mean (SD); n (%) 2 Wilcoxon rank sum test; Pearson's Chi-squared test |
||||
df.matched %>%
lm(Outcome ~ Exposure, data=.) %>%
tidy()## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 59.7 0.154 388. 0.
## 2 Exposure 5.54 0.218 25.5 1.11e-130マッチング後はExposure群間でPSおよびC1~C3の分布が似通っている。
ATE推定値は5.54となり、真のMarginal Effect5.8とは一致しない。
マッチング前のサンプルは全体で40%がC2=1だったが、マッチング後のサンプルでは30%しかC2=1がいない。C2は効果修飾因子なので、その分布が異なればATE推定値も異なる。
IPTW
Unstabilized weightsを使用。
\[ USTW = \frac{1}{Pr[A=a|C1,C2,C3]} \]
正しく設定されたモデル
正しい曝露モデルから推定された傾向スコアを使用。Marginal Effectを推定。理論上は正しいアウトカムモデルを使用したg-formulaと一致するはず。
df.sim <- df.sim %>%
mutate(ustw = ifelse(Exposure == 1,
1/PS.correct, #Exposure =1の人は分母がPS
1/(1-PS.correct))) #Exposure = 1の人は分母が1-PS
# fit MSM
df.sim %>%
geeglm(Outcome ~ Exposure, weights = ustw, data=.,id = ID) %>% #重み付け使用のためGEEを用いてSE補正
tidy()## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 61.6 0.0636 935874. 0
## 2 Exposure 5.82 0.189 943. 0真のMarginal Effect(5.8)と（ほぼ）一致。
モデルの誤設定
誤設定された曝露モデルから推定された傾向スコアを使用。
df.sim <- df.sim %>%
mutate(ustw.misspecified = ifelse(Exposure == 1,
1/PS.misspecified, #Exposure =1の人は分母がPS
1/(1-PS.misspecified))) #Exposure = 1の人は分母が1-PS
# fit MSM
df.sim %>%
geeglm(Outcome ~ Exposure, weights = ustw.misspecified, data=.,id = ID) %>% #重み付け使用のためGEEを用いてSE補正
tidy()## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 61.7 0.0597 1069567. 0
## 2 Exposure 5.12 0.193 702. 0モデルの誤設定によるバイアスが生じている。
Doubly-robust Estimation
Targeted Maximum Likelihood Estimationを使用。Marginal effectを推定。アウトカムモデルと曝露モデルのどちらか少なくとも一方が正しく設定されていれば良い。
Correct outcome model & Misspecified exposure model
正しい推定結果が得られるはず。
TMLE.fit1 <- df.sim %>%
select(C1,C2,C3,Exposure,Outcome) %>%
ltmle(Anodes = "Exposure", Ynodes = "Outcome", abar = list(1,0),
Qform = c(Outcome = "Q.kplus1 ~ C1 + I(C1^2) + C3 + C2*Exposure"), ) #アウトカムモデルを指定 (曝露モデルはデフォルトで掛け算項なし&線形性の仮定)
summary(TMLE.fit1)## Estimator: tmle
## Call:
## ltmle(data = ., Anodes = "Exposure", Ynodes = "Outcome", Qform = c(Outcome = "Q.kplus1 ~ C1 + I(C1^2) + C3 + C2*Exposure"),
## abar = list(1, 0))
##
## Treatment Estimate:
## Parameter Estimate: 67.286
## Estimated Std Err: 0.072615
## p-value: <2e-16
## 95% Conf Interval: (67.143, 67.428)
##
## Control Estimate:
## Parameter Estimate: 61.566
## Estimated Std Err: 0.055996
## p-value: <2e-16
## 95% Conf Interval: (61.457, 61.676)
##
## Additive Treatment Effect:
## Parameter Estimate: 5.7191
## Estimated Std Err: 0.039036
## p-value: <2e-16
## 95% Conf Interval: (5.6426, 5.7956)ATE = 5.72なので微妙なところですが・・・一応真の値（5.8）にニアピンということにしましょう。
Correct exposure model & Misspecified outcome model
これも正しいはず。
TMLE.fit2 <- df.sim %>%
select(C1,C2,C3,Exposure,Outcome) %>%
ltmle(Anodes = "Exposure", Ynodes = "Outcome", abar = list(1,0),
gform = "Exposure ~ C1 + I(C1^2) + C2*C3") #曝露モデルを指定 (アウトカムモデルはデフォルトで掛け算項なし&線形性の仮定)
summary(TMLE.fit2)## Estimator: tmle
## Call:
## ltmle(data = ., Anodes = "Exposure", Ynodes = "Outcome", gform = "Exposure ~ C1 + I(C1^2) + C2*C3",
## abar = list(1, 0))
##
## Treatment Estimate:
## Parameter Estimate: 67.356
## Estimated Std Err: 0.10985
## p-value: <2e-16
## 95% Conf Interval: (67.14, 67.571)
##
## Control Estimate:
## Parameter Estimate: 61.569
## Estimated Std Err: 0.059658
## p-value: <2e-16
## 95% Conf Interval: (61.452, 61.686)
##
## Additive Treatment Effect:
## Parameter Estimate: 5.7864
## Estimated Std Err: 0.10354
## p-value: <2e-16
## 95% Conf Interval: (5.5835, 5.9894)ATE =5.79なのでほぼ真値=5.8と一致。
Misspecified outcome & exposure models
これはモデルの誤設定によるバイアスが生じるはず。
TMLE.fit3 <- df.sim %>%
select(C1,C2,C3,Exposure,Outcome) %>%
ltmle(Anodes = "Exposure", Ynodes = "Outcome", abar = list(1,0)) #アウトカムモデル・曝露モデルともにデフォルト
summary(TMLE.fit3)## Estimator: tmle
## Call:
## ltmle(data = ., Anodes = "Exposure", Ynodes = "Outcome", abar = list(1,
## 0))
##
## Treatment Estimate:
## Parameter Estimate: 66.811
## Estimated Std Err: 0.12335
## p-value: <2e-16
## 95% Conf Interval: (66.57, 67.053)
##
## Control Estimate:
## Parameter Estimate: 61.692
## Estimated Std Err: 0.056397
## p-value: <2e-16
## 95% Conf Interval: (61.582, 61.803)
##
## Additive Treatment Effect:
## Parameter Estimate: 5.1188
## Estimated Std Err: 0.11625
## p-value: <2e-16
## 95% Conf Interval: (4.891, 5.3466)ATE = 5.12となり、かなりズレた。
まとめ
|
手法
|
推定値 | 説明 |
|---|---|---|
| 真のConditional Effect (女性） | 5 | |
| 真のConditional Effect (男性） | 7 | |
| 真のMarginal Effect | 5.8 | |
| 重回帰分析（正しいアウトカムモデル） | Exposure = 5; Exposure*C2 = 2 | 正しく推定できた |
| 重回帰分析（調整なしモデル） | 3.32 | 交絡バイアス |
| 重回帰分析（調整あり＆モデルの誤設定１） | 4.88 | C1の線形性 & Exposur e*C2がないという仮定 |
| 重回帰分析（調整あり＆モデルの誤設定２） | 5.63 | Exposur e*C2がないという仮定 |
| 標 準化/g-formula（正しいアウトカムモデル） | 5.79 | 正しく推定できた |
| 標準化/g-form ula（モデルの誤設定） | 4.88 | C1の線形性 & Exposur e*C2がないという仮定 |
| 傾向スコア：層化 | 4.02 ~ 7.26 | PS層ごとの効果推定 |
| 傾向スコア：回帰モデルでの調整 | 5.63 | PSの線形性＆PSによ る効果修飾なしの仮定 |
| 傾向スコア：回帰モデルでの調整＆g-formula | 5.90 | PSの線形性の仮定 |
| 傾向スコア：マッチング | 5.54 | C2の分布がマ ッチング前後で変化。 効果修飾因子の分布が 異なる集団への推論。 |
| IPTW（正しい曝露モデル） | 5.82 | 正しく推定できた |
| IPTW（モデルの誤設定） | 5.12 | C1の線形性＆C 2*C3がないという仮定 |
| Doubly-robust（アウトカムモデルだけ正しい） | 5.72 | 正し く推定できた（はず） |
| Doubly-robust（曝 露モデルだけ正しい） | 5.79 | 正しく推定できた |
| Doubly-robust（両方誤設定） | 5.12 | モデルの 誤設定によるバイアス |
概ね予想通りの結果になったと言えます。
「どの手法を使ってなにを調整すべきか」だけでなく、モデルの設定が重要ということがわかります。
また傾向スコアマッチングの場合は正しい（曝露）モデルを使っていても推定結果が変わってくるというのが面白いですね。
