pROC

パッケージの概要

pROCパッケージは、ROC曲線（Receiver Operating Characteristic Curve）の描画や分析に特化して作られたパッケージです。ROC曲線は、二値変数である目的変数（アウトカム変数）と、連続変数である説明変数の関係の強さを評価する手法であり、医療分野における診断マーカーや、より一般には二値分類モデルの精度指標として用いられています。

ROC曲線を作成する

データフレームにroc()関数を適用することで、rocオブジェクト（または、rocオブジェクトを要素とするリスト）を作成することができます。rocオブジェクトに対してplot()関数やggroc()関数を適用することで、ROC曲線を描画することができます。

library(ggplot2)
library(pROC)

response <- c(0, 0, 0, 1, 0, 1, 1, 1)
probability <- 1:8 /10
(sample.roc <- roc(response, predictor = probability))

Call:
roc.default(response = response, predictor = probability)

Data: probability in 4 controls (response 0) < 4 cases (response 1).
Area under the curve: 0.9375

ggroc(sample.roc)

なお、ROC曲線は、連続変数を分類に用いるときの「しきい値」を変化させたときの、感度と特異度の変化を平面上にプロットしたものです。

たとえば、予測確率（10%~80%）の値がしきい値 t 以上であるときに response の値を 1 と予測することにすると、t = 10% （すべてを 1 と予測）のとき感度 = 1、特異度 = 0 となり、これはグラフの右上の点に対応しています。t を 20%、30%、…と変化させていったときの点を明示すれば、下のようになります。

df <- data.frame(sensitivity = 1 - cumsum(c(0, response)) / 4,
                 specificity = cumsum(c(0, 1 - response)) / 4,
                 label = paste0("t=", 1:9 * 10, "%"))

# ggplot2 ベースでROCカーブを描画する
ggroc(sample.roc) +
  geom_label(aes(specificity, sensitivity, label = label), df)

AUCを計算する

ROC曲線のグラフにおいて、グラフの下の部分の面積をAUC（Area Under the Curve）と呼びます。AUCは、しばしば分類モデルの性能を評価する指標として用いられます。

cat(auc(sample.roc))

0.9375

# 出力される数値は"auc"クラスのオブジェクトでもあり、表示方法が設定されている。
auc(sample.roc)

Area under the curve: 0.9375

aSAHデータセット

pROCパッケージに収録されているaSAHデータセットは、動脈瘤性くも膜下出血患者113名について、以下の情報を記録したものです。

変数名	概要
gos6	入院6か月後におけるGlasgow Outcome Scaleの値
outcome	予後（gos6が3以下かどうかで分類）
gender	性別
age	年齢
wfns	入院時のWFNS分類（くも膜下出血の重症度分類の一種）
s100b	入院後12時間以内に採取された静脈血中のS100β濃度
ndka	静脈血中のNDKA（ヌクレオシド二リン酸キナーゼA）の濃度

data(aSAH, package = "pROC")
DT::datatable(aSAH)

合計113件のデータのうち、41件が予後良好（Good）、72件が予後不良（Poor）です。これを、データの他の指標から判定ないし予測できるかどうかを、ROC曲線で表現します。

ROCs <- roc(outcome ~ age + s100b + ndka, data = aSAH)

# plot()関数による描画
plot(ROCs$s100b)

# ggroc()関数によるggplot2ベースの描画
ggroc(ROCs, legacy.axes = FALSE)

各 ROC 曲線について AUC（ROC-AUC）の値は以下の通りとなり、くも膜下出血の予後の指標としては s100b の精度が他の指標よりも良いことが読み取れます。

# auc()関数をrocオブジェクトにまとめて適用する
sapply(ROCs, auc)

      age     s100b      ndka 
0.6150068 0.7313686 0.6119580