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 1
multiclass.roc <- function(...)
21 1
  UseMethod("multiclass.roc")
22

23
multiclass.roc.formula <- function(formula, data, ...) {
24 1
	data.missing <- missing(data)
25 1
	call <- match.call()
26 1
	roc.data <- roc.utils.extract.formula(formula, data, ..., 
27 1
										  data.missing = data.missing,
28 1
										  call = call)
29 1
	response <- roc.data$response
30 1
	predictors <- roc.data$predictors
31 1
	if (ncol(predictors) == 1) {
32 1
		predictors <- predictors[, 1]
33
	}
34 1
	multiclass.roc <- multiclass.roc.default(response, predictors, ...)
35 1
	multiclass.roc$call <- call
36 1
	if (! data.missing) {
37 1
		multiclass.roc$data <- data
38
	}
39 1
	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 1
  multiclass.roc <- list(
51 1
                         response = response,
52 1
                         predictor = predictor,
53 1
                         percent = percent)
54 1
  class(multiclass.roc) <- "multiclass.roc"
55 1
  if (is.factor(response) && any(names(table(response))[table(response) == 0] %in% levels)) {
56 1
    missing.levels <- names(table(response))[table(response) == 0]
57 1
    missing.levels.requested <- missing.levels[missing.levels %in% levels]
58 1
    warning(paste("No observation for response level(s):", paste(missing.levels.requested, collapse=", ")))
59 1
    levels <- levels[!(levels %in% missing.levels.requested)]
60
  }
61 1
  multiclass.roc$levels <- levels
62
  
63 1
  rocs <- utils::combn(levels, 2, function(X, response, predictor, percent, ...) {
64 1
    roc(response, predictor, levels=X, percent=percent, auc=FALSE, ci=FALSE, ...)
65 1
  }, simplify=FALSE, response=response, predictor=predictor, percent=percent, direction=direction, ...)
66

67 1
  multiclass.roc$rocs <- rocs
68

69
  # Makes no sense to turn auc off, so remove this option
70
  #if (auc)
71 1
    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 1
  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 1
    pred.i <- pred.matrix[which(ref.outcome == i), i] # p(G = i) assigned to class i observations
84 1
    pred.j <- pred.matrix[which(ref.outcome == j), i] # p(G = i) assigned to class j observations
85 1
    classes <- factor(c(rep(i, length(pred.i)), rep(j, length(pred.j))))
86
    # override levels argument by new levels
87 1
    levels <- unique(classes)
88 1
    predictor <- c(pred.i, pred.j)
89 1
    auc <- roc(classes, predictor, levels = levels, percent = percent, auc = FALSE, ci = FALSE, direction = direction, ...)
90 1
    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 1
    if (!methods::is(predictor, "matrix") && !methods::is(predictor, "data.frame")) {
97 0
        stop("Please provide a matrix or data frame via 'predictor'.")
98
    }
99 1
    if (nrow(predictor) != length(response)) {
100 1
        stop("Number of rows in 'predictor' does not agree with 'response'");
101
    }
102 1
	if (direction == "auto") {
103 1
		stop("'direction=\"auto\"' not available for multivariate multiclass.roc")
104
	}
105 1
	if (is.factor(response) && any(names(table(response))[table(response) == 0] %in% levels)) {
106 1
		missing.levels <- names(table(response))[table(response) == 0]
107 1
		missing.levels.requested <- missing.levels[missing.levels %in% levels]
108 1
		warning(paste("No observation for response level(s):", paste(missing.levels.requested, collapse=", ")))
109 1
		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 1
    m <- match(colnames(predictor), levels)
114 1
    missing.classes <- levels[setdiff(seq_along(levels), m)]
115 1
    levels <- colnames(predictor)[!is.na(m)]
116 1
    if (length(levels) == 1) {
117 1
        stop("For a single decision value, please provide 'predictor' as a vector.")
118 1
    } else if (length(levels) == 0) {
119 1
        stop("The column names of 'predictor' could not be matched to the levels of 'response'.")
120
    }
121 1
    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 1
    additional.classes <- colnames(predictor)[which(is.na(m))]
133 1
    if (length(additional.classes) != 0) {
134 1
        out.classes <- paste0(additional.classes, collapse = ",")
135 1
        warning("The following classes were not found in 'response': ", out.classes, ".")
136
    }
137 1
    multiclass.roc <- list(
138 1
                         response = response,
139 1
                         predictor = predictor,
140 1
                         percent = percent)
141 1
    class(multiclass.roc) <- "mv.multiclass.roc"
142 1
    multiclass.roc$levels <- levels
143 1
    rocs <- utils::combn(levels, 2, function(x, predictor, response, levels, percent, direction, ...) {
144 1
        A1 <- compute.pair.AUC(predictor, x[1], x[2], response, levels, percent, direction, ...)
145 1
        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 1
        A <- list(A1, A2) 
148 1
        return(A)
149 1
    }, simplify = FALSE, predictor = predictor, response = response, 
150 1
        levels = levels, percent = percent, direction, ...)
151 1
    pairs <- unlist(lapply(utils::combn(levels, 2, simplify = FALSE), 
152 1
                           function(x) paste(x, collapse = "/")))
153 1
    names(rocs) <- pairs
154 1
    multiclass.roc$rocs <- rocs
155
    
156 1
    multiclass.roc$auc <- auc.mv.multiclass.roc(multiclass.roc, ...)
157
    
158 1
    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 1
	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 1
    if (methods::is(predictor, "matrix") || methods::is(predictor, "data.frame")) {
173
        # for decision values for multiple classes (e.g. probabilities of individual classes)
174 1
    	if (missing("direction")) {
175
            # need to have uni-directional decision values for consistency
176 1
            direction <- ">" 
177
    	}
178
    	else {
179 1
    		direction <- match.arg(direction)
180
    	}
181 1
        mc.roc <- multiclass.roc.multivariate(response, predictor, levels, percent, direction, ...)
182
    } else {
183
        # for a single decision value for separating the classes
184 1
    	direction <- match.arg(direction)
185 1
    	mc.roc <- multiclass.roc.univariate(response, predictor, levels, percent, direction, ...)
186
    }
187 1
	mc.roc$call <- match.call()
188 1
	return(mc.roc)
189
}

Read our documentation on viewing source code .

Loading