xrobin / pROC
1
# pROC: Tools Receiver operating characteristic (ROC curves) with
2
# (partial) area under the curve, confidence intervals and comparison. 
3
# Copyright (C) 2010-2014 Xavier Robin, Alexandre Hainard, Natacha Turck,
4
# Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez
5
# and Markus Müller
6
#
7
# This program is free software: you can redistribute it and/or modify
8
# it under the terms of the GNU General Public License as published by
9
# the Free Software Foundation, either version 3 of the License, or
10
# (at your option) any later version.
11
#
12
# This program is distributed in the hope that it will be useful,
13
# but WITHOUT ANY WARRANTY; without even the implied warranty of
14
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15
# GNU General Public License for more details.
16
#
17
# You should have received a copy of the GNU General Public License
18
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
19

20
roc.test <- function(...) {
21 8
	UseMethod("roc.test")
22
}
23

24
roc.test.formula <- function (formula, data, ...) {
25 8
	data.missing <- missing(data)
26 8
	call <- match.call()
27 8
	roc.data <- roc.utils.extract.formula(formula, data, ..., 
28 8
										  data.missing = data.missing,
29 8
										  call = call)
30 8
	if (length(roc.data$predictor.name) != 2) {
31 0
		stop("Invalid formula: exactly 2 predictors are required in a formula of type response~predictor1+predictor2.")
32
	}
33 8
	response <- roc.data$response
34 8
	predictors <- roc.data$predictors
35
	
36 8
	testres <- roc.test.default(response, predictors, ...)
37 8
	testres$call <- call
38
	# data.names for pretty print()ing
39 8
	if (data.missing) {
40 8
		testres$data.names <- sprintf("%s and %s by %s (%s, %s)", roc.data$predictor.names[1], roc.data$predictor.names[2], roc.data$response.name, testres$roc1$levels[1], testres$roc1$levels[2])
41
	}
42
	else {
43 8
		testres$data.names <- sprintf("%s and %s in %s by %s (%s, %s)", roc.data$predictor.names[1], roc.data$predictor.names[2], deparse(substitute(data)), roc.data$response.name, testres$roc1$levels[1], testres$roc1$levels[2])
44
	}
45
	
46 8
	return(testres)
47
}
48

49
roc.test.default <- function(response, predictor1, predictor2=NULL, na.rm=TRUE, method=NULL, ...) {
50 8
	if (is.matrix(predictor1) | is.data.frame(predictor1)) {
51 8
		if (!is.null(predictor2))
52 8
			stop("Predictor2 must not be specified if predictor1 is a matrix or a data.frame.")
53 8
		if (dim(predictor1)[2] == 2 & length(response) == dim(predictor1)[1]) {
54 8
			roc1 <- roc(response, predictor1[,1], ...)
55 8
			roc2 <- roc(response, predictor1[,2], ...)
56 8
			if (!is.null(names(predictor1)))
57 8
				data.names <- sprintf("%s and %s in %s by %s (%s, %s)", names(predictor1)[1], names(predictor1)[2], deparse(substitute(predictor1)), deparse(substitute(response)), roc1$levels[1], roc1$levels[2])
58 8
			else if (!is.null(colnames(predictor1)))
59 8
				data.names <- sprintf("%s and %s in %s by %s (%s, %s)", colnames(predictor1)[1], colnames(predictor1)[2], deparse(substitute(predictor1)), deparse(substitute(response)), roc1$levels[1], roc1$levels[2])
60
			else
61 8
				data.names <- sprintf("%s by %s (%s, %s)", deparse(substitute(predictor1)), deparse(substitute(response)), roc1$levels[1], roc1$levels[2])
62
		}
63
		else {
64 0
			stop("Wrong dimension for predictor1 as a matrix or a data.frame.")
65
		}
66
	}
67
	else {
68 8
		if (missing(predictor2))
69 8
			stop("Missing argument predictor2 with predictor1 as a vector.")
70
		# Need to remove NAs
71 8
		if (na.rm) {
72 8
			nas <- is.na(response) | is.na(predictor1) | is.na(predictor2)
73 8
			response <- response[!nas]
74 8
			predictor1 <- predictor1[!nas]
75 8
			predictor2 <- predictor2[!nas]
76
		}
77 8
		roc1 <- roc(response, predictor1, ...)
78 8
		roc2 <- roc(response, predictor2, ...)
79 8
		call <- match.call()
80 8
		data.names <- sprintf("%s and %s by %s (%s, %s)", deparse(call$predictor1), deparse(call$predictor2), deparse(call$response), roc1$levels[1], roc1$levels[2])
81
	}
82 8
	test <- roc.test.roc(roc1, roc2, method=method, ...)
83 8
	test$data.names <- data.names
84 8
	return(test)
85
}
86

