ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

ROC解析で識別能を評価する(pROCパッケージ, ggplot系ggroc関数)[R]

検査を使って真の状態(例えばある疾患に罹患しているかどうか)や予後を予測する際に、その検査によって真の状態(罹患 vs. 未罹患)を区別して言い当てる能力を識別能(discrimination ability)という。
次のような指標で評価する。

  • 感度(sensitivitiy):疾患がある人の中で、検査が陽性になる人の割合
  • 特異度(specificity):疾患がない人の中で、検査が陰性になる人の割合
  • 陽性的中率(positive predictive value):検査が陽性の人の中で、疾患がある人の割合
  • 陰性的中率(negative predictive value):検査が陰性の人の中で、疾患がない人の割合

検査が連続値を返す場合は、どのようなカットオフ値で2値化するかによって上記の指標の値が変わる。
様々なカットオフ値で2値化した感度・特異度をプロットしたものがROC曲線(receiver operating characteristic curve)である。

図1:カットオフ値を変化させてROC曲線を描く

真陽性者の増加に対して偽陽性者の増加が小さいほどROC曲線は左上角に近づく。

パッケージとデータの準備

まず必要なパッケージを読み込む。

library(pROC)
library(ggplot2)

今回はpROCパッケージ内のaSAHデータセットを使う。
これはクモ膜下出血(subarachnoid hemorrhage, SAH)で入院した患者の入院時の血漿マーカーが予後を予測できるかを調べた研究のデータ*1

> data(aSAH)
> str(aSAH)
'data.frame':  113 obs. of  7 variables:
 $ gos6   : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 5 5 5 5 1 1 4 1 5 4 ...
 $ outcome: Factor w/ 2 levels "Good","Poor": 1 1 1 1 2 2 1 2 1 1 ...
 $ gender : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 1 1 1 2 2 ...
 $ age    : int  42 37 42 27 42 48 57 41 49 75 ...
 $ wfns   : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 3 2 5 4 1 2 ...
 $ s100b  : num  0.13 0.14 0.1 0.04 0.13 0.1 0.47 0.16 0.18 0.1 ...
 $ ndka   : num  3.01 8.54 8.09 10.42 17.4 ...
  • outcome:6か月時点のGlasgow Outcome Scale(GOS)をもとに、1-3点を"Good"、4-5点を"Poor"としたもの。
  • s100b:血漿中s100β濃度
  • ndka:血漿中nucleotide diphosphate kinase A濃度

ROC曲線下面積(AUC)を計算する

予測したい状態(2値変数)をresponse、予測に使う検査値(連続変数)をpredictorに指定する。
以下のどちらでもよい。

  • roc(response ~ predictor)
  • roc(response, predictor)

前者の記法では予測に使いたい変数を+で追加していくことでROC解析のリスト化できるので、複数のROC曲線を一括で作成したいときに便利(線形回帰モデルのときと+の意味が違うので注意)。

信頼区間も計算したい場合はci=TRUEとする。
信頼区間の算出に使用する方法は引数ci.methodで指定可。デフォルトは="delong"だが、="bootstrap"でブートストラップ法も使える。

roc_list <- roc(outcome ~ age + s100b + ndka, data=aSAH, ci=TRUE)
> roc_list
$age

Call:
roc.formula(formula = outcome ~ age, data = aSAH, ci = TRUE)

Data: age in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.615
95% CI: 0.5082-0.7219 (DeLong)

$s100b

Call:
roc.formula(formula = outcome ~ s100b, data = aSAH, ci = TRUE)

Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.7314
95% CI: 0.6301-0.8326 (DeLong)

$ndka

Call:
roc.formula(formula = outcome ~ ndka, data = aSAH, ci = TRUE)

Data: ndka in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.612
95% CI: 0.5012-0.7227 (DeLong)

ggroc( )を使ってROC曲線を描く

さっきroc( )で作成したオブジェクトを描画関数に渡す。
グラフについては次のような引数を指定可能。

  • aes:予測変数毎に変更したいROC曲線の色や種類を指定する。共通のものを指定したいときはaesの外で指定。
    • "colour":色
    • "linetype":線種
    • "alpha":透過度
    • "size":線の太さ。予測変数ごとに変えると太くなりすぎるので非推奨。例ではaesの外で指定した。
  • legacy.axes:x軸を偽陽性割合(=1-特異度)としたいときは=TRUE、特異度にしたいときは=FALSEとする。
ggroc(roc_list, 
      aes=c("linetype","color"),
      size=1,
      legacy.axes = TRUE) +
geom_abline(color="dark grey", size=0.5)

図2:ggroc( )を使ってROC曲線を描いた

パイプ演算子%>%を使うことも可能

%>%を使ってggroc( )に渡すことも可能。
ただし、前述の~+には対応していないみたい。

library(tidyverse)
aSAH %>% roc(outcome, s100b) %>% ggroc()

2つのAUCを比較する

リスト化されたROCオブジェクトでは上手くいかないので、2つのROCオブジェクトを別々に作成する

roc1 <- roc(outcome~s100b, data=aSAH)
roc2 <- roc(outcome~ndka, data=aSAH)

methodで手法を変更可能。デフォルトは="delong"だが、="bootstrap"でブートストラップ法も使える。

> roc.test(roc1, roc2, method="bootstrap")
  |==============================================| 100%

    Bootstrap test for two correlated ROC curves

data:  roc1 and roc2
D = 1.4028, boot.n = 2000, boot.stratified = 1, p-value = 0.1607
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2 
  0.7313686   0.6119580 

おわりに

  • キレイな図を描きたいときはggplot系を探せば大抵そろっている。

参考資料

  • CRANの説明書

search.r-project.org

*1:Turck N, Vutskits L, Sanchez-Pena P, Robin X, Hainard A, Gex-Fabry M, Fouda C, Bassem H, Mueller M, Lisacek F, Puybasset L, Sanchez JC. A multiparameter panel method for outcome prediction following aneurysmal subarachnoid hemorrhage. Intensive Care Med. 2010 Jan;36(1):107-15.