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

24
ci.thresholds.formula <- function(formula, data, ...) {
25 33
	data.missing <- missing(data)
26 33
	roc.data <- roc.utils.extract.formula(formula, data, ..., 
27 33
										  data.missing = data.missing,
28 33
										  call = match.call())
29 33
	if (length(roc.data$predictor.name) > 1) {
30 0
		stop("Only one predictor supported in 'ci.thresholds'.")
31
	}
32 33
	response <- roc.data$response
33 33
	predictor <- roc.data$predictors[, 1]
34 33
	ci.thresholds(roc(response, predictor, ci=FALSE, ...), ...)
35
}
36

37
ci.thresholds.default <- function(response, predictor, ...) {
38 33
	if (methods::is(response, "multiclass.roc") || methods::is(response, "multiclass.auc")) {
39 33
		stop("'ci.thresholds' not available for multiclass ROC curves.")
40
	}
41 33
	ci.thresholds(roc.default(response, predictor, ci=FALSE, ...), ...)
42
}
43

44 0
ci.thresholds.smooth.roc <- function(smooth.roc, ...)
45 0
  stop("'ci.thresholds' is not available for smoothed ROC curves.")
46

47
ci.thresholds.roc <- function(roc,
48
                   conf.level = 0.95,
49
                   boot.n = 2000,
50
                   boot.stratified = TRUE,
51
                   thresholds = "local maximas",
52
                   progress = getOption("pROCProgress")$name,
53
                   parallel = FALSE,
54
                   ...
55
                   ) {
56 33
  if (conf.level > 1 | conf.level < 0)
57 0
    stop("'conf.level' must be within the interval [0,1].")
58
  
59 33
  if (roc.utils.is.perfect.curve(roc)) {
60 0
  	warning("ci.thresholds() of a ROC curve with AUC == 1 is always a null interval and can be misleading.")
61
  }
62

63
  # Check and prepare thresholds
64 33
  if (is.character(thresholds)) {
65 33
    if (length(thresholds) != 1)
66 0
      stop("'thresholds' of class character must be of length 1.")
67 33
    thresholds <- match.arg(thresholds, c("all", "best", "local maximas"))
68 33
    thresholds.num <- coords(roc, x=thresholds, input="threshold", ret="threshold", as.matrix = TRUE, transpose = FALSE, ...)[, 1]
69 33
    attr(thresholds.num, "coords") <- thresholds
70
  }
71 33
  else if (is.logical(thresholds)) {
72 0
    thresholds.num <- roc$thresholds[thresholds]
73 0
    attr(thresholds.num, "logical") <- thresholds
74
  }
75 33
  else if (! is.numeric(thresholds)) {
76 0
    stop("'thresholds' is not character, logical or numeric.")
77
  }
78
  else {
79 33
    thresholds.num <- thresholds
80
  }
81

82 33
  if(class(progress) != "list")
83 33
    progress <- roc.utils.get.progress.bar(progress, title="Thresholds confidence interval", label="Bootstrap in progress...", ...)
84

85 33
  if (boot.stratified) {
86 33
    perfs <- laply(1:boot.n, stratified.ci.thresholds, roc=roc, thresholds=thresholds.num, .progress=progress, .parallel=parallel)
87
  }
88
  else {
89 33
    perfs <- laply(1:boot.n, nonstratified.ci.thresholds, roc=roc, thresholds=thresholds.num, .progress=progress, .parallel=parallel)
90
  }
91

92 33
  if (length(thresholds.num) > 1) {
93 33
    if (any(is.na(perfs))) {
94 0
      warning("NA value(s) produced during bootstrap were ignored.")
95 0
      perfs <- perfs[!apply(perfs, 1, function(x) any(is.na(x))),]
96
    }
97
    # laply returns a 3d matrix, with dim 1 = bootstrap replicates, dim 2 = SE/SP and dim 3 = thresholds
98
    # [,1,] = SP and [,2,] = SE
99 33
    sp <- t(apply(perfs[,1,], 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
100 33
    se <- t(apply(perfs[,2,], 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
101
  }
102
  else {
103 33
    if (any(is.na(perfs))) {
104 0
      warning("NaN value(s) in bootstrap ignored in confidence interval.")
105 0
      perfs <- perfs[!apply(perfs, 1, function(x) any(is.na(x))),]
106
    }
107 33
    sp <- as.matrix(t(quantile(perfs[,1], probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2))))
108 33
    se <- as.matrix(t(quantile(perfs[,2], probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2))))
109
  }
110

111 33
  rownames(se) <- rownames(sp) <- thresholds.num
112

113 33
  if (roc$percent) {
114 33
    se <- se * 100
115 33
    sp <- sp * 100
116
  }
117

118 33
  ci <- list(specificity = sp, sensitivity = se)
119 33
  class(ci) <- c("ci.thresholds", "ci", "list")
120 33
  attr(ci, "conf.level") <- conf.level
121 33
  attr(ci, "boot.n") <- boot.n
122 33
  attr(ci, "boot.stratified") <- boot.stratified
123 33
  attr(ci, "thresholds") <- thresholds.num
124 33
  attr(ci, "roc") <- roc
125 33
  return(ci)
126
}

Read our documentation on viewing source code .

Loading