87
roc.test.auc <- function(roc1, roc2, ...) {
88
	# First save the names
89 4
	data.names <- paste(deparse(substitute(roc1)), "and", deparse(substitute(roc2)))
90
	# Change roc1 from an auc to a roc object but keep the auc specifications
91 4
	auc1 <- roc1
92 4
	attr(auc1, "roc") <- NULL
93 4
	roc1 <- attr(roc1, "roc")
94 4
	roc1$auc <- auc1
95
	# Pass to roc.test.roc
96 4
	testres <- roc.test.roc(roc1, roc2, ...)
97 4
	testres$call <- match.call()
98 4
	testres$data.names <- data.names
99 4
	return(testres)
100
}
101

102
roc.test.smooth.roc <- function(roc1, roc2, ...) {
103 4
	testres <- roc.test.roc(roc1, roc2, ...)
104 4
	testres$call <- match.call()
105 4
	testres$data.names <- paste(deparse(substitute(roc1)), "and", deparse(substitute(roc2)))
106 4
	return(testres)
107
}
108

109
roc.test.roc <- function(roc1, roc2,
110
						 method=c("delong", "bootstrap", "venkatraman", "sensitivity", "specificity"),
111
						 sensitivity=NULL, specificity=NULL,
112
						 alternative = c("two.sided", "less", "greater"),
113
						 paired=NULL,
114
						 reuse.auc=TRUE,
115
						 boot.n=2000, boot.stratified=TRUE,
116
						 ties.method="first",
117
						 progress=getOption("pROCProgress")$name,
118
						 parallel=FALSE,
119
						 conf.level=0.95,
120
						 ...) {
121 8
	alternative <- match.arg(alternative)
122 8
	data.names <- paste(deparse(substitute(roc1)), "and", deparse(substitute(roc2)))
123
	# If roc2 is an auc, take the roc but keep the auc specifications
124 8
	if (methods::is(roc2, "auc")) {
125 4
		auc2 <- roc2
126 4
		attr(auc2, "roc") <- NULL
127 4
		roc2 <- attr(roc2, "roc")
128 4
		roc2$auc <- auc2
129
	}
130
	
131 8
	if (roc.utils.is.perfect.curve(roc1) && roc.utils.is.perfect.curve(roc2)) {
132 0
		warning("roc.test() of two ROC curves with AUC == 1 has always p.value = 1 and can be misleading.")
133
	}
134
	
135
	# store which objects are smoothed, and how
136 8
	smoothing.args <- list()
137 8
	if (methods::is(roc1, "smooth.roc")) {
138 4
		smoothing.args$roc1 <- roc1$smoothing.args
139 4
		smoothing.args$roc1$smooth <- TRUE
140 4
		roc1 <- attr(roc1, "roc")
141
	}
142
	else {
143 8
		smoothing.args$roc1 <- list(smooth=FALSE)
144
	}
145 8
	if (methods::is(roc2, "smooth.roc")) {
146 4
		smoothing.args$roc2 <- roc2$smoothing.args
147 4
		smoothing.args$roc2$smooth <- TRUE
148 4
		roc2 <- attr(roc2, "roc")
149
	}
150
	else {
151 8
		smoothing.args$roc2 <- list(smooth=FALSE)
152
	}
153
	
154
	# Check if we do a paired or unpaired roc.test
155 8
	if (is.null(paired)) {
156
		# then determine whether the rocs are paired or not
157 8
		rocs.are.paired <- are.paired(roc1, roc2, return.paired.rocs=TRUE, reuse.auc=TRUE, reuse.ci=FALSE, reuse.smooth=TRUE)
158 8
		if (rocs.are.paired) {
159 8
			paired <- TRUE
160 8
			roc1 <- attr(rocs.are.paired, "roc1")
161 8
			roc2 <- attr(rocs.are.paired, "roc2")
162
		}
163
		else {
164 8
			paired <- FALSE
165 8
			roc1 <- roc1
166 8
			roc2 <- roc2
167
		}
168
	}
169 8
	else if (paired) {
170
		# make sure the rocs are really paired
171 8
		rocs.are.paired <- rocs.are.paired <- are.paired(roc1, roc2, return.paired.rocs=TRUE, reuse.auc=TRUE, reuse.ci=FALSE, reuse.smooth=TRUE)
172 8
		if (! rocs.are.paired) 
173 8
			stop("The paired ROC test cannot be applied to unpaired curves.")
174 4
		roc1 <- attr(rocs.are.paired, "roc1")
175 4
		roc2 <- attr(rocs.are.paired, "roc2")
176
	}
177 8
	else { # assume unpaired
178 8
		rocs.are.paired <- are.paired(roc1, roc2, return.paired.rocs=FALSE)
179 8
		if (rocs.are.paired) 
180 8
			warning("The ROC curves seem to be paired. Consider performing a paired roc.test.")
181 8
		roc1 <- roc1
182 8
		roc2 <- roc2
183
	}
184
	
185
	# check that the AUC was computed, or do it now
186 8
	if (is.null(roc1$auc) | !reuse.auc) {
187 0
		if (smoothing.args$roc1$smooth) {
188 0
			roc1$auc <- auc(smooth.roc=do.call("smooth.roc", c(list(roc=roc1), smoothing.args$roc1)), ...)
189
			# remove partial.auc.* arguments that are now in roc1$auc and that will mess later processing
190
			# (formal argument "partial.auc(.*)" matched by multiple actual arguments)
191
			# This removal should be safe because we always use smoothing.args with roc1 in the following processing,
192
			# however it is a potential source of bugs.
193 0
			smoothing.args$roc1$partial.auc <- NULL
194 0
			smoothing.args$roc1$partial.auc.correct <- NULL
195 0
			smoothing.args$roc1$partial.auc.focus <- NULL
196
		}
197
		else
198 0
			roc1$auc <- auc(roc1, ...)
199
	}
200 8
	if (is.null(roc2$auc) | !reuse.auc) {
201 0
		if (smoothing.args$roc2$smooth) {
202 0
			roc2$auc <- auc(smooth.roc=do.call("smooth.roc", c(list(roc=roc2), smoothing.args$roc2)), ...)
203
			# remove partial.auc.* arguments that are now in roc1$auc and that will mess later processing
204
			# (formal argument "partial.auc(.*)" matched by multiple actual arguments)
205
			# This removal should be safe because we always use smoothing.args with roc2 in the following processing,
206
			# however it is a potential source of bugs.
207 0
			smoothing.args$roc2$partial.auc <- NULL
208 0
			smoothing.args$roc2$partial.auc.correct <- NULL
209 0
			smoothing.args$roc2$partial.auc.focus <- NULL
210
		}
211
		else
212 0
			roc2$auc <- auc(roc2, ...)
213
	}
214
	
215
	# check that the same region was requested in auc. Otherwise, issue a warning
216 8
	if (!identical(attributes(roc1$auc)[names(attributes(roc1$auc))!="roc"], attributes(roc2$auc)[names(attributes(roc2$auc))!="roc"]))
217 8
		warning("Different AUC specifications in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
218
	# check that the same smoothing params were requested in auc. Otherwise, issue a warning
219 8
	if (!identical(smoothing.args$roc1, smoothing.args$roc2))
220 8
		warning("Different smoothing parameters in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
221
	
222
	# Check the method
223 8
	if (missing(method) | is.null(method)) {
224
		# determine method if missing
225 8
		if (has.partial.auc(roc1)) {
226
			# partial auc: go for bootstrap
227 0
			method <- "bootstrap"
228
		}
229 8
		else if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth) {
230
			# smoothing in one or both: bootstrap
231 0
			method <- "bootstrap"
232
		}
233 8
		else if (roc1$direction != roc2$direction) {
234
			# delong doesn't work well with opposite directions (will report high significance if roc1$auc and roc2$auc are similar and high)
235 0
			method <- "bootstrap"
236
		}
237
		else {
238 8
			method <- "delong"
239
		}
240
	}
241
	else {
242 4
		method <- match.arg(method)
243 4
		if (method == "delong") {
244
			# delong NA to pAUC: warn + change
245 0
			if (has.partial.auc(roc1) || has.partial.auc(roc2)) {
246 0
				stop("DeLong's test is not supported for partial AUC. Use method=\"bootstrap\" instead.")
247
			}
248 0
			if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth) {
249 0
				stop("DeLong's test is not supported for smoothed ROCs. Use method=\"bootstrap\" instead.")
250
			}
251 0
			if (roc1$direction != roc2$direction)
252 0
				warning("DeLong's test should not be applied to ROC curves with a different direction.")
253
			
254
			# Check if conf.level is specified correctly. This is currently
255
			# only used for the delong paired method, which is why it lives 
256
			# here for now.
257 0
			if (!is.numeric(conf.level)) {
258 0
				stop("conf.level must be numeric between 0 and 1.")
259 0
			} else if (0 > conf.level | 1 < conf.level) {
260 0
				stop("conf.level must be between 0 and 1.")
261
			}
262
		}
263 4
		else if (method == "venkatraman") {
264 4
			if (has.partial.auc(roc1))
265 4
				stop("Partial AUC is not supported for Venkatraman's test.")
266 4
			if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth)
267 4
				stop("Venkatraman's test is not supported for smoothed ROCs")
268 4
			if (roc1$direction != roc2$direction)
269 4
				warning("Venkatraman's test should not be applied to ROC curves with different directions.")
270 4
			if (alternative != "two.sided") {
271 0
				stop("Only two-sided tests are available for Venkatraman.")
272
			}
273
		}
274
	}
