rrocaucdiagnosticsproc-r-package

Different ROC optimal cutoff obtained from different functions with same methods


I just tried to calculate optimal cutoff from a ROC curve. However, when I tried several functions from different packages, they returned different results. Which one is corrent one if I want to use what is calculated by Youden Statistics.

Thank you very much. Below are the results.

if(!require(install.load)) {install.packages("install.load"); require(install.load)} #load / install+load install.load
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
if (!require("BiocManager", quietly = TRUE)) {install.packages("BiocManager")}
if (!require("limma", quietly = TRUE)) {BiocManager::install("limma", force = TRUE)}

install.load::install_load("openxlsx", "readxl", "shiny", "dplyr", "tidyr", "zoo", "tidyverse", "stringr", "splitstackshape")
install.load::install_load("compiler", "XML")
install.load::install_load("excel.link", "gt")
install.load::install_load("devtools")
install.load::install_load("pROC", "OptimalCutpoints", "verification", "reportROC", "MKmisc")
if(!require(voronoiTreemap)) {devtools::install_github('uRosConf/voronoiTreemap'); library(voronoiTreemap)} #load / install+load install.load
if(!require(multipleROC)) {remotes::install_github("cardiomoon/multipleROC")}; require(multipleROC)


> completedf_01 <- data.frame(
BLDCCFL = c(0,0,1,1,0,0,0,1,0,1,1,0,1,0,1,1,1,0,1,0,0,0,1,0,1,1,0,1,0,0,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,0,1,1,1,0,1,1,0,0,0,1,0,1,1,1,0,1,1,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,0,1,1,1,0,0,0,1,0,1,1,0,1,1,1,0,0,0,1,0,0,1,1,0,1,0,1,0,1,1,1,0,1,1,1,0,0,0,1,1,1,0,1,1,1,1,0,0,1,1,1,0,1,1,1,1,1,0,1,1,1,0,1,1,1,1,0,0,1,1,1,1,0,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,0,1,1,0,1,1,0,0,0,1,1,0,1,1,1,0,1,1,1,1,1,1,1,0,1,1,1,1,1,0,1,0,1,1,1,1,0,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,0,1,1,1,1,0,1,0,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,0,0,1,1,0,0,1,0,1,1,1,0,1,1,0,0,1,1,0,0,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,0,1,1,0,1,1,0,1,1,0,0,1,0,1,1,1,0,0,0,0,1,0,1,1,1,1,1,0,0,1,1,0,1,0,0,1,0,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,0,1,0,0,0,1,0,0,1,1,0,0,0,0,0,0,0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,1,1,1,0,0,0,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,1,1,1,1,1,1,1,0,1,0,1,1,1,1,0,1,0,1,1,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,0,0,1,0,1,1,0,1,0,1,1,1,1,0,0,1,0,1,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,0,1,1,0,0,1,0,0,1,1,1,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,0,0,1,1,0,1,1,0,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,0,1,0,1,1,0,1,0,0,1,0,1,0,0,1,1,1,0,1,0,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,1,0,0,1,0,1,1,1,0,0,1,1,0,0,1,1,1,0,1,0,1,0,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,0,1,0,1,0,1,1,1,1,1,1,0,0,0,1,1,1,0,1,1,1,0,1,1,1,1,0,0,1,0,1,0,0,1,1,0,1,1,0,0,1,1,1,0,0,1,1,0,1,0,1,1,0,0,0,1,1,1,0,1,0,1,0,0,0,1,0,1,1,0,0,0,1,1,1,1,1,0,1,0,1,1,1,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,0,1,1,1,1,1,0,1,1,1,1,1,1,1,0,0,0,1,0,0,1,1,1,0,1,0,0,1,0,1,1,1,1,1,0,1,1,0,1,1,0,0,1,1,1,1,1,1,0,1,0,1,1,1,0,1,1,1,1,0,1,1,1,1,0,1,1,1,1,0,1,1,0,1,0,1),

URIWBC = c(11,312.37,9,1339.82,1.69,16.75,0,14.06,242.62,11.78,74.26,9.42,17.96,72.13,1.76,13.47,30.91,49.24,14.59,
134.06,506.34,9.79,3.37,18.14,0,444.77,53.01,1058.9,1.84,1.8,8.28,0,67.53,10.79,0,11.56,14.49,134.81,11.32,6.52,0,13.47,
444.77,3.11,0,142.33,86.97,511.77,613.05,189.17,23.81,21.58,12.82,5.05,6.73,0,0,145.21,1.5,5.28,0,10.06,0,3.73,0,4.42,116.43,
0,3.55,3.37,5.8,146.88,0,107.07,13.32,6.31,107.95,57.47,1732.96,213.48,3.37,1.53,1.63,26.42,0,102.05,4.29,159.51,58.38,73.9,13.97,
45.27,35.83,178.85,2.22,0,0,6.56,33.77,0,4.19,690.21,45.48,18.36,177.09,20.69,0,1.3,7,0,22.1,20.72,10.36,38.54,16.38,9.37,1.46,8.52,267,
411.84,105.84,64.97,230.47,13.1,0,7.02,1.82,252.7,124.84,526.83,0,0,21.9,290.11,11.51,1.63,166.37,8.63,3.26,305.89,4.89,152.89,5.05,70.15,
140.59,0,1.83,0,78.52,3,29.28,102.63,107.95,3.66,3.35,298.15,58.01,113.39,0,71.55,80.49,150.34,74.34,24,11.04,23.23,256.94,3415.49,246.03,
808.21,10000,57.31,2.63,61,83.38,0,0,61.55,9.81,10.85,5.66,49.07,5.78,20.53,37.17,31.07,0,35.63,0,13.02,1.75,24.19,1.63,310.67,135.93,26.04,
13.98,329.13,205.99,2.22,10000,8.95,25.61,15.71,27.67,28.27,1.29,0,65.56,4.56,73.14,79.38,25.05,1.76,64.21,15.52,81.95,355.18,146.14,0,23.01,
58.45,0,3.7,101.63,5.49,70.65,0,18.52,0,0,75.47,0,0,65.65,6.31,511.42,50.48,0,38.14,13.2,1.68,77.25,0,12.57,0,4.79,20.91,4.04,32.4,3.41,9.15,
0,28.38,1.68,88.08,3.67,45.4,179.3,1296.23,7.92,4.27,67.3,0,24.3,3.45,0,33,0,27.99,0,92.11,6.73,29.14,195.62,120.73,14.06,454.59,25.15,10.06,1299.3,
3.37,26.39,8.74,34.3,333.96,54.88,696.86,28.8,0,0,8.42,156.61,67.16,68.99,24.58,198.42,0,74.93,2.39,2.5,3.55,0,0,8.57,37.74,0,365.41,61.55,51.83,158.47,
68.07,311.03,5.35,120.88,1,49.46,1886.1,75.13,8.72,11.78,42.79,14.68,3.45,9.91,1.68,10.61,12.93,0,16.83,3.67,0,6.85,0,111.08,5.26,0,19.29,53.23,329.09,0,
72.59,18.74,201.99,0,3.45,80.57,41.92,11.75,0,0,17.01,1.7,204.99,0,42.35,898.88,5.37,74.81,127.92,76.31,54.64,16.68,11.16,2.85,61.68,2363.32,39,45.05,
117.87,270.03,128.2,0,13.33,9.65,16.77,106.21,183.91,0,0,46.79,2.69,529.3,78.79,295.35,1.72,0,9.39,9.76,48.2,5.01,1.64,6.97,0,0,0,20.31,7.5,7.3,19.49,
5.12,13.65,57.51,14.3,10.97,47.11,12.59,2.31,0,739.31,3.47,293.62,15.7,25.64,6.01,5.03,1.68,10.94,146.26,0,3.96,122.98,0,93.44,536.63,8.63,3.35,43.88,
0,5.96,0,0,29.91,6.73,8.42,45.9,99.99,10.89,477.67,17.41,38.59,17.57,472.84,0,0,8.42,1.73,1268.07,112.08,0,0,3.44,0,89.6,54.53,6.27,13.97,20.2,21.06,
16.54,5.51,306.3,0,0,624.24,1.67,0,0,38.56,1.68,20.86,4.42,6,6,252.49,0,19.88,0,60.92,0,0,52.85,3.37,91.53,133.83,15.9,71.79,1.68,0,9.47,1.69,8,0,199.24,
1.68,135.66,118.88,8.83,53.17,5.05,99.83,0,69.39,55.96,8.08,3.76,40.89,188.6,8.38,25.34,0,30.94,674.46,0,310.78,6.73,5.05,33.59,33.11,0,1.74,47.96,4.54,
1.68,108.78,0,43.34,78.73,16.71,1.73,1.43,72.52,19.41,0,0,22.51,9.97,8.41,7.48,33.13,0,1189.36,22.52,34.54,76.35,18.94,34.26,68.94,505.66,8.19,0,18.11,
22.38,265.31,90.33,22.66,6.21,125.67,76.22,13.13,15.4,1,15.04,0,14.47,22.09,4.39,68.7,0,45.94,67.09,0,88.65,4.8,21.3,61.61,37.49,211.43,160.61,1.63,0,
6.73,0,188.07,0,1634.34,36.34,3.45,1527.45,18.44,84.84,0,1.7,6.25,1.77,10.17,52.07,58.58,13.47,87.93,0,94.79,0,3.65,364.25,3.7,73.32,1.7,3.26,6.06,0,1725.34,
50.04,193.8,1641.8,0,24.31,171.38,0,3.94,0,0,13.71,106.32,256.27,405.68,12.4,5.51,1.84,15.77,11.91,0,0,97.02,0,11.52,0,92.76,0,1.77,470.11,28.19,95.37,0,0,
16.33,0,14.04,0,17.76,0,105.04,18,0,587.67,469.96,193.48,5.18,122.03,12.09,558.93,90.98,18.7,18.91,0,558.32,244.49,4.14,6.11,0,13.03,103.45,3.26,387.07,
2.99,2.99,22.11,55.95,331.74,4.42,25.77,28.17,517.39,82.12,10.05,5.74,3.62,0,11.62,116.58,20.1,24.94,2.95,4.89,4,46.34,12.37,7.66,0,90.16,38.37,1,1.97,
2.55,0,0,0,620.71,5.85,21.57,150.68,97.79,115.09,940.7,84.85,44.89,0,62.55,74.05,21.27,2.24,6.13,800.23,6.91,62.41,1.79,10.92,30.87,49.34,1,7.87,75.78,
1737.47,16.55,19.4,0,10.7,18.72,2.15,24.12,1.47,0,46.7,2857.1,2.63,366.81,429.46,6.6,24,15.82,0,0,10.87,0,12.42,0,3154,0,2.3,1,33.19,2.67,10.3,7.94,0,
20.99,63.67,4.91,38.44,106.05,0,14.08,4.96,3,19.12,187.05,89.79,95.14,0,1294.4,3.23,1,1607.23,1.96,1,485.35,35.32,35.32,1.52,25.65,1.52,20.94,6.7,6.7,
0,10.1,0,1979.94,344.62,0,1088.05,12.09,1.83,59.91,23.94,0,1.7,1182.92,73.67,29.98,1.7,12.9,3.37,279.14,218.37,0,0,90.38,2.28,72.98,296.75,7.53,1.48,
11.27,55.5,0,18.76,1.47,1,103.13,21.86,66.87,1.99,39.82,0,0,0,5320.51,0,1267.94,65.62,0,29.36,18.52,52.85,107.16,69.01,11.11,0,122.83,565.85,34.02,5.89,
33.57,1.68,2.01,57.08,4.04,3293.64,184.75,317.75,10.17,57.8,16.72,5.05,0,70.78,151.8,1.67,0,14.65,3.61,0,0,25.34,0,0,18.68,95.41,25.87,0,9.86,4.57,38,0,19.03,25,7.96)
)



