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 33
var <- function(...)
21 33
  UseMethod("var")
22

23
var.default <- function(...) {
24 33
  stats::var(...)
25
}
26

27
var.auc <- function(auc, ...) {
28
  # Change roc from an auc to a roc object but keep the auc specifications
29 33
  roc <- auc
30 33
  attr(auc, "roc") <- NULL
31 33
  roc <- attr(roc, "roc")
32 33
  roc$auc <- auc
33
  # Pass to var.roc
34 33
  var.roc(roc, ...)
35
}
36

37
var.smooth.roc <- function(smooth.roc, ...) {
38 33
  var.roc(smooth.roc, ...) # just pass to var.roc that will do the job
39
}
40

41
var.roc <- function(roc,
42
                    method=c("delong", "bootstrap", "obuchowski"),
43
                    boot.n = 2000,
44
                    boot.stratified = TRUE,
45
                    reuse.auc=TRUE,
46
                    progress = getOption("pROCProgress")$name,
47
                    parallel = FALSE,
48
                    ...) {
49
  # We need an auc
50 33
  if (is.null(roc$auc) | !reuse.auc)
51 0
    roc$auc <- auc(roc, ...)
52
  
53 33
  if (roc.utils.is.perfect.curve(roc)) {
54 0
  	warning("var() of a ROC curve with AUC == 1 is always 0 and can be misleading.")
55
  }
56

57
  # do all the computations in fraction, re-transform in percent later
58 33
  percent <- roc$percent
59 33
  if (percent) {
60 33
  	roc <- roc.utils.unpercent(roc)
61
  }
62

63
  # Check the method
64 33
  if (missing(method) | is.null(method)) {
65
    # determine method if missing
66 33
    if (has.partial.auc(roc)) {
67
      # partial auc: go for bootstrap
68 0
      method <- "bootstrap"
69
    }
70 33
    else if (class(roc) == "smooth.roc") {
71
      # smoothing: bootstrap
72 0
      method <- "bootstrap"
73
    }
74
    else {
75 33
      method <- "delong"
76
    }
77
  }
78
  else {
79 33
    method <- match.arg(method)
80
    # delong NA to pAUC: warn + change
81 33
    if (method == "delong") {
82 33
      if (has.partial.auc(roc)) {
83 0
      	stop("DeLong method is not supported for partial AUC. Use method=\"bootstrap\" instead.")
84
      }
85 33
      else if ("smooth.roc" %in% class(roc)) {
86 0
      	stop("DeLong method is not supported for smoothed ROCs. Use method=\"bootstrap\" instead.")
87
      }
88
    }
89

90 33
    else if (method == "obuchowski") {
91 33
      if ("smooth.roc" %in% class(roc)) {
92 0
        stop("Using Obuchowski for smoothed ROCs is not supported. Using bootstrap instead.")
93
      }
94 33
      if (has.partial.auc(roc) && attr(roc$auc, "partial.auc.focus") == "sensitivity") {
95 0
        stop("Using Obuchowski for partial AUC on sensitivity region is not supported. Using bootstrap instead.")
96
      }
97
    }
98
  }
99
  
100 33
  if (method == "delong") {
101 33
    n <- length(roc$controls)
102 33
    m <- length(roc$cases)  
103 33
    V <- delongPlacements(roc)
104 33
    var <- var(V$Y) / n + var(V$X) / m
105
  }
106 33
  else if (method == "obuchowski") {
107 33
    var <- var.roc.obuchowski(roc) / length(roc$cases)
108
  }
109
  else {
110 33
    var <- var.roc.bootstrap(roc, boot.n, boot.stratified, progress, parallel, ...)
111
  }
112
  
113 33
  if (percent) {
114 33
    var <- var * (100^2)
115
  }
116 33
  return(var)
117
}
118

119
var.roc.bootstrap <- function(roc, boot.n, boot.stratified, progress, parallel, ...) {
120 33
  if(class(progress) != "list")
121 33
    progress <- roc.utils.get.progress.bar(progress, title="AUC variance", label="Bootstrap in progress...", ...)
122

123
  ## Smoothed ROC curve variance
124 33
  if (class(roc) == "smooth.roc") {
125 33
    smoothing.args <- roc$smoothing.args
126 33
    smoothing.args$smooth <- TRUE
127 33
    non.smoothed.roc <- attr(roc, "roc")
128 33
    non.smoothed.roc$percent <- FALSE # as we did earlier for the smoothed.roc
129 33
    smooth.roc.call <- as.call(c(utils::getS3method("smooth", "roc"), roc$smoothing.args))
130 33
    auc.args <- attributes(roc$auc)[grep("partial.auc", names(attributes(roc$auc)))]
131 33
    auc.args$allow.invalid.partial.auc.correct <- TRUE
132 33
    auc.call <- as.call(c(utils::getS3method("auc", "smooth.roc"), auc.args))
133

134 33
    if (boot.stratified) {
135 0
      aucs <- unlist(llply(1:boot.n, stratified.ci.smooth.auc, roc=non.smoothed.roc, smooth.roc.call=smooth.roc.call, auc.call=auc.call, .progress=progress, .parallel=parallel))
136
    }
137
    else {
138 33
      aucs <- unlist(llply(1:boot.n, nonstratified.ci.smooth.auc, roc=non.smoothed.roc, smooth.roc.call=smooth.roc.call, auc.call=auc.call, .progress=progress, .parallel=parallel))
139
    }
140
  }
141
  ## Non smoothed ROC curves variance
142
  else {
143 33
    if (boot.stratified) {
144 33
      aucs <- unlist(llply(1:boot.n, stratified.ci.auc, roc=roc, .progress=progress, .parallel=parallel)) # ci.auc: returns aucs just as we need for var, so re-use it!
145
    }
146
    else {
147 33
      aucs <- unlist(llply(1:boot.n, nonstratified.ci.auc, roc=roc, .progress=progress, .parallel=parallel))
148
    }
149
  }
150

151 33
  if ((num.NAs <- sum(is.na(aucs))) > 0) {
152 0
    warning(sprintf("%i NA value(s) produced during bootstrap were ignored.", num.NAs))
153 0
    aucs <- aucs[!is.na(aucs)]
154
  }
155 33
  return(var(aucs))
156
}

Read our documentation on viewing source code .

Loading