ROCカーブを使用して最適なカットオフを見つける方法
概要
ROC曲線を描くことで、トレーニングデータを使って得たモデルがテストデータをどれだけよく説明しているか一目で把握できるため、非常に有用だ。しかし、この曲線はすべてのカットオフに対する分類率を計算して求めたものであるため、結局のところ『どのカットオフで0と1を分類するか』はわからない。これを知るためには、クロスバリデーションの方法論を応用してみよう。
検証データ
トレーニングデータから0と1を最もよく分類する最適カットオフを見つけたとしても、それはトレーニングデータをよく説明するカットオフに過ぎない。当然、テストデータで分類率を計算しても、その最適カットオフはテストデータのみをよく説明するカットオフに過ぎない。したがって、バリデーションデータvalidation Dataというものを別に作成し、どのデータにも偏りのないカットオフを見つける。
- Step 1. トレーニングデータを学習させてモデルを作る。
- Step 2. バリデーションデータに適用して最適カットオフを見つける。
- Step 3. トレーニングデータで得たモデルをバリデーションデータで得たカットオフで分類し、テストデータでどれだけパフォーマンスが高いか確認する。複数のモデル候補がある場合、このパフォーマンスが最も高いモデルを最終モデルとして選択する。
この時、Step 1で得たモデルはStep 2で得たカットオフを含んで初めて完全なモデルとなる。
実習
(ROC曲線の描き方に続いて)
上の図は前のポストでデータをトレーニングデータとテストデータに分けた時のROC曲線だ。最適カットオフを求めるためには、データを三分割する必要がある。この時の比率はデータによって異なるが、特に問題がなければ3:1:1程度が無難だ。
DATANUM<-nrow(DATA)
numeric(DATANUM)
DATANUM*c(0.6,0.2,0.2)
slicing<-sample(1:DATANUM)
slicing[slicing>(DATANUM*0.8)]<-(-1)
slicing[slicing>(DATANUM*0.6)]<-(-2)
slicing[slicing>0]<-(-3)
slicing<-slicing+4
train<-DATA[slicing==1,]; head(train)
valid<-DATA[slicing==2,]; head(valid)
test<-DATA[slicing==3,]; head(test)
データの前処理を終え、上記のコードを実行すると、以下のようにトレーニング、バリデーション、テストの3つのデータに分割される。
Step 1.
トレーニングデータを学習させてモデルを作る。
out0<-glm(vote~.,family=binomial(),data=train); summary(out0)
vif(out0)
out1<-step(out0, direction="both"); summary(out1)
qchisq(0.95,df=1454)
vif(out1)
上記のコードを実行すると、変数選択手順を経てモデルを作成し、適合度検定や多重共線性をチェックする。
モデルを確認すると、特に問題点は見当たらない。
Step 2.
バリデーションデータに適用して最適カットオフを見つける。
p <- predict(out1, newdata=valid, type="response")
pr <- prediction(p, valid$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Validation Data')
上記のコードを実行すると、バリデーションデータをテストデータのように扱い、以下のようにROC曲線を描く。
最適カットオフはデータと目的に応じて異なるが、特に重視する点がなければ、左上の$(0,1)$に最も近い点を探し、その点のカットオフを最適カットオフとする。距離を計算する必要があるため、コードは少し複雑になる。
optid<-(1:length(prf@y.values[[1]][-1]))[((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-11])^2)
==min((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-1])^2)]
points(prf@x.values[[1]][-1][optid],prf@y.values[[1]][-1][optid], col='red', pch=15)
optcut<-prf@alpha.values[[1]][-1][optid]; optcut
上記のコードを実行すると、説明通りに$(0,1)$で最も近い点を赤で示し、その点のカットオフを出力する。
コードは非常に複雑だが、理解する必要はない。複雑なものと難しいものは別の問題だ。上記のコードは長いだけで、概念的には全く難しくない。曲線上のすべての点と$(0,1)$の距離を測り、その距離が最も短い点を選んだだけだ。その点での$alpha.valuesを参照すればカットオフがわかる。このようにして求めたカットオフは、データを最も適切に分類するカットオフとして受け入れることができる。（再度強調するが、これは絶対的な尺度ではない。使用者の目的によって『最適なカットオフ』自体が完全に新たに定義されることがある。）
この例で得た最適カットオフは$0.4564142$であり、これより高ければ1と判定し、低ければ0と判定するのが最も良いと考えて構わない。（すでに3回目の強調だが、最も良いと考えられるものであり、これが最も良いものではない。これに関して合理的な解釈を出すことは、全て分析者に委ねられている。）
Step 3.
テストデータでどれだけ正確に分類できるか確認する。
p <- predict(out1, newdata=test, type="response"); head(p,48)
table(test$vote, p>optcut)
上記のコードを実行すると、テストデータで確率を計算し、最適カットオフに基づくエラーマトリックスを出力する。
上記エラーマトリックスの正分類率は約$81 \%$であり、相当に使え、分析者が満足するならば最終モデルとして受け入れるべきだと考えられる。気づいたかもしれないが、厳密にいえば、最適カットオフを求めるためには必ずしもROC曲線を描く必要はない。計算のためのデータはデータフレームとしてすでに用意されているため、コードをうまく回して値を得るだけでも全く問題ない。
コード
以下は例のコードだ。
install.packages("car")
install.packages("ResourceSelection")
install.packages("ROCR")
library(car)
library(ResourceSelection)
library(ROCR)
set.seed(150421)
?Chile
str(Chile)
nrow(Chile)
DATA<-na.omit(Chile)
DATA$vote[DATA$vote!='Y']<-'N'
DATA$vote<-factor(DATA$vote)
DATANUM<-nrow(DATA)
numeric(DATANUM)
DATANUM*c(0.6,0.2,0.2)
slicing<-sample(1:DATANUM)
slicing[slicing>(DATANUM*0.8)]<-(-1)
slicing[slicing>(DATANUM*0.6)]<-(-2)
slicing[slicing>0]<-(-3)
slicing<-slicing+4
train<-DATA[slicing==1,]; head(train)
valid<-DATA[slicing==2,]; head(valid)
test<-DATA[slicing==3,]; head(test)
out0<-glm(vote~.,family=binomial(),data=train); summary(out0)
vif(out0)
out1<-step(out0, direction="both"); summary(out1)
qchisq(0.95,df=1454)
vif(out1)
p <- predict(out1, newdata=valid, type="response")
pr <- prediction(p, valid$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Validation Data')
optid<-(1:length(prf@y.values[[1]][-1]))[((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-11])^2)
==min((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-1])^2)]
points(prf@x.values[[1]][-1][optid],prf@y.values[[1]][-1][optid], col='red', pch=15)
optcut<-prf@alpha.values[[1]][-1][optid]; optcut
p <- predict(out1, newdata=test, type="response"); head(p,48)
table(test$vote, p>optcut)