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
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