【R】ROC曲線を書く

備忘録的なやつです。

Rを使ってROC曲線を描きます。

使う関数は**roc()**と**plot()**です。

#ROC曲線

ROC曲線については以下のURLを参照。

[URL](http://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub6.html)

#必要なパッケージ

必要なパッケージは**pROC**。

**pROC**パッケージをインストールしておきます。

```R
install.packages("pROC")
```

パッケージを読み込む

```R
library(pROC)
```

#データの読み込みと確認

使用するデータは**pROC**の中にある**aSAH**という113例の動脈瘤破裂のSAHのデータ。

```R
df <- aSAH
```

読み込んだデータを確認する。

```R
summary(df)
```
⬇実行結果

```R
> summary(df)
 gos6   outcome      gender        age       wfns       s100b            ndka       
 1:28   Good:72   Male  :42   Min.   :18.0   1:39   Min.   :0.030   Min.   :  3.01  
 2: 0   Poor:41   Female:71   1st Qu.:42.0   2:32   1st Qu.:0.090   1st Qu.:  9.01  
 3:13                         Median :51.0   3: 4   Median :0.140   Median : 12.22  
 4: 6                         Mean   :51.1   4:16   Mean   :0.247   Mean   : 19.66  
 5:66                         3rd Qu.:61.0   5:22   3rd Qu.:0.330   3rd Qu.: 17.30  
                              Max.   :81.0          Max.   :2.070   Max.   :419.19
```

データ構造を見てみます。

gos6：脳損傷のある患者の回復を5段階で評価し、高いほどよい状態。
outcome:gos6が4,5をGood、それ以外はPoor。
wfns：脳神経外科連合(WFNS)による分類。らしい
s100b：血清マーカー
ndka：血清マーカー

#ROC曲線の描画
ROC曲線の描画には**pROC**パッケージの**roc()**を使います。

まず、**roc()**の結果をオブジェクトに格納しておきます。

``roc(Y ~ X, data, ci)``

``Y``：疾患の有無（2値変数）
``X`` ：検査の変数
``data``：データフレーム
``ci``：信頼区間

今回は``Y``がoutocome、``X``をndkaとします。

```R
ROC <- roc(outcome ~ ndka, data = df, ci = TRUE)
ROC
```
⬇実行結果

```R
> ROC <- roc(outcome ~ ndka, data = df, ci = TRUE)
Setting levels: control = Good, case = Poor
Setting direction: controls < cases

> ROC

Call:
roc.formula(formula = outcome ~ ndka, data = df, 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)
```

Area under the curve (AUC)は0.612（95%信頼区間：0.5012-0.7227）

**plot()**関数で描画します。

```R
plot(ROC)
```
⬇描画結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/0fc319c0-3b94-f312-1380-7dd80c478ecd.png)

ROCオブジェクトを**plot()**関数でただ返すとx軸の左が1、右が0となります。

これは一般的ではないので、次の項で調整していきます。

#カットオフ値の決定
カットオフ値の決め方は、1）Yoden index、2）左上隅との距離が最小となる点を求める方法があります。

詳細はさきほどのURLを参照。

``plot(ROC, identity, print.thres, print.thres.best.method, legacy.axes)``

``ROC``：**roc()**で作成したオブジェクト。
``identity``：**TRUE**でAUC = 0.5となる線を描画。
``print.thres``：**best**でROC上の最適なカットオフ値を表示。
``print.thres.best.method``：カットオフ値の決定方法。**"youden"**でYoden indexを用いて、**"closest.topleft"**で左上隅との距離が最小となる点を求める。
``legacy.axes``：**TRUE**でx軸の目盛りを1-特異度に合わせて左が0、右が1となるように表示。


```R
#ROC曲線の描画（カットオフ値：Youden indexを用いた方法）
plot(ROC,
     identity = TRUE,
     print.thres = "best",
     print.thres.best.method="youden",
     legacy.axes = TRUE
     )　

```
⬇実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/6aba6901-46b9-d0b7-2b90-9f04ec68f0ed.png)

```R
#ROC曲線の描画（カットオフ値：左上隅からの距離を用いた方法）
plot(ROC,
     identity = TRUE,
     print.thres = "best",
     print.thres.best.method="closest.topleft",
     legacy.axes = TRUE
     )
```
⬇実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/a90e2368-90a6-8a65-e355-a06f29657eb9.png)

Youden indexでカットオフ値を求めた場合、11.0801で特異度、感度はそれぞれ0.514、0.707となりました。

左上隅からの距離を用いた方法でカットオフ値を求めた場合、12.730で特異度、感度はそれぞれ0.625、0.585となりました。

#最終スクリプト

```R
#ROC曲線

#Rをきれいにする
rm(list = ls())

#ライブラリの読み込み
library(pROC)

#データの読み込み
df <- aSAH

#データの確認
summary(df)

#ROC曲線で評価
ROC <- roc(outcome ~ ndka, data = df, ci = TRUE)
ROC

#ROC曲線の描画
plot(ROC)

#ROC曲線の描画（カットオフ値：Youden indexを用いた方法）
plot(ROC,　# roc()で作成したオブジェクト
     identity = TRUE,　# TRUE：AUROC=0.5となるラインを描画
     print.thres = "best",　# "best"：ROC上の最適なカットポイントを表示
     print.thres.best.method="youden",　# "youden"：Youden indexを用いてカットオフ値を決定
     legacy.axes = TRUE　# TRUE：x軸の目盛りを1-特異度に合わせて左が0、右が1となるように表示
     )　

#ROC曲線の描画（カットオフ値：左上隅からの距離を用いた方法）
plot(ROC,　# roc()で作成したオブジェクト
     identity = TRUE,　# TRUE：AUROC=0.5となるラインを描画
     print.thres = "best",　# "best"：ROC上の最適なカットポイントを表示
     print.thres.best.method="closest.topleft",　# "closest.topleft"：左上隅からの距離を利用してカットオフ値を決定
     legacy.axes = TRUE　# TRUE：x軸の目盛りを1-特異度に合わせて左が0、右が1となるように表示
     )

```