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

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

37
ci.se.default <- function(response, predictor, ...) {
38 33
	if (methods::is(response, "multiclass.roc") || methods::is(response, "multiclass.auc")) {
39 33
		stop("'ci.sp' not available for multiclass ROC curves.")
40
	}
41 33
	roc <- roc.default(response, predictor, ci = FALSE, ...)
42 33
	if (methods::is(roc, "smooth.roc")) {
43 0
		return(ci.se(smooth.roc = roc, ...))
44
	}
45
	else {
46 33
		return(ci.se(roc = roc, ...))
47
	}
48
}
49

50
ci.se.smooth.roc <- function(smooth.roc,
51
                      specificities = seq(0, 1, .1) * ifelse(smooth.roc$percent, 100, 1),
52
                      conf.level = 0.95,
53
                      boot.n = 2000,
54
                      boot.stratified = TRUE,
55
                      progress = getOption("pROCProgress")$name,
56
                      parallel = FALSE,
57
                      ...
58
                      ) {
59 33
  if (conf.level > 1 | conf.level < 0)
60 0
    stop("'conf.level' must be within the interval [0,1].")
61
  
62 33
  if (roc.utils.is.perfect.curve(smooth.roc)) {
63 0
  	warning("ci.se() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
64
  }
65

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

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

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

77 33
  if(class(progress) != "list")
78 33
    progress <- roc.utils.get.progress.bar(progress, title="SE confidence interval", label="Bootstrap in progress...", ...)
79

80 33
  if (boot.stratified) {
81 33
    perfs <- ldply(1:boot.n, stratified.ci.smooth.se, roc=roc, sp=specificities, smooth.roc.call=smooth.roc.call, .progress=progress, .parallel=parallel)
82
  }
83
  else {
84 33
    perfs <- ldply(1:boot.n, nonstratified.ci.smooth.se, roc=roc, sp=specificities, smooth.roc.call=smooth.roc.call, .progress=progress, .parallel=parallel)
85
  }
86

87 33
  if (any(is.na(perfs))) {
88 0
    warning("NA value(s) produced during bootstrap were ignored.")
89 0
    perfs <- perfs[!apply(perfs, 1, function(x) any(is.na(x))),]
90
  }
91

92 33
  ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
93 33
  rownames(ci) <- paste(specificities, ifelse(roc$percent, "%", ""), sep="")
94

95 33
  class(ci) <- c("ci.se", "ci", class(ci))
96 33
  attr(ci, "conf.level") <- conf.level
97 33
  attr(ci, "boot.n") <- boot.n
98 33
  attr(ci, "boot.stratified") <- boot.stratified
99 33
  attr(ci, "specificities") <- specificities
100 33
  attr(ci, "roc") <- smooth.roc
101 33
  return(ci)
102
}
103

104
ci.se.roc <- function(roc,
105
                      specificities = seq(0, 1, .1) * ifelse(roc$percent, 100, 1),
106
                      conf.level = 0.95,
107
                      boot.n = 2000,
108
                      boot.stratified = TRUE,
109
                      progress = getOption("pROCProgress")$name,
110
                      parallel = FALSE,
111
                      ...
112
                      ) {
113 33
  if (conf.level > 1 | conf.level < 0)
114 0
    stop("'conf.level' must be within the interval [0,1].")
115
  
116 33
  if (roc.utils.is.perfect.curve(roc)) {
117 0
  	warning("ci.se() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
118
  }
119

120 33
  if(class(progress) != "list")
121 33
    progress <- roc.utils.get.progress.bar(progress, title="SE confidence interval", label="Bootstrap in progress...", ...)
122

123 33
  if (boot.stratified) {
124 33
    perfs <- ldply(1:boot.n, stratified.ci.se, roc=roc, sp=specificities, .progress=progress, .parallel=parallel)
125
  }
126
  else {
127 33
    perfs <- ldply(1:boot.n, nonstratified.ci.se, roc=roc, sp=specificities, .progress=progress, .parallel=parallel)
128
  }
129

130 33
  if (any(is.na(perfs))) {
131 0
    warning("NA value(s) produced during bootstrap were ignored.")
132 0
    perfs <- perfs[!apply(perfs, 1, function(x) any(is.na(x))),]
133
  }
134 33
  ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
135 33
  rownames(ci) <- paste(specificities, ifelse(roc$percent, "%", ""), sep="")
136
  
137 33
  class(ci) <- c("ci.se", "ci", class(ci))
138 33
  attr(ci, "conf.level") <- conf.level
139 33
  attr(ci, "boot.n") <- boot.n
140 33
  attr(ci, "boot.stratified") <- boot.stratified
141 33
  attr(ci, "specificities") <- specificities
142 33
  attr(ci, "roc") <- roc
143 33
  return(ci)
144
}

Read our documentation on viewing source code .

Loading