Merge pull request #95 from wz-billings/master

Add confidence interval output to paired DeLong test (issue #73)

R/delong.R

@@ -18,31 +18,15 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # 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) {

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))

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

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", {