不均衡データの分類をクラス重み付けではなくクラス分類事後確率の閾値で補正するとどうなるか

先日ask.fmでこんな質問をいただいたのでやってみました。

不均衡データの分類についてブログを拝見しました。 不躾な質問で恐縮ですが、正例の少ない不均衡データをRandomforestで2値分類を行う際に、ウェイトを使うのであれば、単純にProbで出力される｢正例である確率｣の閾値を下げる、というアイディアでも良いのではないかと思うのですが、何かよい情報ご存知ではないでしょうか？ 諸々の精度指標を元にカットオフを決められるので、問題なければ便利かなぁと思うのですが

ask.fm

もちろん以下の過去記事に関連する話題だと思うのですが、折角なので改めて自分でちろっとRで回してみた次第です。

とりあえずサクッとやった感じがこちら。用いた不均衡データセットはこちらから。ちなみにこの手の計算をする時にRはちょっと不親切なので、極めて不本意ですがfor文を使ってます*1。まずはランダムフォレストから。

> d<-read.delim("xorc_ub.txt")
> d$label<-as.factor(d$label)
> 
> library(randomForest)
> d.rf<-randomForest(label~.,d)
> d.rf.ub<-randomForest(label~.,d,classwt=c(1,30))
> 
> px<-seq(-4,4,0.03)
> py<-seq(-4,4,0.03)
> pgrid<-expand.grid(px,py)
> names(pgrid)<-names(d)[-3]
> 
> out.d.rf<-predict(d.rf,newdata=pgrid,type='response')
> out.d.rf.ub<-predict(d.rf.ub,newdata=pgrid,type='response')
> out.d.rf.prob<-predict(d.rf,newdata=pgrid,type='prob')
> 
> out.d.rf.prob.grid<-rep(0,length(out.d.rf))
> for (i in 1:length(out.d.rf)) {
+     if (out.d.rf.prob[i,1]>30/31) out.d.rf.prob.grid[i]<-0
+     if (out.d.rf.prob[i,2]>1/31) out.d.rf.prob.grid[i]<-1
+ }
> 
> par(mfrow=c(2,2))
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='Without class weights')
> par(new=T)
> contour(px,py,array(out.d.rf,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F)
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='With class weights')
> par(new=T)
> contour(px,py,array(out.d.rf.ub,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F)
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='With prob threshld changed')
> par(new=T)
> contour(px,py,array(out.d.rf.prob.grid,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F)

おー、クラス分類事後確率の閾値をいじったやつ（3番目）は明らかに何か変ですね。同じことをガウシアンカーネルSVMでもやってみましょう。

> library(e1071)
> d.svm<-svm(label~.,d,probability=T)
> d.svm.ub<-svm(label~.,d,class.weights=c('0'=1,'1'=30),probability=T)
> 
> out.d.svm<-predict(d.svm,newdata=pgrid)
> out.d.svm.ub<-predict(d.svm.ub,newdata=pgrid)
> out.d.svm.tmp<-predict(d.svm,newdata=pgrid,probability=T)
> out.d.svm.prob<-attr(out.d.svm.tmp,"probabilities")
> 
> out.d.svm.prob.grid<-rep(0,length(out.d.svm.prob))
> for (i in 1:length(out.d.svm.prob.grid)) {
+     if (out.d.svm.prob[i,1]>30/31) out.d.svm.prob.grid[i]<-0
+     if (out.d.svm.prob[i,2]>1/31) out.d.svm.prob.grid[i]<-1
+ }
> 
> par(mfrow=c(2,2))
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='Without class weights')
> par(new=T)
> contour(px,py,array(out.d.svm,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F)
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='With class weights')
> par(new=T)
> contour(px,py,array(out.d.svm.ub,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F)
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='With prob threshld changed')
> par(new=T)
> contour(px,py,array(out.d.svm.prob.grid,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F)

やっぱり何かが変です。最後にロジスティック回帰。

> d.glm<-glm(label~.,d,family=binomial)
> d.glm.ub<-glm(label~.,d,family=binomial,weights=c(rep(1,600),rep(30,20)))
> 
> out.d.glm<-predict(d.glm,newdata=pgrid,type='response')
> out.d.glm.ub<-predict(d.glm.ub,newdata=pgrid,type='response')
> out.d.glm.prob<-out.d.glm.ub
> out.d.glm.prob<-ifelse(out.d.glm.prob>1/31,1,0)
> 
> par(mfrow=c(2,2))
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='Without class weights')
> par(new=T)
> contour(px,py,array(out.d.glm,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F,levels=0.5)
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='With class weights')
> par(new=T)
> contour(px,py,array(out.d.glm.ub,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F,levels=0.5)
> plot(d[,-3],pch=19,col=c(rep('blue',600),rep('red',20)),cex=2,xlim=c(-4,4),ylim=c(-4,4),main='With prob threshld changed')
> par(new=T)
> contour(px,py,array(out.d.glm.prob,dim=c(length(px),length(py))),xlim=c(-4,4),ylim=c(-4,4),col="purple",lwd=3,drawlabels=F,levels=0.5)

アヒャヒャヒャヒャヒャ、これ全然ダメじゃないですかｗｗ　事後確率の閾値をいじってみたら、思いっきり正例の端っこの方をつかんじゃってますね。

ということで、現段階では「クラス分類事後確率の閾値をいじっても適切な結果にはならないらしい」ことが分かりました。そもそも論として、事後確率の閾値をずらしたことでランダムフォレストとロジスティック回帰ではfalse positiveが物凄く増える結果になってしまっているというのが3つの図からはっきり見て取れると思います。では、何故そんなことになったんでしょうか？というのは理論的な話をしなければ分からないので、ひとまずはじパタの当該箇所を復習中。

はじめてのパターン認識

pp.47-48の辺りを読みつつ、(4.32)式の変形を事後確率の閾値をいじるように改変してやることで何かが分かるはずなんですが、ちょっと今はその余裕がないので続きはまた後日改めてということで。

*1:applyファミリーで普通にいけそうな気もしますが