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
##########  AUC of two ROC curves (roc.test, cov)  ##########
21

22
bootstrap.cov <- function(roc1, roc2, boot.n, boot.stratified, boot.return, smoothing.args, progress, parallel) {
23
  # rename method into smooth.method for roc
24 33
  smoothing.args$roc1$smooth.method <- smoothing.args$roc1$method
25 33
  smoothing.args$roc1$method <- NULL
26 33
  smoothing.args$roc2$smooth.method <- smoothing.args$roc2$method
27 33
  smoothing.args$roc2$method <- NULL
28

29
  # Prepare arguments for later calls to roc
30 33
  auc1skeleton <- attributes(roc1$auc)
31 33
  auc1skeleton$roc <- NULL
32 33
  auc1skeleton$direction <- roc1$direction
33 33
  auc1skeleton$class <- NULL
34 33
  auc1skeleton$fun.sesp <- roc1$fun.sesp
35 33
  auc1skeleton$allow.invalid.partial.auc.correct <- TRUE
36 33
  auc1skeleton <- c(auc1skeleton, smoothing.args$roc1)
37 33
  names(auc1skeleton)[which(names(auc1skeleton) == "n")] <-  "smooth.n"
38 33
  auc2skeleton <- attributes(roc2$auc)
39 33
  auc2skeleton$roc <- NULL
40 33
  auc2skeleton$direction <- roc2$direction
41 33
  auc2skeleton$class <- NULL
42 33
  auc2skeleton$fun.sesp <- roc2$fun.sesp
43 33
  auc2skeleton$allow.invalid.partial.auc.correct <- TRUE
44 33
  auc2skeleton <- c(auc2skeleton, smoothing.args$roc2)
45 33
  names(auc2skeleton)[which(names(auc2skeleton) == "n")] <-  "smooth.n"
46

47 33
  auc1skeleton$auc <- auc2skeleton$auc <- TRUE
48

49
  # Some attributes may be duplicated in AUC skeletons and will mess the boostrap later on when we do.call().
50
  # If this condition happen, it probably means we have a bug elsewhere.
51
  # Rather than making a complicated processing to remove the duplicates,
52
  # just throw an error and let us solve the bug when a user reports it.
53 33
  duplicated.auc1skeleton <- duplicated(names(auc1skeleton))
54 33
  duplicated.auc2skeleton <- duplicated(names(auc2skeleton))
55 33
  if (any(duplicated.auc1skeleton)) {
56 0
  	sessionInfo <- sessionInfo()
57 0
  	save(roc1, roc2, boot.n, boot.stratified, boot.return, smoothing.args, progress, parallel, sessionInfo, file="pROC_bug.RData")
58 0
  	stop(sprintf("pROC: duplicated argument(s) in AUC1 skeleton: \"%s\". Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", paste(names(auc1skeleton)[duplicated(names(auc1skeleton))], collapse=", "), utils::packageDescription("pROC")$BugReports))
59
  	
60
  }
61 33
  if (any(duplicated.auc2skeleton)) {
62 0
  	sessionInfo <- sessionInfo()
63 0
  	save(roc1, roc2, boot.n, boot.stratified, boot.return, smoothing.args, progress, parallel, sessionInfo, file="pROC_bug.RData")
64 0
  	stop(sprintf("duplicated argument(s) in AUC2 skeleton: \"%s\". Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", paste(names(auc2skeleton)[duplicated(names(auc2skeleton))], collapse=", "), utils::packageDescription("pROC")$BugReports))
65
  }
66 33
  if (boot.stratified) { # precompute sorted responses if stratified
67
    #response.roc1 <- factor(c(rep(roc1$levels[1], length(roc1$controls)), rep(roc1$levels[2], length(roc1$cases))), levels=roc1$levels)
68
    #response.roc2 <- factor(c(rep(roc2$levels[1], length(roc2$controls)), rep(roc2$levels[2], length(roc2$cases))), levels=roc2$levels)
69
    #auc1skeleton$response <- response.roc1
70
    #auc2skeleton$response <- response.roc2
71 33
    resampled.values <- laply(1:boot.n, stratified.bootstrap.test, roc1=roc1, roc2=roc2, test="boot", x=NULL, paired=TRUE, auc1skeleton=auc1skeleton, auc2skeleton=auc2skeleton, .progress=progress, .parallel=parallel)
72
  }
73
  else {
74 33
    auc1skeleton$levels <- roc1$levels
75 33
    auc1skeleton$direction <- roc1$direction
76 33
    auc2skeleton$levels <- roc2$levels
77 33
    auc2skeleton$direction <- roc2$direction
78 33
    resampled.values <- laply(1:boot.n, nonstratified.bootstrap.test, roc1=roc1, roc2=roc2, test="boot", x=NULL, paired=TRUE, auc1skeleton=auc1skeleton, auc2skeleton=auc2skeleton, .progress=progress, .parallel=parallel)
79
  }
80

81
  # are there NA values?
82 33
  if ((num.NAs <- sum(apply(resampled.values, 1, is.na))) > 0) {
83 0
    warning(sprintf("%i NA value(s) produced during bootstrap were ignored.", num.NAs))
84 0
    resampled.values <- resampled.values[!apply(resampled.values, 1, function(x) any(is.na(x))),]
85
  }
86

87 33
  cov <- stats::cov(resampled.values[,1], resampled.values[,2])
88 33
  if (boot.return) {
89 0
    attr(cov, "resampled.values") <- resampled.values
90
  }
91 33
  return(cov)
92
}
93

