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
ci.coords <- function(...) {
21 8
  UseMethod("ci.coords")
22
}
23

24
ci.coords.formula <- function(formula, data, ...) {
25 4
	data.missing <- missing(data)
26 4
	roc.data <- roc.utils.extract.formula(formula, data, ..., 
27 4
										  data.missing = data.missing,
28 4
										  call = match.call())
29 4
	if (length(roc.data$predictor.name) > 1) {
30 0
		stop("Only one predictor supported in 'ci.coords'.")
31
	}
32 4
	response <- roc.data$response
33 4
	predictor <- roc.data$predictors[, 1]
34 4
  	ci.coords(roc(response, predictor, ci=FALSE, ...), ...)
35
}
36

37
ci.coords.default <- function(response, predictor, ...) {
38 8
	if (methods::is(response, "multiclass.roc") || methods::is(response, "multiclass.auc")) {
39 8
		stop("'ci.coords' not available for multiclass ROC curves.")
40
	}
41 4
	roc <- roc.default(response, predictor, ci = FALSE, ...)
42 4
	if (methods::is(roc, "smooth.roc")) {
43 0
		return(ci.coords(smooth.roc = roc, ...))
44
	}
45
	else {
46 4
		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 8
  if (conf.level > 1 | conf.level < 0)
62 0
    stop("'conf.level' must be within the interval [0,1].")
63
  
64 8
  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 8
  input <- match.arg(input)
69 8
  ret <- roc.utils.match.coords.ret.args(ret)
70 8
  best.policy <- match.arg(best.policy)
71 8
  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 8
  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 8
  roc <- attr(smooth.roc, "roc")
84 8
  roc$ci <- NULL # remove potential ci in roc to avoid infinite loop with smooth.roc()
85

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

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

92 8
  if (boot.stratified) {
93 8
    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 8
    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 8
  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 8
  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 8
  ci <- sapply(ret, function(x) t(summarized.perfs[,,x]), simplify = FALSE)
111
  
112 8
  class(ci) <- c("ci.coords", "ci", class(ci))
113 8
  attr(ci, "input") <- input
114 8
  attr(ci, "x") <- x
115 8
  attr(ci, "ret") <- ret
116 8
  attr(ci, "conf.level") <- conf.level
117 8
  attr(ci, "boot.n") <- boot.n
118 8
  attr(ci, "boot.stratified") <- boot.stratified
119 8
  attr(ci, "roc") <- smooth.roc
120 8
  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 8
  if (conf.level > 1 | conf.level < 0)
135 0
    stop("'conf.level' must be within the interval [0,1].")
136
  
137 8
  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 8
  input <- match.arg(input)
142
  
143 8
  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 8
  	ret <- roc.utils.match.coords.ret.args(ret)
149
  }
150
  
151 8
  best.policy <- match.arg(best.policy)
152 8
  if (is.character(x)) {
153 8
    x <- match.arg(x, c("all", "local maximas", "best"))
154 8
    if (x == "all" || x == "local maximas") {
155 0
      stop("'all' and 'local maximas' are not available for confidence intervals.")
156
    }
157
  }
158
  
159 8
  if ("threshold" %in% ret && ! (identical(x, "best") || input == "threshold")) {
160 8
  	stop("'threshold' output is only supported for best ROC point ('x = \"best\"') or if \"threshold\" was given as input.")
161
  }
162

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

166 8
  if (boot.stratified) {
167 8
    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 8
    perfs <- raply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights, best.policy), .progress=progress, .drop = FALSE)
171
  }
172

173 8
  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 8
  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 8
  ci <- sapply(ret, function(x) t(summarized.perfs[,,x]), simplify = FALSE)
185

186 8
  class(ci) <- c("ci.coords", "ci", class(ci))
187 8
  attr(ci, "input") <- input
188 8
  attr(ci, "x") <- x
189 8
  attr(ci, "ret") <- ret
190 8
  attr(ci, "conf.level") <- conf.level
191 8
  attr(ci, "boot.n") <- boot.n
192 8
  attr(ci, "boot.stratified") <- boot.stratified
193 8
  attr(ci, "roc") <- roc
194 8
  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

Read our documentation on viewing source code .

Loading