R pROC coords() give suboptimal tresholds for my ROC curve
I have a bunch of predictors that could be useful for a diagnosis. Since I have more than 30 predictors, I selected only the predictors that are more useful in this diagnosis, using the AUC measure (area under the ROC curve).
library(pROC)
library(dplyr)
library(purrr)
#Sorted my variables by AUC in table.AUC (one col with AUCs, one col with variable names)
# Then extracted a list of all variables with an AUC>=0.7
list.predictors=table.auc %>% filter(AUC>=0.7) %>% dplyr::select(Variable) %>% unlist()
# Then did a purrr::map() to calculate the bests thresholds, senibility, specificity and positive and negative predictive values
map(list.predictors, .f=~{
threshold=coords(roc(data$Illness,data[[.x]], levels=c("Healthy","Sick")), "best", ret = "threshold") %>% unlist()
c2=data
c2$NM=ifelse(data[[.x]]>=threshold,"1","0")
TP=c2 %>% filter(NM!="0" & Illness=="Sick") %>% nrow()
TN= c2 %>% filter(NM=="0" & Illness=="Healthy") %>% nrow()
FP=c2 %>% filter(NM!="0" & Illness=="Healthy") %>% nrow()
FN= c2 %>% filter(NM=="0" & Illness=="Sick") %>% nrow()
tribble(~"name",~"threshold",~"Sen",~"Spe",~"PPV",~"PNV",
.x,threshold, TP/(TP+FN), TN/(TN+FP), TP/(TP+FP), TN/(TN+FN)
) %>%
mutate(across((Sen:PNV),~round(.x,2)),
threshold=round(threshold,3)) %>% return()
}) %>% bind_rows()
However, I have some weird sensibility and specificity for one predictor : Sen = 18%, Spe = 41%. Looking at the ROC curve, I think my function has failed to find the best threshold (and found one that doesn't even exist ??) .
Looking at the coords() documentation, it seems that the threshold is calculated by maxing (Sen+Spe). However, here, it is clearely possible to find a threshold that have a higher Sen+Spe value.
Is the error coming from my code ? (doing coords(roc(data$Illness,data$the_predictor, levels=c("Healthy","Sick")), "best", ret = "threshold") %>% unlist() just gives the same threshold. Is it coming from the pROC function ?
Help me, I am so lost here.
Edit : code to copy/paste
data=structure(list(Predictor = c(-30.0265168079189, 1.15099363860717,
-11.9804372907826, -5.77111036023964, -2.12543529096544, 3.92087977078597,
5.34247848261018, 4.22335287554391, -23.6137844687066, -19.9373894314097,
-7.66163269551192, 9.13278170787256, -1.36786740246491, -2.72758047450363,
-12.1955983263269, 4.12642888543807, 0.163977606901045, 1.46547604262282,
-1.8412138361035, -6.14477150770652, -3.17570617995096, -30.064951386853,
10.9975933111011, 4.69525214654113, 7.17071677098559, -42.7958812479347,
-22.0530194817619, 5.22959301386168, -5.37225791749525, -8.41992810521969,
0.214693153960777, -10.2699595330821, 7.5188942409622, -8.35948733431006,
-28.3236974772852, -0.334108822977305, 10.3570377885479, -6.21627576016198,
-1.34577632302221, -0.704958574099198, -8.21298595801565, -0.673991917300366,
10.1477239063002, -13.3154432636137, -24.7929589691277, -13.7662638672144,
-11.8388358393301, 1.2811315099285, -16.6351572937165, -11.2898600124526,
-22.8066837752388, 3.56811422380944, -23.3577893771017, -14.1165843206214,
-3.60022726353809, -3.11503373374385, -12.2597492695888, 7.31226134848711,
-11.5833482372835, -16.7535225022961, -8.56422928031759, -47.9699360179454,
-19.7539978196463, -26.3666055391003, -17.2614678085344, 13.2702046763866,
-13.2264594680634, 15.7449331695833, -5.50032293388125, 2.81599579215209,
1.92271901789842, -1.99266757116988, -27.2703581251132, 0.476318001224534,
-6.65722119385071, 2.35294755630814, 13.2977149882188, -15.6037631636617,
1.38135215321582, -9.71125135113497, -4.11971899794967, 5.21192589771364,
5.63060382748609, 3.93730115204109, 12.9318052438389, 6.16873004521188,
0.987650782375663, -11.352515305607, 13.4704068306721, -8.43764175284198,
18.5651453322091, 2.45966987595847, -42.9176544778171, -3.89773417236148,
-28.3070154341685, -13.6189396001791, -8.10091523989799, 2.45341660282897,
3.57427638074422, 16.1254084424853, -13.234491473745, -37.5823778107753,
3.89815090155948, 11.4130217096843, -26.9867338207, -21.9465772218131,
7.34257329115229, -11.3782968157392, -34.7142009728594, -49.8201162965836,
9.29909876231373, -10.0850944178935, 11.6016303844531, 5.23411199811675,
11.7946167683302, 3.28158141564784, 2.03385869625709, -1.89888171264814,
-23.3807367058501, -4.08156864378285, 9.00363191604415, 12.8526886439116,
-13.0772932176236, -8.00976153286771, 8.65199795385856, 15.0607427725054,
14.1375633988393, -8.66285501996315, 10.6594655613105, -17.5706905139599,
6.12042615000835, -3.31329113354379, -3.99151135257908, 8.04726385642244,
-0.56725238509163, -9.36165609862254, -4.62609609004956, -54.3425894866824,
-31.5916887347727, -22.8418309700193, -5.8349028529916, -24.2643112312585,
10.0928208337933, -10.663258246763, 0.154714176468552, 1.94161020180308,
1.23173520914207, -11.5359195732191, 10.4662441932512, -1.15329871483222,
2.48304873375732, 11.1247749589697, 13.600816862598, 5.4048677111595,
7.75142783090626, 12.407508070537, 14.2936521098022, 8.57492908453424,
7.78937753060301, 15.5758817729099, 12.9689634444848, 19.7716027797074,
17.3722434464835, 25.8398368056633, 22.9286256713902, 12.4332591080305,
20.1646775468313, 13.6359691315215, 22.4075156394212, 2.52052660698105,
-16.1294988287829, 24.2580803169934, 6.53302505803349, 28.0007397454365,
12.5997533567436, 19.8231609146357, 20.4758235394322, -3.14132866132124,
15.3652243575294, 5.51333807966991, 12.6701530511479, 7.76875005452657,
9.50414582122118, 7.18440942260985, 17.0577133339618, 1.75973853154878,
19.3991204516457, 17.2086372124398, 11.4946377474413, 18.1208643720225,
13.4371964205874, -1.41932411310972, -0.217642491330989, 0.559828427316018,
14.9513204550895, 8.40789113945641, 7.32690822343274, -10.6306122848244,
8.78982973978855, 19.9236377854904, -16.146128022031, -4.91319750465569,
1.94096133981224, 19.1978399847426, 5.68816798940584, 16.0286843888299,
10.8532049397931, 14.9332042247407, 7.89131242660244, 4.84817793483408,
-23.5709741092812, -5.0750006627721, 9.82990315423877, 13.445883799078,
-16.6078061062201, 8.09109746282265, 4.03027064373712, 7.04329554573276,
-21.7969169325793, 5.17566711827446, 17.7920853961723, 19.5637542942757,
-2.64115652660286, 15.8666090101656, 19.8301694557142, 6.12840423007928,
10.6956110325132, 20.6264711334891, 13.7039419874073, -5.48902063664391,
-13.5307133212758, -3.57547635643989, 8.28126733549754, -28.1686872403507,
21.1407803859056, -1.65377737760872, 0.831092132680901, 16.9496509311116,
-9.18907250852743, 14.8757097896787, 12.2295286749314), Illness = structure(c(1L,
2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L,
2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), levels = c("Healthy", "Sick"), class = "factor")), row.names = c(NA,
-241L), class = c("data.table", "data.frame"))
list.predictors=list("Predictor")
summary(data)
map(list.predictors, .f=~{
threshold=coords(roc(data$Illness,data[[.x]], levels=c("Healthy","Sick")), "best", ret = "threshold") %>% unlist()
c2=data
c2$NM=ifelse(data[[.x]]>=threshold,"1","0")
TP=c2 %>% filter(NM!="0" & Illness=="Sick") %>% nrow()
TN= c2 %>% filter(NM=="0" & Illness=="Healthy") %>% nrow()
FP=c2 %>% filter(NM!="0" & Illness=="Healthy") %>% nrow()
FN= c2 %>% filter(NM=="0" & Illness=="Sick") %>% nrow()
tribble(~"name",~"threshold",~"Sen",~"Spe",~"PPV",~"PNV",
.x,threshold, TP/(TP+FN), TN/(TN+FP), TP/(TP+FP), TN/(TN+FN)
) %>%
mutate(across((Sen:PNV),~round(.x,2)),
threshold=round(threshold,3)) %>% return()
}) %>% bind_rows()
Solution
I believe pROC gives the correct optimal threshold. The problem lies in your code which calculates the parameters incorrectly.
Notice the following output line when running the code:
Setting direction: controls > cases
This indicates that Predictor takes lower values in Sick cases. As a consequence, you are thresholding the data incorrectly. If you don't know the direction your predictor goes a priori, you should replace:
c2$NM=ifelse(data[[.x]]>=threshold,"1","0")
with:
c2$NM=ifelse(data[[.x]]<threshold,"1","0")
which produces the correct output:
Setting direction: controls > cases
# A tibble: 1 × 6
name threshold Sen Spe PPV PNV
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Predictor 1.33 0.82 0.59 0.31 0.94
In addition you should be able to get the results you want much more easily and without the need to compute the performance values yourself with a pipeline like this:
map(list.predictors, .f = ~{
roc(data$Illness, data[[.x]], levels = c("Healthy","Sick")) %>%
coords("best", ret = c("threshold", "sensitivity", "specificity", "ppv", "npv")) %>%
as_tibble() %>%
mutate(name = .x, .before = "threshold",
across((sensitivity:npv), ~round(.x, 2)),
threshold = round(threshold, 3))
}) %>%
bind_rows()
This pipeline produces the same output as yours, using the performance values computed by the coords function.
- grid.arrange's relative height argument didn't works well in RMarkdown
- How apply a function to every element of a matrix?
- Polar plot in R using ggplot2
- Find total length of time in multiple values of a dataset that overlaps with a given period of time
- Faster than R's `simplify2array`?
- What is the most simple way to label a data.table frequency chart in R?
- How to show dates in a secondary x-axis while the primary x-axis is numerical?
- online compiler for r-markdown document
- Unecessarily large column widths in HTML tables published by R Quarto when the table header contains special characters
- Bootstrapping a vector of results, by group in R
- Speeding Up Monte Carlo Simulations for Weather-Dependent Claims
- Is there an R class which allows to list 2 elements but shows only the first of them by default when called?
- Convert day of week number to weekday name in R
- How do you use "<<-" (scoping assignment) in R?
- Using reactivePoll to monitor and process file changes
- How to find out which graph is most similar to another graph?
- ggplot2: rectangular background behind axis text in circular bar plot
- How to simulate button click using rvest
- CMIP6, automating by month
- Covert dummy variables to single categorical in R?
- Animated reducing bar chart with gganimate
- Create an object class matrix/array from data frame
- How to add a vector of names to `sort_by()` in base R
- Among many subfolders, find CSVs starting with a given string and within ZIPs or not, and merge them while adding their names in a new column with R
- How to add a greek letter to an axis label to a 3D plot in plotly?
- bigrquery + dplyr error when using functions and sourcing from separate R script
- Shiny App Error: /v1/applications/ 400 - Validation Error Execution halted
- Rupee Unicode replacement not working in R
- Mutate multiple variable to create multiple new variables
- How to use the Stata equivalent of a foreach loop in R with mutate()