Rのこと。

R言語のことをまとめる備忘録

Rで因果推論~その13:Consistency of IPW~

Note:スマートフォンでは数式が表示されないのでPC版で開いて見る。

はじめに

この記事は因果推論について学習した内容を自分の備忘録としてまとめたものです。 rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp rlang.hatenablog.jp

逆確率重み付け推定量の一致性

条件付き交換可能性もとで、逆確率重み付け推定量は母集団のATEに対して不偏となるかどうかを調べる。\( Y_{i}(t),T_{i},X_{i} \)は、トリートメントが\(t=1,0\)のもとでの潜在的アウトカムのベクトル。そして、条件付き交換可能性(Conditional Exchangeability)は成り立つと仮定。つまり、\( Y_{i}(t) ⊥ T_{i}|X_{i}, \forall t \)です。\(f(T_{i}|X_{i})\)は、傾向スコア(条件付き確率分布)。これに逆確率重み付けを行う。 $$ \sum_{ i = 1 }^{ n } \frac {T_iY_i}{ T_i } - \sum_{ i = 1 }^{ n } \frac{(1-T_i)Y_i}{ (1- T_i) } $$ 第1項の部分が証明できれば、第2項の部分も同じなので、IPWで重み付けた第1項に焦点をあてる。 $$ \frac{E[ \frac{I(T_{i}=t)}{f(T_{i}|X_{i})}Y_{i}]}{E[ \frac{I(T_{i}=t)}{f(T_{i}|X_{i})}]}=E[Y_{i}(t)] $$ まずは先程の分子に焦点をあてて話を進める。1,4行目は繰り返し期待値の法則、5行目は条件付き交換可能性を利用。 $$ \begin{align} Numerator &= E[\frac{I(T_{i}=t)}{f(T_{i}|X_{i})}Y_{i}] \\ &=E[\frac{I(T_{i}=t)}{f(T_{i}|X_{i})}[Y_{i}(1)T_{i} + Y_{i}(0)(1-T_{i})]] \\ &=E[\frac{I(T_{i}=t)}{f(T_{i}|X_{i})}Y_{i}(t)] \\ &=E_{X_{i}}E_{Y_{i}(t), T_{i}|X_{i}}[\frac{I(T_{i}=t)}{f(T_{i}|X_{i})}Y_{i}(t)|X_{i}] \\ &=E_{X_{i}}[E_{T_{i}|X_{i}}[\frac{I(T_{i}=t)}{f(T_{i}|X_{i})}|X_{i}] E_{Y_{i}(t)|X_{i}}[Y_{i}(t)|X_{i} ]] \\ &=E_{X_{i}}[1\cdot E_{Y_{i}(t)|X_{i}}[Y_{i}(t)|X_{i}]] \\ &=E[Y_{i}(t)] \\ \end{align} $$ 分母に分解して話を進めますが、これは1となる。 $$ \begin{align} Denominator &= E_{X_{i}}[E_{T_{i}|X_{i}}[\frac{I(T_{i}=t)}{f(T_{i}|X_{i})}|X_{i}]] \\ &=1 \\ \end{align} $$ このことから、条件付き交換可能性もとでは、逆確率重み付け推定量は母集団のATEに対して不偏となる。

Rでシュミレーション

データ生成過程を用いて、IPWが真の値を推定できるかどうかをRでシュミレーションする。まずサンプルサイズ1000の共変量を2つ作成。

set.seed(1)
n <- 1000
x1 <- rnorm(n)
x2 <- rnorm(n)

次に傾向スコアの真のモデルと、そこからトリートメントの割付を生成します。x1が高ければ割り付けられやすく、x2が低いと割り付けられにくくモデルを生成する。

sig <- function(x){
  res = 1/(1+exp(-x))
  return(res)
}

## 傾向スコアのモデル
p <- sig(-1*x1 + 1*x2)

## 傾向スコアが高いほうが割り付けられやすい
t <- rbinom(n,1,p)

傾向スコアが大きいと、アウトカムが大きくなるようにする。つまり、x2が大きいと(x1が小さいと)、トリートメントされやすくなるので、モデルはx2の係数をプラス、x1はマイナスに設定。その結果、傾向スコアが高いとアウトカムが大きくなるが、因果効果がマイナスでも、係数によってトリートメント群のほうがyの平均は大きくなるという間違った結果を導ける。そして、こから単純にアウトカムを比較すると、因果効果はアウトカムをプラス1するという誤った結論となる。観察研究の場合、本来はdfの値しかしらないので生成過程はわからない。なので、交絡、セレクションバイアスなどを無視して単純比較しても、因果効果は推定できないし、意味はない。

## アウトカムのモデル
y = 100 + (-10)*t - 5*x1 + 10*x2 + rnorm(n)

df <- data.frame(x1,x2,t,y) 
mean(df$y[t == 1]) - mean(df$y[t == 0])
[1] 1.056606

このような場合に、IPW定量を用いて因果効果を推定する。つまり、傾向スコアを求め、傾向スコアの逆確率で重み付ける。

## 傾向スコア
fit <- glm(t ~ x1 + x2,
           data = df,
           family = binomial(link = "logit"))

df$ps <- fit$fitted

# IPW
df$ipw_ate <- dplyr::if_else(df$t == 1, 1/df$ps, 1/(1-df$ps))

あとは、IPWで重み付けた回帰分析を行う。結果として、トリートメントの効果t=-10.45と正しい値と近い値を推定できている。つまり、傾向スコアの逆数で重み付けることで、ランダムな割り付けを実行したかのようなバランシングが可能になる。

# Wighted Regression
glm(y ~ t,
    data = df,
    family = gaussian(link = "identity"),
    weights = ipw_ate)

Call:  glm(formula = y ~ t, family = gaussian(link = "identity"), data = df, 
    weights = ipw_ate)

Coefficients:
(Intercept)            t  
      98.97       -10.45  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      332200 
Residual Deviance: 277000   AIC: 7923


glm(y ~ t + I(1 - t) + 0,
    data = df,
    family = gaussian(link = "identity"),
    weights = ipw_ate)

Call:  glm(formula = y ~ t + I(1 - t) + 0, family = gaussian(link = "identity"), 
    data = df, weights = ipw_ate)

Coefficients:
       t  I(1 - t)  
   88.52     98.97  

Degrees of Freedom: 1000 Total (i.e. Null);  998 Residual
Null Deviance:      1.8e+07 
Residual Deviance: 277000   AIC: 7923

Data Science by R and PythonThe intuition behind inverse probability weighting in causal inferenceを参考にしました。勉強になりました。ありがとうございます。