この記事は以下の記事の続きです。
機械学習分類器で算出した傾向スコアを調整する話ですが、最後に課題として残ったのがprobability calibrationによる実践。探してみると前回の記事でもやったisotonic regressionとか色々出てくるんですが、もう一つ出てくるのがPlatt's scaling。これはあのSVMのSMOを提案したPlatt*1がprobability calibrationのために提案しているもので、web上にも幾つか資料があるようです。
やってることは極めて単純で、最初の分類確率値の計算をする際に0 / 1の二値への分類モデリングをやるのではなく、
という別の二値に変換した上で(ただしは正例の数では負例の数)、これを回帰したモデルで分類確率値を予測させるというだけです。
ということで、サクッとやってみましょう。ただしPlatt's scalingは流石はと言うべきかSVMの研究者が提案したものゆえ、SVMによる傾向スコア(というか分類確率値)のcalibrationのみを念頭に置いたものということなので、SVMの場合のみここでは扱います。その他の分類器による傾向スコアの調整は皆さんご自身でお願いしますということで。
Platt's scalingをやってみる
非常に簡単で、以下のようにやればおしまいです。
> d.platt <- d > table(d.platt$cm_dummy) 0 1 5856 4144 > d.platt$cm_dummy <- ifelse(d$cm_dummy==1, (4144+1)/(4144+2), 1/(5856+2)) > summary(d.platt$cm_dummy) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0001707 0.0001707 0.0001707 0.4144000 0.9998000 0.9998000
あとは普通に前々回の記事と同じ方法で傾向スコアを算出するだけです。
Platt's scalingされたSVMによる傾向スコアを算出してみる
普通にsvm {e1071}で以下のように回帰します。
> library(e1071) > model.svm.platt <- svm(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.platt) > ps.svm.platt <- predict(model.svm.platt, newdata=d.platt[,-1]) > hist(ps.svm.platt, breaks=100)
このままだと0未満や1より大きい傾向スコアがあるため、ATE / ATTの計算には使えません。そこで、便宜的にPlatt's scalingの下限値と上限値をそれぞれ代入して代用することにします。
> ps.svm.platt[ps.svm.platt<0] <- 1/(5856+2) > ps.svm.platt[ps.svm.platt>1] <- (4144+1)/(4144+2) # Calibration plot > library(caret) > df.platt <- data.frame(obs=factor(d$cm_dummy),svmplatt=1-ps.svm.platt,svmorg=ps.svm[,1]) > xyplot(calibration(obs~svmplatt+svmorg, df.platt))
元の傾向スコアよりも、より理想的な形に近付いていることが分かります。
Platt's scalingされたSVMによる傾向スコアでATE / ATTを算出してみる
後は機械的にATE / ATTを算出するだけです。
# ATE > iestp1.svm.platt <- (ivec1/ps.svm.platt) > iestp0.svm.platt <- (ivec0/(1-ps.svm.platt)) > iestp.svm.platt <- iestp1.svm.platt+iestp0.svm.platt > ipwe_gs.svm.platt <- lm(d$gamesecond ~ ivec+0, weights=iestp.svm.platt) > summary(ipwe_gs.svm.platt) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.svm.platt) Weighted Residuals: Min 1Q Median 3Q Max -472872 -6338 -337 -276 2341185 Coefficients: Estimate Std. Error t value Pr(>|t|) ivecivec1 6178.3 154.3 40.034 <2e-16 *** ivecivec0 269.5 218.7 1.232 0.218 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 79140 on 9998 degrees of freedom Multiple R-squared: 0.1383, Adjusted R-squared: 0.1381 F-statistic: 802.1 on 2 and 9998 DF, p-value: < 2.2e-16 # ATT > iestp1_ATT.svm.platt <- ivec1 > iestp0_ATT.svm.platt <- ivec0*ps.svm.platt/(1-ps.svm.platt) > iestp_ATT.svm.platt <- iestp1_ATT.svm.platt+iestp0_ATT.svm.platt > ipwe_ATT_gs.svm.platt <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.svm.platt) > summary(ipwe_ATT_gs.svm.platt) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.svm.platt) Weighted Residuals: Min 1Q Median 3Q Max -8796 -2478 -99 -30 271832 Coefficients: Estimate Std. Error t value Pr(>|t|) ivecivec1 2478.07 235.07 10.542 < 2e-16 *** ivecivec0 136.62 42.79 3.193 0.00141 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15130 on 9998 degrees of freedom Multiple R-squared: 0.01199, Adjusted R-squared: 0.01179 F-statistic: 60.66 on 2 and 9998 DF, p-value: < 2.2e-16
何じゃこりゃあ。。。。。。確かにCMダミーによる介入効果が正であると出てますが、この平均の推定値はおかしいでしょ(汗)。これまでとは違って想定通りにCMダミーによる介入効果「あり」とはちゃんと出るのですが、肝心の推定値のスケールには問題があるようです。
0未満1より大きい部分のスケールを変えてみる
もしかしたらSVMで推定した傾向スコアの0未満1より大きい部分のスケールのやり方が悪かったかな?と思ったので、単純に縮小させて押し込むような変換をかけてみました。
> ps.svm.platt.scale <- predict(model.svm.platt, newdata=d.platt[,-1]) > max(ps.svm.platt.scale) [1] 1.220675 > min(ps.svm.platt.scale) [1] -0.2071101 > cal.scale <- max(ps.svm.platt.scale) - min(ps.svm.platt.scale) > ps.svm.platt.mod <- (ps.svm.platt.scale - min(ps.svm.platt.scale))/cal.scale > hist(ps.svm.platt.mod, breaks=100) # Calibration plot > df.platt <- data.frame(obs=factor(d$cm_dummy),svmplatt=1-ps.svm.platt.mod,svmorg=ps.svm[,1]) > xyplot(calibration(obs~svmplatt+svmorg, df.platt)) # ATE > iestp1.svm.platt.mod <- (ivec1/ps.svm.platt.mod) > iestp0.svm.platt.mod <- (ivec0/(1-ps.svm.platt.mod)) > iestp.svm.platt.mod <- iestp1.svm.platt.mod+iestp0.svm.platt.mod > ipwe_gs.svm.platt.mod <- lm(d$gamesecond ~ ivec+0, weights=iestp.svm.platt.mod) > summary(ipwe_gs.svm.platt.mod) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.svm.platt.mod) Weighted Residuals: Min 1Q Median 3Q Max -25528 -4607 -3989 -3469 396909 Coefficients: Estimate Std. Error t value Pr(>|t|) ivecivec1 3957.6 281.7 14.05 <2e-16 *** ivecivec0 3147.7 271.0 11.61 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 25530 on 9996 degrees of freedom (2 observations deleted due to missingness) Multiple R-squared: 0.03217, Adjusted R-squared: 0.03198 F-statistic: 166.1 on 2 and 9996 DF, p-value: < 2.2e-16 # ATT > iestp1_ATT.svm.platt.mod <- ivec1 > iestp0_ATT.svm.platt.mod <- ivec0*ps.svm.platt.mod/(1-ps.svm.platt.mod) > iestp_ATT.svm.platt.mod <- iestp1_ATT.svm.platt.mod+iestp0_ATT.svm.platt.mod > ipwe_ATT_gs.svm.platt.mod <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.svm.platt.mod) > summary(ipwe_ATT_gs.svm.platt.mod) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.svm.platt.mod) Weighted Residuals: Min 1Q Median 3Q Max -9710 -2479 -2469 -1494 220617 Coefficients: Estimate Std. Error t value Pr(>|t|) ivecivec1 2478.7 218.9 11.32 <2e-16 *** ivecivec0 3224.3 256.5 12.57 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 14090 on 9996 degrees of freedom (1 observation deleted due to missingness) Multiple R-squared: 0.02783, Adjusted R-squared: 0.02764 F-statistic: 143.1 on 2 and 9996 DF, p-value: < 2.2e-16
今度はATEは良いものの、ATTがCMダミーの効果が負という結果になってしまいました。。。やっぱり難しいですね。キャリブレーションプロットもかなりガクガクしていて、あまり良くなさそうです。
SVMではなくxgboostで無理やりPlatt's scalingをロジットリンク関数にかませてスケールさせてみる
リンク先の資料には「Platt's scalingはそれ自体にはあまり大きな効果はないらしい」というコメントがあるんですが、個人的にはむしろPlatt's scalingした後の回帰のやり方に問題があるのかなと。。。これを無理くり最適化の段階で[0,1]にスケーリングさせるとかした方が良いのかもとは思います。
あと、そもそもこれって「回帰させる目的変数をPlatt's scalingさせる」「その上でロジットの範疇で回帰させる」代物なので、SVMでただ回帰させるのは間違ってるんじゃないかという気も。。。でもSVMで連続値を目的変数にしつつリンク関数をロジットというかソフトマックスにする方法がないんですよね。StackOverFlowを見ても「分類してprobability = TRUEにする」以外のやり方がないようで。
そこで、リンク関数にロジットやソフトマックスを設定できる{xgboost}でやってみようと思います。なお、分類精度が高過ぎると傾向スコアとしては使いづらくなる?ということがこのシリーズ記事で分かってきているので、あえてeta=0.1と汎化を強めにしています。
> library(xgboost) > library(Matrix) > dtrain.platt <- d.platt[,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.p <- sparse.model.matrix(cm_dummy~., data=dtrain.platt) > ddtrain.p <- xgb.DMatrix(d.mx.p, label=d.platt$cm_dummy) > model.gbdt.p <- xgb.train(params=list(objective="binary:logistic",eta=0.1), data=ddtrain.p, nrounds=20) > ps.gbdt.platt <- predict(model.gbdt.p, newdata=ddtrain.p) > hist(ps.gbdt.platt, breaks=100) > df.platt <- data.frame(obs=factor(d$cm_dummy), gbdtp=1-ps.gbdt.platt, gbdtorg=ps.gbdt[,1]) > xyplot(calibration(obs~gbdtp+gbdtorg,df.platt)) > iestp1.gbdt.platt <- (ivec1/ps.gbdt.platt) > iestp0.gbdt.platt <- (ivec0/(1-ps.gbdt.platt)) > iestp.gbdt.platt <- iestp1.gbdt.platt+iestp0.gbdt.platt > ipwe_gs.gbdt.platt <- lm(d$gamesecond ~ ivec+0, weights=iestp.gbdt.platt) > summary(ipwe_gs.gbdt.platt) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.gbdt.platt) Weighted Residuals: Min 1Q Median 3Q Max -7710 -4069 -3507 -3234 405515 Coefficients: Estimate Std. Error t value Pr(>|t|) ivecivec1 3085.4 264.5 11.67 <2e-16 *** ivecivec0 2924.9 240.2 12.18 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 22590 on 9998 degrees of freedom Multiple R-squared: 0.02766, Adjusted R-squared: 0.02746 F-statistic: 142.2 on 2 and 9998 DF, p-value: < 2.2e-16 > iestp1_ATT.gbdt.platt <- ivec1 > iestp0_ATT.gbdt.platt <- ivec0*ps.gbdt.platt/(1-ps.gbdt.platt) > iestp_ATT.gbdt.platt <- iestp1_ATT.gbdt.platt+iestp0_ATT.gbdt.platt > ipwe_ATT_gs.gbdt.platt <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.gbdt.platt) > summary(ipwe_ATT_gs.gbdt.platt) Call: lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.gbdt.platt) Weighted Residuals: Min 1Q Median 3Q Max -6260 -2478 -2252 -1241 220618 Coefficients: Estimate Std. Error t value Pr(>|t|) ivecivec1 2478.1 209.4 11.84 <2e-16 *** ivecivec0 2566.7 246.5 10.41 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 13480 on 9998 degrees of freedom Multiple R-squared: 0.02425, Adjusted R-squared: 0.02405 F-statistic: 124.2 on 2 and 9998 DF, p-value: < 2.2e-16
やっぱりATEはCMダミーに正の効果が出るものの、ATTではそうはなっていません。単にロジットの範囲でスケールさせるだけではうまくいかないようです。