xrobin / pROC
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
# Helper functions for the ROC curves. These functions should not be called directly as they peform very specific tasks and do nearly no argument validity checks. Not documented in RD and not exported.
21

22
# returns a list of sensitivities (se) and specificities (sp) for the given data. Robust algorithm
23
roc.utils.perfs.all.safe <- function(thresholds, controls, cases, direction) {
24 8
  perf.matrix <- sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=direction)
25
  #stopifnot(identical(roc.utils.perfs.all.fast(thresholds, controls, cases, direction), list(se=perf.matrix[2,], sp=perf.matrix[1,])))
26 8
  return(list(se=perf.matrix[2,], sp=perf.matrix[1,]))
27
}
28

29

30
roc.utils.perfs.all.test <- function(thresholds, controls, cases, direction) {
31 8
	perfs.safe <- roc.utils.perfs.all.safe(thresholds=thresholds, controls=controls, cases=cases, direction=direction)
32 8
	perfs.fast <- roc.utils.perfs.all.fast(thresholds=thresholds, controls=controls, cases=cases, direction=direction)
33 8
	perfs.C <- rocUtilsPerfsAllC(thresholds=thresholds, controls=controls, cases=cases, direction=direction)
34 8
	if (! (identical(perfs.safe, perfs.fast) && identical(perfs.safe, perfs.C))) {
35 0
		sessionInfo <- sessionInfo()
36 0
		save(thresholds, controls, cases, direction, sessionInfo, file="pROC_bug.RData")
37 0
		stop(sprintf("pROC: algorithms returned different values. Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", utils::packageDescription("pROC")$BugReports))
38
	}
39 8
	return(perfs.safe)
40
}
41

42

43
# returns a list of sensitivities (se) and specificities (sp) for the given data. Fast algorithm
44
roc.utils.perfs.all.fast <- function(thresholds, controls, cases, direction) {
45 8
  ncontrols <- length(controls)
46 8
  ncases <- length(cases)
47 8
  predictor <- c(controls, cases)
48 8
  response <- c(rep(0, length(controls)), rep(1, length(cases)))
49 8
  decr <- direction=="<"
50 8
  predictor.order <- order(predictor, decreasing=decr)
51 8
  predictor.sorted <- predictor[predictor.order]
52 8
  response.sorted <- response[predictor.order]
53
  
54 8
  tp <- cumsum(response.sorted==1)
55 8
  fp <- cumsum(response.sorted==0)
56 8
  se <- tp / ncases
57 8
  sp <- (ncontrols - fp) / ncontrols
58
  # filter duplicate thresholds
59 8
  dups.pred <- rev(duplicated(rev(predictor.sorted)))
60 8
  dups.sesp <- duplicated(se) & duplicated(sp)
61 8
  dups <- dups.pred | dups.sesp
62
  # Make sure we have the right length
63 8
  if (sum(!dups) != length(thresholds) - 1) {
64 0
  	sessionInfo <- sessionInfo()
65 0
  	save(thresholds, controls, cases, direction, sessionInfo, file="pROC_bug.RData")
66 0
    stop(sprintf("pROC: fast algorithm computed an incorrect number of sensitivities and specificities. Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", utils::packageDescription("pROC")$BugReports))
67
  }
68 8
  if (direction == "<") {
69 8
    se <- rev(c(0, se[!dups]))
70 8
    sp <- rev(c(1, sp[!dups]))
71
  }
72
  else {
73 8
    se <- c(0, se[!dups])
74 8
    sp <- c(1, sp[!dups])
75
  }
76 8
  return(list(se=se, sp=sp))
77
}
78

79
# As roc.utils.perfs.all but returns an "old-style" matrix (pre-fast-algo-compatible)
80
#roc.utils.perfs.all.matrix <- function(...) {
81
#  perfs <- roc.utils.perfs.all(...)
82
#  return(matrix(c(perfs$sp, perfs$se), nrow=2, byrow=TRUE))
83
#}
84

