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

23
coords.smooth.roc <- function(smooth.roc,
24
                              x,
25
                              input,
26
                              ret=c("specificity", "sensitivity"),
27
                              as.list=FALSE,
28
                              drop=TRUE,
29
                              best.method=c("youden", "closest.topleft"),
30
                              best.weights=c(1, 0.5),
31
                              transpose = FALSE,
32
                              as.matrix = FALSE,
33
                              ...) {
34
  # make sure x was provided
35 33
  if (missing(x))
36 0
    stop("'x' must be a numeric or character vector.")
37
  
38
  # Warn about future change in transpose <https://github.com/xrobin/pROC/issues/54>
39 33
  if (missing(transpose) || is.null(transpose)) {
40 33
    transpose <- FALSE
41 33
    warning("The 'transpose' argument is set to FALSE by default since pROC 1.16. Set transpose = FALSE explicitly silence this warning or transpose = TRUE to revert to the previous behavior. Type help(coords_transpose) for additional information.")
42
  }
43
  
44
  # match return
45 33
  ret <- roc.utils.match.coords.ret.args(ret, threshold = FALSE)
46

47 33
  if (is.character(x)) {
48 33
    x <- match.arg(x, c("best")) # no thresholds in smoothed roc: only best is possible
49 33
    partial.auc <- attr(smooth.roc$auc, "partial.auc")
50

51
    # cheat: allow the user to pass "topleft"
52 33
    best.method <- match.arg(best.method[1], c("youden", "closest.topleft", "topleft"))
53 33
    if (best.method == "topleft") {
54 0
    	best.method <- "closest.topleft"
55
    }
56 33
    optim.crit <- roc.utils.optim.crit(smooth.roc$sensitivities, smooth.roc$specificities,
57 33
    								   ifelse(smooth.roc$percent, 100, 1),
58 33
    								   best.weights, best.method)
59
    
60 33
    if (is.null(smooth.roc$auc) || identical(partial.auc, FALSE)) {
61 33
      se <- smooth.roc$sensitivities[optim.crit==max(optim.crit)]
62 33
      sp <- smooth.roc$specificities[optim.crit==max(optim.crit)]
63 33
      optim.crit <- optim.crit[optim.crit==max(optim.crit)]
64
    }
65
    else {
66 33
      if (attr(smooth.roc$auc, "partial.auc.focus") == "sensitivity") {
67 33
        optim.crit.partial <- (optim.crit)[smooth.roc$se <= partial.auc[1] & smooth.roc$se >= partial.auc[2]]
68 33
        se <- smooth.roc$sensitivities[smooth.roc$sensitivities <= partial.auc[1] & smooth.roc$sensitivities >= partial.auc[2]][optim.crit.partial==max(optim.crit.partial)]
69 33
        sp <- smooth.roc$specificities[smooth.roc$sensitivities <= partial.auc[1] & smooth.roc$sensitivities >= partial.auc[2]][optim.crit.partial==max(optim.crit.partial)]
70 33
        optim.crit <- optim.crit[smooth.roc$sensitivities <= partial.auc[1] & smooth.roc$sensitivities >= partial.auc[2]][optim.crit.partial==max(optim.crit.partial)]
71
      }
72
      else {
73 33
        optim.crit.partial <- (optim.crit)[smooth.roc$sp <= partial.auc[1] & smooth.roc$sp >= partial.auc[2]]
74 33
        se <- smooth.roc$sensitivities[smooth.roc$specificities <= partial.auc[1] & smooth.roc$specificities >= partial.auc[2]][optim.crit.partial==max(optim.crit.partial)]
75 33
        sp <- smooth.roc$specificities[smooth.roc$specificities <= partial.auc[1] & smooth.roc$specificities >= partial.auc[2]][optim.crit.partial==max(optim.crit.partial)]
76 33
        optim.crit <- optim.crit[smooth.roc$specificities <= partial.auc[1] & smooth.roc$specificities >= partial.auc[2]][optim.crit.partial==max(optim.crit.partial)]
77
      }
78
    }
79
    
80 33
    if (any(! ret %in% c("specificity", "sensitivity", best.method))) {
81
    	# Deduce additional tn, tp, fn, fp, npv, ppv
82 33
    	res <- roc.utils.calc.coords(smooth.roc, NA, se, sp, best.weights)
83
    }
84
    else {
85 33
    	res <- cbind(
86 33
    		specificity = sp,
87 33
    		sensitivity = se,
88 33
    		best.method = ifelse(best.method == "youden", 1, -1) * optim.crit
89
    	)
90 33
    	colnames(res)[3] <- best.method
91
    }
92

93 33
    if (as.list) {
94 33
      warning("'as.list' is deprecated and will be removed in a future version.")
95 33
    	list <- apply(t(res)[ret, , drop=FALSE], 2, as.list)
96 33
    	if (drop == TRUE && length(x) == 1) {
97 33
    		return(list[[1]])
98
    	}
99 33
    	return(list)
100
    }
101 33
    else if (transpose) {
102 33
      rownames(res) <- NULL
103 33
      return(t(res)[ret,, drop=drop])
104
    }
105
    else {
106 33
      if (missing(drop) ) {
107 33
        drop = FALSE
108
      }
109 33
      if (! as.matrix) {
110 33
        res <- as.data.frame(res)
111
      }
112 33
      return(res[, ret, drop=drop])
113
    }
114
  }
115
  
116
  # Adjust drop for downstream call
117 33
  if (missing(drop) && ! transpose) {
118 33
    drop = FALSE
119
  }
120
  
121
  # match input 
122 33
  input <- roc.utils.match.coords.input.args(input, threshold = FALSE)
123
  
124
  # use coords.roc
125 33
  smooth.roc$thresholds <- rep(NA, length(smooth.roc$specificities))
126 33
  return(coords.roc(smooth.roc, x, input, ret, as.list, drop, 
127 33
                    transpose = transpose, as.matrix = as.matrix, ...))
128
}
129

