目的

生データまたは度数分布表データに基づいて，ROC 曲線を描く。また，ROC 曲線下面積を計算する。

使用法・引数

疾病群，健康群の生データについて
ROC0(disease, normal, lowest=NULL, width=NULL)
　　disease  疾病群の測定値ベクトル
　　normal   健康者群の測定値ベクトル
　　lowest   最も小さい値のキリのよい数値
　　width    階級幅（計算精度）のキリのよい数値
　　         lowest か width がデフォルト値なら，データから適当に決める

度数分布表の形でまとめられているデータについて
ROC(x, disease, normal)
　　x        分割表の下限値のベクトル（例題参照）
　　disease  疾病群の度数分布ベクトル
　　normal   健康者群の度数分布ベクトル

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/ROC.R", encoding="euc-jp")

# 生データまたは度数分布表データに基づいて，ROC 曲線を描く。また，ROC 曲線下面積を計算する
# 生データがあるとき（後述する ROC 関数も必要なので注意）
ROC0 <- function(    disease,                                # 疾病群の測定値ベクトル
                        normal,                                 # 健康者群の測定値ベクトル
                        lowest=NULL,                            # データの最小値より小さいキリのよい数値
                        width=NULL)                             # 度数分布を作成する階級幅のキリのよい数値
{
        my.hist <- function(x, brks)                         # R の hist 関数は，級限界の扱いがイヤラシイので自前の関数を書く
        {
                k <- length(brks)
                freq <- numeric(k)
                for (i in 1:(k-1)) {
                        freq[i] <- sum(brks[i] <= x & x < brks[i+1])
                }
                freq[k] <- sum(x >= brks[k])
                freq
        }

        x <- c(disease, normal)                                      # データをプールする
        min.x <- min(x)                                              # 最小値
        max.x <- max(x)                                              # 最小値
        cat("最小値 x = ", min.x, "\n")
        cat("最大値 x = ", max.x, "\n\n")
        if (is.null(lowest) || is.null(width)) {
                temp <- pretty(c(disease, normal), n=min(length(disease)+length(normal), 50))
                lowest <- temp[1]
                width <- diff(temp)[1]
                cat("最小値より小さいキリのよい数値 = ", lowest, "\n")
                cat("度数分布を作成する階級幅の切りのよい数値 = ", width, "\n\n")
        }
        
        brks <- seq(lowest, max.x+width, by=width)
        ROC(brks, my.hist(disease, brks), my.hist(normal, brks))
        
}
# 度数分布表しかないとき（生データから計算されるときにも，下請けとして使う）    
ROC <- function(     x,                                      # 分割表の下限値のベクトル
                        disease,                                # 疾病群の度数分布ベクトル
                        normal)                                 # 健康者群の度数分布ベクトル
{
        k <- length(x)                                               # 度数分布表の長さ
        stopifnot(k == length(disease) && k == length(normal))    # データのチェック
        Sensitivity <- c(rev(cumsum(rev(disease)))/sum(disease), 0)
        False.Positive.Rate <- c(rev(cumsum(rev(normal)))/sum(normal), 0)
        Specificity <- 1-False.Positive.Rate
        plot(False.Positive.Rate, Sensitivity, type="b")
        abline(h=c(0, 1), v=c(0, 1))
        c.index <- sum(sapply(1:k, function(i) (False.Positive.Rate[i]-False.Positive.Rate[i+1])*(Sensitivity[i+1]+Sensitivity[i])/2)) # area under ROC curve
        result <- cbind(x, disease, normal, Sensitivity[-k-1], Specificity[-k-1], False.Positive.Rate[-k-1])
        rownames(result) <- as.character(1:k)
        colnames(result) <- c("Value", "Disease", "Normal", "Sensitivity", "Specificity", "F.P. rate")
        list(result=result, c.index=c.index)
}


使用例

参考にしたページにある例。

FRA <- c(100, 220, 230, 240, 250, 260, 270, 280, 290, 300, 320, 340, 360, 400)
FRA.disease <- c(3, 2, 1, 4, 7, 4, 16, 5, 3, 9, 10, 5, 10, 21)
FRA.normal <- c(25, 7, 19, 17, 7, 8, 7, 6, 2, 2, 0, 0, 0, 0)

ROC(FRA, FRA.disease, FRA.normal)

生データがあるときの例。

disease.x <- as.integer(rnorm(200, mean=110, sd=10))	# 平均値=110，標準偏差=10 のテストデータ（200ケース）
normal.x <- as.integer(rnorm(300, mean=100, sd=10))	# 平均値=100，標準偏差=10 のテストデータ（300ケース）

ROC0(disease.x, normal.x)

出力結果例

> FRA <- c(100, 220, 230, 240, 250, 260, 270, 280, 290, 300, 320, 340, 360, 400)
> FRA.disease <- c(3, 2, 1, 4, 7, 4, 16, 5, 3, 9, 10, 5, 10, 21)
> FRA.normal <- c(25, 7, 19, 17, 7, 8, 7, 6, 2, 2, 0, 0, 0, 0)

> ROC(FRA, FRA.disease, FRA.normal)
$result
   Value Disease Normal Sensitivity Specificity F.P. rate	# Value は階級の下限値（区間に含まれる）
1    100       3     25        1.00        0.00      1.00	# つまり，最初の階級は，100 以上，220 未満と考える
2    220       2      7        0.97        0.25      0.75	# F.P. rate は偽陽性率
3    230       1     19        0.95        0.32      0.68
4    240       4     17        0.94        0.51      0.49
5    250       7      7        0.90        0.68      0.32
6    260       4      8        0.83        0.75      0.25
7    270      16      7        0.79        0.83      0.17
8    280       5      6        0.63        0.90      0.10
9    290       3      2        0.58        0.96      0.04
10   300       9      2        0.55        0.98      0.02
11   320      10      0        0.46        1.00      0.00
12   340       5      0        0.36        1.00      0.00
13   360      10      0        0.31        1.00      0.00
14   400      21      0        0.21        1.00      0.00

$c.index	# ROC 曲線下面積
[1] 0.88215

以下のような図が描かれる



> disease.x <- as.integer(rnorm(200, mean=110, sd=10))
> normal.x <- as.integer(rnorm(300, mean=100, sd=10))

> ROC0(disease.x, normal.x)
最小値 x =  68 
最大値 x =  132 

最小値より小さいキリのよい数値 =  68 
度数分布を作成する階級幅の切りよい数値 =  1 

$result
   Value Disease Normal Sensitivity Specificity   F.P. rate
1     68       0      1       1.000 0.000000000 1.000000000
2     69       0      0       1.000 0.003333333 0.996666667
3     70       0      0       1.000 0.003333333 0.996666667
4     71       0      0       1.000 0.003333333 0.996666667
5     72       0      0       1.000 0.003333333 0.996666667
6     73       0      1       1.000 0.003333333 0.996666667
7     74       0      0       1.000 0.006666667 0.993333333
8     75       0      2       1.000 0.006666667 0.993333333
9     76       0      1       1.000 0.013333333 0.986666667
　　途中省略
59   126       2      0       0.035 1.000000000 0.000000000
60   127       1      0       0.025 1.000000000 0.000000000
61   128       1      0       0.020 1.000000000 0.000000000
62   129       1      0       0.015 1.000000000 0.000000000
63   130       0      0       0.010 1.000000000 0.000000000
64   131       0      0       0.010 1.000000000 0.000000000
65   132       2      0       0.010 1.000000000 0.000000000
66   133       0      0       0.000 1.000000000 0.000000000

$c.index
[1] 0.7734917

 参考にしたページ