ci.coords.R at 9240f16 ⋅ xrobin/pROC

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

ci.coords <- function(...) {
  UseMethod("ci.coords")
}

ci.coords.formula <- function(formula, data, ...) {
	data.missing <- missing(data)
	roc.data <- roc.utils.extract.formula(formula, data, ..., 
										  data.missing = data.missing,
										  call = match.call())
	if (length(roc.data$predictor.name) > 1) {
		stop("Only one predictor supported in 'ci.coords'.")
	}
	response <- roc.data$response
	predictor <- roc.data$predictors[, 1]
  	ci.coords(roc(response, predictor, ci=FALSE, ...), ...)
}

ci.coords.default <- function(response, predictor, ...) {
	if (methods::is(response, "multiclass.roc") || methods::is(response, "multiclass.auc")) {
		stop("'ci.coords' not available for multiclass ROC curves.")
	}
	roc <- roc.default(response, predictor, ci = FALSE, ...)
	if (methods::is(roc, "smooth.roc")) {
		return(ci.coords(smooth.roc = roc, ...))
	}
	else {
		return(ci.coords(roc = roc, ...))
	}
}

ci.coords.smooth.roc <- function(smooth.roc,
											x, 
											input=c("specificity", "sensitivity"), ret=c("specificity", "sensitivity"), 
											best.method=c("youden", "closest.topleft"), best.weights=c(1, 0.5),
											best.policy = c("stop", "omit", "random"),
											conf.level = 0.95,
											boot.n = 2000,
											boot.stratified = TRUE,
											progress = getOption("pROCProgress")$name,
											...
                      ) {
  if (conf.level > 1 | conf.level < 0)
    stop("'conf.level' must be within the interval [0,1].")
  
  if (roc.utils.is.perfect.curve(smooth.roc)) {
  	warning("ci.coords() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
  }
 
  input <- match.arg(input)
  ret <- roc.utils.match.coords.ret.args(ret)
  best.policy <- match.arg(best.policy)
  if (is.character(x)) {
    x <- match.arg(x, c("all", "local maximas", "best"))
    if (x == "all" || x == "local maximas") {
      stop("'all' and 'local maximas' are not available for confidence intervals.")
    }
  }

  # Check if called with density.cases or density.controls
  if (is.null(smooth.roc$smoothing.args) || is.numeric(smooth.roc$smoothing.args$density.cases) || is.numeric(smooth.roc$smoothing.args$density.controls))
    stop("Cannot compute CI of ROC curves smoothed with numeric density.controls and density.cases.")

  # Get the non smoothed roc.
  roc <- attr(smooth.roc, "roc")
  roc$ci <- NULL # remove potential ci in roc to avoid infinite loop with smooth.roc()

  # prepare the calls
  smooth.roc.call <- as.call(c(utils::getS3method("smooth", "roc"), smooth.roc$smoothing.args))

  if(class(progress) != "list")
    progress <- roc.utils.get.progress.bar(progress, title="Coords confidence interval", label="Bootstrap in progress...", ...)

  if (boot.stratified) {
    perfs <- raply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call, best.policy), .progress=progress, .drop=FALSE)
  }
  else {
    perfs <- raply(boot.n, nonstratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights,smooth.roc.call, best.policy), .progress=progress, .drop=FALSE)
  }

  if (any(which.ones <- apply(perfs, 1, function(x) all(is.na(x))))) {
  	if (all(which.ones)) {
  		warning("All bootstrap iterations produced NA values only.")
  	}
  	else {
  		how.many <- sum(which.ones)
  		warning(sprintf("%s NA value(s) produced during bootstrap were ignored.", how.many))
  	}
  }

  summarized.perfs <- apply(perfs, c(2, 3), quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2), na.rm=TRUE)
  ci <- sapply(ret, function(x) t(summarized.perfs[,,x]), simplify = FALSE)
  
  class(ci) <- c("ci.coords", "ci", class(ci))
  attr(ci, "input") <- input
  attr(ci, "x") <- x
  attr(ci, "ret") <- ret
  attr(ci, "conf.level") <- conf.level
  attr(ci, "boot.n") <- boot.n
  attr(ci, "boot.stratified") <- boot.stratified
  attr(ci, "roc") <- smooth.roc
  return(ci)
}