85
# returns a vector with two elements, sensitivity and specificity, given the threshold at which to evaluate the performance, the values of controls and cases and the direction of the comparison, a character '>' or '<' as controls CMP cases
86
# sp <- roc.utils.perfs(...)[1,]
87
# se <- roc.utils.perfs(...)[2,]
88
roc.utils.perfs <- function(threshold, controls, cases, direction) {
89 8
  if (direction == '>') {
90 8
    tp <- sum(cases <= threshold)
91 8
    tn <- sum(controls > threshold)
92
  }
93 8
  else if (direction == '<') {
94 8
    tp <- sum(cases >= threshold)
95 8
    tn <- sum(controls < threshold)
96
  }
97
  # return(c(sp, se))
98 8
  return(c(sp=tn/length(controls), se=tp/length(cases)))
99
}
100

101
# as roc.utils.perfs, but for densities
102
roc.utils.perfs.dens <- function(threshold, x, dens.controls, dens.cases, direction) {
103 8
  if (direction == '>') {
104 8
    tp <- sum(dens.cases[x <= threshold])
105 8
    tn <- sum(dens.controls[x > threshold])
106
  }
107 8
  else if (direction == '<') {
108 8
    tp <- sum(dens.cases[x >= threshold])
109 8
    tn <- sum(dens.controls[x < threshold])
110
  }
111
  # return(c(sp, se))
112 8
  return(c(sp=tn/sum(dens.controls), se=tp/sum(dens.cases)))
113
}
114

115
# return the thresholds to evaluate in the ROC curve, given the 'predictor' values. Returns all unique values of 'predictor' plus 2 extreme values
116
roc.utils.thresholds <- function(predictor, direction) {
117 8
  unique.candidates <- sort(unique(predictor))
118 8
  thresholds1 <- (c(-Inf, unique.candidates) + c(unique.candidates, +Inf))/2
119 8
  thresholds2 <- (c(-Inf, unique.candidates)/2 + c(unique.candidates, +Inf)/2)
120 8
  thresholds <- ifelse(abs(thresholds1) > 1e100, thresholds2, thresholds1)
121 8
  if (any(ties <- thresholds %in% predictor)) {
122
  	# If we get here, some thresholds are identical to the predictor
123
  	# This is caused by near numeric ties that caused the mean to equal
124
  	# one of the candidate
125
  	# We need to make sure we select the right threshold more carefully
126 8
  	if (direction == '>') {
127
  		# We have:
128
  		# tp <- sum(cases <= threshold)
129
  		# tn <- sum(controls > threshold)
130
  		# We need to make sure the selected threshold
131
  		# Corresponds to the lowest observation of the predictor
132
  		# Identify problematic thresholds
133
  		# rows <- which(ties)
134 8
  		for (tie.idx in which(ties)) {
135 8
  			if (thresholds[tie.idx] == unique.candidates[tie.idx - 1]) {
136
  				# We're already good, nothing to do
137
  			}
138 8
  			else if (thresholds[tie.idx] == unique.candidates[tie.idx]) {
139 8
  				thresholds[tie.idx] <- unique.candidates[tie.idx - 1]
140
  			}
141
  			else {
142 0
  				sessionInfo <- sessionInfo()
143 0
  				save(predictor, direction, sessionInfo, file="pROC_bug.RData")
144 0
  				stop(sprintf("Couldn't fix near ties in thresholds: %s, %s, %s, %s. Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", thresholds[tie.idx], unique.candidates[tie.idx - 1], unique.candidates[tie.idx], direction, utils::packageDescription("pROC")$BugReports))
145
  			}
146
  		}
147
  	}
148 8
  	else if (direction == '<') {
149
  		# We have:
150
  		# tp <- sum(cases >= threshold)
151
  		# tn <- sum(controls < threshold)
152
  		# We need to make sure the selected threshold
153
  		# Corresponds to the highest observation of the predictor
154
  		# Identify the problematic thresholds:
155
  		# rows <- which(apply(o, 1, any))
156 8
  		for (tie.idx in which(ties)) {
157 8
  			if (thresholds[tie.idx] == unique.candidates[tie.idx - 1]) {
158
  				# Easy to fix: should be unique.candidates[tie.idx]
159 8
  				thresholds[tie.idx] <- unique.candidates[tie.idx]
160 8
  			} else if (thresholds[tie.idx] == unique.candidates[tie.idx]) {
161
  				# We're already good, nothing to do
162
  			}
163
  			else {
164 0
  				sessionInfo <- sessionInfo()
165 0
  				save(predictor, direction, sessionInfo, file="pROC_bug.RData")
166 0
  				stop(sprintf("Couldn't fix near ties in thresholds: %s, %s, %s, %s. Diagnostic data saved in pROC_bug.RData. Please report this bug to <%s>.", thresholds[tie.idx], unique.candidates[tie.idx - 1], unique.candidates[tie.idx], direction, utils::packageDescription("pROC")$BugReports))
167
  			}
168
  		}
169
  	}
170
  }
171 8
  return(thresholds)
172
}
173

