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
plot.roc <- function(x, ...) {
21 14
  UseMethod("plot.roc")
22
}
23

24
plot.roc.formula <- function(x, data, subset, na.action, ...) {
25 14
	data.missing <- missing(data)
26 14
	call <- match.call()
27 14
	names(call)[2] <- "formula" # forced to be x by definition of plot
28 14
	roc.data <- roc.utils.extract.formula(formula=x, data, subset, na.action, ..., 
29 14
										  data.missing = data.missing,
30 14
										  call = call)
31 14
	if (length(roc.data$predictor.name) > 1) {
32 0
		stop("Only one predictor supported in 'plot.roc'.")
33
	}
34 14
	response <- roc.data$response
35 14
	predictor <- roc.data$predictors[, 1]
36
	
37 14
	roc <- roc(response, predictor, plot=TRUE, ...)
38 14
	roc$call <- match.call()
39 14
	invisible(roc)
40
}
41

42
plot.roc.default <- function(x, predictor, ...) {
43 14
  roc <- roc(x, predictor, plot=TRUE, ...)
44 14
  roc$call <- match.call()
45 14
  invisible(roc)
46
}
47

48
plot.roc.smooth.roc <- plot.smooth.roc <- function(x, ...) {
49 0
  invisible(plot.roc.roc(x, ...)) # force usage of plot.roc.roc: only print.thres not working
50
}
51