ci.coords.roc <- function(roc,
								  x,
								  input=c("threshold", "specificity", "sensitivity"), ret=c("threshold", "specificity", "sensitivity"),
								  best.method=c("youden", "closest.topleft"), best.weights=c(1, 0.5),
								  best.policy = c("stop", "omit", "random"),
								  conf.level = 0.95,
								  boot.n = 2000,
								  boot.stratified = TRUE,
								  progress = getOption("pROCProgress")$name,
                      ...
                      ) {
  if (conf.level > 1 | conf.level < 0)
    stop("'conf.level' must be within the interval [0,1].")
  
  if (roc.utils.is.perfect.curve(roc)) {
  	warning("ci.coords() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
  }
 
  input <- match.arg(input)
  
  if (missing(ret) && input != "threshold") {
  	# Don't show NA thresholds by default
  	ret <- roc.utils.match.coords.ret.args(ret, threshold = FALSE)
  }
  else {
  	ret <- roc.utils.match.coords.ret.args(ret)
  }
  
  best.policy <- match.arg(best.policy)
  if (is.character(x)) {
    x <- match.arg(x, c("all", "local maximas", "best"))
    if (x == "all" || x == "local maximas") {
      stop("'all' and 'local maximas' are not available for confidence intervals.")
    }
  }
  
  if ("threshold" %in% ret && ! (identical(x, "best") || input == "threshold")) {
  	stop("'threshold' output is only supported for best ROC point ('x = \"best\"') or if \"threshold\" was given as input.")
  }

  if(class(progress) != "list")
    progress <- roc.utils.get.progress.bar(progress, title="Coords confidence interval", label="Bootstrap in progress...", ...)

  if (boot.stratified) {
    perfs <- raply(boot.n, stratified.ci.coords(roc, x, input, ret, best.method, best.weights, best.policy), .progress=progress, .drop = FALSE)
  }
  else {
    perfs <- raply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights, best.policy), .progress=progress, .drop = FALSE)
  }

  if (any(which.ones <- apply(perfs, 1, function(x) all(is.na(x))))) {
  	if (all(which.ones)) {
  		warning("All bootstrap iterations produced NA values only.")
  	}
  	else {
  		how.many <- sum(which.ones)
  		warning(sprintf("%s NA value(s) produced during bootstrap were ignored.", how.many))
  	}
  }
  
  summarized.perfs <- apply(perfs, c(2, 3), quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2), na.rm=TRUE)
  ci <- sapply(ret, function(x) t(summarized.perfs[,,x]), simplify = FALSE)

  class(ci) <- c("ci.coords", "ci", class(ci))
  attr(ci, "input") <- input
  attr(ci, "x") <- x
  attr(ci, "ret") <- ret
  attr(ci, "conf.level") <- conf.level
  attr(ci, "boot.n") <- boot.n
  attr(ci, "boot.stratified") <- boot.stratified
  attr(ci, "roc") <- roc
  return(ci)
}

# Function to be called when "best" threshold returned more than 1 column
# Will follow the action defined by best.policy
# For instance:
#   if (x == "best" && nrow(res) != 1) {
# return(enforce.best.policy(res, best.policy))
# }
enforce.best.policy <- function(res, best.policy) {
	if (best.policy == "stop") {
		stop("More than one \"best\" threshold was found, aborting. Change 'best.policy' to alter this behavior.")
	}
	else if (best.policy == "omit") {
		res[1, ] <- NA
		return(res[1, drop = FALSE])
	}
	else {
		return(res[sample(seq_len(nrow(res)), size = 1), , drop = FALSE])
	}
}


