1
# pROC: Tools Receiver operating characteristic (ROC curves) with
2
# (partial) area under the curve, confidence intervals and comparison. 
3
# Copyright (C) 2010-2019 Xavier Robin, Matthias Doering, 
4
# Alexandre Hainard, Natacha Turck, Natalia Tiberti, 
5
# Frédérique Lisacek, Jean-Charles Sanchez 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 3
multiclass.roc <- function(...)
21 3
  UseMethod("multiclass.roc")
22

23
multiclass.roc.formula <- function(formula, data, ...) {
24 3
	data.missing <- missing(data)
25 3
	call <- match.call()
26 3
	roc.data <- roc.utils.extract.formula(formula, data, ..., 
27 3
										  data.missing = data.missing,
28 3
										  call = call)
29 3
	response <- roc.data$response
30 3
	predictors <- roc.data$predictors
31 3
	if (ncol(predictors) == 1) {
32 3
		predictors <- predictors[, 1]
33
	}
34 3
	multiclass.roc <- multiclass.roc.default(response, predictors, ...)
35 3
	multiclass.roc$call <- call
36 3
	if (! data.missing) {
37 3
		multiclass.roc$data <- data
38
	}
39 3
	return(multiclass.roc)
40
}
41

42
multiclass.roc.univariate <- function(response, predictor,
43
                                   levels=base::levels(as.factor(response)),
44
                                   percent=FALSE, # Must sensitivities, specificities and AUC be reported in percent? Note that if TRUE, and you want a partial area, you must pass it in percent also (partial.area=c(100, 80))
45
                                   direction,
46
                                   # what computation must be done
47
                                   #auc=TRUE, # call auc.roc on the current object
48
                                   #ci=FALSE, # call ci.roc on the current object
49
                                   ...) {
50 3
  multiclass.roc <- list(
51 3
                         response = response,
52 3
                         predictor = predictor,
53 3
                         percent = percent)
54 3
  class(multiclass.roc) <- "multiclass.roc"
55 3
  if (is.factor(response) && any(names(table(response))[table(response) == 0] %in% levels)) {
56 3
    missing.levels <- names(table(response))[table(response) == 0]
57 3
    missing.levels.requested <- missing.levels[missing.levels %in% levels]
58 3
    warning(paste("No observation for response level(s):", paste(missing.levels.requested, collapse=", ")))
59 3
    levels <- levels[!(levels %in% missing.levels.requested)]
60
  }
61 3
  multiclass.roc$levels <- levels
62
  
63 3
  rocs <- utils::combn(levels, 2, function(X, response, predictor, percent, ...) {
64 3
    roc(response, predictor, levels=X, percent=percent, auc=FALSE, ci=FALSE, ...)
65 3
  }, simplify=FALSE, response=response, predictor=predictor, percent=percent, direction=direction, ...)
66

67 3
  multiclass.roc$rocs <- rocs
68

69
  # Makes no sense to turn auc off, so remove this option
70
  #if (auc)
71 3
    multiclass.roc$auc <- auc.multiclass.roc(multiclass.roc, ...)
72
  # CI is not implemented yet.
73
  #if (ci)
74
  #  multiclass.roc$ci <- ci.multiclass.roc(multiclass.roc, ...)
75

76 3
  return(multiclass.roc)
77
}
78

79
compute.pair.AUC <- function(pred.matrix, i, j, ref.outcome, levels, percent, direction, ... ) {
80
    # computes A(i|j), the probability that a randomly 
81
    # chosen member of class j has a lower estimated probability (or score) 
82
    # of belonging to class i than a randomly chosen member of class i
83 3
    pred.i <- pred.matrix[which(ref.outcome == i), i] # p(G = i) assigned to class i observations
84 3
    pred.j <- pred.matrix[which(ref.outcome == j), i] # p(G = i) assigned to class j observations
85 3
    classes <- factor(c(rep(i, length(pred.i)), rep(j, length(pred.j))))
86
    # override levels argument by new levels
87 3
    levels <- unique(classes)
88 3
    predictor <- c(pred.i, pred.j)
89 3
    auc <- roc(classes, predictor, levels = levels, percent = percent, auc = FALSE, ci = FALSE, direction = direction, ...)
90 3
    return(auc)
91
}
92

