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を参考にしました。勉強になりました。ありがとうございます。

Rで因果推論~その12: 傾向スコアの可視化~

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

傾向スコアの可視化

因果推論を行う上で傾向スコアを求め、IPW法の計算で使用し、ATEやATT推定量の計算時の重みとして利用したり、回帰分析の調整変数にしたり、マッチングのキーにしたり、傾向スコアは様々な方法で使われる。今回は、その傾向スコアを使って様々な可視化を通して、傾向スコアへの理解を深める。

サンプルデータは、岩波データサイエンスvol3の加藤・星野先生によって書かれた「因果効果推定の応用~CM接触の因果効果と調整効果~」という章で使用されているダミーのデータ。

準備

ここでは傾向スコアを計算し、データフレームに戻し、傾向スコアの判別精度を確認する指標であるc統計量を可視化する。c統計量といってもROCのAUCのことである。{pROC}パッケージを利用し、1000回のブートストラップによって、信頼区間を計算している。

library(RCurl)
library(tidyverse)
library(scales)
library(tableone)
library(pROC)
library(Matching)

url_txt <- "https://raw.githubusercontent.com/iwanami-datascience/vol3/master/kato%26hoshino/q_data_x.csv"
url <- getURL(url_txt)
df <- read.csv(text = url, header = TRUE)

logi <- glm(cm_dummy ~ TVwatch_day + age + sex + marry_dummy + child_dummy + inc + pmoney +
              area_kanto +area_tokai + area_keihanshin + 
              job_dummy1 + job_dummy2 + job_dummy3 + job_dummy4 + job_dummy5 + job_dummy6 + job_dummy7  + 
              fam_str_dummy1 + fam_str_dummy2 + fam_str_dummy3 + fam_str_dummy4, 
            family=binomial(link="logit") , data = df)

df2 <- df %>% bind_cols(.,  tibble(ps = logi$fitted.values)) 

fit_roc <- roc(cm_dummy ~ ps, data = df2, ci = TRUE)
cis_boot <- ci(fit_roc, of = "se",sp = seq(0, 1, 0.01), boot.n=1000)

plot(fit_roc, print.auc=TRUE, legacy.axes = TRUE)
plot(cis_boot, type="shape", col="#00860022", no.roc=TRUE)

f:id:AZUMINO:20190505153338p:plain

傾向スコアのオーバラップ

実際に得られた傾向スコアのオーバーラップ具合を計算する。%で計算している。

ggplot() +
  geom_histogram(data = df2 %>% filter(cm_dummy == 0), 
                 aes(ps, y = (..count../sum(..count..))),
                 position = "identity", alpha = 0.4, binwidth = 0.02, fill = "darkred", col = "black") + 
  geom_histogram(data = df2 %>% filter(cm_dummy == 1), 
                 aes(ps, y = (..count../sum(..count..))),
                 position = "identity", alpha = 0.4, binwidth = 0.02, fill = "royalblue", col = "black") + 
  scale_y_continuous(labels = percent) +
  xlab("Propensity Score") + ylab("Percent") + theme_classic()

f:id:AZUMINO:20190505151053p:plain

医学系の論文でよく見るパターンの可視化。こっちはカウントで計算。

df_tmp <- df2 %>%
  dplyr::select(cm_dummy, ps) %>% 
  dplyr::mutate(id = row_number()) %>% 
  tidyr::spread(cm_dummy, ps, sep = "_")

ggplot(df_tmp) + 
  geom_histogram(aes(x = cm_dummy_0, y =  ..count..), alpha = 0.4, binwidth = 0.02, fill = "darkred", col = "black") + 
  geom_histogram(aes(x = cm_dummy_1, y = -..count..), alpha = 0.4, binwidth = 0.02, fill = "royalblue", col = "black") + # `-` sign needs to flip
  geom_hline(yintercept = 0, lwd = 0.5) +
  scale_y_continuous(label = abs, limits = c(-400, 400)) +  
  ylab("Count") + xlab("Propensity Score")  + theme_classic() + 
  coord_flip()

f:id:AZUMINO:20190505151130p:plain

バランスの確認

{tableone}パッケージのCreateTableOne()を使って、ベースライン特性を把握する。ここでは、マッチング前のデータに対してテーブル表を作成しているので、ベースラインがどの程度ずれているのかを確認できる。

table <- CreateTableOne(
  vars = c("age","sex","marry_dummy","inc","pmoney" ,"TVwatch_day",
           "area_kanto","area_keihan","area_tokai","area_keihanshin",
           "job_dummy1","job_dummy2","job_dummy3","job_dummy4","job_dummy5","job_dummy6","job_dummy7","job_dummy8",
           "fam_str_dummy1","fam_str_dummy2","fam_str_dummy3","fam_str_dummy4","fam_str_dummy5","child_dummy",
           "T","F1","F2","F3","M1","M2","M3"),
  strata = "cm_dummy",
  data = df2,
  factorVars = c("area_kanto","area_keihan","area_tokai","area_keihanshin",
                 "job_dummy1","job_dummy2","job_dummy3","job_dummy4","job_dummy5","job_dummy6","job_dummy7","job_dummy8",
                 "fam_str_dummy1","fam_str_dummy2","fam_str_dummy3","fam_str_dummy4","fam_str_dummy5","child_dummy",
                 "T","F1","F2","F3","M1","M2","M3"))

