RでのROC解析:ROCRパッケージを使ったROC曲線とAUCの求め方

研究でROC解析を行う必要があり、Rでどうやったらできるのか調べてみました。
そうしたところ、ROCRというパッケージが公開されており、比較的簡単にROC解析を行い、グラフを作成できることがわかりました。

  • ROCRパッケージのインストール(Ubuntu)
  • 既にRはインストールされているとします。Ubuntuの場合、ROCRパッケージはapt経由で簡単に入手できます。

    $ sudo apt-get install r-cran-rocr
    

    これでインストール完了です。

  • ROCRを使うための準備
  • ROC解析に必要なものは、何らかの指標と、それが属するグループの一覧です。具体例を挙げると、以下のようになります。
    第1列に指標、第2列に属するグループ(0か1)が記載されています。
    これをroc_data.txtとという名前で保存することとします。保存したディレクトリをRのワーキングディレクトリとします。

    0.9706	1
    0.9572	1
    0.4854	1
    0.8003	1
    0.1419	1
    0.4218	1
    0.9157	1
    0.7922	1
    0.9595	1
    0.6557	1
    0.0357	1
    0.8491	1
    0.934	1
    0.6787	1
    0.7577	1
    0.7431	1
    0.3922	1
    0.6555	1
    0.1712	1
    0.706	1
    0.4797	0
    0.4551	0
    0.0374	0
    0.081	0
    0.2984	0
    0.7597	0
    0.1404	0
    0.3853	0
    0.0238	0
    0.5513	0
    0.0551	0
    0.306	0
    0.4991	0
    0.6909	0
    0.7593	0
    0.3472	0
    0.0614	0
    0.0507	0
    0.0575	0
    0.6407	0
    
  • ROCRの起動
  • ROCRはRを立ち上げた後に、library(ROCR)で起動できます。

    $ R
    > library(ROCR)
    
  • データの読み込み
  • 先ほどのroc_data.txtをrocdataという変数に読み込みます。変数名は何でもいいのですが、ここではそうします。
    read.tableという関数で表を読み込めるので、それを使います。

    rocdata <- read.table("roc_data.txt")
    
  • ROCR ステップ1: prediction
  • ROCRのステップ1はpredictionで、値と属するグループを指定します。
    今、rocdataは20行2列の行列になっています。1列目はrocdata[,1]で、2列目はrocdata[,2]であらわすことができますので、以下のように記載します。

    pred <- prediction(rocdata[,1], rocdata[,2])
    

    ここで、変数predはpredictionの頭文字4文字です。もちろん、別の名前でもかまいません。
    私はよくカンマを忘れるので、カンマも忘れないようにしましょう。

    ここで、何が行われているかというと、データを大きい順にソートし、真陽性(TP)、偽陽性(FP)、偽陰性(FN)、真陰性(TN) の数を算出します。

  • ROCR ステップ2: performance
  • ROCRのステップ2はperformanceです。ここでは、感度、すなわち真陽性率 (TP/(TP+FN)で定義)と、1-特異度、すなわち偽陽性率(FP/(FP+TN)で定義)を求めます。以下のようにタイプします。

    perf <- performance(pred, "tpr", "fpr")
    

    ここでtprはtrue positive rateを、fprはfalse positive rateを意味します。

  • ROCR ステップ3: グラフの描画
  • それでは、ROC曲線を描きます。非常に簡単です。

    plot(perf)
    

    そうすると、次のようなグラフが現れるはずです。

    roc-curve

  • グラフの保存
  • グラフをPNG形式で保存するには、次のように行うことで、roc-curve.pngという名前でワーキングディレクトリに保存されます。

    png("roc-curve.png")
    plot(perf)
    dev.off()
    
  • AUCの算出
  • 先ほどのperformanceの際に”auc”と指定するとAUCも計算されます。ただ、1クッション入れる必要があります。具体的な方法はこちらのサイトに記載されていましたが、それを転載します。

    auc.tmp <- performance(pred,"auc")
    auc <- as.numeric(auc.tmp@y.values)
    

    performance(pred,”auc”)の結果をauc.tmpという変数に代入し、
    auc.tmpの中からy.valuesの値を取り出して、その値を変数aucに代入します。

    最後に変数aucを表示させてみます。

    auc
    [1] 0.8
    

    これでAUCが0.8だということがわかります。

  • 正診率の算出
  • ここからさらに一歩踏み込んで、正診率を求めたいと思います。

    正診率、感度、特異度は以下で定義されます。

    正診率=(TP+TN)/総数
    感度(真陽性率)=TP/(TP+FN)
    特異度=TN/(FP+TN)

    今、総数は、rocdataの行数を求めればよいですから、nrow(rocdata)で求められます

    カットオフ値を少しずつずらした時に、TP, FP, FN, TNは変わっていきますので、その一覧を表に出力しましょう。

    表の列は以下のようにしたいと思います

    Cutoff TP FP FN TN Sensitivity Specificity Accuracy

    table <- data.frame(Cutoff=unlist(pred@cutoffs),
      TP=unlist(pred@tp), FP=unlist(pred@fp),
      FN=unlist(pred@fn), TN=unlist(pred@tn),
      Sensitivity=unlist(pred@tp)/(unlist(pred@tp)+unlist(pred@fn)),
      Specificity=unlist(pred@tn)/(unlist(pred@fp)+unlist(pred@tn)),
      Accuracy=((unlist(pred@tp)+unlist(pred@tn))/nrow(rocdata))
      )
    

    これで、tableを表示させると、以下のように表示されます。

    > table
       Cutoff TP FP FN TN Sensitivity Specificity Accuracy
    1     Inf  0  0 20 20        0.00        1.00    0.500
    2  0.9706  1  0 19 20        0.05        1.00    0.525
    3  0.9595  2  0 18 20        0.10        1.00    0.550
    4  0.9572  3  0 17 20        0.15        1.00    0.575
    5  0.9340  4  0 16 20        0.20        1.00    0.600
    6  0.9157  5  0 15 20        0.25        1.00    0.625
    7  0.8491  6  0 14 20        0.30        1.00    0.650
    8  0.8003  7  0 13 20        0.35        1.00    0.675
    9  0.7922  8  0 12 20        0.40        1.00    0.700
    10 0.7597  8  1 12 19        0.40        0.95    0.675
    11 0.7593  8  2 12 18        0.40        0.90    0.650
    12 0.7577  9  2 11 18        0.45        0.90    0.675
    13 0.7431 10  2 10 18        0.50        0.90    0.700
    14 0.7060 11  2  9 18        0.55        0.90    0.725
    15 0.6909 11  3  9 17        0.55        0.85    0.700
    16 0.6787 12  3  8 17        0.60        0.85    0.725
    17 0.6557 13  3  7 17        0.65        0.85    0.750
    18 0.6555 14  3  6 17        0.70        0.85    0.775
    19 0.6407 14  4  6 16        0.70        0.80    0.750
    20 0.5513 14  5  6 15        0.70        0.75    0.725
    21 0.4991 14  6  6 14        0.70        0.70    0.700
    22 0.4854 15  6  5 14        0.75        0.70    0.725
    23 0.4797 15  7  5 13        0.75        0.65    0.700
    24 0.4551 15  8  5 12        0.75        0.60    0.675
    25 0.4218 16  8  4 12        0.80        0.60    0.700
    26 0.3922 17  8  3 12        0.85        0.60    0.725
    27 0.3853 17  9  3 11        0.85        0.55    0.700
    28 0.3472 17 10  3 10        0.85        0.50    0.675
    29 0.3060 17 11  3  9        0.85        0.45    0.650
    30 0.2984 17 12  3  8        0.85        0.40    0.625
    31 0.1712 18 12  2  8        0.90        0.40    0.650
    32 0.1419 19 12  1  8        0.95        0.40    0.675
    33 0.1404 19 13  1  7        0.95        0.35    0.650
    34 0.0810 19 14  1  6        0.95        0.30    0.625
    35 0.0614 19 15  1  5        0.95        0.25    0.600
    36 0.0575 19 16  1  4        0.95        0.20    0.575
    37 0.0551 19 17  1  3        0.95        0.15    0.550
    38 0.0507 19 18  1  2        0.95        0.10    0.525
    39 0.0374 19 19  1  1        0.95        0.05    0.500
    40 0.0357 20 19  0  1        1.00        0.05    0.525
    41 0.0238 20 20  0  0        1.00        0.00    0.500
    

    Accuracyがもっとも高いところを見つけるには、

    max(table$Accuracy)
    

    とします。そうすると、今は、

    > max(table$Accuracy)
    [1] 0.775
    

    となりますので、該当するところをみると、感度70%、特異度85%、正診率77.5%達成できるということがわかりました。

  • おまけ:感度特異度曲線
  • 感度と特異度の曲線も簡単に書けます。predictionまで行った後に、次のようにします。

    perf <- performance(pred, "sens", "spec")
    png("sens-spec-curve.png")
    plot(perf)
    dev.off()
    

    これで、下図のような感度、特異度の曲線がsens-spec-curve.pngという名前で保存されます。

    sens-spec-curve