94
# Bootstrap test, used by roc.test.roc
95
bootstrap.test <- function(roc1, roc2, test, x, paired, boot.n, boot.stratified, smoothing.args, progress, parallel) {
96
  # rename method into smooth.method for roc
97 33
  smoothing.args$roc1$smooth.method <- smoothing.args$roc1$method
98 33
  smoothing.args$roc1$method <- NULL
99 33
  smoothing.args$roc2$smooth.method <- smoothing.args$roc2$method
100 33
  smoothing.args$roc2$method <- NULL
101

102
  # Prepare arguments for later calls to roc
103 33
  auc1skeleton <- attributes(roc1$auc)
104 33
  auc1skeleton$roc <- NULL
105 33
  auc1skeleton$direction <- roc1$direction
106 33
  auc1skeleton$class <- NULL
107 33
  auc1skeleton$fun.sesp <- roc1$fun.sesp
108 33
  auc1skeleton$allow.invalid.partial.auc.correct <- TRUE
109 33
  auc1skeleton <- c(auc1skeleton, smoothing.args$roc1)
110 33
  names(auc1skeleton)[which(names(auc1skeleton) == "n")] <-  "smooth.n"
111 33
  auc2skeleton <- attributes(roc2$auc)
112 33
  auc2skeleton$roc <- NULL
113 33
  auc2skeleton$direction <- roc2$direction
114 33
  auc2skeleton$class <- NULL
115 33
  auc2skeleton$fun.sesp <- roc2$fun.sesp
116 33
  auc2skeleton$allow.invalid.partial.auc.correct <- TRUE
117 33
  auc2skeleton <- c(auc2skeleton, smoothing.args$roc2)
118 33
  names(auc2skeleton)[which(names(auc2skeleton) == "n")] <-  "smooth.n"
119

120 33
  auc1skeleton$auc <- auc2skeleton$auc <- test == "boot"
121

122
  # Some attributes may be duplicated in AUC skeletons and will mess the boostrap later on when we do.call().
123
  # If this condition happen, it probably means we have a bug elsewhere.
124
  # Rather than making a complicated processing to remove the duplicates,
125
  # just throw an error and let us solve the bug when a user reports it.
126 33
  duplicated.auc1skeleton <- duplicated(names(auc1skeleton))
127 33
  duplicated.auc2skeleton <- duplicated(names(auc2skeleton))
128 33
  if (any(duplicated.auc1skeleton)) {
129 0
  	sessionInfo <- sessionInfo()
130 0
  	save(roc1, roc2, test, x, paired, boot.n, boot.stratified, smoothing.args, progress, parallel, sessionInfo, file="pROC_bug.RData")
131 0
  	stop(sprintf("pROC: duplicated argument(s) in AUC1 skeleton: \"%s\". Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", paste(names(auc1skeleton)[duplicated(names(auc1skeleton))], collapse=", "), utils:: packageDescription("pROC")$BugReports))
132
  	
133
  }
134 33
  if (any(duplicated.auc2skeleton)) {
135 0
  	sessionInfo <- sessionInfo()
136 0
  	save(roc1, roc2, test, x, paired, boot.n, boot.stratified, smoothing.args, progress, parallel, sessionInfo, file="pROC_bug.RData")
137 0
  	stop(sprintf("duplicated argument(s) in AUC2 skeleton: \"%s\". Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", paste(names(auc2skeleton)[duplicated(names(auc2skeleton))], collapse=", "), utils:: packageDescription("pROC")$BugReports))
138
  }
139

140 33
  if (boot.stratified) { # precompute sorted responses if stratified
141
    #response.roc1 <- factor(c(rep(roc1$levels[1], length(roc1$controls)), rep(roc1$levels[2], length(roc1$cases))), levels=roc1$levels)
142
    #response.roc2 <- factor(c(rep(roc2$levels[1], length(roc2$controls)), rep(roc2$levels[2], length(roc2$cases))), levels=roc2$levels)
143
    #auc1skeleton$response <- response.roc1
144
    #auc2skeleton$response <- response.roc2
145 33
    resampled.values <- laply(1:boot.n, stratified.bootstrap.test, roc1=roc1, roc2=roc2, test=test, x=x, paired=paired, auc1skeleton=auc1skeleton, auc2skeleton=auc2skeleton, .progress=progress, .parallel=parallel)
146
  }
147
  else {
148 33
    auc1skeleton$levels <- roc1$levels
149 33
    auc1skeleton$direction <- roc1$direction
150 33
    auc2skeleton$levels <- roc2$levels
151 33
    auc2skeleton$direction <- roc2$direction
152 33
    resampled.values <- laply(1:boot.n, nonstratified.bootstrap.test, roc1=roc1, roc2=roc2, test=test, x=x, paired=paired, auc1skeleton=auc1skeleton, auc2skeleton=auc2skeleton, .progress=progress, .parallel=parallel)
153
  }
154

155
  # compute the statistics
156 33
  diffs <- resampled.values[,1] - resampled.values[,2]
157

158
  # are there NA values?
159 33
  if ((num.NAs <- sum(is.na(diffs))) > 0) {
160 0
    warning(sprintf("%i NA value(s) produced during bootstrap were ignored.", num.NAs))
161 0
    diffs <- diffs[!is.na(diffs)]
162
  }
163

164
  # Restore smoothing if necessary
165 33
  if (smoothing.args$roc1$smooth) {
166 33
    smoothing.args$roc1$method <- smoothing.args$roc1$smooth.method
167 33
    roc1 <- do.call("smooth.roc", c(list(roc=roc1), smoothing.args$roc1))
168
  }
169 33
  if (smoothing.args$roc2$smooth) {
170 33
    smoothing.args$roc2$method <- smoothing.args$roc2$smooth.method
171 33
    roc2 <- do.call("smooth.roc", c(list(roc=roc2), smoothing.args$roc2))
172
  }
173

174 33
  if (test == "sp") {
175 33
    coord1 <- coords(roc1, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
176 33
    coord2 <- coords(roc2, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
177 33
    D <- (coord1 - coord2) / sd(diffs)
178
  }
179 33
  else if (test == "se") {
180 33
    coord1 <- coords(roc1, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
181 33
    coord2 <- coords(roc2, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
182 33
    D <- (coord1 - coord2) / sd(diffs)
183
  }
184
  else {
185 33
    D <- (roc1$auc - roc2$auc) / sd(diffs)
186
  }
187 33
  if (is.nan(D) && all(diffs == 0) && roc1$auc == roc2$auc)
188 0
    D <- 0 # special case: no difference between AUCs produces a NaN
189

190 33
  return(D)
191
}
192

193
stratified.bootstrap.test <- function(n, roc1, roc2, test, x, paired, auc1skeleton, auc2skeleton) {
194
  # sample control and cases separately for a stratified bootstrap
195 33
  idx.controls.roc1 <- sample(1:length(roc1$controls), replace=TRUE)
196 33
  idx.cases.roc1 <- sample(1:length(roc1$cases), replace=TRUE)
197
  # finish roc skeletons
198 33
  auc1skeleton$controls <- roc1$controls[idx.controls.roc1]
199 33
  auc1skeleton$cases <- roc1$cases[idx.cases.roc1]
200

201 33
  if (paired) {
202 33
    auc2skeleton$controls <- roc2$controls[idx.controls.roc1]
203 33
    auc2skeleton$cases <- roc2$cases[idx.cases.roc1]
204
  }
205
  else { # for unpaired, resample roc2 separately
206 33
    idx.controls.roc2 <- sample(1:length(roc2$controls), replace=TRUE)
207 33
    idx.cases.roc2 <- sample(1:length(roc2$cases), replace=TRUE)
208 33
    auc2skeleton$controls <- roc2$controls[idx.controls.roc2]
209 33
    auc2skeleton$cases <- roc2$cases[idx.cases.roc2]
210
  }
211

212
  # re-compute the resampled ROC curves
213 33
  roc1 <- try(do.call("roc.cc.nochecks", auc1skeleton), silent=TRUE)
214 33
  roc2 <- try(do.call("roc.cc.nochecks", auc2skeleton), silent=TRUE)
215

216
  # resampled ROCs might not be smoothable: return NA
217 33
  if (methods::is(roc1, "try-error") || methods::is(roc2, "try-error")) {
218 0
    return(c(NA, NA))
219
  }
220
  else {
221 33
    if (test == "sp") {
222 33
      coord1 <- coords(roc1, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
223 33
      coord2 <- coords(roc2, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
224 33
      return(c(coord1, coord2))
225
    }
226 33
    else if (test == "se") {
227 33
      coord1 <- coords(roc1, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
228 33
      coord2 <- coords(roc2, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
229 33
      return(c(coord1, coord2))
230
    }
231
    else {
232 33
      return(c(roc1$auc, roc2$auc))
233
    }
234
  }
235
}
236

237
nonstratified.bootstrap.test <- function(n, roc1, roc2, test, x, paired, auc1skeleton, auc2skeleton) {
238
  # sample all patients
239 33
  idx.all.roc1 <- sample(1:length(roc1$response), replace=TRUE)
240
  # finish roc skeletons
241 33
  auc1skeleton$response <- roc1$response[idx.all.roc1]
242 33
  auc1skeleton$predictor <- roc1$predictor[idx.all.roc1]
243 33
  if (paired) { # if paired, resample roc2 as roc1
244 33
    auc2skeleton$response <- roc2$response[idx.all.roc1]
245 33
    auc2skeleton$predictor <- roc2$predictor[idx.all.roc1]
246
  }
247
  else { # if unpaired, resample roc2 separately
248 33
    idx.all.roc2 <- sample(1:length(roc2$response), replace=TRUE)
249 33
    auc2skeleton$response <- roc2$response[idx.all.roc2]
250 33
    auc2skeleton$predictor <- roc2$predictor[idx.all.roc2]
251
  }
252

253
  # re-compute the resampled ROC curves
254 33
  roc1 <- try(do.call("roc.rp.nochecks", auc1skeleton), silent=TRUE)
255 33
  roc2 <- try(do.call("roc.rp.nochecks", auc2skeleton), silent=TRUE)
256
  # resampled ROCs might not be smoothable: return NA
257 33
  if (methods::is(roc1, "try-error") || methods::is(roc2, "try-error")) {
258 0
    return(c(NA, NA))
259
  }
260
  else {
261 33
    if (test == "sp") {
262 33
      coord1 <- coords(roc1, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
263 33
      coord2 <- coords(roc2, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
264 33
      return(c(coord1, coord2))
265
    }
266 33
    else if (test == "se") {
267 33
      coord1 <- coords(roc1, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
268 33
      coord2 <- coords(roc2, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
269 33
      return(c(coord1, coord2))
270
    }
271
    else {
272 33
      return(c(roc1$auc, roc2$auc))
273
    }
274
  }
275
}
276

277
##########  AUC of one ROC curves (ci.auc, var)  ##########
278

279
ci.auc.bootstrap <- function(roc, conf.level, boot.n, boot.stratified, progress, parallel, ...) {
280 33
  if(class(progress) != "list")
281 33
    progress <- roc.utils.get.progress.bar(progress, title="AUC confidence interval", label="Bootstrap in progress...", ...)
282

283 33
  if (boot.stratified) {
284 33
    aucs <- unlist(llply(1:boot.n, .fun=stratified.ci.auc, roc=roc, .progress=progress, .parallel=parallel))
285
  }
286
  else {
287 33
    aucs <- unlist(llply(1:boot.n, .fun=nonstratified.ci.auc, roc=roc, .progress=progress, .parallel=parallel))
288
  }
289

290 33
  if (sum(is.na(aucs)) > 0) {
291 0
    warning("NA value(s) produced during bootstrap were ignored.")
292 0
    aucs <- aucs[!is.na(aucs)]
293
  }
294
  # TODO: Maybe apply a correction (it's in the Tibshirani?) What do Carpenter-Bithell say about that?
295
  # Prepare the return value
296 33
  return(quantile(aucs, c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
297
}
298

299
stratified.ci.auc <- function(n, roc) {
300 33
  controls <- sample(roc$controls, replace=TRUE)
301 33
  cases <- sample(roc$cases, replace=TRUE)
302 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
303
  
304 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
305 33
  roc$sensitivities <- perfs$se
306 33
  roc$specificities <- perfs$sp
307

308 33
  auc.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"), allow.invalid.partial.auc.correct = TRUE)
309
}
310

311
nonstratified.ci.auc <- function(n, roc) {
312 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
313 33
  predictor <- roc$predictor[tmp.idx]
314 33
  response <- roc$response[tmp.idx]
315 33
  splitted <- split(predictor, response)
316 33
  controls <- splitted[[as.character(roc$levels[1])]]
317 33
  cases <- splitted[[as.character(roc$levels[2])]]
318 33
  thresholds <- roc.utils.thresholds(c(controls, cases), roc$direction)
319

320 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
321 33
  roc$sensitivities <- perfs$se
322 33
  roc$specificities <- perfs$sp
323
  
324 33
  auc.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"), allow.invalid.partial.auc.correct = TRUE)
325
}
326

327
##########  AUC of a smooth ROC curve (ci.smooth.auc)  ##########
328

329
# Returns a smoothed auc in a stratified manner
330
stratified.ci.smooth.auc <- function(n, roc, smooth.roc.call, auc.call) {
331 33
  controls <- sample(roc$controls, replace=TRUE)
332 33
  cases <- sample(roc$cases, replace=TRUE)
333
  # need to rebuild a ROC and smooth it
334 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
335
  
336 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
337

338
  # update ROC
339 33
  roc$sensitivities <- perfs$se
340 33
  roc$specificities <- perfs$sp
341 33
  roc$cases <- cases
342 33
  roc$controls <- controls
343 33
  roc$predictor <- c(controls, cases)
344 33
  roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
345 33
  roc$thresholds <- thresholds
346

347
  # call smooth.roc and auc.smooth.roc
348 33
  smooth.roc.call$roc <- roc
349 33
  auc.call$smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
350 33
  if (methods::is(auc.call$smooth.roc, "try-error")) {
351 0
    return(NA)
352
  }
353 33
  return(eval(auc.call))
354
}
355

356
# Returns a smoothed auc in a non stratified manner
357
nonstratified.ci.smooth.auc <- function(n, roc, smooth.roc.call, auc.call) {
358 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
359 33
  predictor <- roc$predictor[tmp.idx]
360 33
  response <- roc$response[tmp.idx]
361 33
  splitted <- split(predictor, response)
362 33
  controls <- splitted[[as.character(roc$levels[1])]]
363 33
  cases <- splitted[[as.character(roc$levels[2])]]
364 33
  thresholds <- roc.utils.thresholds(c(controls, cases), roc$direction)
365

366 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
367

368
  # update ROC
369 33
  roc$sensitivities <- perfs$se
370 33
  roc$specificities <- perfs$sp
371 33
  roc$cases <- cases
372 33
  roc$controls <- controls
373 33
  roc$predictor <- predictor
374 33
  roc$response <- response
375 33
  roc$thresholds <- thresholds
376

377
  # call smooth.roc and auc.smooth.roc
378 33
  smooth.roc.call$roc <- roc
379 33
  auc.call$smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
380 33
  if (methods::is(auc.call$smooth.roc, "try-error")) {
381 0
    return(NA)
382
  }
383 33
  return(eval(auc.call))
384
}
385

386
# ##########  AUC of a multiclass ROC (ci.multiclass.auc)  ##########
387
#   
388
# ci.multiclass.auc.bootstrap <- function(roc, conf.level, boot.n, boot.stratified, progress, parallel, ...) {
389
#   if(class(progress) != "list")
390
#     progress <- roc.utils.get.progress.bar(progress, title="Multi-class AUC confidence interval", label="Bootstrap in progress...", ...)
391
# 
392
#   if (boot.stratified) {
393
#     aucs <- unlist(llply(1:boot.n, stratified.ci.multiclass.auc, roc=roc, .progress=progress, .parallel=parallel))
394
#   }
395
#   else {
396
#     aucs <- unlist(llply(1:boot.n, nonstratified.ci.multiclass.auc, roc=roc, .progress=progress, .parallel=parallel))
397
#   }
398
# 
399
#   if (sum(is.na(aucs)) > 0) {
400
#     warning("NA value(s) produced during bootstrap were ignored.")
401
#     aucs <- aucs[!is.na(aucs)]
402
#   }
403
#   # TODO: Maybe apply a correction (it's in the Tibshirani?) What do Carpenter-Bithell say about that?
404
#   # Prepare the return value
405
#   return(quantile(aucs, c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
406
# 
407
# }
408
# 
409
# # Returns an auc in a stratified manner
410
# stratified.ci.multiclass.auc <- function(n, roc) {
411
#   controls <- sample(roc$controls, replace=TRUE)
412
#   cases <- sample(roc$cases, replace=TRUE)
413
#   thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
414
#   
415
#   perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
416
#   roc$sensitivities <- perfs$se
417
#   roc$specificities <- perfs$sp
418
# 
419
#   auc.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"))
420
# }
421
# 
422
# 
423
# # Returns an auc in a non stratified manner
424
# nonstratified.ci.multiclass.auc <- function(n, roc) {
425
#   tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
426
#   predictor <- roc$predictor[tmp.idx]
427
#   response <- roc$response[tmp.idx]
428
#   splitted <- split(predictor, response)
429
#   controls <- splitted[[as.character(roc$levels[1])]]
430
#   cases <- splitted[[as.character(roc$levels[2])]]
431
#   thresholds <- roc.utils.thresholds(c(controls, cases), roc$direction)
432
# 
433
#   perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
434
#   roc$sensitivities <- perfs$se
435
#   roc$specificities <- perfs$sp
436
#   
437
#   auc.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"))
438
# }
439

440
##########  SE of a ROC curve (ci.se)  ##########
441

442
stratified.ci.se <- function(n, roc, sp) {
443 33
  controls <- sample(roc$controls, replace=TRUE)
444 33
  cases <- sample(roc$cases, replace=TRUE)
445 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
446
  
447 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
448 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
449 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
450 33
  roc$thresholds <- thresholds
451

452 33
  return(coords.roc(roc, sp, input = "specificity", ret = "sensitivity", transpose = FALSE, as.matrix = TRUE)[,1])
453
}
454

455
nonstratified.ci.se <- function(n, roc, sp) {
456 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
457 33
  predictor <- roc$predictor[tmp.idx]
458 33
  response <- roc$response[tmp.idx]
459 33
  splitted <- split(predictor, response)
460 33
  controls <- splitted[[as.character(roc$levels[1])]]
461 33
  cases <- splitted[[as.character(roc$levels[2])]]
462 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
463
  
464 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
465 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
466 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
467 33
  roc$thresholds <- thresholds
468

469 33
  return(coords.roc(roc, sp, input = "specificity", ret = "sensitivity", transpose = FALSE, as.matrix = TRUE)[, 1])
470
}
471

472
##########  SE of a smooth ROC curve (ci.se)  ##########
473

474
stratified.ci.smooth.se <- function(n, roc, sp, smooth.roc.call) {
475 33
  controls <- sample(roc$controls, replace=TRUE)
476 33
  cases <- sample(roc$cases, replace=TRUE)
477 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
478
  
479 33
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
480

481
  # update ROC
482 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
483 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
484 33
  roc$cases <- cases
485 33
  roc$controls <- controls
486 33
  roc$predictor <- c(controls, cases)
487 33
  roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
488 33
  roc$thresholds <- thresholds
489

490
  # call smooth.roc and auc.smooth.roc
491 33
  smooth.roc.call$roc <- roc
492 33
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
493 33
  if (methods::is(smooth.roc, "try-error"))
494 0
    return(NA)
495 33
  return(coords.smooth.roc(smooth.roc, sp, input = "specificity", ret = "sensitivity", transpose = FALSE, as.matrix = TRUE)[, 1])
496
}
497

498
nonstratified.ci.smooth.se <- function(n, roc, sp, smooth.roc.call) {
499 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
500 33
  predictor <- roc$predictor[tmp.idx]
501 33
  response <- roc$response[tmp.idx]
502 33
  splitted <- split(predictor, response)
503 33
  controls <- splitted[[as.character(roc$levels[1])]]
504 33
  cases <- splitted[[as.character(roc$levels[2])]]
505 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
506
  
507 33
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
508

509
  # update ROC
510 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
511 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
512 33
  roc$cases <- cases
513 33
  roc$controls <- controls
514 33
  roc$predictor <- predictor
515 33
  roc$response <- response
516 33
  roc$thresholds <- thresholds
517

518
  # call smooth.roc and auc.smooth.roc
519 33
  smooth.roc.call$roc <- roc
520 33
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
521 33
  if (methods::is(smooth.roc, "try-error"))
522 0
    return(NA)
523 33
  return(coords.smooth.roc(smooth.roc, sp, input = "specificity", ret = "sensitivity", transpose = FALSE, as.matrix = TRUE)[, 1])
524
}
525

526
##########  SP of a ROC curve (ci.sp)  ##########
527

528
stratified.ci.sp <- function(n, roc, se) {
529 33
  controls <- sample(roc$controls, replace=TRUE)
530 33
  cases <- sample(roc$cases, replace=TRUE)
531 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
532
  
533 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
534 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
535 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
536 33
  roc$thresholds <- thresholds
537

538 33
  return(coords.roc(roc, se, input = "sensitivity", ret = "specificity", transpose = FALSE, as.matrix = TRUE)[, 1])
539
}
540

541
nonstratified.ci.sp <- function(n, roc, se) {
542 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
543 33
  predictor <- roc$predictor[tmp.idx]
544 33
  response <- roc$response[tmp.idx]
545 33
  splitted <- split(predictor, response)
546 33
  controls <- splitted[[as.character(roc$levels[1])]]
547 33
  cases <- splitted[[as.character(roc$levels[2])]]
548 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
549
  
550 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
551 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
552 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
553 33
  roc$thresholds <- thresholds
554

555 33
  return(coords.roc(roc, se, input = "sensitivity", ret = "specificity", transpose = FALSE, as.matrix = TRUE)[, 1])
556
}
557

558
##########  SP of a smooth ROC curve (ci.sp)  ##########
559

560
stratified.ci.smooth.sp <- function(n, roc, se, smooth.roc.call) {
561 33
  controls <- sample(roc$controls, replace=TRUE)
562 33
  cases <- sample(roc$cases, replace=TRUE)
563 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
564
  
565 33
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
566

567
  # update ROC
568 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
569 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
570 33
  roc$cases <- cases
571 33
  roc$controls <- controls
572 33
  roc$predictor <- c(controls, cases)
573 33
  roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
574 33
  roc$thresholds <- thresholds
575

576
  # call smooth.roc and auc.smooth.roc
577 33
  smooth.roc.call$roc <- roc
578 33
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
579 33
  if (methods::is(smooth.roc, "try-error"))
580 0
    return(NA)
581 33
  return(coords.smooth.roc(smooth.roc, se, input = "sensitivity", ret = "specificity", transpose = FALSE, as.matrix = TRUE)[, 1])
582
}
583

584
nonstratified.ci.smooth.sp <- function(n, roc, se, smooth.roc.call) {
585 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
586 33
  predictor <- roc$predictor[tmp.idx]
587 33
  response <- roc$response[tmp.idx]
588 33
  splitted <- split(predictor, response)
589 33
  controls <- splitted[[as.character(roc$levels[1])]]
590 33
  cases <- splitted[[as.character(roc$levels[2])]]
591 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
592
  
593 33
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
594

595
  # update ROC
596 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
597 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
598 33
  roc$cases <- cases
599 33
  roc$controls <- controls
600 33
  roc$predictor <- predictor
601 33
  roc$response <- response
602 33
  roc$thresholds <- thresholds
603

604
  # call smooth.roc and auc.smooth.roc
605 33
  smooth.roc.call$roc <- roc
606 33
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
607 33
  if (methods::is(smooth.roc, "try-error"))
608 0
    return(NA)
609 33
  return(coords.smooth.roc(smooth.roc, se, input = "sensitivity", ret = "specificity", transpose = FALSE, as.matrix = TRUE)[, 1])
610
}
611

612
##########  Threshold of a ROC curve (ci.thresholds)  ##########
613

614
stratified.ci.thresholds <- function(n, roc, thresholds) {
615 33
  controls <- sample(roc$controls, replace=TRUE)
616 33
  cases <- sample(roc$cases, replace=TRUE)
617
  
618 33
  return(sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction))
619
}
620

621
# Returns an auc in a non stratified manner
622
nonstratified.ci.thresholds <- function(n, roc, thresholds) {
623 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
624 33
  predictor <- roc$predictor[tmp.idx]
625 33
  response <- roc$response[tmp.idx]
626 33
  splitted <- split(predictor, response)
627 33
  controls <- splitted[[as.character(roc$levels[1])]]
628 33
  cases <- splitted[[as.character(roc$levels[2])]]
629

630 33
  return(sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction))
631
}
632

633

634
##########  Coords of one ROC curves (ci.coords)  ##########
635
stratified.ci.coords <- function(roc, x, input, ret, best.method, best.weights, best.policy) {
636 33
  controls <- sample(roc$controls, replace=TRUE)
637 33
  cases <- sample(roc$cases, replace=TRUE)
638 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
639
  
640 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
641

642
  # update ROC
643 33
  roc$sensitivities <- perfs$se
644 33
  roc$specificities <- perfs$sp
645 33
  roc$cases <- cases
646 33
  roc$controls <- controls
647 33
  roc$predictor <- c(controls, cases)
648 33
  roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
649 33
  roc$thresholds <- thresholds
650

651 33
  res <- coords.roc(roc, x = x, input = input, ret = ret, 
652 33
                    best.method = best.method, best.weights = best.weights,
653 33
                    drop = FALSE, transpose = FALSE, as.matrix = TRUE)
654
  # Return a random column with "best"
655 33
  if (length(x) == 1 && x == "best" && nrow(res) != 1) {
656 0
  	return(enforce.best.policy(res, best.policy))
657
  }
658
  else {
659 33
  	return(res)
660
  }
661
}
662

663
nonstratified.ci.coords <- function(roc, x, input, ret, best.method, best.weights, best.policy) {
664 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
665 33
  predictor <- roc$predictor[tmp.idx]
666 33
  response <- roc$response[tmp.idx]
667 33
  splitted <- split(predictor, response)
668 33
  controls <- splitted[[as.character(roc$levels[1])]]
669 33
  cases <- splitted[[as.character(roc$levels[2])]]
670 33
  thresholds <- roc.utils.thresholds(c(controls, cases), roc$direction)
671

672

673 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
674

675
  # update ROC
676 33
  roc$sensitivities <- perfs$se
677 33
  roc$specificities <- perfs$sp
678 33
  roc$cases <- cases
679 33
  roc$controls <- controls
680 33
  roc$predictor <- c(controls, cases)
681 33
  roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
682 33
  roc$thresholds <- thresholds
683
  
684 33
  res <- coords.roc(roc, x = x, input = input, ret = ret,
685 33
                    best.method = best.method, best.weights = best.weights,
686 33
                    drop = FALSE, transpose = FALSE, as.matrix = TRUE)
687
  # Return a random column with "best"
688 33
  if (length(x) == 1 && x == "best" && nrow(res) != 1) {
689 0
  	return(enforce.best.policy(res, best.policy))
690
  }
691
  else {
692 33
  	return(res)
693
  }
694
}
695

696
##########  Coords of a smooth ROC curve (ci.coords)  ##########
697

698
stratified.ci.smooth.coords <- function(roc, x, input, ret, best.method, best.weights, smooth.roc.call, best.policy) {
699 33
  controls <- sample(roc$controls, replace=TRUE)
700 33
  cases <- sample(roc$cases, replace=TRUE)
701 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
702
  
703

704 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
705

706
  # update ROC
707 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
708 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
709 33
  roc$cases <- cases
710 33
  roc$controls <- controls
711 33
  roc$predictor <- c(controls, cases)
712 33
  roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
713 33
  roc$thresholds <- thresholds
714

715
  # call smooth.roc and auc.smooth.roc
716 33
  smooth.roc.call$roc <- roc
717 33
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
718 33
  if (methods::is(smooth.roc, "try-error"))
719 0
    return(NA)
720 33
  res <- coords.roc(smooth.roc, x = x, input = input, ret = ret,
721 33
                    best.method = best.method, best.weights = best.weights,
722 33
                    drop = FALSE, transpose = FALSE, as.matrix = TRUE)
723
  # Return a random column with "best"
724 33
  if (length(x) == 1 && x == "best" && nrow(res) != 1) {
725 0
  	return(enforce.best.policy(res, best.policy))
726
  }
727
  else {
728 33
  	return(res)
729
  }
730
}
731

732
nonstratified.ci.smooth.coords <- function(roc, x, input, ret, best.method, best.weights, smooth.roc.call, best.policy) {
733 33
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
734 33
  predictor <- roc$predictor[tmp.idx]
735 33
  response <- roc$response[tmp.idx]
736 33
  splitted <- split(predictor, response)
737 33
  controls <- splitted[[as.character(roc$levels[1])]]
738 33
  cases <- splitted[[as.character(roc$levels[2])]]
739 33
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
740
  
741 33
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
742

743
  # update ROC
744 33
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
745 33
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
746 33
  roc$cases <- cases
747 33
  roc$controls <- controls
748 33
  roc$predictor <- predictor
749 33
  roc$response <- response
750 33
  roc$thresholds <- thresholds
751

752
  # call smooth.roc and auc.smooth.roc
753 33
  smooth.roc.call$roc <- roc
754 33
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
755 33
  if (methods::is(smooth.roc, "try-error"))
756 0
    return(NA)
757 33
  res <- coords.roc(smooth.roc, x = x, input = input, ret = ret,
758 33
                    best.method = best.method, best.weights = best.weights,
759 33
                    drop = FALSE, transpose = FALSE, as.matrix = TRUE)
760
  # Return a random column with "best"
761 33
  if (length(x) == 1 && x == "best" && nrow(res) != 1) {
762 0
  	return(enforce.best.policy(res, best.policy))
763
  }
764
  else {
765 33
  	return(res)
766
  }
767
}

Read our documentation on viewing source code .

Loading