> completedf_01 <- data_raw[complete.cases(data_raw[, c("BLDCCFL", "URIWBC")]), c("BLDCCFL", "URIWBC")]
> rocobj_01 <- pROC::roc(completedf_01$BLDCCFL, completedf_01$URIWBC)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> rocauc_01 <- pROC::auc(rocobj_01)
> rocaucci_01 <- pROC::ci.auc(rocauc_01, conf.level=0.95, method=c("delong"))
> # reportROC(gold=completedf_01$BLDCCFL,predictor=completedf_01$URIWBC, important="se",plot=FALSE, exact = TRUE)
> rocpval_01 <- verification::roc.area(rocobj_01[["response"]],rocobj_01[["predictor"]])$p.value
> rocpval_01 <- ifelse(rocpval_01 < 0.001, "italic(p) < 0.001", paste0("italic(p) == ", round(rocpval_01, digits = 3)))
> rocoptcut_01 <- cutoff::roc(score = rocobj_01[["predictor"]], class = rocobj_01[["response"]])
> rocoptcut_01
                     type  auc cutoff sensitivity specificity
1 positive classification 0.57   9.91   0.6263158   0.5345912
> # OptimalCutpoints::optimalCutoff(actuals=rocobj_01[["response"]], predictedScores=rocobj_01[["predictor"]], optimiseFor="misclasserror", returnDiagnostics=TRUE)
> rocoptcutoff_01 <- MKmisc::optCutoff(pred = rocobj_01[["predictor"]], truth = rocobj_01[["response"]], namePos = 1)
> rocoptcutoff_01
     Optimal Cut-off Youden's J statistic 
               9.910                0.161 