275
	
276
	# Prepare the return value htest
277 8
	if (smoothing.args$roc1$smooth)
278 8
		estimate <- do.call("smooth.roc", c(list(roc=roc1), smoothing.args$roc1))$auc
279
	else
280 8
		estimate <- roc1$auc
281 8
	if (smoothing.args$roc2$smooth)
282 8
		estimate <- c(estimate, do.call("smooth.roc", c(list(roc=roc2), smoothing.args$roc2))$auc)
283
	else
284 8
		estimate <- c(estimate, roc2$auc)
285 8
	if (identical(attr(roc1$auc, "partial.auc"), FALSE)) {
286 8
		nest <- paste(ifelse(smoothing.args$roc1$smooth, "Smoothed ", ""), "AUC of roc1", sep="")
287
	}
288
	else {
289 4
		nest <- paste(ifelse (attr(roc1$auc, "partial.auc.correct"), "Corrected ", ""),
290 4
					  ifelse (smoothing.args$roc1$smooth, "Smoothed ", ""),
291 4
					  "pAUC (", attr(roc1$auc, "partial.auc")[1], "-", attr(roc1$auc, "partial.auc")[2], " ", attr(roc1$auc, "partial.auc.focus"),
292 4
					  ") of roc1", sep="")
293
	}
294 8
	if (identical(attr(roc2$auc, "partial.auc"), FALSE)) {
295 8
		nest <- c(nest, paste(ifelse(smoothing.args$roc2$smooth, "Smoothed ", ""), "AUC of roc2", sep=""))
296
	}
