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 1
  UseMethod("cov")
22
}
23

24
cov.default <- function(...) {
25 1
  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 1
  auc1 <- roc1
31 1
  attr(auc1, "roc") <- NULL
32 1
  roc1 <- attr(roc1, "roc")
33 1
  roc1$auc <- auc1
34
  # Pass to cov.roc
35 1
  return(cov.roc(roc1, roc2, ...))
36
}
37

38
cov.smooth.roc <- function(roc1, roc2, ...) {
39 1
  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 1
  if (methods::is(roc2, "auc")) {
51 1
    auc2 <- roc2
52 1
    attr(auc2, "roc") <- NULL
53 1
    roc2 <- attr(roc2, "roc")
54 1
    roc2$auc <- auc2
55
  }
56
  
57 1
  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 1
  smoothing.args <- list()
63 1
  if ("smooth.roc" %in% class(roc1)) {
64 1
    smoothing.args$roc1 <- roc1$smoothing.args
65 1
    smoothing.args$roc1$smooth <- TRUE
66 1
    roc1 <- attr(roc1, "roc")
67
    #oroc1$auc <- roc1$auc
68
  }
69
  else {
70 1
    smoothing.args$roc1 <- list(smooth=FALSE)
71
  }
72 1
  if ("smooth.roc" %in% class(roc2)) {
73 1
    smoothing.args$roc2 <- roc2$smoothing.args
74 1
    smoothing.args$roc2$smooth <- TRUE
75 1
    roc2 <- attr(roc2, "roc")
76
    #oroc2$auc <- roc2$auc
77
  }
78
  else {
79 1
    smoothing.args$roc2 <- list(smooth=FALSE)
80
  }
81

82
  # then determine whether the rocs are paired or not
83 1
  rocs.are.paired <- are.paired(roc1, roc2, return.paired.rocs=FALSE, reuse.auc=TRUE, reuse.ci=FALSE, reuse.smooth=TRUE)
84 1
  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 1
  if (is.null(roc1$auc) | !reuse.auc) {
91 1
    if (smoothing.args$roc1$smooth) {
92 1
      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 1
      smoothing.args$roc1$partial.auc <- NULL
98 1
      smoothing.args$roc1$partial.auc.correct <- NULL
99 1
      smoothing.args$roc1$partial.auc.focus <- NULL
100
    }
101
    else
102 0
      roc1$auc <- auc(roc1, ...)
103
  }
104 1
  if (is.null(roc2$auc) | !reuse.auc) {
105 1
    if (smoothing.args$roc2$smooth) {
106 1
      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 1
      smoothing.args$roc2$partial.auc <- NULL
112 1
      smoothing.args$roc2$partial.auc.correct <- NULL
113 1
      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 1
  if (!identical(attributes(roc1$auc)[names(attributes(roc1$auc))!="roc"], attributes(roc2$auc)[names(attributes(roc2$auc))!="roc"]))
121 1
    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 1
  if (!identical(smoothing.args$roc1, smoothing.args$roc2))
124 1
    warning("Different smoothing parameters in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
125

126
  # Check the method
127 1
  if (missing(method) | is.null(method)) {
128
    # determine method if missing
129 1
    if (has.partial.auc(roc1)) {
130
      # partial auc: go for bootstrap
131 1
      method <- "bootstrap"
132
    }
133 1
    else if (smoothing.args$roc1$smooth || smoothing.args$roc2$smooth) {
134
      # smoothing in one or both: bootstrap
135 0
      method <- "bootstrap"
136
    }
137 1
    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 1
      method <- "delong"
143
    }
144
  }
145
  else {
146 1
    method <- match.arg(method)
147 1
    if (method == "delong") {
148
      # delong NA to pAUC: warn + change
149 1
      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 1
      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 1
      if (roc1$direction != roc2$direction)
156 0
        warning("DeLong method should not be applied to ROC curves with a different direction.")
157
    }
158 1
    else if (method == "obuchowski") {
159 1
      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 1
      if ((has.partial.auc(roc1) && attr(roc1$auc, "partial.auc.focus") == "sensitivity")
163 1
          || (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 1
      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 1
  if (method == "delong") {
172 1
    n <- length(roc1$controls)
173 1
    m <- length(roc1$cases)
174

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

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

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

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

191 1
    if (roc1$percent) {
192 1
      cov <- cov * (100^2)
193
    }
194
  }
195
  else { # method == "bootstrap"
196
    # Check if called with density.cases or density.controls
197 1
    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 1
    if(class(progress) != "list") {
201 1
      progress <- roc.utils.get.progress.bar(progress, title="Bootstrap covariance", label="Bootstrap in progress...", ...)
202
    }
203
    
204 1
    cov <- bootstrap.cov(roc1, roc2, boot.n, boot.stratified, boot.return, smoothing.args, progress, parallel)
205
  }
206

207 1
  return(cov)
208
}

Read our documentation on viewing source code .

Loading