Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement the approach from Hand & Till for generalization the AUC to the multi-class setting #37

Merged
merged 5 commits into from
Jan 27, 2019
Merged

Implement the approach from Hand & Till for generalization the AUC to the multi-class setting #37

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1107c09
implement the approach from Hand & Till for generalization the AUC to…
Dec 5, 2018
bf899b6
use of roc function from pROC for multivariate calculations; remove r…
Dec 5, 2018
9d5f082
fix inconsistent AUC evaluation for multi-class by disallowing direct…
Dec 6, 2018
e7c9998
add documentation for direction setting
Dec 12, 2018
ffd5345
add documentation for direction setting
Dec 12, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 123 additions & 5 deletions R/multiclass.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,19 +30,22 @@ multiclass.roc.formula <- function(formula, data, ...) {
m <- eval(m, parent.frame())
term.labels <- attr(attr(m, "terms"), "term.labels")
response <- model.extract(m, "response")

if (length(response) == 0) {
stop("Error in the formula: a response is required in a formula of type response~predictor.")
}

multiclass.roc <- multiclass.roc.default(response, m[[term.labels]], ...)
predictor <- m[term.labels]
if (ncol(predictor) == 1) {
predictor <- as.vector(unlist(predictor))
}
multiclass.roc <- multiclass.roc.default(response, predictor, ...)
multiclass.roc$call <- match.call()
return(multiclass.roc)
}

multiclass.roc.default <- function(response, predictor,
multiclass.roc.univariate <- function(response, predictor,
levels=base::levels(as.factor(response)),
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))
direction,
# what computation must be done
#auc=TRUE, # call auc.roc on the current object
#ci=FALSE, # call ci.roc on the current object
Expand All @@ -53,7 +56,6 @@ multiclass.roc.default <- function(response, predictor,
percent = percent,
call=match.call())
class(multiclass.roc) <- "multiclass.roc"