297
	else {
298 4
		nest <- c(nest, paste(ifelse (attr(roc2$auc, "partial.auc.correct"), "Corrected ", ""),
299 4
							  ifelse (smoothing.args$roc2$smooth, "Smoothed ", ""),
300 4
							  "pAUC (", attr(roc2$auc, "partial.auc")[1], "-", attr(roc2$auc, "partial.auc")[2], " ", attr(roc2$auc, "partial.auc.focus"),
301 4
							  ") of roc2", sep=""))
302
	}
303 8
	nest <- sub("Corrected Smoothed", "Corrected smoothed", nest) # no upper on smoothed if corrected.
304 8
	names(estimate) <- nest
305 8
	null.value <- 0
306 8
	names(null.value) <- "difference in AUC"
307 8
	htest <- list(
308 8
		alternative = alternative,
309 8
		data.names = data.names,
310 8
		estimate = estimate,
311 8
		null.value = null.value
312
	)
313 8
	class(htest) <- "htest"
314
	
315 8
	if (method == "delong") {
316 8
		if (paired) {
317 8
			delong.calcs <- delong.paired.calculations(roc1, roc2)
318 8
			stat <- delong.paired.test(delong.calcs)
319 8
			stat.ci <- ci.delong.paired(delong.calcs, conf.level)
320 8
			names(stat) <- "Z"
321 8
			htest$statistic <- stat
322 8
			htest$method <- "DeLong's test for two correlated ROC curves"
323 8
			htest$conf.int <- c(stat.ci$lower, stat.ci$upper)
324 8
			attr(htest$conf.int, "conf.level") <- stat.ci$level
325
			
326 8
			if (alternative == "two.sided")
327 8
				pval <- 2*pnorm(-abs(stat))
328 8
			else if (alternative == "greater")
329 8
				pval <- pnorm(-stat)
330
			else
331 8
				pval <- pnorm(stat)
332 8
			htest$p.value <- pval
333
		}
334
		else {
335 8
			stats <- delong.unpaired.test(roc1, roc2)
336 8
			stat <- stats[1]
337 8
			df <- stats[2]
338 8
			htest$statistic <- c("D"=stat)
339 8
			htest$parameter <- c("df"=df)
340 8
			htest$method <- "DeLong's test for two ROC curves"
341
			
342 8
			if (alternative == "two.sided")
343 8
				pval <- 2*pt(-abs(stat), df=df)
344 8
			else if (alternative == "greater")
345 8
				pval <- pt(-stat, df=df)
346
			else
347 8
				pval <- pt(stat, df=df)
348 8
			htest$p.value <- pval
349
		}
350
	}
