- 因果推論したいとき、交絡をいかにコントロールするかは大切
- 交絡のコントロールの方法はたくさんある
- 層別分析
- 多変量解析（重回帰分析など）
- 傾向スコア（マッチング、逆確率重みづけ、など）
などなど
- 交絡が生じているときにコントロールしないとどうなるか実験
交絡が生じていない状況では？
# 処置あり(x = 1)のとき、処置なし(x = 0)と比べて、yが平均50だけ上がるモデル
set.seed(1234)
d <- tibble(
x = rbinom(n = sample_size, size = 1, prob = 0.5),
y = rnorm(n = sample_size, mean = 10 + 50 * x, sd = 10)
)
# xとyの関係を可視化
ggplot(data = d %>% mutate(x = as.character(x))) +
geom_histogram(aes(x = y, fill = x), alpha = 2/6)
# 単純な線形回帰で処置効果を推定
model_lm <- lm(y ~ x, data = d)
broom::tidy(model_lm)
# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 10.1 0.140 72.2 0
2 x 50.0 0.198 252. 0
- xのestimateが真値の50付近なので、うまく推定できとる
- 交絡が生じていないのだから当然
交絡が生じている状況では？
# 処置あり(x = 1)のとき、処置なし(x = 0)と比べてyが平均50だけ上がるモデル
# ただし、交絡因子cがある状況
set.seed(1234)
d <- tibble(
c = rnorm(n = sample_size, mean = 10, sd = 5),
x = rbinom(n = sample_size, size = 1, prob = 1 / (1 + exp(-1 * (- 10 + c)))),
y = rnorm(n = sample_size, mean = 10 + 50 * x + 20 * c, sd = 10)
)
# 処置効果を推定
model_lm <- lm(y ~ x, data = d)
broom::tidy(model_lm)
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 136. 0.941 145. 0
2 x 198. 1.33 149. 0
- 真値の50からめちゃ逸脱。過大評価ワロタ
- 交絡が生じているときにコントロールしてないから当然
重回帰モデルでコントロール
# 処置効果を推定
model_lm_with_covariate <- lm(y ~ x + c, data = d)
broom::tidy(model_lm_with_covariate)
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 9.73 0.240 40.5 0
2 x 49.9 0.303 165. 0
3 c 20.0 0.0307 653. 0
- 交絡因子を共変量として組み込めばよい
- 真値の50とほぼ一致。ええ感じ。
傾向スコアをつかったマッチング
m_near <- MatchIt::matchit(
formula = x ~ c,
data = d,
method = 'nearest',
replace = T
)
match_data <- MatchIt::match.data(m_near)
model_lm_with_ps_score_matching <- lm(y ~ x, data = match_data)
broom::tidy(model_lm_with_ps_score_matching)
love.plot(m_near, threshold = 0.1, abs = T)
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 210. 2.20 95.7 0
2 x 124. 2.37 52.2 0
- 傾向スコアを求めて最近隣法でマッチング
- 処置ありだけでの層別分析と枠組みは同じ（らしい）
- 過大評価のときよりは推定は多少マシになっている。ただし不完全。
- ATTだから？？う～む勉強不足
傾向スコアを使った逆確率重みづけ
weight <- WeightIt::weightit(
formula = x ~ c,
data = d,
method = 'ps',
estimand = 'ATE'
)
model_lm_with_ps_score <- lm(y ~ x + c, data = d, weights = weight$weights)
broom::tidy(model_lm_with_ps_score)
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 6.18 0.268 23.0 2.00e-114
2 x 46.4 0.207 224. 0
3 c 20.5 0.0233 878. 0
- 傾向スコアが小さいデータは欠損しやすいと考えられる
- 傾向スコアの逆数で欠損値補完をしている（らしい）
- ATEなので推定値が真値がかなり近づいた