if ("ordered" %in% class(response) && any(names(table(response))[table(response) == 0] %in% levels)) {
missing.levels <- names(table(response))[table(response) == 0]
missing.levels.requested <- missing.levels[missing.levels %in% levels]
Expand All @@ -77,3 +79,119 @@ multiclass.roc.default <- function(response, predictor,

return(multiclass.roc)
}

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

compute.new.AUC <- function(pred.matrix, i, j, ref.outcome, levels, percent, ... ) {
# computes A(i|j), the probability that a randomly
# chosen member of class j has a lower estimated probability (or score)
# of belonging to class i than a randomly chosen member of class i

pred.i <- pred.matrix[which(ref.outcome == i), i] # p(G = i) assigned to class i observations
pred.j <- pred.matrix[which(ref.outcome == j), i] # p(G = i) assigned to class j observations
classes <- factor(c(rep(i, length(pred.i)), rep(j, length(pred.j))))
# override levels argument by new levels
levels <- unique(classes)
predictor <- c(pred.i, pred.j)
auc <- roc(classes, predictor, levels = levels, percent = percent, auc = FALSE, ci = FALSE, ...)
return(auc)
}

multiclass.roc.multivariate <- function(response, predictor, levels, percent, direction, ...) {
# Reference: "A Simple Generalisation of the Area Under the ROC
# Curve for Multiple Class Classification Problems" (Hand and Till, 2001)
if (!is(predictor, "matrix") && !is(predictor, "data.frame")) {
stop("Please provide a matrix or data frame via 'predictor'.")
}
if (nrow(predictor) != length(response)) {
stop("Number of rows in 'predictor' does not agree with 'response'");
}
# check whether the columns of the prediction matrix agree with the factors in 'response'
m <- match(colnames(predictor), levels)
missing.classes <- levels[setdiff(seq_along(levels), m)]
levels <- colnames(predictor)[!is.na(m)]
if (length(levels) == 1) {
stop("For a single decision value, please provide 'predictor' as a vector.")
} else if (length(levels) == 0) {
stop("The column names of 'predictor' could not be matched to the levels of 'response'.")
}
if (length(missing.classes) != 0) {
out.classes <- paste0(missing.classes, collapse = ",")
if (length(missing.classes) == length(levels)) {
# no decision values found
stop(paste0("Could not find any decision values in 'predictor' matching the 'response' levels.",
" Could not find the following classes: ", out.classes, ". Check your column names!"))
} else {
# some decision values not found
warning("You did not provide decision values for the following classes: ", out.classes, ".")
}
}
additional.classes <- colnames(predictor)[which(is.na(m))]
if (length(additional.classes) != 0) {
out.classes <- paste0(additional.classes, collapse = ",")
warning("The following classes were not found in 'response': ", out.classes, ".")
}
multiclass.roc <- list(
response = response,
predictor = predictor,
percent = percent,
call=match.call())
class(multiclass.roc) <- "multiclass.roc"
multiclass.roc$levels <- levels
rocs <- utils::combn(levels, 2, function(x, predictor, response, levels, percent, direction, ...) {
A1 <- compute.pair.AUC(predictor, x[1], x[2], response, levels, percent, direction, ...)
A2 <- compute.pair.AUC(predictor, x[2], x[1], response, levels, percent, direction, ...)
# merging A1 and A2 is infeasible as auc() would not be well-defined
A <- list(A1, A2)
return(A)
}, simplify = FALSE, predictor = predictor, response = response,
levels = levels, percent = percent, direction, ...)
pairs <- unlist(lapply(combn(levels, 2, simplify = FALSE),
function(x) paste(x, collapse = "/")))
names(rocs) <- pairs
multiclass.roc$rocs <- rocs
# can't use 'auc.multiclass.roc' because mean of ROCs needs to be computed before ...
#A.total <- auc.multiclass.roc(multiclass.roc)
mean.AUCs <- unlist(lapply(rocs, function(x) mean(auc(x[[1]]), auc(x[[2]]))))
A.ij.total <- sum(mean.AUCs)
c <- length(levels)
M <- 2 / (c * (c-1)) * A.ij.total
multiclass.roc$auc <- M # manually store auc; TODO: improve functioning of auc function
multiclass.roc$pair_aucs <- mean.AUCs # manually store pairwise AUCs
return(multiclass.roc)
}

multiclass.roc.default <- function(response, predictor,
levels = base::levels(as.factor(response)),
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)),
direction = ">",
# what computation must be done
#auc=TRUE, # call auc.roc on the current object
#ci=FALSE, # call ci.roc on the current object
...) {
# implements the approach from Hand & Till (2001)
if (is(predictor, "matrix") || is(predictor, "data.frame")) {
# for decision values for multiple classes (e.g. probabilities of individual classes)
if (direction == "auto") {
# need to have uni-directional decision values for consistency
direction <- ">"
message("direction = 'auto' is not allowed for multivariate input. Setting 'direction' to '>'.")
}
return(multiclass.roc.multivariate(response, predictor, levels, percent, direction, ...))
} else {
# for a single decision value for separating the classes
return(multiclass.roc.univariate(response, predictor, levels, percent, direction, ...))
}
}
78 changes: 74 additions & 4 deletions man/multiclass.Rd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ multiclass.roc(...)
\S3method{multiclass.roc}{formula}(formula, data, ...)
\S3method{multiclass.roc}{default}(response, predictor,
levels=base::levels(as.factor(response)),
percent=FALSE, ...)
percent=FALSE, direction = ">", ...)

}

Expand All @@ -24,8 +24,9 @@ percent=FALSE, ...)
responses, typically encoded with 0 (controls) and 1 (cases), as in
\code{\link{roc}}.
}
\item{predictor}{a numeric vector, containing the value of each
observation, as in \code{\link{roc}}.
\item{predictor}{either a numeric vector, containing the value of each
observation, as in \code{\link{roc}}, or, a matrix giving the decision value
(e.g. probability) for each class.
}
\item{formula}{a formula of the type \code{response~predictor}.}
\item{data}{a matrix or data.frame containing the variables in the
Expand All @@ -38,6 +39,15 @@ percent=FALSE, ...)
\item{percent}{if the sensitivities, specificities and AUC must be
given in percent (\code{TRUE}) or in fraction (\code{FALSE}, default).
}
\item{direction}{in which direction to make the comparison?
\dQuote{>} (default): if the predictor values for the control group are
higher than the values of the case group (controls > t >= cases).
\dQuote{<}: if the predictor values for the control group are lower
or equal than the values of the case group (controls < t <= cases).
Since setting the parameter to \dQuote{auto} would lead to inconsistent
AUC values, a setting of \dQuote{auto} is automatically adjusted to
the default setting.
}
\item{...}{further arguments passed to \code{\link{roc}}.
}
}
Expand All @@ -52,7 +62,7 @@ intervals and comparison tests are not implemented yet.
}

