pROC_threshhold.R

rm(list = ls())
## Load pROC
library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## 载入程辑包：'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

## Load data
data(aSAH)
## Print a roc object:
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 ...

rocobj <- roc(aSAH$outcome, aSAH$s100b)

## Setting levels: control = Good, case = Poor

## Setting direction: controls < cases

rocobj

## 
## Call:
## roc.default(response = aSAH$outcome, predictor = aSAH$s100b)
## 
## Data: aSAH$s100b in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
## Area under the curve: 0.7314

names(rocobj)

##  [1] "percent"            "sensitivities"      "specificities"     
##  [4] "thresholds"         "direction"          "cases"             
##  [7] "controls"           "fun.sesp"           "auc"               
## [10] "call"               "original.predictor" "original.response" 
## [13] "predictor"          "response"           "levels"

## Plot deciles of possible thresholds

plot(rocobj, print.thres = quantile(rocobj$thresholds, seq(0,1,0.1)))

coords(rocobj, "b", ret = "t", best.method = "youden")

##   threshold
## 1     0.205

plot(rocobj, print.thres="best", print.thres.best.method="closest.topleft",
     print.thres.best.weights=c(5, 0.1), print.auc = TRUE)

###################################
plot(smooth(rocobj))

##################################
plot.roc(aSAH$outcome, aSAH$s100b,
         main="Confidence interval of a threshold", percent=F,
         ci=TRUE, of="thresholds", # compute AUC (of threshold)
         thresholds="best", # select the (best) threshold
         print.thres="best") # also highlight this threshold on the plot

## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases

########################################
#https://rpubs.com/Wangzf/pROC
#ref https://rpubs.com/kaz_yos/pROC-coords