読者です 読者をやめる 読者になる 読者になる

六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

統計的因果推論(3): 傾向スコア算出を機械学習に置き換えてみると

R 統計学 統計的因果推論 機械学習

この記事は以下の記事の続きです。

前回の記事では普通にロジスティック回帰で傾向スコアを求めたのですが、傾向スコアというのは元はと言えば「共変量に基づいてそれぞれの群に割り付けられる確率値を求めたもの」なので、やろうと思えば機械学習分類器で代替しても良いわけです。実際、岩波DS3にもそのように書かれています。

岩波データサイエンス Vol.3

岩波データサイエンス Vol.3

ということで、前回記事のCM接触データセットに対して任意の機械学習分類器を用いて傾向スコアを算出した歳の、各種効果指標の違いを見ていこうかと思います。なおデータセットは前回から引き続きdというデータフレームに入っているとします。またAUCを求めるに当たっては{ROCR}パッケージを用いています。下準備として以下のようにインデックスを用意しておきます。

> ivec1 <- d$cm_dummy # Treated group
> ivec0 <- rep(1, nrow(d))-ivec1 # Untreated group
> ivec <- cbind(ivec1, ivec0)

あとは機械学習分類器を使っていくだけです。


SVMによる場合


まず{e1071}のsvmでやってみます。どうやっても出せるとは思いますが、一番簡単なのは以下のやり方です。

> library(e1071)
> model.svm <- svm(as.factor(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, d, probability = TRUE)
> ps.svm.source <- predict(model.svm, newdata=d[,-1], probability=TRUE)
> ps.svm <- attr(ps.svm.source, "probabilities")

# ROCカーブからAUC(即ちc統計量)を求める
> library(ROCR)
> roc.svm <- cbind(ps.svm[,2],d$cm_dummy)
> roc.pred.svm <- prediction(roc.svm[,1],roc.svm[,2])
> auc.svm.tmp <- performance(roc.pred.svm, "auc")
> as.numeric(auc.svm.tmp@y.values)
[1] 0.9035095


c統計量は大幅にロジスティック回帰で求めた場合を上回っています。このSVM版傾向スコアを用いてATEとATTを算出してみましょう。

# ATE

> iestp1.svm <- (ivec1/ps.svm[,2])
> iestp0.svm <- (ivec0/ps.svm[,1])
> iestp.svm <- iestp1.svm+iestp0.svm
> ipwe_gs.svm <- lm(d$gamesecond ~ ivec+0, weights=iestp.svm)
> summary(ipwe_gs.svm)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.svm)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-23385  -4654  -3507  -3040 389299 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   4318.1      282.7   15.28   <2e-16 ***
ivecivec0   2856.7      263.9   10.82   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26220 on 9998 degrees of freedom
Multiple R-squared:  0.03388,	Adjusted R-squared:  0.03368 
F-statistic: 175.3 on 2 and 9998 DF,  p-value: < 2.2e-16

# ATT

> iestp1_ATT.svm <- ivec1
> iestp0_ATT.svm <- ivec0*ps.svm[,2]/ps.svm[,1]
> iestp_ATT.svm <- iestp1_ATT.svm+iestp0_ATT.svm
> ipwe_ATT_gs.svm <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.svm)
> summary(ipwe_ATT_gs.svm)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.svm)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-16272  -2478  -1771   -917 220618 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   2478.1      212.5   11.66   <2e-16 ***
ivecivec0   2490.6      215.9   11.54   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13680 on 9998 degrees of freedom
Multiple R-squared:  0.02621,	Adjusted R-squared:  0.02602 
F-statistic: 134.6 on 2 and 9998 DF,  p-value: < 2.2e-16

あれれれれれ??? ATEではともかく、何故かATTではCM接触群の方が非接触群よりゲームプレー時間が短いことになっています。。。


ランダムフォレストによる場合


もちろん{randomForest}でやります。

> library(randomForest)
> > library(e1071)
> model.rf <- randomForest(as.factor(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, d)
> ps.rf <- predict(model.rf, newdata=d[,-1], type='prob')

# ROCカーブからAUC(即ちc統計量)を求める
> roc.rf <- cbind(ps.rf[,2],d$cm_dummy)
> roc.pred.rf <- prediction(roc.rf[,1],roc.rf[,2])
> auc.rf.tmp <- performance(roc.pred.rf, "auc")
> as.numeric(auc.rf.tmp@y.values)
[1] 0.9967271


c統計量はロジスティック回帰の時を大幅に上回っています。このランダムフォレスト版傾向スコアを用いて、ATEとATTを求めてみましょう。

# ATE

> iestp1.rf <- (ivec1/ps.rf[,2])
> iestp0.rf <- (ivec0/ps.rf[,1])
> iestp.rf <- iestp1.rf+iestp0.rf
> ipwe_gs.rf <- lm(d$gamesecond ~ ivec+0, weights=iestp.rf)
> summary(ipwe_gs.rf)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.rf)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
 -7734  -3550  -3336  -3254 364137 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   3151.0      282.4   11.16   <2e-16 ***
ivecivec0   3234.6      255.6   12.65   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20670 on 9786 degrees of freedom
  (212 observations deleted due to missingness)
Multiple R-squared:  0.02826,	Adjusted R-squared:  0.02807 
F-statistic: 142.3 on 2 and 9786 DF,  p-value: < 2.2e-16

# ATT