174
# Find all the local maximas of the ROC curve. Returns a logical vector
175
roc.utils.max.thresholds.idx <- function(thresholds, sp, se) {
176 8
  reversed <- FALSE
177 8
  if (is.unsorted(sp)) {
178
    # make sure SP is sorted increasingly, and sort thresholds accordingly
179 0
    thresholds <- rev(thresholds)
180 0
    sp <- rev(sp)
181 0
    se <- rev(se)
182 0
    reversed <- TRUE
183
  }
184
  # TODO: find whether the duplicate check is still needed.
185
  # Should have been fixed by passing only c(controls, cases)
186
  # instead of whole 'predictor' to roc.utils.thresholds in roc.default
187
  # but are there other potential issues like that?
188 8
  dup <- duplicated(data.frame(sp, se))
189 8
  thresholds <- thresholds[!dup]
190 8
  sp <- sp[!dup]
191 8
  se <- se[!dup]
192
  # Find the local maximas
193 8
  if (length(thresholds) == 1) {
194 0
    local.maximas <- TRUE # let's consider that if there is only 1 point, we should print it.
195
  }
196 8
  else if (length(thresholds) == 2) {
197 0
    local.maximas <- c(se[1] > se[2], sp[2] > sp[1])
198
  }
199
  else {
200 8
    local.maximas <- se[1] > se[2]
201 8
    for (i in 2:(length(thresholds)-1)) {
202 8
      if (sp[i] > sp[i-1] & se[i] > se[i+1])
203 8
        local.maximas <- c(local.maximas, TRUE)
204 8
      else if (sp[i] > sp[i-1] & se[i] == 0)
205 0
        local.maximas <- c(local.maximas, TRUE)
206 8
      else if (se[i] > se[i-1] & sp[i] == 1)
207 0
        local.maximas <- c(local.maximas, TRUE)
208
      else
209 8
        local.maximas <- c(local.maximas, FALSE)
210
    }
211 8
    local.maximas <- c(local.maximas, sp[length(thresholds)] > sp[length(thresholds)-1])
212
  }
213 8
  if (any(dup)) {
214 0
    lms <- rep(FALSE, length(dup))
215 0
    lms[!dup] <- local.maximas
216 0
    local.maximas <- lms
217
  }
218 8
  if (reversed)
219 0
    rev(local.maximas)
220

221
  # Remove +-Inf at the limits of the curve
222
  #local.maximas <- local.maximas & is.finite(thresholds)
223
  # Question: why did I do that? It breaks coords.roc when partial.auc contains only the extreme point
224

225 8
  return(local.maximas)
226
}
227

228
# Define which progress bar to use
229
roc.utils.get.progress.bar <- function(name = getOption("pROCProgress")$name, title = "Bootstrap", label = "", width = getOption("pROCProgress")$width, char = getOption("pROCProgress")$char, style = getOption("pROCProgress")$style, ...) {
230 8
  if (name == "tk") { # load tcltk if possible
231 0
    if (!requireNamespace("tcltk")) {
232
      # If tcltk cannot be loaded fall back to default text progress bar
233 0
      name <- "text"
234 0
      style <- 3
235 0
      char <- "="
236 0
      width <- NA
237 0
      warning("Package tcltk required with progress='tk' but could not be loaded. Falling back to text progress bar.")
238
    }
239
  }
240 8
  if (name == "none")
241 8
    progress_none()
242 0
  else if (name == "text") {
243
    # Put some default values if user only passed a name
244 0
    if (missing(style) && missing(char) && missing(width) && getOption("pROCProgress")$name != "text") {
245 0
      style <- 3
246 0
      char <- "="
247 0
      width <- NA
248
    }
249 0
    progress_text(char=char, style=style, width=width)
250
  }
251 0
  else if (name == "tk" || name == "win")
252 0
    match.fun(paste("progress", name, sep = "_"))(title=title, label=label, width=width)
253
  else # in the special case someone made a progress_other function
254 0
    match.fun(paste("progress", name, sep = "_"))(title=title, label=label, width=width, char=char, style=style)
255
}
256

