How should two cross-validated logistic regression models be compared?

Question

I'm using 100 times 10-fold repeated cross-validation to assess the ROC-AUC performance improvement of adding a biomarker to an existing model: Model_A : pred1 + pred2 Model_B :pred1 + pred2 + pred3

I've seen advice before to use the Wilcoxon rank test to compare the AUCs between each fold. Averaging ROC curves over folds in cross-validation

Should I pull the median p-value from this? Is it acceptable to use the diff.resamples function in the Caret package and use the Wilcoxin rank instead of the default t-test? Does it need Bonferonni correction if only looking at AUC? https://www.rdocumentation.org/packages/caret/versions/6.0-86/topics/diff.resamples

Lastly, any thoughts on using DeLong or the likelihood ratio test. instead?

Answer 1

As these are nested logistic regression models there is no doubt that Frank Harrell's comment shows how to proceed: do the standard likelihood ratio test on the 2 models,* based on all of the data, to determine whether adding the third predictor improves performance. That has a well established theoretical basis, is more sensitive for detecting model differences than AUC, and it doesn't inherently require cross validation.

Cross validation or bootstrapping to evaluate model optimism and calibration would certainly help bolster your case that your modeling approach is correct, but the emphasis shouldn't be on AUC. There's no harm in showing how much the AUC changes, but that should be a secondary consideration. The validate function in Harrell's rms package provides several measures of model quality based on bootstrapping or cross validation, including a Dxy rank-correlation value (both original and optimism-corrected) that can be transformed into an AUC value.

---

*I'm a bit worried that you seem to be including so few predictors in your model. Logistic regression can have an omitted-variable bias if a predictor associated with outcome is left out of the model. Unlike linear regression, the omitted predictor doesn't even need to be correlated with the included predictors to get biased estimates. That's not to say you should be overfitting, but there are usually so many clinical variables associated with some condition or outcome that only including 2 or 3 would tend to be risky.

Answer 2

Instead of averaging AUCs per fold you can calculate two ROC curve per iteration for Model_A and Model_B (since every instance is exactly predicted once in k-fold CV). To calculate whether the addition of a biomarker results in a model with significantly different AUC you can use DeLong's test. Here, I wouldn't use the median of the p-values - a simple count will do (e.g: around 5 significant p-values out of 100 times 10-fold CV can be explained by chance and indicate no improvement in model performance).

Different approaches to "combine" your p-values are mentioned in "Statistical Methods for Meta-Analysis" by Larry V. Hedges and Ingram Olkin.

If you are using Python and want to use DeLong's test, this blog post might be helpful (altough still in draft): https://biasedml.com/roc-comparison/

Answer 3

DeLong’s test may not be suitable for nested model comparisons (e.g. https://pubmed.ncbi.nlm.nih.gov/22415937/). You may consider R^2 based test (https://www.sciencedirect.com/science/article/pii/S0002929723000046), which is available in CRAN (https://cran.r-project.org/web/packages/r2redux/index.html).

For your model comparisons, there are two ways.

Model_A : y = pred1 + pred2 + e Model_B : y = pred1 + pred2 + pred3 + e

mod = lm (y ~ pred1 + pred2 + pred3) merged_predictor = cbind(pred1, pred2) %*% mod$coefficients[2:3]

Then, the comparison can be equivalently expressed as Model_A : y = merged_predictor + e Model_B : y = merged_predictor + pred3 + e

Then, you can use r2redux as r2_diff(dat,c(1,2),1,nv) (please see manual in https://cran.r-project.org/web/packages/r2redux/index.html).

Note that whether using lm or glm with logit link wouldn’t be really matter, i.e. the result would be very similar.

Secondly, you can use preadjusted y as Model_B : y* = pred3 + e where y* is adjusted y for pred1 and pred2, i.e. y* = mod2$residuals with mod2 = lm(y ~ pred1 + pred2).

Then, you can use r2redux as r2_var(dat,1,nv) with properly specified inputs (please see manual in https://cran.r-project.org/web/packages/r2redux/index.html).

Additional note. There is one-to-one mapping between R^2 and AUC (https://cran.r-project.org/web/packages/R2ROC/index.html), which can also be used in such model comparisons. Please see https://www.biorxiv.org/content/10.1101/2023.08.01.551571v1.full.