93
multiclass.roc.multivariate <- function(response, predictor, levels, percent, direction, ...) {
94
    # Reference: "A Simple Generalisation of the Area Under the ROC 
95
    # Curve for Multiple Class Classification Problems" (Hand and Till, 2001)
96 3
    if (!methods::is(predictor, "matrix") && !methods::is(predictor, "data.frame")) {
97 0
        stop("Please provide a matrix or data frame via 'predictor'.")
98
    }
99 3
    if (nrow(predictor) != length(response)) {
100 3
        stop("Number of rows in 'predictor' does not agree with 'response'");
101
    }
102 3
	if (direction == "auto") {
103 3
		stop("'direction=\"auto\"' not available for multivariate multiclass.roc")
104
	}
105 3
	if (is.factor(response) && any(names(table(response))[table(response) == 0] %in% levels)) {
106 3
		missing.levels <- names(table(response))[table(response) == 0]
107 3
		missing.levels.requested <- missing.levels[missing.levels %in% levels]
108 3
		warning(paste("No observation for response level(s):", paste(missing.levels.requested, collapse=", ")))
109 3
		levels <- levels[!(levels %in% missing.levels.requested)]
110
	}
111
	
112
    # check whether the columns of the prediction matrix agree with the factors in 'response'
113 3
    m <- match(colnames(predictor), levels)
114 3
    missing.classes <- levels[setdiff(seq_along(levels), m)]
115 3
    levels <- colnames(predictor)[!is.na(m)]
116 3
    if (length(levels) == 1) {
117 3
        stop("For a single decision value, please provide 'predictor' as a vector.")
118 3
    } else if (length(levels) == 0) {
119 3
        stop("The column names of 'predictor' could not be matched to the levels of 'response'.")
120
    }
121 3
    if (length(missing.classes) != 0) {
122 0
        out.classes <- paste0(missing.classes, collapse = ",")
123 0
        if (length(missing.classes) == length(levels)) {
124
            # no decision values found
125 0
            stop(paste0("Could not find any decision values in 'predictor' matching the 'response' levels.",
126 0
                 " Could not find the following classes: ", out.classes, ". Check your column names!"))
127
        } else {
128
            # some decision values not found
129 0
            warning("You did not provide decision values for the following classes: ", out.classes, ".")
130
        }
131
    }
132 3
    additional.classes <- colnames(predictor)[which(is.na(m))]
133 3
    if (length(additional.classes) != 0) {
134 3
        out.classes <- paste0(additional.classes, collapse = ",")
135 3
        warning("The following classes were not found in 'response': ", out.classes, ".")
136
    }
137 3
    multiclass.roc <- list(
138 3
                         response = response,
139 3
                         predictor = predictor,
140 3
                         percent = percent)
141 3
    class(multiclass.roc) <- "mv.multiclass.roc"
142 3
    multiclass.roc$levels <- levels
143 3
    rocs <- utils::combn(levels, 2, function(x, predictor, response, levels, percent, direction, ...) {
144 3
        A1 <- compute.pair.AUC(predictor, x[1], x[2], response, levels, percent, direction, ...)
145 3
        A2 <- compute.pair.AUC(predictor, x[2], x[1], response, levels, percent, direction, ...)
146
        # merging A1 and A2 is infeasible as auc() would not be well-defined
147 3
        A <- list(A1, A2) 
148 3
        return(A)
149 3
    }, simplify = FALSE, predictor = predictor, response = response, 
150 3
        levels = levels, percent = percent, direction, ...)
151 3
    pairs <- unlist(lapply(utils::combn(levels, 2, simplify = FALSE), 
152 3
                           function(x) paste(x, collapse = "/")))
153 3
    names(rocs) <- pairs
154 3
    multiclass.roc$rocs <- rocs
155
    
156 3
    multiclass.roc$auc <- auc.mv.multiclass.roc(multiclass.roc, ...)
157
    
158 3
    return(multiclass.roc)
159
}
160

161
multiclass.roc.default <- function(response, predictor,
162
                                   levels = base::levels(as.factor(response)),
163
                                   percent = FALSE, # Must sensitivities, specificities and AUC be reported in percent? Note that if TRUE, and you want a partial area, you must pass it in percent also (partial.area=c(100, 80)),
164
                                   direction = c("auto", "<", ">"),
165
                                   ...) {
166
	# We need at least two levels in response
167 3
	if (length(unique(response)) < 2) {
168 0
		stop("'response' must have at least two levels")
169
	}
170
	
171
    # implements the approach from Hand & Till (2001)
172 3
    if (methods::is(predictor, "matrix") || methods::is(predictor, "data.frame")) {
173
        # for decision values for multiple classes (e.g. probabilities of individual classes)
174 3
    	if (missing("direction")) {
175
            # need to have uni-directional decision values for consistency
176 3
            direction <- ">" 
177
    	}
178
    	else {
179 3
    		direction <- match.arg(direction)
180
    	}
181 3
        mc.roc <- multiclass.roc.multivariate(response, predictor, levels, percent, direction, ...)
182
    } else {
183
        # for a single decision value for separating the classes
184 3
    	direction <- match.arg(direction)
185 3
    	mc.roc <- multiclass.roc.univariate(response, predictor, levels, percent, direction, ...)
186
    }
187 3
	mc.roc$call <- match.call()
188 3
	return(mc.roc)
189
}

Read our documentation on viewing source code .

Loading