この記事は以下の記事の続きです。
前回の記事では普通にロジスティック回帰で傾向スコアを求めたのですが、傾向スコアというのは元はと言えば「共変量に基づいてそれぞれの群に割り付けられる確率値を求めたもの」なので、やろうと思えば機械学習分類器で代替しても良いわけです。実際、岩波DS3にもそのように書かれています。
ということで、前回記事の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.gbdt <- xgb.train(params=list(objective="binary:logistic",eta=0.5), data=ddtrain, nrounds=20) > pred.gbdt <- predict(model.gbdt, newdata=ddtrain) > ps.gbdt <- cbind(1-pred.gbdt, pred.gbdt) # ROCカーブからAUC(即ちc統計量)を求める > roc.gbdt <- cbind(ps.gbdt[,2],d$cm_dummy) > roc.pred.gbdt <- prediction(roc.gbdt[,1],roc.gbdt[,2]) > auc.gbdt.tmp <- performance(roc.pred.gbdt, "auc") > as.numeric(auc.gbdt.tmp@y.values) [1] 0.9926357
ランダムフォレストの時と同様、c統計量が0.99を上回っています。これでATEとATTを求めてみると。。。
# ATE > iestp1.gbdt <- (ivec1/ps.gbdt[,2]) > iestp0.gbdt <- (ivec0/ps.gbdt[,1]) > iestp.gbdt <- iestp1.gbdt+iestp0.gbdt > ipwe_gs.gbdt <- lm(d$gamesecond ~ ivec+0, weights=iestp.gbdt) > summary(ipwe_gs.gbdt) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.gbdt) 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.gbdt <- ivec1 > iestp0_ATT.gbdt <- ivec0*ps.gbdt[,2]/ps.gbdt[,1] > iestp_ATT.gbdt <- iestp1_ATT.gbdt+iestp0_ATT.gbdt > ipwe_ATT_gs.gbdt <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.gbdt) > summary(ipwe_ATT_gs.gbdt) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.gbdt) 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.gbdt,breaks=50, main='Xgboost')
ランダムフォレストやxgboostといった精度が高く出やすい分類器だと両端にピークが来る分布になるのは当然だと思うんですが、これがオリジナルのロジスティック回帰による傾向スコアだとかなり平べったい感じのする分布になっていることが分かります。
そして元来c統計量=AUCが0.8以上であれば傾向スコアとして信頼できるという話になっていたわけですが、機械学習分類器3種による傾向スコアはいずれも0.9を上回っているにもかかわらず、IPW推定量であるATEとATTは奇妙な値を示しています。
もしかしたら、この「分布の形状」こそが傾向スコアとしての信頼性に影響するのかな?とも思ったのですが(そしてもしかしたらその辺の話が近年の機械学習を用いた傾向スコアに関する研究で指摘されているのかもしれませんが)、ちょっと今回はこれ以上は分かりませんでした。また何か新しい材料が出たら取り組んでみようと思います。
そうそう、僕のコードが間違っている可能性もいつも通り否定できませんので、もし見つけた方はコメントなりブコメなりでお知らせ下さいm(_ _)m