Unboundedly

統計的因果推論・疫学についてのお話

モデルに基づく因果推論の各種手法をRで実装&結果を比較してみた

因果推論のための分析手法は様々ありますが、回帰モデルを使った主なアプローチのRでの実装方法とその推定結果の比較をします。

モチベーション的な部分は以下をご参照ください。 

シミュレーションデータを使って、各手法がどのような(主にモデリングに関する)仮定に基づいているのか、それが結果の違いにどのように影響しているのかをみていきます。

なお、Rマークダウンで書いたものをコピペしただけなのであまり見た目は美しくないですが、予め ご了承ください。

そのうち書籍化予定ですので、そのときにはもっと見やすく&詳しく説明します・・・

Person Encoding in Laptop

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.62

Exposureに対する推定値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.79

ATEの推定値=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.88

ATEの推定値=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.        0

Exposureの推定値=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.90

ATE推定値 = 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 交絡バイアス
重回帰分析(調整あり&モデルの誤設定1) 4.88 C1の線形性 & Exposur e*C2がないという仮定
重回帰分析(調整あり&モデルの誤設定2) 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 モデルの 誤設定によるバイアス

概ね予想通りの結果になったと言えます。

「どの手法を使ってなにを調整すべきか」だけでなく、モデルの設定が重要ということがわかります。

また傾向スコアマッチングの場合は正しい(曝露)モデルを使っていても推定結果が変わってくるというのが面白いですね。