> iestp1_ATT.rf <- ivec1
> iestp0_ATT.rf <- ivec0*ps.rf[,2]/ps.rf[,1]
> iestp_ATT.rf <- iestp1_ATT.rf+iestp0_ATT.rf
> ipwe_ATT_gs.rf <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.rf)
> summary(ipwe_ATT_gs.rf)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.rf)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
 -5329  -2334  -1614   -622 220762 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   2334.2      176.1  13.253   <2e-16 ***
ivecivec0   3422.3      384.2   8.906   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11300 on 9786 degrees of freedom
  (31 observations deleted due to missingness)
Multiple R-squared:  0.02539,	Adjusted R-squared:  0.02519 
F-statistic: 127.5 on 2 and 9786 DF,  p-value: < 2.2e-16

あれれれれれれれれれれれ????? 今度はATEもATTもCM接触群より非接触群の方がゲームプレー時間が長いことになっています。。。これは何なんでしょうか。


Xgboostによる場合


言うまでもなく{xgboost}でやります。

> library(xgboost)
> library(Matrix)
> dtrain <- d[,c(1,3,5,6,7,8,9,10,11,12,13,14,15,16,18,19,20,21,22,23,25,33)]
> d.mx <- sparse.model.matrix(cm_dummy~., data=dtrain)
> ddtrain <- xgb.DMatrix(d.mx, label=d$cm_dummy)
> model.gdbt <- xgb.train(params=list(objective="binary:logistic",eta=0.5), data=ddtrain, nrounds=20)
> pred.gdbt <- predict(model.gdbt, newdata=ddtrain)
> ps.gdbt <- cbind(1-pred.gdbt, pred.gdbt)

# ROCカーブからAUC(即ちc統計量)を求める
> roc.gdbt <- cbind(ps.gdbt[,2],d$cm_dummy)
> roc.pred.gdbt <- prediction(roc.gdbt[,1],roc.gdbt[,2])
> auc.gdbt.tmp <- performance(roc.pred.gdbt, "auc")
> as.numeric(auc.gdbt.tmp@y.values)
[1] 0.9926357


ランダムフォレストの時と同様、c統計量が0.99を上回っています。これでATEとATTを求めてみると。。。

# ATE

> iestp1.gdbt <- (ivec1/ps.gdbt[,2])
> iestp0.gdbt <- (ivec0/ps.gdbt[,1])
> iestp.gdbt <- iestp1.gdbt+iestp0.gdbt
> ipwe_gs.gdbt <- lm(d$gamesecond ~ ivec+0, weights=iestp.gdbt)
> summary(ipwe_gs.gdbt)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.gdbt)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
 -7706  -3482  -3196  -3060 368640 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   2984.2      276.3   10.80   <2e-16 ***
ivecivec0   3033.6      242.4   12.52   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20550 on 9998 degrees of freedom
Multiple R-squared:  0.02661,	Adjusted R-squared:  0.02641 
F-statistic: 136.7 on 2 and 9998 DF,  p-value: < 2.2e-16

# ATT

> iestp1_ATT.gdbt <- ivec1
> iestp0_ATT.gdbt <- ivec0*ps.gdbt[,2]/ps.gdbt[,1]
> iestp_ATT.gdbt <- iestp1_ATT.gdbt+iestp0_ATT.gdbt
> ipwe_ATT_gs.gdbt <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.gdbt)
> summary(ipwe_ATT_gs.gdbt)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.gdbt)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
 -6326  -2478  -1649   -585 220618 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   2478.1      182.8  13.557   <2e-16 ***
ivecivec0   2709.0      321.8   8.419   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11770 on 9998 degrees of freedom
Multiple R-squared:  0.02484,	Adjusted R-squared:  0.02464 
F-statistic: 127.3 on 2 and 9998 DF,  p-value: < 2.2e-16

ランダムフォレストに引き続き、xgboostでもやはりATE, ATTともCM接触群よりも非接触群の方がゲームプレー時間が長いという逆転した結果になってしまいました。


これは一体どういうことなのか


もちろん何かがおかしいからこういうことになるのでしょうが、ちょっと気になってオリジナルのロジスティック回帰によるものと合わせて4つの傾向スコアのヒストグラムを比べてみました。

> fit <- glm(cm_dummy~., data=dtrain, family=binomial)
> ps <- fit$fitted.values # オリジナルの傾向スコア
> par(mfrow=c(2,2))
> hist(ps,breaks=50, main='Original PS: Logistic')
> hist(ps.svm,breaks=50, main='SVM')
> hist(ps.rf,breaks=50, main='Random Forest')
> hist(ps.gdbt,breaks=50, main='Xgboost')

f:id:TJO:20161005182601p:plain


ランダムフォレストやxgboostといった精度が高く出やすい分類器だと両端にピークが来る分布になるのは当然だと思うんですが、これがオリジナルのロジスティック回帰による傾向スコアだとかなり平べったい感じのする分布になっていることが分かります。


そして元来c統計量=AUCが0.8以上であれば傾向スコアとして信頼できるという話になっていたわけですが、機械学習分類器3種による傾向スコアはいずれも0.9を上回っているにもかかわらず、IPW推定量であるATEとATTは奇妙な値を示しています。


もしかしたら、この「分布の形状」こそが傾向スコアとしての信頼性に影響するのかな?とも思ったのですが(そしてもしかしたらその辺の話が近年の機械学習を用いた傾向スコアに関する研究で指摘されているのかもしれませんが)、ちょっと今回はこれ以上は分かりませんでした。また何か新しい材料が出たら取り組んでみようと思います。


そうそう、僕のコードが間違っている可能性もいつも通り否定できませんので、もし見つけた方はコメントなりブコメなりでお知らせ下さいm(_ _)m