257
# sort roc curve. Make sure specificities are increasing.
258
sort.roc <- function(roc) {
259 8
  if (is.unsorted(roc$specificities)) {
260 8
    roc$sensitivities <- rev(roc$sensitivities)
261 8
    roc$specificities <- rev(roc$specificities)
262 8
    roc$thresholds <- rev(roc$thresholds)
263
  }
264 8
  return(roc)
265
}
266

267
# sort smoothed roc curve. Make sure specificities are increasing.
268
sort.smooth.roc <- function(roc) {
269 8
  if (is.unsorted(roc$specificities)) {
270 8
    roc$sensitivities <- rev(roc$sensitivities)
271 8
    roc$specificities <- rev(roc$specificities)
272
  }
273 8
  return(roc)
274
}
275

276
# The list of valid coordinate arguments, without 'thresholds'
277
roc.utils.valid.coords <- c("specificity", "sensitivity", "accuracy",
278
	"tn", "tp", "fn", "fp",
279
	"npv", "ppv", "fdr",
280
	"fpr", "tpr", "tnr", "fnr", 
281
	"1-specificity", "1-sensitivity", "1-accuracy", "1-npv", "1-ppv",
282
	"precision", "recall",
283
	"youden", "closest.topleft")
284

285
# Arguments which can be returned by coords
286
# @param threshold: FALSE for smooth.roc where threshold isn't valid
287
roc.utils.match.coords.ret.args <- function(x, threshold = TRUE) {
288 8
	valid.ret.args <- roc.utils.valid.coords
289 8
	if (threshold) {
290 8
		valid.ret.args <- c("threshold", valid.ret.args)
291
	}
292

293 8
	if ("all" %in% x) {
294 8
		if (length(x) > 1) {
295 8
			stop("ret='all' can't be used with other 'ret' options.")
296
		}
297 8
		x <- valid.ret.args
298
	}
299 8
	x <- replace(x, x=="topleft", "closest.topleft")
300 8
	x <- replace(x, x=="t", "threshold")
301 8
	x <- replace(x, x=="npe", "1-npv")
302 8
	x <- replace(x, x=="ppe", "1-ppv")
303 8
	return(match.arg(x, valid.ret.args, several.ok=TRUE))
304
}
305

306
# Arguments which can be used as input for coords
307
# @param threshold: FALSE for smooth.roc where threshold isn't valid
308
roc.utils.match.coords.input.args <- function(x, threshold = TRUE) {
309 8
	valid.args <- roc.utils.valid.coords
310 8
	if (threshold) {
311 8
		valid.args <- c("threshold", valid.args)
312
	}
313 8
	x <- replace(x, x=="topleft", "closest.topleft")
314 8
	x <- replace(x, x=="t", "threshold")
315 8
	x <- replace(x, x=="npe", "1-npv")
316 8
	x <- replace(x, x=="ppe", "1-ppv")
317 8
	matched <- match.arg(x, valid.args, several.ok=FALSE)
318
	# We only handle monotone coords
319 8
	if (! coord.is.monotone[matched]) {
320 8
		stop(sprintf("Coordinate '%s' is not monotone and cannot be used as input.", matched))
321
	}
322 8
	return(matched)
323
}
324

325
# Compute the min/max for partial AUC
326
# ... with an auc
327
roc.utils.min.partial.auc.auc <- function(auc) {
328 8
  roc.utils.min.partial.auc(attr(auc, "partial.auc"), attr(auc, "percent"))
329
}
330

