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

24
cov.default <- function(...) {
25 33
  stats::cov(...)
26
}
27

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

38
cov.smooth.roc <- function(roc1, roc2, ...) {
39 33
  cov.roc(roc1, roc2, ...)
40
}
41

42
cov.roc <- function(roc1, roc2,
43
                         method=c("delong", "bootstrap", "obuchowski"),
44
                         reuse.auc=TRUE,
45
                         boot.n=2000, boot.stratified=TRUE, boot.return=FALSE,
46
                         progress=getOption("pROCProgress")$name,
47
                         parallel = FALSE,
48
                         ...) {
49
  # If roc2 is an auc, take the roc but keep the auc specifications
50 33
  if (methods::is(roc2, "auc")) {
51 33
    auc2 <- roc2
52 33
    attr(auc2, "roc") <- NULL
53 33
    roc2 <- attr(roc2, "roc")
54 33
    roc2$auc <- auc2
55
  }
56
  
57 33
  if (roc.utils.is.perfect.curve(roc1) && roc.utils.is.perfect.curve(roc2)) {
58 0
  	warning("cov() of two ROC curves with AUC == 1 is always 0 and can be misleading.")
59
  }
60

61
  # store which objects are smoothed, and how
62 33
  smoothing.args <- list()
63 33
  if ("smooth.roc" %in% class(roc1)) {
64 33
    smoothing.args$roc1 <- roc1$smoothing.args
65 33
    smoothing.args$roc1$smooth <- TRUE
66 33
    roc1 <- attr(roc1, "roc")
67
    #oroc1$auc <- roc1$auc
68
  }
69
  else {
70 33
    smoothing.args$roc1 <- list(smooth=FALSE)
71
  }
72 33
  if ("smooth.roc" %in% class(roc2)) {
73 33
    smoothing.args$roc2 <- roc2$smoothing.args
74 33
    smoothing.args$roc2$smooth <- TRUE
75 33
    roc2 <- attr(roc2, "roc")
76
    #oroc2$auc <- roc2$auc
77
  }
78
  else {
79 33
    smoothing.args$roc2 <- list(smooth=FALSE)
80
  }
81

82
  # then determine whether the rocs are paired or not
83 33
  rocs.are.paired <- are.paired(roc1, roc2, return.paired.rocs=FALSE, reuse.auc=TRUE, reuse.ci=FALSE, reuse.smooth=TRUE)
84 33
  if (! rocs.are.paired) {
85 0
    message("ROC curves are unpaired.")
86 0
    return(0)
87
  }    
88

89
  # check that the AUC was computed, or do it now
90 33
  if (is.null(roc1$auc) | !reuse.auc) {
91 33
    if (smoothing.args$roc1$smooth) {
92 33
      roc1$auc <- auc(smooth.roc=do.call("smooth.roc", c(list(roc=roc1), smoothing.args$roc1)), ...)
93
      # remove partial.auc.* arguments that are now in roc1$auc and that will mess later processing
94
      # (formal argument "partial.auc(.*)" matched by multiple actual arguments)
95
      # This removal should be safe because we always use smoothing.args with roc1 in the following processing,
96
      # however it is a potential source of bugs.
97 33
      smoothing.args$roc1$partial.auc <- NULL
98 33
      smoothing.args$roc1$partial.auc.correct <- NULL
99 33
      smoothing.args$roc1$partial.auc.focus <- NULL
100
    }
101
    else
102 0
      roc1$auc <- auc(roc1, ...)
103
  }
104 33
  if (is.null(roc2$auc) | !reuse.auc) {
105 33
    if (smoothing.args$roc2$smooth) {
106 33
      roc2$auc <- auc(smooth.roc=do.call("smooth.roc", c(list(roc=roc2), smoothing.args$roc2)), ...)
107
      # remove partial.auc.* arguments that are now in roc1$auc and that will mess later processing
108
      # (formal argument "partial.auc(.*)" matched by multiple actual arguments)
109
      # This removal should be safe because we always use smoothing.args with roc2 in the following processing,
110
      # however it is a potential source of bugs.
111 33
      smoothing.args$roc2$partial.auc <- NULL
112 33
      smoothing.args$roc2$partial.auc.correct <- NULL
113 33
      smoothing.args$roc2$partial.auc.focus <- NULL
114
    }
115
    else
116 0
      roc2$auc <- auc(roc2, ...)
117
  }
118
    
119
  # check that the same region was requested in auc. Otherwise, issue a warning
120 33
  if (!identical(attributes(roc1$auc)[names(attributes(roc1$auc))!="roc"], attributes(roc2$auc)[names(attributes(roc2$auc))!="roc"]))
121 33
    warning("Different AUC specifications in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
122
  # check that the same smoothing params were requested in auc. Otherwise, issue a warning
123 33
  if (!identical(smoothing.args$roc1, smoothing.args$roc2))
124 33
    warning("Different smoothing parameters in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
125

126
  # Check the method
127 33
  if (missing(method) | is.null(method)) {
128
    # determine method if missing
129 33
    if (has.partial.auc(roc1)) {
130
      # partial auc: go for bootstrap
131 33
      method <- "bootstrap"
132
    }
133 33
    else if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth) {
134
      # smoothing in one or both: bootstrap
135 0
      method <- "bootstrap"
136
    }
137 33
    else if (roc1$direction != roc2$direction) {
138
      # delong doesn't work well with opposite directions (will report high significance if roc1$auc and roc2$auc are similar and high)
139 0
      method <- "bootstrap"
140
    }
141
    else {
142 33
      method <- "delong"
143
    }
144
  }
145
  else {
146 33
    method <- match.arg(method)
147 33
    if (method == "delong") {
148
      # delong NA to pAUC: warn + change
149 33
      if (has.partial.auc(roc1) || has.partial.auc(roc2)) {
150 0
      	stop("DeLong method is not supported for partial AUC. Use method=\"bootstrap\" instead.")
151
      }
152 33
      if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth) {
153 0
      	stop("DeLong method is not supported for smoothed ROCs. Use method=\"bootstrap\" instead.")
154
      }
155 33
      if (roc1$direction != roc2$direction)
156 0
        warning("DeLong method should not be applied to ROC curves with a different direction.")
157
    }
158 33
    else if (method == "obuchowski") {
159 33
      if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth) {
160 0
        stop("Obuchowski method is not supported for smoothed ROCs. Use method=\"bootstrap\" instead.")
161
      }
162 33
      if ((has.partial.auc(roc1) && attr(roc1$auc, "partial.auc.focus") == "sensitivity")
163 33
          || (has.partial.auc(roc2) && attr(roc2$auc, "partial.auc.focus") == "sensitivity")) {
164 0
        stop("Obuchowski method is not supported for partial AUC on sensitivity region. Use method=\"bootstrap\" instead.")
165
      }
166 33
      if (roc1$direction != roc2$direction)
167 0
        warning("Obuchowski method should not be applied to ROC curves with a different direction.")
168
    }
169
  }
