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 1
  smoothing.args$roc1$smooth.method <- smoothing.args$roc1$method
25 1
  smoothing.args$roc1$method <- NULL
26 1
  smoothing.args$roc2$smooth.method <- smoothing.args$roc2$method
27 1
  smoothing.args$roc2$method <- NULL
28

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

47 1
  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 1
  duplicated.auc1skeleton <- duplicated(names(auc1skeleton))
54 1
  duplicated.auc2skeleton <- duplicated(names(auc2skeleton))
55 1
  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 1
  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 1
  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 1
    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 1
    auc1skeleton$levels <- roc1$levels
75 1
    auc1skeleton$direction <- roc1$direction
76 1
    auc2skeleton$levels <- roc2$levels
77 1
    auc2skeleton$direction <- roc2$direction
78 1
    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 1
  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 1
  cov <- stats::cov(resampled.values[,1], resampled.values[,2])
88 1
  if (boot.return) {
89 0
    attr(cov, "resampled.values") <- resampled.values
90
  }
91 1
  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 1
  smoothing.args$roc1$smooth.method <- smoothing.args$roc1$method
98 1
  smoothing.args$roc1$method <- NULL
99 1
  smoothing.args$roc2$smooth.method <- smoothing.args$roc2$method
100 1
  smoothing.args$roc2$method <- NULL
101

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

120 1
  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 1
  duplicated.auc1skeleton <- duplicated(names(auc1skeleton))
127 1
  duplicated.auc2skeleton <- duplicated(names(auc2skeleton))
128 1
  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 1
  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 1
  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 1
    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 1
    auc1skeleton$levels <- roc1$levels
149 1
    auc1skeleton$direction <- roc1$direction
150 1
    auc2skeleton$levels <- roc2$levels
151 1
    auc2skeleton$direction <- roc2$direction
152 1
    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 1
  diffs <- resampled.values[,1] - resampled.values[,2]
157

158
  # are there NA values?
159 1
  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 1
  if (smoothing.args$roc1$smooth) {
166 1
    smoothing.args$roc1$method <- smoothing.args$roc1$smooth.method
167 1
    roc1 <- do.call("smooth.roc", c(list(roc=roc1), smoothing.args$roc1))
168
  }
169 1
  if (smoothing.args$roc2$smooth) {
170 1
    smoothing.args$roc2$method <- smoothing.args$roc2$smooth.method
171 1
    roc2 <- do.call("smooth.roc", c(list(roc=roc2), smoothing.args$roc2))
172
  }
173

