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

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

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

統計的因果推論(5): Platt's scalingで機械学習分類器による傾向スコアを調整してみる

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

機械学習分類器で算出した傾向スコアを調整する話ですが、最後に課題として残ったのがprobability calibrationによる実践。探してみると前回の記事でもやったisotonic regressionとか色々出てくるんですが、もう一つ出てくるのがPlatt's scaling。これはあのSVMのSMOを提案したPlatt*1がprobability calibrationのために提案しているもので、web上にも幾つか資料があるようです。


やってることは極めて単純で、最初の分類確率値の計算をする際に0 / 1の二値への分類モデリングをやるのではなく、

y_+ = \frac{N_+ + 1}{N_+ + 2}, ~ y_- = \frac{1}{N_- + 2}

という別の二値に変換した上で(ただしN_+は正例の数でN_-は負例の数)、これを回帰したモデルで分類確率値を予測させるというだけです。


ということで、サクッとやってみましょう。ただし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)

f:id:TJO:20161212130108p:plain


このままだと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))

f:id:TJO:20161212125853p:plain

元の傾向スコアよりも、より理想的な形に近付いていることが分かります。


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

f:id:TJO:20161212143453p:plain

f:id:TJO:20161212143510p:plain

今度は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.gdbt.p <- xgb.train(params=list(objective="binary:logistic",eta=0.1), data=ddtrain.p, nrounds=20)
> ps.gdbt.platt <- predict(model.gdbt.p, newdata=ddtrain.p)
> hist(ps.gdbt.platt, breaks=100)
> df.platt <- data.frame(obs=factor(d$cm_dummy), gdbtp=1-ps.gdbt.platt, gdbtorg=ps.gdbt[,1])
> xyplot(calibration(obs~gdbtp+gdbtorg,df.platt))
> iestp1.gdbt.platt <- (ivec1/ps.gdbt.platt)
> iestp0.gdbt.platt <- (ivec0/(1-ps.gdbt.platt))
> iestp.gdbt.platt <- iestp1.gdbt.platt+iestp0.gdbt.platt
> ipwe_gs.gdbt.platt <- lm(d$gamesecond ~ ivec+0, weights=iestp.gdbt.platt)
> summary(ipwe_gs.gdbt.platt)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp.gdbt.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.gdbt.platt <- ivec1
> iestp0_ATT.gdbt.platt <- ivec0*ps.gdbt.platt/(1-ps.gdbt.platt)
> iestp_ATT.gdbt.platt <- iestp1_ATT.gdbt.platt+iestp0_ATT.gdbt.platt
> ipwe_ATT_gs.gdbt.platt <- lm(d$gamesecond ~ ivec+0, weights=iestp_ATT.gdbt.platt)
> summary(ipwe_ATT_gs.gdbt.platt)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT.gdbt.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

f:id:TJO:20161212152100p:plain

f:id:TJO:20161212152018p:plain

やっぱりATEはCMダミーに正の効果が出るものの、ATTではそうはなっていません。単にロジットの範囲でスケールさせるだけではうまくいかないようです。


感想など


これって結局普通のロジスティック回帰による傾向スコアとの対比で言うと「ちゃんとロジットを使ってロジットっぽい分布になるように制約をかけているものでないとダメ」ということなんじゃないかな、と思いました。当たり前の話なんですが、機械学習分類器だとその辺はうまく制約をかけられない*2ので、傾向スコアを算出させたらうまくいかなくて当然なんじゃないかと。


と、色々言ってますが僕に別にアイデアがあるわけではありませんので、そのうちまた気が向いたら調べてみますということで。。。お後がよろしいようで。

*1:昨年移籍して現在同じ会社にいるのだと知りました笑

*2:というかそれだとそもそも機械学習分類器としてのパフォーマンスが出ない