G-computation and IPW
Hiroyuki Kuroda
2023-04-28
0. はじめに
Implementation of G-Computation on a Simulated Data Set:
Demonstration of a Causal Inference Technique
https://doi.org/10.1093/aje/kwq472
本ページは、上記論文中で扱われたシミュレーションデータを用いたG-computation
(parametric g-formula)と、IPWによるATEの推定を実践するものです。
1. 準備
1.1. パッケージの準備
下記のパッケージを使用します。
pacman::p_loadを使用すると、複数のパッケージを一度に読み込めます。
また、インストールしていないものは自動でインストールしてくれます。
pacman::p_load(tidyverse, # 主にtidyverse記法を用います。
boot, # 信頼区間の推定にbootstrap法を用います。
broom, # broom::tidyで回帰分析の結果を簡潔にまとめられます。
cobalt, # 傾向スコアを用いた共変量分布の提示に用います。
survey) # IPW法による重みづけ分析に用います1.2. シミュレーションデータの生成
乱数シードは論文のsupplimentaryと同じく、285としました。
rbinomの引数に線形予測子をそのまま使用していますが、本来はlogistic関数を通してPr(A)にするべきと思います。今回はサプリのコードのまま話を進めます。
n <- 300
set.seed(285)
simdata <- data.frame(W1 = rbinom(n, 1, 0.4),
W2 = rbinom(n, 1, 0.5))
simdata <- simdata %>%
mutate(A = rbinom(n, 1, (0.5 + 0.2*W1 - 0.3*W2))) %>%
mutate(Y = rnorm(n, (3 - 0.5*A + W1 + 0.3*A*W2), 0.4))
# 補足：Pr(A)に戻すなら
# # ロジスティック関数を定義
# logistic <- function(x) {
# return(1 / (1 + exp(-x)))
# }
#
# simdata_botu <- simdata %>%
# mutate(A = rbinom(n, 1, logistic(0.5 + 0.2*W1 - 0.3*W2))) %>%
# mutate(Y = rnorm(n, (3 - 0.5*A + W1 + 0.3*A*W2), 0.4))2. G-computation
2.1. Q-modelの作成
論文のsupplimentaryにならい、product termを含まないmodel 2と、真のモデルであるmodel 4について検討する。
# model 2
reg2 <- glm(Y ~ A + W1 + W2, data = simdata, family = gaussian)
# model 4
reg4 <- glm(Y ~ A + W1 + A:W2, data = simdata, family = gaussian)2.2. Y^a=1, Y^a=0の推定
全てのAを0もしくは1にした反事実データを作成する。
simdata_A0 <- simdata %>% mutate(A = 0)
simdata_A1 <- simdata %>% mutate(A = 1)
モデルreg2、reg4それぞれでsimdata_A0/A1それぞれに対してYをpredictし、2つのモデルに基づくY^a=1,
Y^a=0を推定する。
# model 2
Y0_m2 <- predict(reg2, newdata = simdata_A0)
Y1_m2 <- predict(reg2, newdata = simdata_A1)
# model 4
Y0_m4 <- predict(reg4, newdata = simdata_A0)
Y1_m4 <- predict(reg4, newdata = simdata_A1)
表にまとめて、確認してみる。その際、Y^a=1 -
Y^a=0も算出しておく。
summary <- simdata %>%
mutate(Y0_m2 = Y0_m2,
Y1_m2 = Y1_m2,
Y1_Y0_m2 = Y1_m2 - Y0_m2,
Y0_m4 = Y0_m4,
Y1_m4 = Y1_m4,
Y1_Y0_m4 = Y1_m4 - Y0_m4)
head(summary)## W1 W2 A Y Y0_m2 Y1_m2 Y1_Y0_m2 Y0_m4 Y1_m4 Y1_Y0_m4
## 1 1 1 0 2.852289 4.053675 3.697305 -0.3563703 3.998715 3.819475 -0.1792394
## 2 0 1 0 2.783062 2.993058 2.636688 -0.3563703 2.947648 2.768409 -0.1792394
## 3 0 1 0 2.790881 2.993058 2.636688 -0.3563703 2.947648 2.768409 -0.1792394
## 4 1 0 0 3.760526 3.925975 3.569604 -0.3563703 3.998715 3.509826 -0.4888889
## 5 1 1 0 4.538255 4.053675 3.697305 -0.3563703 3.998715 3.819475 -0.1792394
## 6 1 0 1 3.506819 3.925975 3.569604 -0.3563703 3.998715 3.509826 -0.48888892.3. marginal effectの推定
E[Y^a=1] - E[Y^a=0] = E[Y^a=1 -
Y^a=0]なので、下記のように推定できる。
詳細な根拠は発表内容を参照。ここでは点推定値のみ。信頼区間は下記へ続く。
# model 2
ME_m2 <- mean(summary$Y1_Y0_m2)
cat("The marginal effect based on model 2 is", round(ME_m2, 3), "\n")## The marginal effect based on model 2 is -0.356# model 4
ME_m4 <- mean(summary$Y1_Y0_m4)
cat("The marginal effect based on model 4 is", round(ME_m4, 3), "\n")## The marginal effect based on model 4 is -0.3362.4 bootstrap法による信頼区間の算出
ここではmodel 4のみ扱う。
まず、bootパッケージのboot関数で動かすためのbootstrap関数を定義する。
# bootstrap関数の定義
boot_glm_m4 <- function(data, indices) {
data_boot <- data[indices, ]
model_boot_m4 <- glm(Y ~ A + W1 + A:W2,
data = data_boot, family = gaussian)
# Y^a=0の算出
data_boot_A0_m4 <- data_boot %>% mutate(A = 0)
boot_Y0_m4 <- predict(model_boot_m4, newdata = data_boot_A0_m4)
# Y^a=1の算出
data_boot_A1_m4 <- data_boot %>% mutate(A = 1)
boot_Y1_m4 <- predict(model_boot_m4, newdata = data_boot_A1_m4)
# re-sampling dataでのmarginal effect推定
# これが定義した関数の返り値
boot_Y1_Y0_m4 = boot_Y1_m4 - boot_Y0_m4
return(mean(boot_Y1_Y0_m4))
}
bootstrapを実行する。サンプリング回数は1000とする。
boot_result_m4 <- boot(simdata, boot_glm_m4, R = 1000)
点推定値を表示する。
cat("The point estimate is", round(boot_result_m4$t0, 3), "\n")## The point estimate is -0.336
bootstrap信頼区間を表示する。
CI_boot <- boot.ci(boot_result_m4, conf = 0.95, type = "norm")
CI_boot## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_result_m4, conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% (-0.4474, -0.2236 )
## Calculations and Intervals on Original Scale3. IPW
3.1. 共変量(W1, W2)分布の確認
cobaltパッケージのlove.plotが便利。 treatment(A)ごとに分布はそろっていないことが分かる。
love_crude <- love.plot(A ~ W1 + W2,
data = simdata,
abs = TRUE,
stats = "mean.diffs",
thresholds = c(m = 0.1),
s.d.denom = "pooled")
love_crude
3.2. PSの推定
色々パッケージはあるが、今回は単純なケースなので用いない。
後でstabilized weightも試すため、Pr(A)も保存しておく。
# PSをロジスティック回帰で推定
psmodel <- glm(A ~ W1 + W2, data = simdata, family = binomial)
# PSを変数psに保存
ps <- psmodel$fitted.values
# stabilized weightの分子に用いるPr(A)
prA <- mean(simdata$A)
# PS、IP weight、stabilized weightもとのデータフレームに保存
simdata_ipw <- simdata %>%
mutate(ps = ps,
weight = if_else(A == 1 , 1/ps, 1/(1-ps)),
sw = if_else(A == 1 , prA/ps, (1-prA)/(1-ps)))
weightとstabilized weightを確認してみる。
weightはmean
2, stabilized weightはmean 1に（ほぼ）なっているのが分かる。
# non-stabilized weight
summary(simdata_ipw$weight)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.219 1.491 1.888 1.992 2.126 5.568# stabilized weight
summary(simdata_ipw$sw)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5960 0.7354 0.8435 0.9981 1.1389 2.20883.3. PSの分布と、重みづけ後の共変量のバランス
# PSの分布
plot <- bal.plot(simdata_ipw, treat = simdata_ipw$A, var.name = "ps")
plot
# 重みづけ後のW1, W2（共変量）のバランス
love_adj <- love.plot(A ~ W1 + W2,
data = simdata_ipw,
weights = simdata_ipw$weight,
abs = TRUE,
stats = "mean.diffs",
thresholds = c(m = 0.1),
s.d.denom = "pooled")
love_adj
共変量がバランシングされたことが分かる。
3.4. glmで重みづけ解析
point estimateはこのまま採用してもOK
信頼区間については、別途ロバスト分散やブートストラップで算出する必要あり。
ipw_model <- glm(Y ~ A,
data = simdata_ipw, family = gaussian,
weights = weight)
result_ipw_model <- ipw_model %>% tidy(conf.int = TRUE)
result_ipw_model## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.38 0.0566 59.7 7.20e-168 3.27 3.49
## 2 A -0.332 0.0805 -4.12 4.85e- 5 -0.490 -0.1743.5. surveyパッケージを使用したIPW
surveyパッケージではロバスト信頼区間が算出される。
design <- svydesign(id = ~1, weights = ~weight, data = simdata_ipw)
ipw <- svyglm(Y ~ A, design = design, family = gaussian)
結果を表示する。
3.4.の結果と比べて広い信頼区間が算出されているのが分かる。
result <- ipw %>% tidy(conf.int = TRUE)
result## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.38 0.0496 68.2 7.34e-184 3.28 3.48
## 2 A -0.332 0.0916 -3.62 3.41e- 4 -0.512 -0.152confint(ipw_model)[2, ] #再掲## Waiting for profiling to be done...## 2.5 % 97.5 %
## -0.4895755 -0.1741052
stabilized weightを用いたIPWもやってみる。
design_s <- svydesign(id = ~1, weights = ~sw, data = simdata_ipw)
ipw_s <- svyglm(Y ~ A, design = design, family = gaussian)
stabilized IPWによる結果の表示。
今回はsaturated
modelなので、普通のweightと結果は変わらない。
time-varyingやcontinuous treatmentの場合、CIがnarrowになる。
result_s <- ipw_s %>% tidy(conf.int = TRUE)
result_s## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.38 0.0496 68.2 7.34e-184 3.28 3.48
## 2 A -0.332 0.0916 -3.62 3.41e- 4 -0.512 -0.1523.6. bootstrapを用いたIPW estimateのCI推定
bootstrap関数を定義する。
boot_ipw <- function(data, indices) {
data_boot_ipw <- data[indices, ]
result_boot_ipw <- glm(Y ~ A, data = data_boot_ipw,
family = gaussian, weights = weight)
# interceptはいらないので、[2]を指定した
return(coef(result_boot_ipw)[2])
}
bootstrapを実行する。サンプリング回数は1000とする。
boot_result_ipw <- boot(simdata_ipw, boot_ipw, R = 1000)
CIを計算する。
robust
CIとさほど変わらない結果が得られた。
CI_ipw <- boot.ci(boot_result_ipw, conf = 0.95, type = "norm")
CI_ipw## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_result_ipw, conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% (-0.5167, -0.1569 )
## Calculations and Intervals on Original Scale4. doubly robust estimation
4.1. b(W)の推定
A=0のサンプルと、A=1のサンプルをそれぞれ抽出する。
# A=0のサンプルだけ抽出
simdata_filterA0 <- simdata %>% filter(A == 0)
# A=1
simdata_filterA1 <- simdata %>% filter(A == 1)
それぞれの分割データでアウトカム回帰を行う。
ここで、アウトカム回帰について誤特定モデル(model 2と正しいモデル(model
4)を考えてみる。
model 4はAとのproduct
termがあるので、A=1の回帰の時だけW2を含める。
（A=0の時は必然的にA*W2の項がゼロになるので。
# model 2でアウトカム回帰
reg_A0_m2 <- glm(Y ~ W1 + W2, data = simdata_filterA0, family = gaussian)
reg_A1_m2 <- glm(Y ~ W1 + W2, data = simdata_filterA1, family = gaussian)
# model 4でアウトカム回帰、A1の時だけW2を含める
reg_A0_m4 <- glm(Y ~ W1, data = simdata_filterA0, family = gaussian)
reg_A1_m4 <- glm(Y ~ W1 + W2, data = simdata_filterA1, family = gaussian)
それぞれ全データにフィットさせて、b(W)を推定する。
# model 2ベースのb^a=0(W)とb^a=1(W)
b_A0_m2 <- predict(reg_A0_m2, newdata = simdata)
b_A1_m2 <- predict(reg_A1_m2, newdata = simdata)
# model 4
b_A0_m4 <- predict(reg_A0_m4, newdata = simdata)
b_A1_m4 <- predict(reg_A1_m4, newdata = simdata)4.2. doubly robust estimation(DR)
PSを推定しておく（前のチャプターと一緒）。
psmodel <- glm(A ~ W1 + W2, data = simdata, family = binomial)
ps <- psmodel$fitted.values
定義通りに推定する。
平均をとればATEになるよう、DR^a=1 -
DR^a=0の変数も作っておく。
result <- simdata %>%
mutate(ps = ps) %>%
mutate(b_A0_m2 = b_A0_m2,
b_A1_m2 = b_A1_m2,
b_A0_m4 = b_A0_m4,
b_A1_m4 = b_A1_m4) %>%
mutate(DR0_m2 = b_A0_m2 + (((1-A)*(Y-b_A0_m2))/(1-ps)),
DR1_m2 = b_A1_m2 + ((A*(Y - b_A1_m2))/ps),
DR0_m4 = b_A0_m4 + (((1-A)*(Y-b_A0_m4))/(1-ps)),
DR1_m4 = b_A1_m4 + ((A*(Y - b_A1_m4))/ps)) %>%
mutate(DR_m2 = DR1_m2 - DR0_m2,
DR_m4 = DR1_m4 - DR0_m4)
最後にmodel 2, 4それぞれの推定値を比較してみる。
ATE_m2 <- mean(result$DR_m2)
ATE_m4 <- mean(result$DR_m4)
cat("The DR-ATE based on outcome reg model 2 is", round(ATE_m2, 3), "\n")## The DR-ATE based on outcome reg model 2 is -0.338cat("The DR-ATE based on outcome reg model 4 is", round(ATE_m4, 3), "\n")## The DR-ATE based on outcome reg model 4 is -0.338
アウトカムモデルの誤特定があっても、推定値に差がないことが分かる。
傾向スコアモデルも誤特定したものを試してもよいかも。
4.3. おわりに
doubly robust
estimationの信頼区間もブートストラップで求めます。
基本は上述のコードと変わらないので、ここでは省略します。。。
お疲れさまでした！！！