RでROC曲線を描く方法
定義
エラーマトリックスのFalse Positive RateとTrue Positive Rateをそれぞれ軸にとって描かれた図をROC曲線Receiver Operating Characteristic curveという。
説明
ROC曲線はモデルのパフォーマンスを一目で示すだけでなく、最適なカットオフを見つけるためや、モデル間の比較にも使われるなど、多くの場面で役立ちます。以下の例を通じて、RでROC曲線を描き、その意味を理解しましょう。重要なパッケージとしてROCRが使われます。
実践
Carパッケージの組み込みデータであるChileデータを読み込んでみましょう。
Chileはチリの独裁者アウグスト・ピノチェトが1988年に自身の政権延長の是非を問う投票に対する調査データで、region（地域）、population（回答者コミュニティの人口）、sex（性別）、age（年齢）、education（教育水準）、income（収入）、statusquo（現状支持度）に基づくvote（投票傾向）を把握できます。
投票傾向はA（棄権）、U（保留）、N（反対）、Y（賛成）の4つのクラスを持ち、ロジスティック回帰分析を使用するために、Y以外はすべてNに変えました。このポストではYを「賛成」、Nを「反対」とします。
モデルのパフォーマンスを確認するために、トレーニングデータとテストデータに分けました。
トレーニングデータを使って得たロジスティック回帰モデルは上記の通りです。次に、このモデルにテストデータを入力して有権者が賛成する確率を計算してみましょう。predict()関数にnewdata=testオプションを指定すると、モデルにテストデータを入力して自動的に計算してくれます。この時、type="response"を指定する必要があることに注意してください。これにより、確率の形で出力されます。
out0<-glm(vote~.,family=binomial(),data=train); summary(out0)
p <- predict(out0, newdata=test, type="response"); p
出力結果は次の通りです。
問題は、この確率をどのように解釈し、使用するかです。実際に計算したいのは確率そのものではなく、どの人が「賛成」したのかを当てることです。したがって、一定の確率以上の場合に「賛成」と予測し、低い場合に「反対」と予測したいのです。この「一定の確率」をカットオフcutoffまたは閾値thresholdといい、一般的にはカットオフという表現がよく使われます。
最適カットオフ
このカットオフの設定によって正分類率は変わってきます。極端に言えば、0.99以下をすべて「反対」と予測しても、「反対」をした人だけは正解するので、正分類率はある程度高くなります。逆に、0.01以上をすべて「賛成」と予測しても同様にある程度正解します。例えば、1年365日のうち120日が雨の日がある国では、1年中次の日の天気を「晴れ」と予測しても正解率は約$67\%$になります。
我々はより良い分析と予測を必要としているため、良いカットオフを見つけたいのです。単純に考えると、0.5を基準にして高ければ「賛成」、低ければ「反対」と予測すれば良いように思えます。しかし、計算してみると必ずしもそうではないことが多いです。（また、常に正分類率だけが基準になるわけではありません。状況によっては、他の重要な指標があるかもしれません。）そこで、この良いカットオフを見つけるために行うのが、上記で得たpのすべての確率をカットオフとして試してみることです。0.89572040でエラーマトリックスを作成し、False Positive RateとTrue Positive Rateを計算し、次に0.81810785で計算し、次に0.70215561で計算し、… これを繰り返して、横軸をFalse Positive Rate、縦軸をTrue Positive RateとしてプロットしたグラフがROCカーブです。
言葉では難しいですが、コード上では簡単なので心配しないでください。エラーマトリックスを作成し、必要な数値を計算するためには、次のようなコードを実行します。実際に出力結果を見る必要はないので、各関数の役割だけを簡単に説明します。prediction()関数は上記で計算した確率pと実際のテストデータのtest$voteを比較して分類率を計算してくれます。performance()関数は上記で計算したエラーマトリックスの数値prから必要なデータを抽出し、plot()関数に渡すことでROCカーブを描くためのデータを返してくれます。
pr <- prediction(p, test$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Test Data')
上記のコードを実行すると、次のようにROCカーブが描かれます。
このカーブは、すべてのカットオフについてFPRとTPRを計算し、それをそれぞれ$x$軸の座標、$y$軸の座標とするカーブです。カーブは可能な限りその下の面積が広い方が良いです。軌道としては、カーブが左上の点$(0,1)$に近づくような形が望ましいです。カットオフ自体は図に表示されませんが、データ構造を調べるとalpha.valuesを通じて参照できることがわかります。
str()関数で調べると、カットオフが0.939の時、TPRが0、FPRが0.00556であることがわかります。FPRが0.00556ということは、「反対」を誤って予測した割合が$0.5%$しかないことを意味します。ここまでは良さそうに見えますが、TPRが0なので「賛成」を正確に予測したケースが一件もないことになります。直感的に考えると、「賛成」の確率が0.939以下の人をすべて「反対」と予測したため、「賛成」を「賛成」と予測するハードルが高すぎたのです。このため、「反対」はすべて正解しても「賛成」は正解できないのです。
グラフが角ばって見える
自分で描いたROCカーブが角ばっていておかしいと感じる人がいます。これは単にテストデータが少ないために起こることです。データが少ないこと自体は心配ですが、角ばった形になることは全く心配いりません。例えば、この例でもテストデータの数を20に減らすと次のような形になります。
コード
以下は、例の全体のコードです。
install.packages("car")
install.packages("ROCR")
library(car)
library(ROCR)
set.seed(150421)
?Chile
str(Chile)
nrow(Chile)
head(Chile); tail(Chile)
DATA<-na.omit(Chile)
DATA$vote[DATA
$vote!='Y']<-'N'
DATA$vote<-factor(DATA$vote)
head(DATA); tail(DATA)
DATANUM<-nrow(DATA)
train<-sample(1:DATANUM)<(DATANUM*0.8)
test<-DATA[!train,]; head(test)
train<-DATA[train,]; head(train)
out0<-glm(vote~.,family=binomial(),data=train); summary(out0)
p <- predict(out0, newdata=test, type="response"); p
pr <- prediction(p, test$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Test Data')
str(prf)
smalltest<-test[sample(nrow(test),20),]
p <- predict(out0, newdata=smalltest, type="response"); p
pr <- prediction(p, smalltest$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Small Test Data')