Cross-validation AUC x Logistic Regression

Dear all,

I would like to ask how I can perform cross-validation AUC? I saw the following R code online and am not sure if this is the correct way to do it. If not, can someone pls share the correct R code for 10-fold cvAUC?

Thank you very much :slight_smile:

(cv_logit <-  train(life_exp_cat ~ total_expenditure + 
             schooling + adult_mortality + infant_deaths + thinness_1_19_years +
            alcohol + status, data=life_exp_new, 
            method="glm", 
            trControl = trainControl(method="cv", number=10)))

Hi @Rlearner123,

I believe the above code is based on the {caret} package. The tidymodels family of packages sort of supersedes that package, so here is how you could do that with the modern and more modular API.

We will use the Ames Housing data, since I don't have access to the data in the original problem. Let's see if we can predict whether a house was built before or after 1960 based on just a couple predictors. We will measure model performance with cross-validated AUC.

library(tidymodels)
library(AmesHousing)

# Prepare the data
data <- 
  make_ames() %>% 
  mutate(
    before_1960 = factor(if_else(Year_Built<1960, 1L, 0L))
  ) %>% 
  select(-Year_Built)

# Get the resampling folds
folds <- vfold_cv(data, v = 10, strata = 'before_1960')

# Define a model
model <- logistic_reg(engine = 'glm')

# Define a workflow
wflow <- 
  workflow() %>% 
  add_formula(before_1960 ~ Sale_Price + Lot_Area + Roof_Style) %>% 
  add_model(model)

# Fit the model to the resamples, compute metrics
results <-
  fit_resamples(
    wflow,
    resamples = folds,
    metrics = metric_set(roc_auc),
    control = control_resamples(event_level = 'second')
  )

# Compute CV-AUC
collect_metrics(results)
#> # A tibble: 1 × 6
#>   .metric .estimator  mean     n std_err .config             
#>   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
#> 1 roc_auc binary     0.830    10 0.00751 Preprocessor1_Model1
1 Like

Hi Matt,
Thank you very much for the help. Can you pls explain the output from the last code? Is the mean of 0.83 the cvAUC? Also, can you pls teach me how to get the 95% CI and p-value of this cvAUC model?

#> # A tibble: 1 × 6
#>   .metric .estimator  mean     n std_err .config             
#>   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
#> 1 roc_auc binary     0.830    10 0.00751 Preprocessor1_Model1

Since we fit a model to each of the 10 training sets and then measure the AUC in the 10 testing sets, we have 10 test AUCs in total. We can then average the 10 AUCs to get the overall cross-validated AUC, which is presented in the mean column. The n column shows how many values were used in computing the average, and this number may change if you use more/less resamples, such as with bootstrapping, LOO-CV, or just a different number of folds in vfold_cv. Hope that helps.

As for getting the 95% CI, this is not straightforward to do in the tidymodels framework. You can definitely do this yourself by using control_resamples(save_pred = TRUE) and then computing the AUCs yourself using the pROC package, or whichever package you prefer. It is a bit longwinded, but something like this should work...

# Fit the model to the resamples, save predictions this time
results <-
  fit_resamples(
    wflow,
    resamples = folds,
    metrics = metric_set(roc_auc),
    control = control_resamples(
      event_level = 'second',
      save_pred = TRUE
    )
  )

library(pROC)

collect_predictions(results) %>% 
  group_nest(id) %>% 
  mutate(
    aucs = map(data, function(x) {
      roc = roc(before_1960 ~ .pred_1, data = x, 
          direction = '<', ci = TRUE)
      as.numeric(roc$ci)
    })
  ) %>% 
  unnest_wider(aucs) %>% 
  rename(auc = ...2, conf.low = ...1, conf.high = ...3) %>% 
  summarise(
    n = n(),
    across(c(auc, conf.low, conf.high), mean)
  )

#> # A tibble: 1 × 4
#>       n   auc conf.low conf.high
#>   <int> <dbl>    <dbl>     <dbl>
#> 1    10 0.829    0.780     0.878
1 Like

Hi Matt
Thank you very much for the detailed response. I have now understood it better!
Would you mind sharing how I can obtain the P-value, sensitivity, and specificity for this cvAUC model?

Sorry for the multiple responses, Matt! But I just tried the second code you share and it is not working for my dataset. Can you pls explain where do "x" and ".pred_1" come from?

Would you mind sharing how I can obtain the P-value, sensitivity, and specificity for this cvAUC model?

I think this question needs to be refined a bit. The cross-validated AUC is a metric we use to evaluate our model(s) performance. Sensitivity and specificity are other metrics we can use. However, AUC is based on measuring sensitivity and specificity at every probability threshold in the range of model-predicted probabilities, where as for sens/spec we need to assign class labels, not probabilities. Second, what p-value are you looking for? For what test/comparison? This needs to be well defined.

Can you pls explain where do "x" and ".pred_1" come from?

.pred_1 is the named of the variable that contains the predicted probabilities for belong to class 1. This is a naming convention used by the tidymodels functions. x is just the name of the argument passed to the anonymous function I used in the map() call.

1 Like

Hi Matt,

I am under the impression that every AUC is associated with sensitivity & specificity. Does this apply to the cvAUC as well (maybe there are 10 sens & 10 specs and we have to average them just like the AUC)? Same as p-value. In the unadjusted AUC, I know the p-value is against the AUC=0.5. Can we do that for the cvAUC as well?

You can obtain the Sensitivity and Specificity alongside the AUC via:

results <-
  fit_resamples(
    wflow,
    resamples = folds,
    metrics = metric_set(roc_auc, sensitivity, specificity),
    control = control_resamples(
      event_level = 'second',
      save_pred = TRUE
    )
  )

# Compute CV-AUC, CV-Sens, CV-Spec
collect_metrics(results)
#> # A tibble: 3 × 6
#>   .metric .estimator  mean     n std_err .config             
#>   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
#> 1 roc_auc binary     0.831    10 0.00713 Preprocessor1_Model1
#> 2 sens    binary     0.550    10 0.0175  Preprocessor1_Model1
#> 3 spec    binary     0.870    10 0.0111  Preprocessor1_Model1

If you want to test the AUC against an AUC of 0.5, this should work:

avg_preds <- 
  collect_predictions(results) %>% 
  group_by(.row) %>% 
  summarise(
    y = max(before_1960==1),
    p = mean(.pred_1)
  )

verification::roc.area(avg_preds$y, avg_preds$p)$p.value
#> Registered S3 method overwritten by 'verification':
#>   method    from
#>   lines.roc pROC
#> [1] 5.521164e-186
1 Like

Thank you Matt

Re: the 95% CI, I somehow got the below error message. Do you know what this means?

Error: Problem with mutate() column aucs.
i aucs = map(...).
x object 'disease' not found

(disease is my outcome, I simply replace the before_1990 with it, but it wont work sadly!)

Hi Matt - I am following up regarding the 95% CI question in my previous post. I am wondering if you know how I can fix this issue? Thank you!