130
coords.roc <- function(roc,
131
                       x,
132
                       input="threshold",
133
                       ret=c("threshold", "specificity", "sensitivity"),
134
                       as.list=FALSE,
135
                       drop=TRUE,
136
                       best.method=c("youden", "closest.topleft"),
137
                       best.weights=c(1, 0.5), 
138
                       transpose = FALSE,
139
                       as.matrix = FALSE,
140
                       ...) {
141
  # make sure x was provided
142 33
  if (missing(x) || is.null(x) || (length(x) == 0 && !is.numeric(x))) {
143 33
    x <- "all"
144
  }
145 33
  else if (length(x) == 0 && is.numeric(x)) {
146 33
    stop("Numeric 'x' has length 0")
147
  }
148
  
149
  # Warn about future change in transpose <https://github.com/xrobin/pROC/issues/54>
150 33
  if (missing(transpose) || is.null(transpose)) {
151 33
    transpose <- FALSE
152 33
    warning("The 'transpose' argument to FALSE by default since pROC 1.16. Set transpose = TRUE explicitly to revert to the previous behavior, or transpose = TRUE to silence this warning. Type help(coords_transpose) for additional information.")
153
  }
154
  
155
  # match input 
156 33
  input <- roc.utils.match.coords.input.args(input)
157
  # match return
158 33
  ret <- roc.utils.match.coords.ret.args(ret)
159
  # make sure the sort of roc is correct
160 33
  roc <- sort(roc)
161

162 33
  if (is.character(x)) {
163 33
    x <- match.arg(x, c("all", "local maximas", "best"))
164 33
    partial.auc <- attr(roc$auc, "partial.auc")
165 33
    if (x == "all") {
166
      # Pre-filter thresholds based on partial.auc
167 33
      if (is.null(roc$auc) || identical(partial.auc, FALSE)) {
168 33
        se <- roc$se
169 33
        sp <- roc$sp
170 33
        thres <- roc$thresholds
171
      }
172
      else {
173 33
        if (attr(roc$auc, "partial.auc.focus") == "sensitivity") {
174 33
          se <- roc$se[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
175 33
          sp <- roc$sp[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
176 33
          thres <- roc$thresholds[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
177
        }
178
        else {
179 33
          se <- roc$se[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
180 33
          sp <- roc$sp[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
181 33
          thres <- roc$thresholds[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
182
        }
183
      }
184 33
      if (length(thres) == 0) {
185 0
      	warning("No coordinates found, returning NULL. This is possibly cased by a too small partial AUC interval.")
186 0
      	return(NULL)
187
      }
188 33
      res <- cbind(
189 33
      	threshold = thres,
190 33
      	specificity = sp,
191 33
      	sensitivity = se
192
      )
193
    }
194 33
    else if (x == "local maximas") {
195
      # Pre-filter thresholds based on partial.auc
196 33
      if (is.null(roc$auc) || identical(partial.auc, FALSE)) {
197 33
        se <- roc$se
198 33
        sp <- roc$sp
199 33
        thres <- roc$thresholds
200
      }
201
      else {
202 33
        if (attr(roc$auc, "partial.auc.focus") == "sensitivity") {
203 33
          se <- roc$se[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
204 33
          sp <- roc$sp[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
205 33
          thres <- roc$thresholds[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
206
        }
207
        else {
208 33
          se <- roc$se[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
209 33
          sp <- roc$sp[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
210 33
          thres <- roc$thresholds[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
211
        }
212
      }
213 33
      if (length(thres) == 0) {
214 0
      	warning("No coordinates found, returning NULL. This is possibly cased by a too small partial AUC interval.")
215 0
      	return(NULL)
216
      }
217 33
      lm.idx <- roc.utils.max.thresholds.idx(thres, sp=sp, se=se)
218 33
      res <- cbind(
219 33
      	threshold = thres[lm.idx],
220 33
      	specificity = sp[lm.idx],
221 33
      	sensitivity = se[lm.idx]
222
      )
223
    }
224
    else { # x == "best"
225
      # cheat: allow the user to pass "topleft"
226 33
      best.method <- match.arg(best.method[1], c("youden", "closest.topleft", "topleft"))
227 33
      if (best.method == "topleft") {
228 0
        best.method <- "closest.topleft"
229
      }
230 33
      optim.crit <- roc.utils.optim.crit(roc$sensitivities, roc$specificities,
231 33
      								   ifelse(roc$percent, 100, 1),
232 33
      								   best.weights, best.method)
233

234
      # Filter thresholds based on partial.auc
235 33
      if (is.null(roc$auc) || identical(partial.auc, FALSE)) {
236 33
      	se <- roc$se[optim.crit==max(optim.crit)]
237 33
      	sp <- roc$sp[optim.crit==max(optim.crit)]
238 33
        thres <- roc$thresholds[optim.crit==max(optim.crit)]
239 33
        optim.crit <- optim.crit[optim.crit==max(optim.crit)]
240
      }
241
      else {
242 33
        if (attr(roc$auc, "partial.auc.focus") == "sensitivity") {
243 33
          optim.crit <- (optim.crit)[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]]
244 33
          se <- roc$se[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]][optim.crit==max(optim.crit)]
245 33
          sp <- roc$sp[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]][optim.crit==max(optim.crit)]
246 33
          thres <- roc$thresholds[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]][optim.crit==max(optim.crit)]
247 33
          optim.crit <- optim.crit[roc$se <= partial.auc[1] & roc$se >= partial.auc[2]][optim.crit==max(optim.crit)]
248
        }
249
        else {
250 33
          optim.crit <- (optim.crit)[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]]
251 33
          se <- roc$se[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]][optim.crit==max(optim.crit)]
252 33
          sp <- roc$sp[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]][optim.crit==max(optim.crit)]
253 33
          thres <- roc$thresholds[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]][optim.crit==max(optim.crit)]
254 33
          optim.crit <- optim.crit[roc$sp <= partial.auc[1] & roc$sp >= partial.auc[2]][optim.crit==max(optim.crit)]
255
        }
256
      }
257 33
      if (length(thres) == 0) {
258 0
      	warning("No coordinates found, returning NULL. This is possibly cased by a too small partial AUC interval.")
259 0
      	return(NULL)
260
      }
261 33
      res <- cbind(
262 33
      	threshold = thres,
263 33
      	specificity = sp,
264 33
      	sensitivity = se,
265 33
      	best.method = ifelse(best.method == "youden", 1, -1) * optim.crit
266
      )
267 33
      colnames(res)[4] <- best.method
268
    }
269
  }
270 33
  else if (is.numeric(x)) {
271 33
    if (input == "threshold") {
272 33
    	thr_idx <- roc.utils.thr.idx(roc, x)
273 33
    	res <- cbind(
274 33
    		threshold = x, # roc$thresholds[thr_idx], # user-supplied vs ours.
275 33
    		specificity = roc$specificities[thr_idx],
276 33
    		sensitivity = roc$sensitivities[thr_idx]
277
    	)
278
    }
279
    else {
280
      # Arbitrary coord given in input.
281
      # We could be tempted to use all_coords directly.
282
      # However any non monotone coordinate in ret will be inaccurate
283
      # when interpolated. Therefore it is safer to only interpolate
284
      # se and sp and re-calculate the remaining coords later.
285 33
      res <- cbind(threshold = rep(NA, length(x)),
286 33
                        sensitivity = rep(NA, length(x)),
287 33
                        specificity = rep(NA, length(x))
288
                   )
289 33
      if (input %in% c("sensitivity", "specificity")) {
290
        # Shortcut slow roc.utils.calc.coords
291 33
        se <- roc$sensitivities
292 33
        sp <- roc$specificities
293 33
        if (methods::is(roc, "smooth.roc")) {
294 33
          thr <- rep(NA, length(roc$sensitivities))
295
        }
296
        else {
297 33
          thr <- roc$thresholds
298
        }
299 33
        if (input == "sensitivity") {
300 33
          input_values <- se
301
        }
302
        else {
303 33
          input_values <- sp
304
        }
305
      }
306
      else {
307 0
        all_coords <- roc.utils.calc.coords(roc, rep(NA, length(roc$sensitivities)), roc$sensitivities, roc$specificities, best.weights)
308 0
        input_values <- all_coords[, input]
309 0
        se <- all_coords[, "sensitivity"]
310 0
        sp <- all_coords[, "specificity"]
311 0
        thr <- all_coords[, "threshold"]
312
      }
313 33
      for (i in seq_along(x)) {
314 33
        value <- x[i]
315 33
        if (value < min(input_values) || value > max(input_values)) {
316 33
          stop(sprintf("Input %s (%s) not in range (%s-%s)", input, value, 
317 33
                       min(input_values), max(input_values)))
318
        }
319
        
320 33
        idx <- which(input_values == value)
321 33
        if (length(idx) > 1) {
322
          # More than one to pick from. Need to take best
323
          # according to sorting
324 33
          if (coord.is.decreasing[input]) {
325 33
            idx <- idx[length(idx)] # last
326
          }
327
          else {
328 33
            idx <- idx[1] # first
329
          }
330
        }
331 33
        if (length(idx) == 1) {
332
          # Exactly one to pick from
333 33
          res[i,] <- c(thr[idx], se[idx], sp[idx])
334
        }
335
        else {
336
          # Need to interpolate
337 33
          if (coord.is.decreasing[input]) {
338 33
            idx.next <- match(TRUE, input_values < value)
339
          }
340
          else {
341 33
            idx.next <- match(TRUE, input_values > value)
342
          }
343 33
          proportion <- (value - input_values[idx.next]) / (input_values[idx.next - 1] - input_values[idx.next])
344 33
          int.se <- se[idx.next] + proportion * (se[idx.next - 1] - se[idx.next])
345 33
          int.sp <- sp[idx.next] + proportion * (sp[idx.next - 1] - sp[idx.next])
346 33
          res[i, 2:3] <- c(int.se, int.sp)
347
        }
348
      }
349
    }
350
  }
351
  else {
352 33
    stop("'x' must be a numeric or character vector.")
353
  }
354
  
355 33
  if (any(! ret %in% colnames(res))) {
356
  	# Deduce additional tn, tp, fn, fp, npv, ppv
357 33
  	res <- roc.utils.calc.coords(roc, res[, "threshold"], res[, "sensitivity"], res[, "specificity"], best.weights)
358
  }
359

360 33
  if (as.list) {
361 33
    warning("'as.list' is deprecated and will be removed in a future version.")
362 33
  	list <- apply(t(res)[ret, , drop=FALSE], 2, as.list)
363 33
  	if (drop == TRUE && length(x) == 1) {
364 33
  		return(list[[1]])
365
  	}
366 33
  	return(list)
367
  }
368 33
  else if (transpose) {
369 33
    rownames(res) <- NULL
370 33
    return(t(res)[ret,, drop=drop])
371
  }
372
  else {
373
  	# HACK:
374
  	# We need an exception for r4lineups that will keep the old drop = TRUE 
375
  	# behavior, until r4lineups gets updated. This is an ugly hack but allows
376
  	# us to switch to a better drop = FALSE for everyone else
377 33
    if (missing(drop)) {
378 33
    	if (sys.nframe() > 2 && 
379 33
    	    length(deparse(sys.call(-2))) == 1 && 
380 33
    	    deparse(sys.call(-2)) == 'make_rocdata(df_confacc)' && 
381 33
    	    length(deparse(sys.call(-1))) == 1 && (
382 33
    	      deparse(sys.call(-1)) == 'coords(rocobj, "all", ret = c("tp"))' ||
383 33
    	      deparse(sys.call(-1)) == 'coords(rocobj, "all", ret = "fp")'
384
    	    )) {
385
    		# We're in r4lineups
386 0
    		drop = TRUE
387
    	}
388
    	else {
389 33
    		drop = FALSE
390
    	}
391
    }
392 33
    if (! as.matrix) {
393 33
      res <- as.data.frame(res)
394
    }
395 33
    return(res[, ret, drop=drop])
396
  	
397
  }
398
}

Read our documentation on viewing source code .

Loading