170
  
171 33
  if (method == "delong") {
172 33
    n <- length(roc1$controls)
173 33
    m <- length(roc1$cases)
174

175 33
    V1 <- delongPlacements(roc1)
176 33
    var1 <- var(V1$Y) / n + var(V1$X) / m
177

178 33
    V2 <- delongPlacements(roc2)
179 33
    var2 <- var(V2$Y) / n + var(V2$X) / m
180

181 33
    cov <- cov(V2$X, V1$X) / m + cov(V2$Y, V1$Y) / n
182

183 33
    if (roc1$percent) {
184 33
      cov <- cov * (100^2)
185
    }
186
  }
187
  
188 33
  else if (method == "obuchowski") {
189 33
    cov <- cov.roc.obuchowski(roc1, roc2) / length(roc1$cases)
190

191 33
    if (roc1$percent) {
192 33
      cov <- cov * (100^2)
193
    }
194
  }
195
  else { # method == "bootstrap"
196
    # Check if called with density.cases or density.controls
197 33
    if (is.null(smoothing.args) || is.numeric(smoothing.args$density.cases) || is.numeric(smoothing.args$density.controls))
198 0
      stop("Cannot compute the covariance of ROC curves smoothed with numeric density.controls and density.cases.")
199

200 33
    if(class(progress) != "list") {
201 33
      progress <- roc.utils.get.progress.bar(progress, title="Bootstrap covariance", label="Bootstrap in progress...", ...)
202
    }
203
    
204 33
    cov <- bootstrap.cov(roc1, roc2, boot.n, boot.stratified, boot.return, smoothing.args, progress, parallel)
205
  }
206

207 33
  return(cov)
208
}

Read our documentation on viewing source code .

Loading