Read our documentation on viewing source code .

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

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		ci.coords <- function(...) {
21	1	UseMethod("ci.coords")
22		}
23
24		ci.coords.formula <- function(formula, data, ...) {
25	1	data.missing <- missing(data)
26	1	roc.data <- roc.utils.extract.formula(formula, data, ...,
27	1	data.missing = data.missing,
28	1	call = match.call())
29	1	if (length(roc.data$predictor.name) > 1) {
30	0	stop("Only one predictor supported in 'ci.coords'.")
31		}
32	1	response <- roc.data$response
33	1	predictor <- roc.data$predictors[, 1]
34	1	ci.coords(roc(response, predictor, ci=FALSE, ...), ...)
35		}
36
37		ci.coords.default <- function(response, predictor, ...) {
38	1	if (methods::is(response, "multiclass.roc") \|\| methods::is(response, "multiclass.auc")) {
39	1	stop("'ci.coords' not available for multiclass ROC curves.")
40		}
41	1	roc <- roc.default(response, predictor, ci = FALSE, ...)
42	1	if (methods::is(roc, "smooth.roc")) {
43	0	return(ci.coords(smooth.roc = roc, ...))
44		}
45		else {
46	1	return(ci.coords(roc = roc, ...))
47		}
48		}
49
50		ci.coords.smooth.roc <- function(smooth.roc,
51		x,
52		input=c("specificity", "sensitivity"), ret=c("specificity", "sensitivity"),
53		best.method=c("youden", "closest.topleft"), best.weights=c(1, 0.5),
54		best.policy = c("stop", "omit", "random"),
55		conf.level = 0.95,
56		boot.n = 2000,
57		boot.stratified = TRUE,
58		progress = getOption("pROCProgress")$name,
59		...
60		) {
61	1	if (conf.level > 1 \| conf.level < 0)
62	0	stop("'conf.level' must be within the interval [0,1].")
63
64	1	if (roc.utils.is.perfect.curve(smooth.roc)) {
65	0	warning("ci.coords() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
66		}
67
68	1	input <- match.arg(input)
69	1	ret <- roc.utils.match.coords.ret.args(ret)
70	1	best.policy <- match.arg(best.policy)
71	1	if (is.character(x)) {
72	0	x <- match.arg(x, c("all", "local maximas", "best"))
73	0	if (x == "all" \|\| x == "local maximas") {
74	0	stop("'all' and 'local maximas' are not available for confidence intervals.")
75		}
76		}
77
78		# Check if called with density.cases or density.controls
79	1	if (is.null(smooth.roc$smoothing.args) \|\| is.numeric(smooth.roc$smoothing.args$density.cases) \|\| is.numeric(smooth.roc$smoothing.args$density.controls))
80	0	stop("Cannot compute CI of ROC curves smoothed with numeric density.controls and density.cases.")
81
82		# Get the non smoothed roc.
83	1	roc <- attr(smooth.roc, "roc")
84	1	roc$ci <- NULL # remove potential ci in roc to avoid infinite loop with smooth.roc()
85
86		# prepare the calls
87	1	smooth.roc.call <- as.call(c(utils::getS3method("smooth", "roc"), smooth.roc$smoothing.args))
88
89	1	if(class(progress) != "list")
90	1	progress <- roc.utils.get.progress.bar(progress, title="Coords confidence interval", label="Bootstrap in progress...", ...)
91
92	1	if (boot.stratified) {
93	1	perfs <- raply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call, best.policy), .progress=progress, .drop=FALSE)
94		}
95		else {
96	1	perfs <- raply(boot.n, nonstratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights,smooth.roc.call, best.policy), .progress=progress, .drop=FALSE)
97		}
98
99	1	if (any(which.ones <- apply(perfs, 1, function(x) all(is.na(x))))) {
100	0	if (all(which.ones)) {
101	0	warning("All bootstrap iterations produced NA values only.")
102		}
103		else {
104	0	how.many <- sum(which.ones)
105	0	warning(sprintf("%s NA value(s) produced during bootstrap were ignored.", how.many))
106		}
107		}
108
109	1	summarized.perfs <- apply(perfs, c(2, 3), quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2), na.rm=TRUE)
110	1	ci <- sapply(ret, function(x) t(summarized.perfs[,,x]), simplify = FALSE)
111
112	1	class(ci) <- c("ci.coords", "ci", class(ci))
113	1	attr(ci, "input") <- input
114	1	attr(ci, "x") <- x
115	1	attr(ci, "ret") <- ret
116	1	attr(ci, "conf.level") <- conf.level
117	1	attr(ci, "boot.n") <- boot.n
118	1	attr(ci, "boot.stratified") <- boot.stratified
119	1	attr(ci, "roc") <- smooth.roc
120	1	return(ci)
121		}
122
123		ci.coords.roc <- function(roc,
124		x,
125		input=c("threshold", "specificity", "sensitivity"), ret=c("threshold", "specificity", "sensitivity"),
126		best.method=c("youden", "closest.topleft"), best.weights=c(1, 0.5),
127		best.policy = c("stop", "omit", "random"),
128		conf.level = 0.95,
129		boot.n = 2000,
130		boot.stratified = TRUE,
131		progress = getOption("pROCProgress")$name,
132		...
133		) {
134	1	if (conf.level > 1 \| conf.level < 0)
135	0	stop("'conf.level' must be within the interval [0,1].")
136
137	1	if (roc.utils.is.perfect.curve(roc)) {
138	0	warning("ci.coords() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
139		}
140
141	1	input <- match.arg(input)
142
143	1	if (missing(ret) && input != "threshold") {
144		# Don't show NA thresholds by default
145	0	ret <- roc.utils.match.coords.ret.args(ret, threshold = FALSE)
146		}
147		else {
148	1	ret <- roc.utils.match.coords.ret.args(ret)
149		}
150
151	1	best.policy <- match.arg(best.policy)
152	1	if (is.character(x)) {
153	1	x <- match.arg(x, c("all", "local maximas", "best"))
154	1	if (x == "all" \|\| x == "local maximas") {
155	0	stop("'all' and 'local maximas' are not available for confidence intervals.")
156		}
157		}
158
159	1	if ("threshold" %in% ret && ! (identical(x, "best") \|\| input == "threshold")) {
160	1	stop("'threshold' output is only supported for best ROC point ('x = \"best\"') or if \"threshold\" was given as input.")
161		}
162
163	1	if(class(progress) != "list")
164	1	progress <- roc.utils.get.progress.bar(progress, title="Coords confidence interval", label="Bootstrap in progress...", ...)
165
166	1	if (boot.stratified) {
167	1	perfs <- raply(boot.n, stratified.ci.coords(roc, x, input, ret, best.method, best.weights, best.policy), .progress=progress, .drop = FALSE)
168		}
169		else {
170	1	perfs <- raply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights, best.policy), .progress=progress, .drop = FALSE)
171		}
172
173	1	if (any(which.ones <- apply(perfs, 1, function(x) all(is.na(x))))) {
174	0	if (all(which.ones)) {
175	0	warning("All bootstrap iterations produced NA values only.")
176		}
177		else {
178	0	how.many <- sum(which.ones)
179	0	warning(sprintf("%s NA value(s) produced during bootstrap were ignored.", how.many))
180		}
181		}
182
183	1	summarized.perfs <- apply(perfs, c(2, 3), quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2), na.rm=TRUE)
184	1	ci <- sapply(ret, function(x) t(summarized.perfs[,,x]), simplify = FALSE)
185
186	1	class(ci) <- c("ci.coords", "ci", class(ci))
187	1	attr(ci, "input") <- input
188	1	attr(ci, "x") <- x
189	1	attr(ci, "ret") <- ret
190	1	attr(ci, "conf.level") <- conf.level
191	1	attr(ci, "boot.n") <- boot.n
192	1	attr(ci, "boot.stratified") <- boot.stratified
193	1	attr(ci, "roc") <- roc
194	1	return(ci)
195		}
196
197		# Function to be called when "best" threshold returned more than 1 column
198		# Will follow the action defined by best.policy
199		# For instance:
200		# if (x == "best" && nrow(res) != 1) {
201		# return(enforce.best.policy(res, best.policy))
202		# }
203		enforce.best.policy <- function(res, best.policy) {
204	0	if (best.policy == "stop") {
205	0	stop("More than one \"best\" threshold was found, aborting. Change 'best.policy' to alter this behavior.")
206		}
207	0	else if (best.policy == "omit") {
208	0	res[1, ] <- NA
209	0	return(res[1, drop = FALSE])
210		}
211		else {
212	0	return(res[sample(seq_len(nrow(res)), size = 1), , drop = FALSE])
213		}
214		}
215