検査を使って真の状態(例えばある疾患に罹患しているかどうか)や予後を予測する際に、その検査によって真の状態(罹患 vs. 未罹患)を区別して言い当てる能力を識別能(discrimination ability)という。
次のような指標で評価する。
- 感度(sensitivitiy):疾患がある人の中で、検査が陽性になる人の割合
- 特異度(specificity):疾患がない人の中で、検査が陰性になる人の割合
- 陽性的中率(positive predictive value):検査が陽性の人の中で、疾患がある人の割合
- 陰性的中率(negative predictive value):検査が陰性の人の中で、疾患がない人の割合
検査が連続値を返す場合は、どのようなカットオフ値で2値化するかによって上記の指標の値が変わる。
様々なカットオフ値で2値化した感度・特異度をプロットしたものがROC曲線(receiver operating characteristic curve)である。
真陽性者の増加に対して偽陽性者の増加が小さいほど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)
パイプ演算子%>%を使うことも可能
%>%
を使って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の説明書
*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.