diff --git a/R/delong.R b/R/delong.R index 20ab494..13573b9 100644 --- a/R/delong.R +++ b/R/delong.R @@ -18,31 +18,15 @@ # along with this program. If not, see . # Delong's test paired, used by roc.test.roc -delong.paired.test <- function(roc1, roc2) { - n <- length(roc1$controls) - m <- length(roc1$cases) - - VR <- delongPlacements(roc1) - VS <- delongPlacements(roc2) - - SX <- matrix(NA, ncol=2, nrow=2) - SX[1,1] <- sum((VR$X - VR$theta) * (VR$X - VR$theta))/(m-1) - SX[1,2] <- sum((VR$X - VR$theta) * (VS$X - VS$theta))/(m-1) - SX[2,1] <- sum((VS$X - VS$theta) * (VR$X - VR$theta))/(m-1) - SX[2,2] <- sum((VS$X - VS$theta) * (VS$X - VS$theta))/(m-1) - - SY <- matrix(NA, ncol=2, nrow=2) - SY[1,1] <- sum((VR$Y - VR$theta) * (VR$Y - VR$theta))/(n-1) - SY[1,2] <- sum((VR$Y - VR$theta) * (VS$Y - VS$theta))/(n-1) - SY[2,1] <- sum((VS$Y - VS$theta) * (VR$Y - VR$theta))/(n-1) - SY[2,2] <- sum((VS$Y - VS$theta) * (VS$Y - VS$theta))/(n-1) - - S <- SX/m + SY/n - L <- c(1,-1) - sig <- sqrt(L%*%S%*%L) - zscore <- (VR$theta-VS$theta)/sig[1] - if (is.nan(zscore) && VR$theta == VR$theta && sig[1] == 0) - zscore <- 0 # special case: no difference between theta's produces a NaN +delong.paired.test <- function(calcs) { + + # Input calcs is a list returned by delong.paired.calculations(). + + zscore <- with(calcs, d/sig) + + if (is.nan(zscore) && calc$d == 0 && calcs$sig == 0) + zscore <- 0 # special case: no difference between theta's produces a NaN + return(zscore) } @@ -111,6 +95,49 @@ ci.auc.delong <- function(roc, conf.level) { return(ci) } +# function to calculate the CI +ci.delong.paired <- function(calcs, conf.level) { + + # Input calcs is a list generated by delong.paired.calculations(). + # CI is calculated using the normally distributed pivot given in + # DeLong's 1988 paper. + + crit_z <- qnorm(1 - ((1 - conf.level)/2)) + out <- list() + out$upper <- with(calcs, d + crit_z * sig) + out$lower <- with(calcs, d - crit_z * sig) + out$level <- conf.level + return(out) +} + +# Runs the placements and main calculations for the paired DeLong's test +# so that they can be easily used by both the test and CI functions. +delong.paired.calculations <- function(roc1, roc2) { + n <- length(roc1$controls) + m <- length(roc1$cases) + + VR <- delongPlacements(roc1) + VS <- delongPlacements(roc2) + + SX <- matrix(NA, ncol=2, nrow=2) + SX[1,1] <- sum((VR$X - VR$theta) * (VR$X - VR$theta))/(m-1) + SX[1,2] <- sum((VR$X - VR$theta) * (VS$X - VS$theta))/(m-1) + SX[2,1] <- sum((VS$X - VS$theta) * (VR$X - VR$theta))/(m-1) + SX[2,2] <- sum((VS$X - VS$theta) * (VS$X - VS$theta))/(m-1) + + SY <- matrix(NA, ncol=2, nrow=2) + SY[1,1] <- sum((VR$Y - VR$theta) * (VR$Y - VR$theta))/(n-1) + SY[1,2] <- sum((VR$Y - VR$theta) * (VS$Y - VS$theta))/(n-1) + SY[2,1] <- sum((VS$Y - VS$theta) * (VR$Y - VR$theta))/(n-1) + SY[2,2] <- sum((VS$Y - VS$theta) * (VS$Y - VS$theta))/(n-1) + + S <- SX/m + SY/n + L <- c(1,-1) + sig <- sqrt(L%*%S%*%L) + d <- VR$theta - VS$theta + return(list("d" = d, "sig" = sig[[1]])) +} + # Calls delongPlacementsCpp safely # Ensures that the theta value calculated is correct delongPlacements <- function(roc) { diff --git a/R/roc.test.R b/R/roc.test.R index 30550bf..9df2d37 100644 --- a/R/roc.test.R +++ b/R/roc.test.R @@ -116,6 +116,7 @@ roc.test.roc <- function(roc1, roc2, ties.method="first", progress=getOption("pROCProgress")$name, parallel=FALSE, + conf.level=0.95, ...) { alternative <- match.arg(alternative) data.names <- paste(deparse(substitute(roc1)), "and", deparse(substitute(roc2))) @@ -249,6 +250,15 @@ roc.test.roc <- function(roc1, roc2, } if (roc1$direction != roc2$direction) warning("DeLong's test should not be applied to ROC curves with a different direction.") + + # Check if conf.level is specified correctly. This is currently + # only used for the delong paired method, which is why it lives + # here for now. + if (!is.numeric(conf.level)) { + stop("conf.level must be numeric between 0 and 1.") + } else if (0 > conf.level | 1 < conf.level) { + stop("conf.level must be between 0 and 1.") + } } else if (method == "venkatraman") { if (has.partial.auc(roc1)) @@ -304,10 +314,14 @@ roc.test.roc <- function(roc1, roc2, if (method == "delong") { if (paired) { - stat <- delong.paired.test(roc1, roc2) + delong.calcs <- delong.paired.calculations(roc1, roc2) + stat <- delong.paired.test(delong.calcs) + stat.ci <- ci.delong.paired(delong.calcs, conf.level) names(stat) <- "Z" htest$statistic <- stat htest$method <- "DeLong's test for two correlated ROC curves" + htest$conf.int <- c(stat.ci$lower, stat.ci$upper) + attr(htest$conf.int, "conf.level") <- stat.ci$level if (alternative == "two.sided") pval <- 2*pnorm(-abs(stat)) diff --git a/man/roc.test.Rd b/man/roc.test.Rd index fbcaee4..80443e4 100644 --- a/man/roc.test.Rd +++ b/man/roc.test.Rd @@ -25,7 +25,7 @@ specificity = NULL, alternative = c("two.sided", "less", "greater"), paired=NULL, reuse.auc=TRUE, boot.n=2000, boot.stratified=TRUE, ties.method="first", progress=getOption("pROCProgress")$name, -parallel=FALSE, ...) +parallel=FALSE, conf.level=0.95, ...) \S3method{roc.test}{auc}(roc1, roc2, ...) \S3method{roc.test}{smooth.roc}(roc1, roc2, ...) \S3method{roc.test}{formula}(formula, data, ...) @@ -97,6 +97,8 @@ na.rm=TRUE, method=NULL, ...) \item{parallel}{if TRUE, the bootstrap is processed in parallel, using parallel backend provided by plyr (foreach). } + \item{conf.level}{a numeric scalar between 0 and 1 (non-inclusive) which + species the confidence level to use for any calculated CI's.} \item{\dots}{further arguments passed to or from other methods, especially arguments for \code{\link{roc}} and \code{roc.test.roc} when calling \code{roc.test.default} or \code{roc.test.formula}. @@ -300,6 +302,7 @@ na.rm=TRUE, method=NULL, ...) \item{statistic}{the value of the Z (\code{method="delong"}) or D (\code{method="bootstrap"}) statistics. } + \item{conf.int}{the confidence interval of the test (currently only returned for the paired DeLong's test). Has an attribute \code{conf.level} specifiying the level of the test.} \item{alternative}{the alternative hypothesis.} \item{method}{the character string \dQuote{DeLong's test for two correlated ROC curves} (if \code{method="delong"}) or diff --git a/tests/testthat/test-roc.test.R b/tests/testthat/test-roc.test.R index cf89377..83d5b82 100644 --- a/tests/testthat/test-roc.test.R +++ b/tests/testthat/test-roc.test.R @@ -16,21 +16,30 @@ test_that("roc.test works", { test_that("roc.test statistic and p are as expected with defaults", { expect_equal(t1$statistic, c(Z=2.20898359144091)) expect_equal(t1$p.value, 0.0271757822291882) + expect_equal(t1$conf.int[[1]], 0.0104061769564846) + expect_equal(t1$conf.int[[2]], 0.174214419249478) expect_match(t1$method, "DeLong") expect_match(t1$method, "correlated") expect_identical(t1$alternative, "two.sided") + expect_identical(attr(t1$conf.int, "conf.level"), 0.95) expect_equal(t2$statistic, c(Z=2.79777591868904)) expect_equal(t2$p.value, 0.00514557970691098) + expect_equal(t2$conf.int[[1]], 0.0634011709339876) + expect_equal(t2$conf.int[[2]], 0.3600405634833566) expect_match(t2$method, "DeLong") expect_match(t2$method, "correlated") expect_identical(t2$alternative, "two.sided") + expect_identical(attr(t2$conf.int, "conf.level"), 0.95) expect_equal(t3$statistic, c(Z=-1.39077002573558)) expect_equal(t3$p.value, 0.164295175223054) + expect_equal(t3$conf.int[[1]], -0.2876917446341914) + expect_equal(t3$conf.int[[2]], 0.0488706064228094) expect_match(t3$method, "DeLong") expect_match(t3$method, "correlated") expect_identical(t3$alternative, "two.sided") + expect_identical(attr(t3$conf.int, "conf.level"), 0.95) }) test_that("two.sided roc.test produces identical p values when roc curves are reversed", {