174 1
  if (test == "sp") {
175 1
    coord1 <- coords(roc1, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
176 1
    coord2 <- coords(roc2, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
177 1
    D <- (coord1 - coord2) / sd(diffs)
178
  }
179 1
  else if (test == "se") {
180 1
    coord1 <- coords(roc1, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
181 1
    coord2 <- coords(roc2, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
182 1
    D <- (coord1 - coord2) / sd(diffs)
183
  }
184
  else {
185 1
    D <- (roc1$auc - roc2$auc) / sd(diffs)
186
  }
187 1
  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 1
  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 1
  idx.controls.roc1 <- sample(1:length(roc1$controls), replace=TRUE)
196 1
  idx.cases.roc1 <- sample(1:length(roc1$cases), replace=TRUE)
197
  # finish roc skeletons
198 1
  auc1skeleton$controls <- roc1$controls[idx.controls.roc1]
199 1
  auc1skeleton$cases <- roc1$cases[idx.cases.roc1]
200

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

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

216
  # resampled ROCs might not be smoothable: return NA
217 1
  if (methods::is(roc1, "try-error") || methods::is(roc2, "try-error")) {
218 0
    return(c(NA, NA))
219
  }
220
  else {
221 1
    if (test == "sp") {
222 1
      coord1 <- coords(roc1, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
223 1
      coord2 <- coords(roc2, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
224 1
      return(c(coord1, coord2))
225
    }
226 1
    else if (test == "se") {
227 1
      coord1 <- coords(roc1, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
228 1
      coord2 <- coords(roc2, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
229 1
      return(c(coord1, coord2))
230
    }
231
    else {
232 1
      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 1
  idx.all.roc1 <- sample(1:length(roc1$response), replace=TRUE)
240
  # finish roc skeletons
241 1
  auc1skeleton$response <- roc1$response[idx.all.roc1]
242 1
  auc1skeleton$predictor <- roc1$predictor[idx.all.roc1]
243 1
  if (paired) { # if paired, resample roc2 as roc1
244 1
    auc2skeleton$response <- roc2$response[idx.all.roc1]
245 1
    auc2skeleton$predictor <- roc2$predictor[idx.all.roc1]
246
  }
247
  else { # if unpaired, resample roc2 separately
248 1
    idx.all.roc2 <- sample(1:length(roc2$response), replace=TRUE)
249 1
    auc2skeleton$response <- roc2$response[idx.all.roc2]
250 1
    auc2skeleton$predictor <- roc2$predictor[idx.all.roc2]
251
  }
252

253
  # re-compute the resampled ROC curves
254 1
  roc1 <- try(do.call("roc.rp.nochecks", auc1skeleton), silent=TRUE)
255 1
  roc2 <- try(do.call("roc.rp.nochecks", auc2skeleton), silent=TRUE)
256
  # resampled ROCs might not be smoothable: return NA
257 1
  if (methods::is(roc1, "try-error") || methods::is(roc2, "try-error")) {
258 0
    return(c(NA, NA))
259
  }
260
  else {
261 1
    if (test == "sp") {
262 1
      coord1 <- coords(roc1, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
263 1
      coord2 <- coords(roc2, x=x, input=c("specificity"), ret=c("sensitivity"), as.matrix=TRUE, transpose=FALSE)[1]
264 1
      return(c(coord1, coord2))
265
    }
266 1
    else if (test == "se") {
267 1
      coord1 <- coords(roc1, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
268 1
      coord2 <- coords(roc2, x=x, input=c("sensitivity"), ret=c("specificity"), as.matrix=TRUE, transpose=FALSE)[1]
269 1
      return(c(coord1, coord2))
270
    }
271
    else {
272 1
      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 1
  if(class(progress) != "list")
281 1
    progress <- roc.utils.get.progress.bar(progress, title="AUC confidence interval", label="Bootstrap in progress...", ...)
282

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

290 1
  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 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
301 1
  cases <- sample(roc$cases, replace=TRUE)
302 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
303
  
304 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
305 1
  roc$sensitivities <- perfs$se
306 1
  roc$specificities <- perfs$sp
307

308 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
313 1
  predictor <- roc$predictor[tmp.idx]
314 1
  response <- roc$response[tmp.idx]
315 1
  splitted <- split(predictor, response)
316 1
  controls <- splitted[[as.character(roc$levels[1])]]
317 1
  cases <- splitted[[as.character(roc$levels[2])]]
318 1
  thresholds <- roc.utils.thresholds(c(controls, cases), roc$direction)
319

320 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
321 1
  roc$sensitivities <- perfs$se
322 1
  roc$specificities <- perfs$sp
323
  
324 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
332 1
  cases <- sample(roc$cases, replace=TRUE)
333
  # need to rebuild a ROC and smooth it
334 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
335
  
336 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
337

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

347
  # call smooth.roc and auc.smooth.roc
348 1
  smooth.roc.call$roc <- roc
349 1
  auc.call$smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
350 1
  if (methods::is(auc.call$smooth.roc, "try-error")) {
351 0
    return(NA)
352
  }
353 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
359 1
  predictor <- roc$predictor[tmp.idx]
360 1
  response <- roc$response[tmp.idx]
361 1
  splitted <- split(predictor, response)
362 1
  controls <- splitted[[as.character(roc$levels[1])]]
363 1
  cases <- splitted[[as.character(roc$levels[2])]]
364 1
  thresholds <- roc.utils.thresholds(c(controls, cases), roc$direction)
365

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

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

377
  # call smooth.roc and auc.smooth.roc
378 1
  smooth.roc.call$roc <- roc
379 1
  auc.call$smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
380 1
  if (methods::is(auc.call$smooth.roc, "try-error")) {
381 0
    return(NA)
382
  }
383 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
444 1
  cases <- sample(roc$cases, replace=TRUE)
445 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
446
  
447 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
448 1
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
449 1
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
450 1
  roc$thresholds <- thresholds
451

452 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
457 1
  predictor <- roc$predictor[tmp.idx]
458 1
  response <- roc$response[tmp.idx]
459 1
  splitted <- split(predictor, response)
460 1
  controls <- splitted[[as.character(roc$levels[1])]]
461 1
  cases <- splitted[[as.character(roc$levels[2])]]
462 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
463
  
464 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
465 1
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
466 1
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
467 1
  roc$thresholds <- thresholds
468

469 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
476 1
  cases <- sample(roc$cases, replace=TRUE)
477 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
478
  
479 1
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
480

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

490
  # call smooth.roc and auc.smooth.roc
491 1
  smooth.roc.call$roc <- roc
492 1
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
493 1
  if (methods::is(smooth.roc, "try-error"))
494 0
    return(NA)
495 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
500 1
  predictor <- roc$predictor[tmp.idx]
501 1
  response <- roc$response[tmp.idx]
502 1
  splitted <- split(predictor, response)
503 1
  controls <- splitted[[as.character(roc$levels[1])]]
504 1
  cases <- splitted[[as.character(roc$levels[2])]]
505 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
506
  
507 1
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
508

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

518
  # call smooth.roc and auc.smooth.roc
519 1
  smooth.roc.call$roc <- roc
520 1
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
521 1
  if (methods::is(smooth.roc, "try-error"))
522 0
    return(NA)
523 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
530 1
  cases <- sample(roc$cases, replace=TRUE)
531 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
532
  
533 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
534 1
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
535 1
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
536 1
  roc$thresholds <- thresholds
537

538 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
543 1
  predictor <- roc$predictor[tmp.idx]
544 1
  response <- roc$response[tmp.idx]
545 1
  splitted <- split(predictor, response)
546 1
  controls <- splitted[[as.character(roc$levels[1])]]
547 1
  cases <- splitted[[as.character(roc$levels[2])]]
548 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
549
  
550 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
551 1
  roc$sensitivities <- perfs$se * ifelse(roc$percent, 100, 1)
552 1
  roc$specificities <- perfs$sp * ifelse(roc$percent, 100, 1)
553 1
  roc$thresholds <- thresholds
554

555 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
562 1
  cases <- sample(roc$cases, replace=TRUE)
563 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
564
  
565 1
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
566

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

576
  # call smooth.roc and auc.smooth.roc
577 1
  smooth.roc.call$roc <- roc
578 1
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
579 1
  if (methods::is(smooth.roc, "try-error"))
580 0
    return(NA)
581 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
586 1
  predictor <- roc$predictor[tmp.idx]
587 1
  response <- roc$response[tmp.idx]
588 1
  splitted <- split(predictor, response)
589 1
  controls <- splitted[[as.character(roc$levels[1])]]
590 1
  cases <- splitted[[as.character(roc$levels[2])]]
591 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
592
  
593 1
    perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
594

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

604
  # call smooth.roc and auc.smooth.roc
605 1
  smooth.roc.call$roc <- roc
606 1
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
607 1
  if (methods::is(smooth.roc, "try-error"))
608 0
    return(NA)
609 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
616 1
  cases <- sample(roc$cases, replace=TRUE)
617
  
618 1
  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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
624 1
  predictor <- roc$predictor[tmp.idx]
625 1
  response <- roc$response[tmp.idx]
626 1
  splitted <- split(predictor, response)
627 1
  controls <- splitted[[as.character(roc$levels[1])]]
628 1
  cases <- splitted[[as.character(roc$levels[2])]]
629

630 1
  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 1
  controls <- sample(roc$controls, replace=TRUE)
637 1
  cases <- sample(roc$cases, replace=TRUE)
638 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
639
  
640 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
641

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

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

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

672

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

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

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

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

715
  # call smooth.roc and auc.smooth.roc
716 1
  smooth.roc.call$roc <- roc
717 1
  smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
718 1
  if (methods::is(smooth.roc, "try-error"))
719 0
    return(NA)
720 1
  res <- coords.roc(smooth.roc, x = x, input = input, ret = ret,
721 1
                    best.method = best.method, best.weights = best.weights,
722 1
                    drop = FALSE, transpose = FALSE, as.matrix = TRUE)
723
  # Return a random column with "best"
724 1
  if (length(x) == 1 && x == "best" && nrow(res) != 1) {
725 0
  	return(enforce.best.policy(res, best.policy))
726
  }
727
  else {
728 1
  	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 1
  tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
734 1
  predictor <- roc$predictor[tmp.idx]
735 1
  response <- roc$response[tmp.idx]
736 1
  splitted <- split(predictor, response)
737 1
  controls <- splitted[[as.character(roc$levels[1])]]
738 1
  cases <- splitted[[as.character(roc$levels[2])]]
739 1
  thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
740
  
741 1
  perfs <- roc$fun.sesp(thresholds=thresholds, controls=controls, cases=cases, direction=roc$direction)
742

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

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

Read our documentation on viewing source code .

Loading