比較的簡単に求められるので便利です。

Installing R and JGR on Ubuntu Lucid Lynx

The other day I wanted to install R on Ubuntu. Googling led me to this site, but I found some of them are not up-to-date. In addition, I don’t use Sun-Java but open-JDK. So I needed to change things a little bit.

Below is what I did. Note that I already installed openJDK.

  1. Add an R repository to /etc/apt/sources.list
  2. First, you need to decide a CRAN mirror close to where you live. The list can be found here. Since I live in Tsukuba Japan, I added the following line to /etc/apt/sources.list

    deb http://cran.md.tsukuba.ac.jp/bin/linux/ubuntu lucid/

    Below is the format.

    deb http://my.favorite.cran.mirror/bin/linux/ubuntu lucid/

  3. Obtain the public key of Ubuntu repository
  4. The Ubuntu archives on CRAN are signed with the key of “Michael Rutter
    ” with key ID E084DAB9. You can fetch this key then feed it to apt-key with the following lines.

    $ gpg --keyserver keyserver.ubuntu.com --recv-key E084DAB9
    $ gpg -a --export E084DAB9 | sudo apt-key add -

  5. Update the respository
  6. $ sudo apt-get update

  7. Install r-base
  8. sudo apt-get install r-base r-base-dev

  9. Upgrade the repository
  10. By upgrading the repository, r-cran-* will be installed.

    $ sudo apt-get upgrade

    The following should be displayed.

    Reading package lists… Done
    Building dependency tree
    Reading state information… Done
    The following packages will be upgraded:
    r-base-html r-cran-boot r-cran-class r-cran-cluster r-cran-codetools
    r-cran-foreign r-cran-kernsmooth r-cran-lattice r-cran-mass r-cran-matrix
    r-cran-mgcv r-cran-nlme r-cran-nnet r-cran-rpart r-cran-spatial
    r-cran-survival r-doc-html
    17 upgraded, 0 newly installed, 0 to remove and 0 not upgraded.
    Need to get 12.2MB of archives.
    After this operation, 827kB disk space will be freed.
    Do you want to continue [Y/n]?

  11. Set java configuration
  12. $ sudo R CMD javareconf

  13. Start R environment
  14. $ sudo R

  15. Install JGR inside R environment
  16. install.packages('JGR')

    A GUI will popup asking you to select a mirror to download from. Select the same mirror as before, though it doesn’t matter.

    After installation is completed. Add the library and install the ggplot2 library.

    library(JGR)
    install.packages('ggplot2', dep = TRUE)

    This will take some time.

    After that, you can start JGR by typing

    JGR()

    This will bring up JGR.

    But this is not the end. Please look at the shell carefully, which says,

    Starting JGR run script. This can be done from the shell as well, just run
    /usr/local/lib/R/site-library/JGR/scripts/run

    That means that you don’t have to run R first before running JGR.

    Now I will make a launcher of JGR.

    Right-click of GNOME menu –> “Edit Menu”

    That will pop up Main Menu Dialogue.

    In this case, I want to put the JGR launcher under Science, so click “New Item”

    Then you will see JRG in the menu. Note that I got the made the icon of JGR from the JGR site. I also customized the user interface. The screenshot is the from Lin4Neuro.