52
plot.roc.roc <- function(x,
53
                         add=FALSE,
54
                         reuse.auc=TRUE,
55
                         axes=TRUE,
56
                         legacy.axes=FALSE,
57 14
                         xlim=if(x$percent){c(100, 0)} else{c(1, 0)},
58 14
                         ylim=if(x$percent){c(0, 100)} else{c(0, 1)},
59
                         xlab=ifelse(x$percent, ifelse(legacy.axes, "100 - Specificity (%)", "Specificity (%)"), ifelse(legacy.axes, "1 - Specificity", "Specificity")),
60
                         ylab=ifelse(x$percent, "Sensitivity (%)", "Sensitivity"),
61
                         asp=1,
62
                         mar=c(4, 4, 2, 2)+.1,
63
                         mgp=c(2.5, 1, 0),
64
                         # col, lty and lwd for the ROC line only
65
                         col=par("col"),
66
                         lty=par("lty"),
67
                         lwd=2,
68
                         type="l",
69
                         # Identity line
70
                         identity=!add,
71
                         identity.col="darkgrey",
72
                         identity.lty=1,
73
                         identity.lwd=1,
74
                         # Print the thresholds on the plot
75
                         print.thres=FALSE,
76
                         print.thres.pch=20,
77
                         print.thres.adj=c(-.05,1.25),
78
                         print.thres.col="black",
79
                         print.thres.pattern=ifelse(x$percent, "%.1f (%.1f%%, %.1f%%)", "%.3f (%.3f, %.3f)"),
80
                         print.thres.cex=par("cex"),
81
                         print.thres.pattern.cex=print.thres.cex,
82
                         print.thres.best.method=NULL,
83
                         print.thres.best.weights=c(1, 0.5),
84
                         # Print the AUC on the plot
85
                         print.auc=FALSE,
86
                         print.auc.pattern=NULL,
87
                         print.auc.x=ifelse(x$percent, 50, .5), 
88
                         print.auc.y=ifelse(x$percent, 50, .5),
89
                         print.auc.adj=c(0,1),
90
                         print.auc.col=col,
91
                         print.auc.cex=par("cex"),
92
                         # Grid
93
                         grid=FALSE,
94
                         grid.v={
95 0
                           if(is.logical(grid) && grid[1]==TRUE){seq(0, 1, 0.1) * ifelse(x$percent, 100, 1)}
96 0
                           else if(is.numeric(grid)) {seq(0, ifelse(x$percent, 100, 1), grid[1])}
97 14
                           else {NULL}
98
                         },
99
                         grid.h={
100 14
                           if (length(grid) == 1) {grid.v}
101 0
                           else if (is.logical(grid) && grid[2]==TRUE){seq(0, 1, 0.1) * ifelse(x$percent, 100, 1)}
102 0
                           else if(is.numeric(grid)) {seq(0, ifelse(x$percent, 100, 1), grid[2])}
103 0
                           else {NULL}
104
                         },
105
                         # for grid.lty, grid.lwd and grid.col, a length 2 value specifies both values for vertical (1) and horizontal (2) grid
106
                         grid.lty=3,
107
                         grid.lwd=1,
108
                         grid.col="#DDDDDD",
109
                         # Polygon for the auc
110
                         auc.polygon=FALSE,
111
                         auc.polygon.col="gainsboro", # Other arguments can be passed to polygon() using "..." (for these two we cannot)
112
                         auc.polygon.lty=par("lty"),
113
                         auc.polygon.density=NULL,
114
                         auc.polygon.angle=45,
115
                         auc.polygon.border=NULL,
116
                         # Should we show the maximum possible area as another polygon?
117
                         max.auc.polygon=FALSE,
118
                         max.auc.polygon.col="#EEEEEE", # Other arguments can be passed to polygon() using "..." (for these two we cannot)
119
                         max.auc.polygon.lty=par("lty"),
120
                         max.auc.polygon.density=NULL,
121
                         max.auc.polygon.angle=45,
122
                         max.auc.polygon.border=NULL,
123
                         # Confidence interval
124
                         ci=!is.null(x$ci),
125
                         ci.type=c("bars", "shape", "no"),
126
                         ci.col=ifelse(ci.type=="bars", par("fg"), "gainsboro"),
127
                         ...
128
                         ) {
129 14
  percent <- x$percent
130
  
131 14
  if (max.auc.polygon | auc.polygon | print.auc) {# we need the auc here
132 14
    if (is.null(x$auc) | !reuse.auc)
133 14
      x$auc <- auc(x, ...)
134 14
    partial.auc <- attr(x$auc, "partial.auc")
135 14
    partial.auc.focus <- attr(x$auc, "partial.auc.focus")
136
  }
137

138
  # compute a reasonable default for print.auc.pattern if required
139 14
  if (print.auc & is.null(print.auc.pattern)) {
140 14
    print.auc.pattern <- ifelse(identical(partial.auc, FALSE), "AUC: ", "Partial AUC: ")
141 14
    print.auc.pattern <- paste(print.auc.pattern, ifelse(percent, "%.1f%%", "%.3f"), sep="")
142 14
    if (ci && methods::is(x$ci, "ci.auc"))
143 14
      print.auc.pattern <- paste(print.auc.pattern, " (", ifelse(percent, "%.1f%%", "%.3f"), "\u2013", ifelse(percent, "%.1f%%", "%.3f"), ")",sep="")
144
  }
145
    
146
  # get and sort the sensitivities and specificities
147 14
  se <- sort(x$se, decreasing=TRUE)
148 14
  sp <- sort(x$sp, decreasing=FALSE) 
149 14
  if (!add) {
150 14
    opar <- par(mar=mar, mgp=mgp)
151 14
    on.exit(par(opar))
152
    # type="n" to plot background lines and polygon shapes first. We will add the line later. axes=FALSE, we'll add them later according to legacy.axis
153 14
    suppressWarnings(plot(x$sp, x$se, xlab=xlab, ylab=ylab, type="n", axes=FALSE, xlim=xlim, ylim=ylim, lwd=lwd, asp=asp, ...))
154

155
    # As we had axes=FALSE we need to add them again unless axes=FALSE
156 14
    if (axes) {
157 14
      box()
158
      # axis behave differently when at and labels are passed (no decimals on 1 and 0), 
159
      # so handle each case separately and consistently across axes
160 14
      if (legacy.axes) {
161 14
      	lab.at <- axTicks(side=1)
162 14
      	lab.labels <- format(ifelse(x$percent, 100, 1) - lab.at)
163 14
      	suppressWarnings(axis(side=1, at=lab.at, labels=lab.labels, ...))
164 14
      	lab.at <- axTicks(side=2)
165 14
      	suppressWarnings(axis(side=2, at=lab.at, labels=format(lab.at), ...))
166
      } 
167
      else {
168 14
      	suppressWarnings(axis(side=1, ...))
169 14
      	suppressWarnings(axis(side=2, ...))
170
      }
171
    }
172
  }
173

174
  # Plot the grid
175
  # make sure grid.lty, grid.lwd and grid.col are at least of length 2
176 14
  grid.lty <- rep(grid.lty, length.out=2)
177 14
  grid.lwd <- rep(grid.lwd, length.out=2)
178 14
  grid.col <- rep(grid.col, length.out=2)
179 14
  if (!is.null(grid.v)) {
180 0
    suppressWarnings(abline(v=grid.v, lty=grid.lty[1], col=grid.col[1], lwd=grid.lwd[1], ...))
181
  }
182 14
  if (!is.null(grid.h)) {
183 0
    suppressWarnings(abline(h=grid.h, lty=grid.lty[2], col=grid.col[2], lwd=grid.lwd[2], ...))
184
  }
185

186
  # Plot the polygon displaying the maximal area
187 14
  if (max.auc.polygon) {
188 14
    if (identical(partial.auc, FALSE)) {
189 0
      map.y <- c(0, 1, 1, 0) * ifelse(percent, 100, 1)
190 0
      map.x <- c(1, 1, 0, 0) * ifelse(percent, 100, 1)
191
    }
192
    else {
193 14
      if (partial.auc.focus=="sensitivity") {
194 14
        map.y <- c(partial.auc[2], partial.auc[2], partial.auc[1], partial.auc[1]) 
195 14
        map.x <- c(0, 1, 1, 0) * ifelse(percent, 100, 1) 
196
      }
197
      else {
198 14
        map.y <- c(0, 1, 1, 0) * ifelse(percent, 100, 1) 
199 14
        map.x <- c(partial.auc[2], partial.auc[2], partial.auc[1], partial.auc[1])
200
      }
201
    }
202 14
    suppressWarnings(polygon(map.x, map.y, col=max.auc.polygon.col, lty=max.auc.polygon.lty, border=max.auc.polygon.border, density=max.auc.polygon.density, angle=max.auc.polygon.angle, ...))
203
  }
204
  # Plot the ci shape
205 14
  if (ci && ! methods::is(x$ci, "ci.auc")) {
206 14
    ci.type <- match.arg(ci.type)
207 14
    if (ci.type=="shape")
208 14
      plot(x$ci, type="shape", col=ci.col, no.roc=TRUE, ...)
209
  }
210
  # Plot the polygon displaying the actual area
211 14
  if (auc.polygon) {
212 14
    if (identical(partial.auc, FALSE)) {
213 0
      suppressWarnings(polygon(c(sp, 0), c(se, 0), col=auc.polygon.col, lty=auc.polygon.lty, border=auc.polygon.border, density=auc.polygon.density, angle=auc.polygon.angle, ...))
214
    }
215
    else {
216 14
      if (partial.auc.focus == "sensitivity") {
217 14
        x.all <- rev(se)
218 14
        y.all <- rev(sp)
219
      }
220
      else {
221 14
        x.all <- sp
222 14
        y.all <- se
223
      }
224
      # find the SEs and SPs in the interval
225 14
      x.int <- x.all[x.all <= partial.auc[1] & x.all >= partial.auc[2]]
226 14
      y.int <- y.all[x.all <= partial.auc[1] & x.all >= partial.auc[2]]
227
      # if the upper limit is not exactly present in SPs, interpolate
228 14
      if (!(partial.auc[1] %in% x.int)) {
229 0
        x.int <- c(x.int, partial.auc[1])
230
        # find the limit indices
231 0
        idx.out <- match(FALSE, x.all < partial.auc[1])
232 0
        idx.in <- idx.out - 1
233
        # interpolate y
234 0
        proportion.start <- (partial.auc[1] - x.all[idx.out]) / (x.all[idx.in] - x.all[idx.out])
235 0
        y.start <- y.all[idx.out] - proportion.start * (y.all[idx.out] - y.all[idx.in])
236 0
        y.int <- c(y.int, y.start)
237
      }
238
      # if the lower limit is not exactly present in SPs, interpolate
239 14
      if (!(partial.auc[2] %in% x.int)) {
240 14
        x.int <- c(partial.auc[2], x.int)
241
        # find the limit indices
242 14
        idx.out <- length(x.all) - match(TRUE, rev(x.all) < partial.auc[2]) + 1
243 14
        idx.in <- idx.out + 1
244
        # interpolate y
245 14
        proportion.end <- (x.all[idx.in] - partial.auc[2]) / (x.all[idx.in] - x.all[idx.out])
246 14
        y.end <- y.all[idx.in] + proportion.end * (y.all[idx.out] - y.all[idx.in])
247 14
        y.int <- c(y.end, y.int)
248
      }
249
      # anchor to baseline
250 14
      x.int <- c(partial.auc[2], x.int, partial.auc[1])
251 14
      y.int <- c(0, y.int, 0)
252 14
      if (partial.auc.focus == "sensitivity") {
253
        # for SE, invert x and y again
254 14
        suppressWarnings(polygon(y.int, x.int, col=auc.polygon.col, lty=auc.polygon.lty, border=auc.polygon.border, density=auc.polygon.density, angle=auc.polygon.angle, ...))
255
      }
256
      else {
257 14
        suppressWarnings(polygon(x.int, y.int, col=auc.polygon.col, lty=auc.polygon.lty, border=auc.polygon.border, density=auc.polygon.density, angle=auc.polygon.angle, ...))
258
      }
259
    }
260
  }
261
  # Identity line
262 14
  if (identity) suppressWarnings(abline(ifelse(percent, 100, 1), -1, col=identity.col, lwd=identity.lwd, lty=identity.lty, ...))
263
  # Actually plot the ROC curve
264 14
  suppressWarnings(lines(sp, se, type=type, lwd=lwd, col=col, lty=lty, ...))
265
  # Plot the ci bars
266 14
  if (ci && !methods::is(x$ci, "ci.auc")) {
267 14
    if (ci.type=="bars")
268 14
      plot(x$ci, type="bars", col=ci.col, ...)
269
  }
270
  # Print the thresholds on the curve if print.thres is TRUE
271 14
  if (isTRUE(print.thres))
272 0
    print.thres <- "best"
273 14
  if (is.character(print.thres))
274 14
    print.thres <- match.arg(print.thres, c("no", "all", "local maximas", "best"))
275 14
  if (methods::is(x, "smooth.roc")) {
276 0
    if (is.numeric(print.thres))
277 0
      stop("Numeric 'print.thres' unsupported on a smoothed ROC plot.")
278 0
    else if (print.thres == "all" || print.thres == "local maximas")
279 0
      stop("'all' and 'local maximas' 'print.thres' unsupported on a smoothed ROC plot.") 
280 0
    else if (print.thres == "best") {
281 0
      co <- coords(x, print.thres, best.method=print.thres.best.method, best.weights=print.thres.best.weights, transpose=FALSE)
282 0
        suppressWarnings(points(co$specificity, co$sensitivity, pch=print.thres.pch, cex=print.thres.cex, col=print.thres.col, ...))
283 0
        suppressWarnings(text(co$specificity, co$sensitivity, sprintf(print.thres.pattern, NA, co$specificity, co$sensitivity), adj=print.thres.adj, cex=print.thres.pattern.cex, col=print.thres.col, ...))
284
     
285
    } # else print.thres == no > do nothing
286
  }
287 14
  else if (is.numeric(print.thres) || is.character(print.thres)) {
288 14
    if (is.character(print.thres) && print.thres == "no") {} # do nothing
289
    else {
290 14
      co <- coords(x, print.thres, best.method=print.thres.best.method, best.weights=print.thres.best.weights, transpose=FALSE)
291 14
      suppressWarnings(points(co$specificity, co$sensitivity, pch=print.thres.pch, cex=print.thres.cex, col=print.thres.col, ...))
292 14
      suppressWarnings(text(co$specificity, co$sensitivity, sprintf(print.thres.pattern, co$threshold, co$specificity, co$sensitivity), adj=print.thres.adj, cex=print.thres.pattern.cex, col=print.thres.col, ...))
293
    }
294
  }
295

296
  # Print the AUC on the plot
297 14
  if (print.auc) {
298 14
    if (ci && methods::is(x$ci, "ci.auc")) {
299 14
      labels <- sprintf(print.auc.pattern, x$auc, x$ci[1], x$ci[3])
300 14
      suppressWarnings(text(print.auc.x, print.auc.y, labels, adj=print.auc.adj, cex=print.auc.cex, col=print.auc.col, ...))
301
    }
302
    else
303 14
      labels <- sprintf(print.auc.pattern, x$auc)
304 14
    suppressWarnings(text(print.auc.x, print.auc.y, labels, adj=print.auc.adj, cex=print.auc.cex, col=print.auc.col, ...))
305
  }
306
  
307 14
  invisible(x)
308
}

Read our documentation on viewing source code .

Loading