\value{
A list of class \dQuote{multiclass.roc} with the following fields:
If \code{predictor} is a vector, a list of class \dQuote{multiclass.roc} with the following fields:
\item{auc}{if called with \code{auc=TRUE}, a numeric of class \dQuote{auc} as
defined in \code{\link{auc}}. Note that this is not the standard AUC
but the multi-class AUC as defined by Hand and Till.
Expand All @@ -75,6 +85,8 @@ intervals and comparison tests are not implemented yet.
\item{call}{how the function was called. See \code{\link{match.call}} for
more details.
}
If \code{predictor} is a matrix, a numeric giving the overall AUC. The attribute
\code{pair_AUCs} is a vector indicating the pairwise AUCs, A(i,j), for individual pairs of classes (i,j).
}

\section{Warnings}{
Expand All @@ -95,6 +107,9 @@ intervals and comparison tests are not implemented yet.
}

\examples{
####
# Examples for a univariate decision value
####
data(aSAH)

# Basic example
Expand All @@ -107,7 +122,62 @@ multiclass.roc(aSAH$gos6, aSAH$s100b, levels=c(3, 4, 5))
# Give the result in percent
multiclass.roc(aSAH$gos6, aSAH$s100b, percent=TRUE)

####
# Examples for multivariate decision values (e.g. class probabilities)
####
n <- c(100, 80, 150)
responses <- factor(c(rep("X1", n[1]), rep("X2", n[2]), rep("X3", n[3])))
# construct prediction matrix: one column per class
set.seed(42)

# Perfect separation
preds <- lapply(n, function(x) runif(x, 0.8, 1))
predictor <- as.matrix(data.frame("X1" = c(preds[[1]], runif(n[2] + n[3], 0, 0.3)),
"X2" = c(runif(n[1], 0.1, 0.4), preds[[2]], runif(n[3], 0, 0.2)),
"X3" = c(runif(n[1] + n[2], 0, 0.5), preds[[3]])))
multiclass.roc(responses, predictor)

# Imperfect separation

preds <- lapply(n, function(x) runif(x, 0.4, 0.6))
predictor <- as.matrix(data.frame(
"X1" = c(preds[[1]], runif(n[2] + n[3], 0, 0.7)),
"X2" = c(runif(n[1], 0.1, 0.4), preds[[2]], runif(n[3], 0.2, 0.8)),
"X3" = c(runif(n[1] + n[2], 0.3, 0.7), preds[[3]])
))
multiclass.roc(responses, predictor)
multiclass.roc(responses, predictor, direction = ">")

# AUC of random classifier

preds <- lapply(n, function(x) runif(x, 0, 1))
predictor <- as.matrix(data.frame("X1" = c(preds[[1]], runif(n[2] + n[3], 0, 1)),
"X2" = c(runif(n[1], 0, 1), preds[[2]], runif(n[3], 0, 1)),
"X3" = c(runif(n[1] + n[2], 0, 1), preds[[3]])))
multiclass.roc(responses, predictor)

# AUC of worst-case classifier

preds <- lapply(n, function(x) rep(0, x))
predictor <- as.matrix(data.frame("X1" = c(preds[[1]], runif(n[2] + n[3], 0.1, 1)),
"X2" = c(runif(n[1], 0, 1), preds[[2]], runif(n[3], 0.1, 1)),
"X3" = c(runif(n[1] + n[2], 0.1, 1), preds[[3]])))
multiclass.roc(responses, predictor)
multiclass.roc(formula = X3 ~ X2, data = predictor)


# Minimum example for preventing the use of direction = "auto"
responses <- factor(c("X1", "X2", "X3"))
predictor <- data.frame("X1" = c(0.9, 0.0, 0.1),
"X2" = c(0.05, 0.1, 0.7),
"X3" = c(0.05, 0.9, 0.2))
multiclass.roc(responses, predictor)
multiclass.roc(responses, predictor, percent=TRUE)
# Limit set of levels
multiclass.roc(responses, predictor, levels = c("X1", "X2"))
# Use with formula
data <- cbind(predictor, "response" = responses)
multiclass.roc(response ~ X1+X3, data)

}

Expand Down