### test-auc.R --- ##---------------------------------------------------------------------- ## Author: Brice Ozenne ## Created: dec 2 2019 (16:55) ## Version: ## Last-Updated: Jan 4 2022 (10:25) ## By: Brice Ozenne ## Update #: 60 ##---------------------------------------------------------------------- ## ### Commentary: ## ### Change Log: ##---------------------------------------------------------------------- ## ### Code: if(FALSE){ library(testthat) library(BuyseTest) library(data.table) library(pROC) library(cvAUC) } context("Check auc calculation vs. cvAUC") ## * Compare AUC and CI n <- 200 set.seed(10) X <- rnorm(n) dt <- data.table(Y = as.factor(rbinom(n, size = 1, prob = 1/(1+exp(1/2-X)))), X = X, fold = unlist(lapply(1:10,function(iL){rep(iL,n/10)}))) ## boxplot(X~Y, data = dt) ## ** no CV test_that("AUC - BuyseTest vs pROC",{ ## example from BuyseTest test <- BuyseTest::auc(labels = dt$Y, predictions = dt$X, direction = ">", transformation = FALSE, pooling = "mean") test2 <- BuyseTest::auc(labels = as.character(dt$Y), predictions = dt$X, direction = ">", transformation = FALSE, pooling = "mean") ## GS2 <- cvAUC::cvAUC(predictions = dt$X, labels = dt$Y) GS <- pROC::roc(dt$Y, -dt$X, ci = TRUE, direction = ">", quiet = TRUE) GS2 <- cvAUC::ci.cvAUC(predictions = dt$X, labels = dt$Y) expect_equal(GS$auc[1], as.double(test[test$fold == "global","estimate"]), tol = 1e-6) expect_equal(GS$auc[1], as.double(test2[test$fold == "global","estimate"]), tol = 1e-6) expect_equal(GS$ci[1], as.double(test[test$fold == "global","lower"]), tol = 1e-3) expect_equal(GS$ci[3], as.double(test[test$fold == "global","upper"]), tol = 1e-3) expect_equal(GS2$cvAUC[1], as.double(test[test$fold == "global","estimate"]), tol = 1e-6) expect_equal(GS2$se, as.double(test[test$fold == "global","se"]), tol = 1e-3) expect_equal(GS2$ci[1], as.double(test[test$fold == "global","lower"]), tol = 1e-3) expect_equal(GS2$ci[2], as.double(test[test$fold == "global","upper"]), tol = 1e-3) ## butils::object2script(test, digit = 6) expect_equal(test$estimate, c(0.7054435), tol = 1e-6) expect_equal(test$se, c(0.03639812), tol = 1e-6) ## example from pROC data(aSAH, package = "pROC") roc.s100b <- pROC::roc(aSAH$outcome, aSAH$s100b, quiet = TRUE) e.BTauc <- BuyseTest::auc(labels = aSAH$outcome, prediction = aSAH$s100b, transformation = TRUE, pooling = "mean") expect_equal(e.BTauc[1,"estimate"],as.numeric(pROC::auc(roc.s100b)), tol = 1e-6) }) ## ** with CV test_that("AUC after CV - BuyseTest vs cvAUC",{ ## cvAUC incorrectly estimate standard error on unequally sized folds if(FALSE){ dtS <- rbind(cbind(dt[,.(Y,X)], fold3 = "1", fold2 = "1", id = 1:NROW(dt)), cbind(dt[,.(Y,X)], fold3 = "2", fold2 = "2", id = 1:NROW(dt)), cbind(dt[,.(Y,X)], fold3 = "3", fold2 = "2", id = 1:NROW(dt))) testS3 <- cvAUC::ci.cvAUC(predictions = dtS$X, labels = dtS$Y, folds = dtS$fold3) testS2 <- cvAUC::ci.cvAUC(predictions = dtS$X, labels = dtS$Y, folds = dtS$fold2) GS <- cvAUC::ci.cvAUC(predictions = dt$X, labels = dt$Y) testBT2 <- auc(labels = dtS$Y, prediction = dtS$X, fold = dtS$fold2, observation = 1:NROW(dtS), pooling = "mean") testBT3 <- auc(labels = dtS$Y, prediction = dtS$X, fold = dtS$fold3, observation = 1:NROW(dtS), pooling = "mean") ## testS3$cvAUC - GS$cvAUC ## testS2$cvAUC - GS$cvAUC ## testS3$se - GS$se/sqrt(3) testS2$se - testS3$se ## this does not seems right testS2$se - testBT2$se[3] weight <- c(1/2,1/2) sum(weight * testBT2$estimate[1:2]) - testBT2$estimate[3] sqrt(sum(weight^2 * testBT2$se[1:2]^2)) - testBT2$se[3] } ## GS0 <- cvAUC::cvAUC(predictions = dt$X, labels = dt$Y, folds = dt$fold) ## gives wrong results as it ignores fold argument GS1 <- cvAUC::ci.cvAUC(predictions = dt$X, labels = dt$Y, folds = dt$fold) test <- BuyseTest::auc(labels = dt$Y, prediction = dt$X, fold = dt$fold, observation = 1:NROW(dt), pooling = "mean") expect_equal(test[test$fold=="global", "estimate"], GS1$cvAUC, tol = 1e-6) expect_equal(test[test$fold=="global", "estimate"], 0.703265, tol = 1e-6) ## expect_equal(test[test$fold=="global", "se"], GS1$se, tol = 1e-6) ## differs (cf issue in cvAUC above) expect_equal(test[test$fold=="global", "se"], 0.0357925, tol = 1e-6) ## another example ## remove warning for uncertainty estimation (P-value/confidence intervals will not be valid with only one observation.) test.auc <- suppressWarnings(BuyseTest::auc(label = c(0,0,0,0,1,1,1,1), prediction = c(0.3,0.4,0.5,0.6,0.55,0.5,0.6,0.4,0.7,0.8), observation = c(1,2,3,4,7,1,2,6,7,8), fold = c(1,1,1,1,1,2,2,2,2,2), direction = ">")) test.auc }) test_that("AUC after CV - ties",{ y <- c(1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1) fit <- c(0.53947368, 0.53947368, 0.53947368, 0.53947368, 0.51315789, 0.51315789, 0.51315789, 0.51315789, 0.52631579, 0.52631579, 0.52631579, 0.52631579, 0.52631579, 0.52631579, 0.52631579, 0.52631579, 0.51315789, 0.51315789, 0.51315789, 0.51315789, 0.53947368, 0.53947368, 0.53947368, 0.53947368, 0.51315789, 0.51315789, 0.51315789, 0.51315789, 0.52631579, 0.52631579, 0.52631579, 0.52631579, 0.53947368, 0.53947368, 0.53947368, 0.53947368, 0.51315789, 0.51315789, 0.51315789, 0.51315789) fold <- c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10) observation <- c(62, 77, 50, 52, 64, 58, 74, 2, 26, 58, 2, 52, 12, 28, 17, 52, 68, 69, 14, 65, 22, 58, 26, 68, 9, 13, 32, 8, 21, 30, 20, 55, 36, 30, 22, 9, 31, 18, 65, 11) test <- suppressWarnings(BuyseTest::auc(labels = y, predictions = fit, fold = fold, observation = observation, add.halfNeutral = TRUE)) expect_true(all(abs(test$estimate-0.5)<1e-6)) }) ###################################################################### ### test-auc.R ends here