> pROC::coords(rocobj_01, "best", ret = c("all"), best.method = "youden")
          threshold specificity sensitivity  accuracy  tn  tp  fn  fp       npv       ppv       fdr
threshold     9.885   0.5345912   0.6263158 0.5934685 170 357 213 148 0.4438642 0.7069307 0.2930693
                fpr       tpr       tnr       fnr 1-specificity 1-sensitivity 1-accuracy     1-npv
threshold 0.4654088 0.6263158 0.5345912 0.3736842     0.4654088     0.3736842  0.4065315 0.5561358
              1-ppv precision    recall   youden closest.topleft
threshold 0.2930693 0.7069307 0.6263158 1.160907       0.3562452

Basically, function cutoff::roc and MKmisc::optCutoff got same optimals: 9.91; while coords from pROC packages reported 9.885. Why are these methods generate different results. Which one shall I use?

Hope there is a professional expalain the reason and provide best solution to choose correct results from these outputs


Solution

  • 9.91 and 9.885 are essentially equivalent with your dataset. As the data is not continuous, any cutoff value t that satisfies 9.86 < t <= 9.91 will give you identical performance values.

    The difference stems from the way the cutoffs are chosen: pROC takes the average between 9.86 and 9.91, while the others take the upper value of 9.91. You can check that the values are identical directly with the pROC::coords function:

    > pROC::coords(rocobj_01, c(9.885, 9.91), best.method = "youden")
      threshold specificity sensitivity
    2     9.885   0.5345912   0.6263158
    3     9.910   0.5345912   0.6263158