331
roc.utils.max.partial.auc.auc <- function(roc) {
332 0
  roc.utils.max.partial.auc(attr(auc, "partial.auc"), attr(auc, "percent"))
333
}
334

335
# ... with partial.auc/percent
336
roc.utils.min.partial.auc <- function(partial.auc, percent) {
337 8
  if (!identical(partial.auc, FALSE)) {
338 8
    min <- sum(ifelse(percent, 100, 1)-partial.auc)*abs(diff(partial.auc))/2/ifelse(percent, 100, 1)
339
  }
340
  else {
341 8
    min <- 0.5 * ifelse(percent, 100, 1)
342
  }
343 8
  return(min)
344
}
345

346
roc.utils.max.partial.auc <- function(partial.auc, percent) {
347 8
  if (!identical(partial.auc, FALSE)) {
348 8
    max <- abs(diff(partial.auc))
349
  }
350
  else {
351 0
    max <- 1 * ifelse(percent, 100, 1)
352
  }
353 8
  return(max)
354
}
355

356
# Checks if the 
357
# Input: roc object
358
# Output: boolean, true the curve reaches 100%/100%, false otherwise
359
roc.utils.is.perfect.curve <- function(roc) {
360 8
	best.point <- max(roc$sensitivities + roc$specificities) / ifelse(roc$percent, 100, 1)
361 8
	return(abs(best.point - 2) < .Machine$double.eps ^ 0.5) # or best.point == 2, with numerical tolerance
362
}
363

364
# Load package namespace 'pkg'. 
365
# Input: package name
366
# Returns: TRUE upon success (or if the package was already loaded)
367
# Stops if the package can't be loaded
368
load.suggested.package <- function(pkg) {
369 8
	if (requireNamespace(pkg)) {
370 8
		return(TRUE)
371
	}
372 8
	else if (interactive()) {
373 0
		if (getRversion() < "3.5.0") { # utils::askYesNo not available
374 0
			message(sprintf("Package %s not available, do you want to install it now?", pkg))
375 0
			auto.install <- utils::menu(c("Yes", "No")) == 1
376
		}
377
		else {
378 0
			auto.install <- utils::askYesNo(sprintf("Package %s not available, do you want to install it now?", pkg))
379
		}
380 0
		if (isTRUE(auto.install)) {
381 0
			utils::install.packages(pkg)
382 0
			if (requireNamespace(pkg)) {
383 0
				return(TRUE)
384
			}
385
			else {
386 0
				stop(sprintf("Installation of %s failed!", pkg))
387
			}
388
		}
389
	}
390 0
	stop(sprintf("Package '%s' not available.", pkg))
391
}
392

393

394
# Calculate coordinates
395
# @param roc: the roc curve, used to guess if data is in percent and number of cases and controls.
396
# @param thr, se, sp
397
# @param best.weights: see coords 
398
# @return data.frame
399
roc.utils.calc.coords <- function(roc, thr, se, sp, best.weights) {
400 8
	ncases <- ifelse(methods::is(roc, "smooth.roc"), length(attr(roc, "roc")$cases), length(roc$cases))
401 8
	ncontrols <- ifelse(methods::is(roc, "smooth.roc"), length(attr(roc, "roc")$controls), length(roc$controls))
402 8
	substr.percent <- ifelse(roc$percent, 100, 1)
403
	
404 8
	tp <- se * ncases / substr.percent
405 8
	fn <- ncases - tp
406 8
	tn <- sp * ncontrols / substr.percent
407 8
	fp <- ncontrols - tn
408 8
	npv <- substr.percent * tn / (tn + fn)
409 8
	ppv <- substr.percent * tp / (tp + fp)
410
	#res <- matrix(NA, nrow = length(ret), ncol = length(se))
411
	#if ("tp" %in% ret) {}
412 8
	accuracy <- substr.percent * (tp + tn) / (ncases + ncontrols)
413 8
	precision <- ppv
414 8
	recall <- tpr <- se
415 8
	fpr <- substr.percent - sp
416 8
	tnr <- sp
417 8
	fnr <- substr.percent * fn / (tp + fn)
418 8
	fdr <- substr.percent * fp / (tp + fp)
419 8
	youden <- roc.utils.optim.crit(se, sp, substr.percent, best.weights, "youden")
420 8
	closest.topleft <- - roc.utils.optim.crit(se, sp, substr.percent, best.weights, "closest.topleft") / substr.percent
421
	
422 8
	return(cbind(
423 8
		threshold=thr,
424 8
		sensitivity=se, 
425 8
		specificity=sp, 
426 8
		accuracy=accuracy, 
427 8
		tn=tn, 
428 8
		tp=tp,
429 8
		fn=fn,
430 8
		fp=fp,
431 8
		npv=npv,
432 8
		ppv=ppv,
433 8
		tpr=tpr,
434 8
		tnr=tnr,
435 8
		fpr=fpr,
436 8
		fnr=fnr,
437 8
		fdr=fdr,
438 8
		"1-specificity"=substr.percent-sp,
439 8
		"1-sensitivity"=substr.percent-se,
440 8
		"1-accuracy"=substr.percent-accuracy,
441 8
		"1-npv"=substr.percent-npv,
442 8
		"1-ppv"=substr.percent-ppv,
443 8
		precision=precision,
444 8
		recall=recall,
445 8
		youden=youden,
446 8
		closest.topleft=closest.topleft
447
	))
448
}
449