351 8
	else if (method == "venkatraman") {
352 4
		if(class(progress) != "list")
353 4
			progress <- roc.utils.get.progress.bar(progress, title="Venkatraman ROC test", label="Permutations in progress...", ...)
354 4
		if (paired) {
355 4
			stats <- venkatraman.paired.test(roc1, roc2, boot.n, ties.method, progress, parallel)
356 4
			htest$method <- "Venkatraman's test for two paired ROC curves"
357
		}
358
		else {
359 4
			stats <- venkatraman.unpaired.test(roc1, roc2, boot.n, ties.method, progress, parallel)
360 4
			htest$method <- "Venkatraman's test for two unpaired ROC curves"
361
		}
362 4
		stat <- stats[[1]]
363 4
		names(stat) <- "E"
364 4
		htest$statistic <- stat
365 4
		parameter <- c(boot.n)
366 4
		names(parameter) <- "boot.n"
367 4
		htest$parameter <- parameter
368 4
		pval <- sum(stats[[2]]>=stats[[1]])/boot.n
369 4
		htest$p.value <- pval
370 4
		names(htest$null.value) <- "difference in ROC operating points"
371 4
		htest$estimate <- NULL # AUC not relevant in venkatraman
372
	}
373 8
	else { # method == "bootstrap" or "sensitivity" or "specificity"
374
		# Check if called with density.cases or density.controls
375 4
		if (is.null(smoothing.args) || is.numeric(smoothing.args$density.cases) || is.numeric(smoothing.args$density.controls))
376 4
			stop("Cannot compute the statistic on ROC curves smoothed with numeric density.controls and density.cases.")
377
		
378 4
		if(class(progress) != "list")
379 4
			progress <- roc.utils.get.progress.bar(progress, title="Bootstrap ROC test", label="Bootstrap in progress...", ...)
380
		
381 4
		if (method == "specificity") {
382 4
			if (! is.numeric(specificity) || length(specificity) != 1) {
383 0
				stop("Argument 'specificity' must be numeric of length 1 for a specificity test.")
384
			}
385 4
			stat <- bootstrap.test(roc1, roc2, "sp", specificity, paired, boot.n, boot.stratified, smoothing.args, progress, parallel)
386 4
			if (paired)
387 4
				htest$method <- "Specificity test for two correlated ROC curves"
388
			else
389 4
				htest$method <- "Specificity test for two ROC curves"
390 4
			names(htest$null.value) <- sprintf("difference in sensitivity at %s specificity",
391 4
											   specificity)
392
		}
393 4
		else if (method == "sensitivity") {
394 4
			if (! is.numeric(sensitivity) || length(sensitivity) != 1) {
395 0
				stop("Argument 'sensitivity' must be numeric of length 1 for a sensitivity test.")
396
			}
397 4
			stat <- bootstrap.test(roc1, roc2, "se", sensitivity, paired, boot.n, boot.stratified, smoothing.args, progress, parallel)
398 4
			if (paired)
399 4
				htest$method <- "Sensitivity test for two correlated ROC curves"
400
			else
401 4
				htest$method <- "Sensitivity test for two ROC curves"
402
			
403 4
			names(htest$null.value) <- sprintf("difference in specificity at %s sensitivity",
404 4
											   sensitivity)
405
		}
406
		else {
407 4
			stat <- bootstrap.test(roc1, roc2, "boot", NULL, paired, boot.n, boot.stratified, smoothing.args, progress, parallel)
408 4
			if (paired)
409 4
				htest$method <- "Bootstrap test for two correlated ROC curves"
410
			else
411 4
				htest$method <- "Bootstrap test for two ROC curves"
412
		}
413 4
		stat <- as.vector(stat) # remove auc attributes
414 4
		names(stat) <- "D"
415 4
		htest$statistic <- stat
416 4
		parameter <- c(boot.n, boot.stratified)
417 4
		names(parameter) <- c("boot.n", "boot.stratified")
418 4
		htest$parameter <- parameter
419
		
420 4
		if (alternative == "two.sided")
421 4
			pval <- 2*pnorm(-abs(stat))
422 4
		else if (alternative == "greater")
423 4
			pval <- pnorm(-stat)
424
		else
425 4
			pval <- pnorm(stat)
426 4
		htest$p.value <- pval
427
	}
428
	
429 8
	htest$roc1 <- roc1
430 8
	htest$roc2 <- roc2
431
	# Remove name from p value
432 8
	htest$p.value <- unname(htest$p.value)
433
	# Restore smoothing if necessary
434 8
	if (smoothing.args$roc1$smooth)
435 8
		htest$roc1 <- do.call("smooth.roc", c(list(roc=roc1), smoothing.args$roc1))
436 8
	if (smoothing.args$roc2$smooth)
437 8
		htest$roc2 <- do.call("smooth.roc", c(list(roc=roc2), smoothing.args$roc2))
438 8
	return(htest)
439
}

Read our documentation on viewing source code .

Loading