table
                        Stratified by cm_dummy
                          0                 1                  p      test
  n                          5856               4144                      
  age (mean (SD))           40.19 (10.64)      41.77 (10.15)   <0.001     
  sex (mean (SD))            0.67 (0.47)        0.60 (0.49)    <0.001     
  marry_dummy (mean (SD))    0.63 (0.48)        0.67 (0.47)    <0.001     
  inc (mean (SD))          369.25 (264.49)    341.70 (270.70)  <0.001     
  pmoney (mean (SD))         3.56 (3.36)        3.54 (3.40)     0.841     
  TVwatch_day (mean (SD)) 5714.98 (5690.37) 11461.88 (8851.09) <0.001     
  area_kanto = 1 (%)          369 ( 6.3)         543 (13.1)    <0.001     
  area_keihan = 1 (%)        2981 (50.9)        2906 (70.1)    <0.001     
  area_tokai = 1 (%)          729 (12.4)         386 ( 9.3)    <0.001     
  area_keihanshin = 1 (%)    1777 (30.3)         309 ( 7.5)    <0.001     
  job_dummy1 = 1 (%)         3517 (60.1)        2145 (51.8)    <0.001     
  job_dummy2 = 1 (%)          334 ( 5.7)         208 ( 5.0)     0.149     
  job_dummy3 = 1 (%)          403 ( 6.9)         356 ( 8.6)     0.002     
  job_dummy4 = 1 (%)           67 ( 1.1)          56 ( 1.4)     0.404     
  job_dummy5 = 1 (%)          531 ( 9.1)         646 (15.6)    <0.001     
  job_dummy6 = 1 (%)          539 ( 9.2)         460 (11.1)     0.002     
  job_dummy7 = 1 (%)          272 ( 4.6)         127 ( 3.1)    <0.001     
  job_dummy8 = 1 (%)          193 ( 3.3)         146 ( 3.5)     0.573     
  fam_str_dummy1 = 1 (%)      881 (15.0)         599 (14.5)     0.430     
  fam_str_dummy2 = 1 (%)      750 (12.8)         698 (16.8)    <0.001     
  fam_str_dummy3 = 1 (%)     3625 (61.9)        2579 (62.2)     0.752     
  fam_str_dummy4 = 1 (%)      489 ( 8.4)         210 ( 5.1)    <0.001     
  fam_str_dummy5 = 1 (%)      111 ( 1.9)          58 ( 1.4)     0.069     
  child_dummy = 1 (%)        2471 (42.2)        1759 (42.4)     0.818     
  T = 1 (%)                    79 ( 1.3)          53 ( 1.3)     0.831     
  F1 = 1 (%)                  810 (13.8)         468 (11.3)    <0.001     
  F2 = 1 (%)                  797 (13.6)         936 (22.6)    <0.001     
  F3 = 1 (%)                  293 ( 5.0)         229 ( 5.5)     0.266     
  M1 = 1 (%)                  977 (16.7)         426 (10.3)    <0.001     
  M2 = 1 (%)                 2001 (34.2)        1287 (31.1)     0.001     
  M3 = 1 (%)                  899 (15.4)         745 (18.0)     0.001 

マッチング後のオーバラップ

実際にマッチングした後の傾向スコアの重なり具合を確認。

set.seed(1234)
match_res <- Match(Y = df2$gamesecond, Tr = df2$cm_dummy, X = df2$ps, replace = TRUE, M = 1)
summary(match_res)


df3 <- rbind(df2 %>% slice(match_res$index.treated) %>% mutate(flg = "treat"),
             df2 %>% slice(match_res$index.control) %>% mutate(flg = "control"))

ggplot() +
  geom_histogram(data = df3 %>% filter(flg == "treat"), 
                 aes(ps, y = (..count../sum(..count..))),
                 position = "identity", alpha = 0.4, binwidth = 0.02, fill = "darkred", col = "black") + 
  geom_histogram(data = df3 %>% filter(flg == "control"),
                 aes(ps, y = (..count../sum(..count..))),
                 position = "identity", alpha = 0.4, binwidth = 0.02, fill = "royalblue", col = "black") + 
  scale_y_continuous(labels = percent) +
  xlab("Matching Distribution") + ylab("Percent") + theme_classic()

f:id:AZUMINO:20190505151241p:plain

マッチング後の共変量

実際にマッチングした後の共変量の重なり具合を確認。loessで計算している線の重なり具合でマッチングの具合を確認できる。

ggplot(df3 %>% sample_n(10000), aes(ps, sex, col = factor(cm_dummy))) + 
  geom_point(alpha = 0.1, size = 1) +
  geom_smooth(method = "loess", se = FALSE) +
  xlab("Propensity score") + ylab("Sex") + theme_bw() + 
  scale_color_hue(name = "cm_dummy", labels = c(`0` = "Ctrl",`1` = "Treat"))

f:id:AZUMINO:20190505151410p:plain

ggplot(df3 %>% sample_n(10000), aes(ps, age, col = factor(cm_dummy))) + 
  geom_point(alpha = 0.1, size = 1) +
  geom_smooth(method = "loess", se = FALSE) +
  xlab("Propensity score") + ylab("Age") + theme_bw() + 
  scale_color_hue(name = "cm_dummy", labels = c(`0` = "Ctrl",`1` = "Treat"))

f:id:AZUMINO:20190505151420p:plain

ggplot(df3 %>% sample_n(10000), aes(ps, pmoney, col = factor(cm_dummy))) + 
  geom_point(alpha = 0.1, size = 1) +
  geom_smooth(method = "loess", se = FALSE) +
  xlab("Propensity score") + ylab("Pmoney") + theme_bw() + 
  scale_color_hue(name = "cm_dummy", labels = c(`0` = "Ctrl",`1` = "Treat"))

f:id:AZUMINO:20190505151443p:plain