450
# Match arbitrary user-supplied thresholds to the threshold of the ROC curve.
451
# We need to be careful to assign x to the right thresholds around exact data point
452
# values. This means this function cannot look at the ROC thresholds themselves
453
# but must instead use the predictor values to assess the thresholds exactly.
454
# Returns the indices of the thresholds x.
455
# @param roc: the roc curve
456
# @param x: the threshold to determine indices
457
# @return integer vector of indices along roc$thresholds/roc$se/roc$sp.
458
roc.utils.thr.idx <- function(roc, x) {
459 8
	cut_points <- sort(unique(roc$predictor))
460 8
	thr_idx <- rep(NA_integer_, length(x))
461 8
	if (roc$direction == "<") {
462 8
		cut_points <- c(cut_points, Inf)
463 8
		j <- 1
464 8
		o <- order(x)
465 8
		for (i in seq_along(x)) {
466 8
			t <- x[o[i]]
467 8
			while (cut_points[j] < t) {
468 8
				j <- j + 1
469
			}
470 8
			thr_idx[o[i]] <- j
471
		}
472
	}
473
	else {
474 8
		cut_points <- c(rev(cut_points), -Inf)
475 8
		j <- 1
476 8
		o <- order(x, decreasing = TRUE)
477 8
		for (i in seq_along(x)) {
478 8
			t <- x[o[i]]
479 8
			while (cut_points[j] > t) {
480 8
				j <- j + 1
481
			}
482 8
			thr_idx[o[i]] <- j
483
		}
484
	}
485 8
	return(thr_idx)
486
}
487

488

489
# Get optimal criteria Youden or Closest Topleft
490
# @param se, sp: the roc curve
491
# @param max: the maximum value, 1 or 100, based on percent. Namely ifelse(percent, 100, 1)
492
# @param weights: see coords(best.weights=)
493
# @param method: youden or closest.topleft coords(best.method=)
494
# @return numeric vector along roc$thresholds/roc$se/roc$sp.
495
roc.utils.optim.crit <- function(se, sp, max, weights, method) {
496 8
	if (is.numeric(weights) && length(weights) == 2) {
497 8
		r <- (1 - weights[2]) / (weights[1] * weights[2]) # r should be 1 by default
498
	}
499
	else {
500 8
		stop("'best.weights' must be a numeric vector of length 2")
501
	}
502 8
	if (weights[2] <= 0 || weights[2] >= 1) {
503 8
		stop("prevalence ('best.weights[2]') must be in the interval ]0,1[.")
504
	}
505
	
506 8
	if (method == "youden") {
507 8
		optim.crit <- se + r * sp
508
	}
509 8
	else if (method == "closest.topleft" || method == "topleft") {
510 8
		optim.crit <- - ((max - se)^2 + r * (max - sp)^2)
511
	}
512 8
	return(optim.crit)
513
}
514

515
coord.is.monotone <- c(
516
	"threshold"=TRUE,
517
	"sensitivity"=TRUE,
518
	"specificity"=TRUE,
519
	"accuracy"=FALSE,
520
	"tn"=TRUE, 
521
	"tp"=TRUE,
522
	"fn"=TRUE,
523
	"fp"=TRUE,       
524
	"npv"=FALSE,
525
	"ppv"=FALSE,
526
	"tpr"=TRUE,
527
	"tnr"=TRUE,
528
	"fpr"=TRUE,
529
	"fnr"=TRUE,
530
	"fdr"=FALSE,
531
	"1-specificity"=TRUE,
532
	"1-sensitivity"=TRUE,
533
	"1-accuracy"=FALSE,
534
	"1-npv"=FALSE,
535
	"1-ppv"=FALSE,
536
	"precision"=FALSE,
537
	"recall"=TRUE,
538
	"youden"=FALSE,
539
	"closest.topleft"=FALSE
540
)
541

542
coord.is.decreasing <- c(
543
	"threshold"=NA, # Depends on direction
544
	"sensitivity"=TRUE,
545
	"specificity"=FALSE,
546
	"accuracy"=NA,
547
	"tn"=FALSE, 
548
	"tp"=TRUE,
549
	"fn"=FALSE,
550
	"fp"=TRUE,       
551
	"npv"=NA,
552
	"ppv"=NA,
553
	"tpr"=TRUE,
554
	"tnr"=FALSE,
555
	"fpr"=TRUE,
556
	"fnr"=FALSE,
557
	"fdr"=NA,
558
	"1-specificity"=TRUE,
559
	"1-sensitivity"=FALSE,
560
	"1-accuracy"=NA,
561
	"1-npv"=NA,
562
	"1-ppv"=NA,
563
	"precision"=NA,
564
	"recall"=TRUE,
565
	"youden"=NA,
566
	"closest.topleft"=NA
567
)
568

569
# Get response and predictor(s) from a formula.
570
# This function takes care of all the logic to handle
571
# weights, subset, na.action etc. It handles formulas with 
572
# and without data. It rejects weights and certain na.actions.
573
# @param formula
574
# @param data
575
# @param data.missing 
576
# @param call the call from the parent function
577
# @param ... the ... from the parent function
578
# @return a list with 3 elements: response (vector), predictor.names (character),
579
#         predictors (data.frame).
580
roc.utils.extract.formula <- function(formula, data, data.missing, call, ...) {
581
	# Get predictors (easy)
582 8
	if (data.missing) {
583 8
		predictors <- attr(terms(formula), "term.labels")
584
	}
585
	else {
586 8
		predictors <- attr(terms(formula, data = data), "term.labels")
587
	}
588
	
589 8
	indx <- match(c("formula", "data", "weights", "subset", "na.action"), names(call), nomatch=0)
590 8
	if (indx[1] == 0) {
591 0
		stop("A formula argument is required")
592
	}
593
	# Keep the standard arguments and run them in model.frame
594 8
	temp <- call[c(1,indx)]  
595 8
	temp[[1]] <- as.name("model.frame")
596
	# Only na.pass and na.fail should be used
597 8
	if (indx[5] != 0) {
598 8
		na.action.value = as.character(call[indx[5]])
599 8
		if (! na.action.value %in% c("na.pass", "na.fail")) {
600 8
			warning(paste0(sprintf("Value %s of na.action is not supported ", na.action.value),
601 8
						   "and will break pairing in roc.test and are.paired. ",
602 8
						   "Please use 'na.rm = TRUE' instead."))
603
		}
604
	}
605
	else {
606 8
		temp$na.action = "na.pass"
607
	}
608
	# Adjust call with data from caller
609 8
	if (data.missing) {
610 8
		temp$data <- NULL
611
	}
612
	
613
	# Run model.frame in the parent
614 8
	m <- eval.parent(temp, n = 2)
615
	
616 8
	if (!is.null(model.weights(m))) stop("weights are not supported")
617
	
618 8
	return(list(response.name = names(m)[1],
619 8
				response = model.response(m),
620 8
				predictor.names = predictors,
621 8
				predictors = m[predictors]))
622
}

Read our documentation on viewing source code .

Loading