Skip to content

Commit

Permalink
Merge branch 'master' into parallel
Browse files Browse the repository at this point in the history
  • Loading branch information
xrobin committed Mar 15, 2012
2 parents 8bcefb3 + 89e8ca8 commit deda1b3
Show file tree
Hide file tree
Showing 15 changed files with 312 additions and 30 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Expand Up @@ -4,7 +4,7 @@ Title: display and analyze ROC curves
Version: 1.5.9
Date: 2011-??-??
Encoding: UTF-8
Depends: plyr, R (>= 2.10), utils
Depends: plyr, R (>= 2.10), utils, methods
Suggests: tcltk, MASS, logcondens
Author: Xavier Robin, Natacha Turck, Alexandre Hainard, Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez and Markus Müller
Maintainer: Xavier Robin <Xavier.Robin@unige.ch>
Expand Down
3 changes: 3 additions & 0 deletions NAMESPACE
Expand Up @@ -2,6 +2,7 @@ export(
are.paired, are.paired.auc, are.paired.roc, are.paired.smooth.roc,
auc, auc.formula, auc.roc, auc.smooth.roc, auc.default, auc.multiclass.roc,
ci, ci.formula, ci.default, ci.roc, ci.smooth.roc,
ci.coords, ci.coords.formula, ci.coords.default, ci.coords.roc, ci.coords.smooth.roc,
ci.thresholds, ci.thresholds.formula, ci.thresholds.default, ci.thresholds.roc, ci.thresholds.smooth.roc,
ci.sp, ci.sp.formula, ci.sp.default, ci.sp.roc, ci.sp.smooth.roc,
ci.se, ci.se.formula, ci.se.default, ci.se.roc, ci.se.smooth.roc,
Expand All @@ -20,3 +21,5 @@ export(
smooth, smooth.roc, smooth.default, smooth.smooth.roc,
var, var.default, var.auc, var.smooth.roc, var.roc
)

import(methods, utils)
4 changes: 4 additions & 0 deletions NEWS
@@ -1,8 +1,12 @@
1.6 (201?-??-??)
* New 'power.roc.test' function for sample size and power computations
* New 'cov' and 'var' functions supports new "obuchowski" method
* New 'ci.coords' function to compute CI of arbitrary coords
* 'coords' accepts new 'ret' value "1-accuracy"

1.5.1 (2012-03-09)
* Faster loading of the package (thanks to Prof Brian Ripley and Glenn Lawyer for the report)

1.5 (2011-12-11)
* New 'cov' and 'var' functions
* 'coords' accepts new 'ret' values: "accuracy", "tn", "tp", "fn", "fp", "npv", "ppv", "1-specificity", "1-sensitivity", "1-npv", "1-ppv", "npe" and "ppe"
Expand Down
96 changes: 96 additions & 0 deletions R/bootstrap.R
Expand Up @@ -591,3 +591,99 @@ nonstratified.ci.thresholds <- function(n, roc, thresholds) {

return(sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction))
}


########## Coords of one ROC curves (ci.coords) ##########
stratified.ci.coords <- function(roc, x, input, ret, best.method, best.weights) {
controls <- sample(roc$controls, replace=TRUE)
cases <- sample(roc$cases, replace=TRUE)
thresholds <- roc.utils.thresholds(c(cases, controls))

perfs <- sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction)
# update ROC
roc$sensitivities <- perfs[2,]
roc$specificities <- perfs[1,]
roc$cases <- cases
roc$controls <- controls
roc$predictor <- c(controls, cases)
roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
roc$thresholds <- thresholds

as.numeric(coords.roc(roc, x=x, input=input, ret=ret, best.method=best.method, best.weights=best.weights))
}

nonstratified.ci.coords <- function(roc, x, input, ret, best.method, best.weights) {
tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
predictor <- roc$predictor[tmp.idx]
response <- roc$response[tmp.idx]
splitted <- split(predictor, response)
controls <- splitted[[as.character(roc$levels[1])]]
cases <- splitted[[as.character(roc$levels[2])]]
thresholds <- roc.utils.thresholds(c(controls, cases))

perfs <- sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction)
# update ROC
roc$sensitivities <- perfs[2,]
roc$specificities <- perfs[1,]
roc$cases <- cases
roc$controls <- controls
roc$predictor <- c(controls, cases)
roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
roc$thresholds <- thresholds

as.numeric(coords.roc(roc, x=x, input=input, ret=ret, best.method=best.method, best.weights=best.weights))
}

########## Coords of a smooth ROC curve (ci.coords) ##########

stratified.ci.smooth.coords <- function(roc, x, input, ret, best.method, best.weights, smooth.roc.call) {
controls <- sample(roc$controls, replace=TRUE)
cases <- sample(roc$cases, replace=TRUE)
thresholds <- roc.utils.thresholds(c(cases, controls))

perfs <- sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction) * ifelse(roc$percent, 100, 1)

# update ROC
roc$sensitivities <- perfs[2,]
roc$specificities <- perfs[1,]
roc$cases <- cases
roc$controls <- controls
roc$predictor <- c(controls, cases)
roc$response <- c(rep(roc$levels[1], length(controls)), rep(roc$levels[2], length(cases)))
roc$thresholds <- thresholds

# call smooth.roc and auc.smooth.roc
smooth.roc.call$roc <- roc
smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
if (is(smooth.roc, "try-error"))
return(NA)
as.numeric(coords.roc(smooth.roc, x=x, input=input, ret=ret, best.method=best.method, best.weights=best.weights))
}

nonstratified.ci.smooth.coords <- function(roc, x, input, ret, best.method, best.weights, smooth.roc.call) {
tmp.idx <- sample(1:length(roc$predictor), replace=TRUE)
predictor <- roc$predictor[tmp.idx]
response <- roc$response[tmp.idx]
splitted <- split(predictor, response)
controls <- splitted[[as.character(roc$levels[1])]]
cases <- splitted[[as.character(roc$levels[2])]]
thresholds <- roc.utils.thresholds(c(cases, controls))

perfs <- sapply(thresholds, roc.utils.perfs, controls=controls, cases=cases, direction=roc$direction) * ifelse(roc$percent, 100, 1)

# update ROC
roc$sensitivities <- perfs[2,]
roc$specificities <- perfs[1,]
roc$cases <- cases
roc$controls <- controls
roc$predictor <- predictor
roc$response <- response
roc$thresholds <- thresholds

# call smooth.roc and auc.smooth.roc
smooth.roc.call$roc <- roc
smooth.roc <- try(eval(smooth.roc.call), silent=TRUE)
if (is(smooth.roc, "try-error"))
return(NA)
as.numeric(coords.roc(smooth.roc, x=x, input=input, ret=ret, best.method=best.method, best.weights=best.weights))
}
4 changes: 2 additions & 2 deletions R/ci.R
Expand Up @@ -29,7 +29,7 @@ ci.default <- function(response, predictor, ...) {
ci.roc(roc.default(response, predictor, ...), ...)
}

ci.smooth.roc <- function(smooth.roc, of = c("auc", "sp", "se"), ...) {
ci.smooth.roc <- function(smooth.roc, of = c("auc", "sp", "se", "coords"), ...) {
of <- match.arg(of)

if (of == "auc")
Expand All @@ -42,7 +42,7 @@ ci.smooth.roc <- function(smooth.roc, of = c("auc", "sp", "se"), ...) {
return(ci)
}

ci.roc <- function(roc, of = c("auc", "thresholds", "sp", "se"), ...) {
ci.roc <- function(roc, of = c("auc", "thresholds", "sp", "se", "coords"), ...) {
of <- match.arg(of)

if (of == "auc")
Expand Down
151 changes: 151 additions & 0 deletions R/ci.coords.R
@@ -0,0 +1,151 @@
# pROC: Tools Receiver operating characteristic (ROC curves) with
# (partial) area under the curve, confidence intervals and comparison.
# Copyright (C) 2010, 2011 Xavier Robin, Alexandre Hainard, Natacha Turck,
# Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez
# and Markus Müller
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

ci.coords <- function(...) {
UseMethod("ci.coords")
}

ci.coords.formula <- function(formula, data, ...) {
ci.coords(roc.formula(formula, data, ci=FALSE, ...), ...)
}

ci.coords.default <- function(response, predictor, ...) {
ci.coords(roc.default(response, predictor, ci=FALSE, ...), ...)
}

ci.coords.smooth.roc <- function(smooth.roc,
x,
input=c("specificity", "sensitivity"), ret=c("specificity", "sensitivity"),
best.method=c("youden", "closest.topleft"), best.weights=c(1, 0.5),
conf.level = 0.95,
boot.n = 2000,
boot.stratified = TRUE,
progress = getOption("pROCProgress")$name,
...
) {
if (conf.level > 1 | conf.level < 0)
stop("'conf.level' must be within the interval [0,1].")
input <- match.arg(input)
ret <- roc.utils.match.coords.ret.args(ret)
if (is.character(x)) {
x <- match.arg(x, c("all", "local maximas", "best"))
if (x == "all" || x == "local maximas") {
stop("'all' and 'local maximas' are not available for confidence intervals.")
}
}

# Check if called with density.cases or density.controls
if (is.null(smooth.roc$smoothing.args) || is.numeric(smooth.roc$smoothing.args$density.cases) || is.numeric(smooth.roc$smoothing.args$density.controls))
stop("Cannot compute CI of ROC curves smoothed with numeric density.controls and density.cases.")

# Get the non smoothed roc.
roc <- attr(smooth.roc, "roc")
roc$ci <- NULL # remove potential ci in roc to avoid infinite loop with smooth.roc()

# prepare the calls
smooth.roc.call <- as.call(c(match.fun("smooth.roc"), smooth.roc$smoothing.args))

if(class(progress) != "list")
progress <- roc.utils.get.progress.bar(progress, title="Coords confidence interval", label="Bootstrap in progress...", ...)

if (boot.stratified) {
perfs <- raply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call), .progress=progress)
}
else {
perfs <- raply(boot.n, nonstratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights,smooth.roc.call), .progress=progress)
}

if (any(is.na(perfs))) {
warning("NA value(s) produced during bootstrap were ignored.")
perfs <- perfs[!apply(perfs, 1, function(x) any(is.na(x))),]
}

if (length(x) > 1) {
inputs <- paste(input, x)
rownames.ret <- paste(rep(inputs, each=length(ret)), ret, sep=": ")
}
else {
rownames.ret <- ret
}

ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
rownames(ci) <- rownames.ret

class(ci) <- "ci.coords"
attr(ci, "conf.level") <- conf.level
attr(ci, "boot.n") <- boot.n
attr(ci, "boot.stratified") <- boot.stratified
attr(ci, "roc") <- smooth.roc
return(ci)
}

ci.coords.roc <- function(roc,
x,
input=c("threshold", "specificity", "sensitivity"), ret=c("threshold", "specificity", "sensitivity"),
best.method=c("youden", "closest.topleft"), best.weights=c(1, 0.5),
conf.level = 0.95,
boot.n = 2000,
boot.stratified = TRUE,
progress = getOption("pROCProgress")$name,
...
) {
if (conf.level > 1 | conf.level < 0)
stop("'conf.level' must be within the interval [0,1].")
input <- match.arg(input)
ret <- roc.utils.match.coords.ret.args(ret)
if (is.character(x)) {
x <- match.arg(x, c("all", "local maximas", "best"))
if (x == "all" || x == "local maximas") {
stop("'all' and 'local maximas' are not available for confidence intervals.")
}
}

if(class(progress) != "list")
progress <- roc.utils.get.progress.bar(progress, title="Coords confidence interval", label="Bootstrap in progress...", ...)

if (boot.stratified) {
perfs <- raply(boot.n, stratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress, .drop=FALSE)
}
else {
perfs <- raply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress, .drop=FALSE)
}

if (any(is.na(perfs))) {
warning("NA value(s) produced during bootstrap were ignored.")
perfs <- perfs[!apply(perfs, 1, function(x) any(is.na(x))),]
}

if (length(x) > 1) {
inputs <- paste(input, x)
rownames.ret <- paste(rep(inputs, each=length(ret)), ret, sep=": ")
}
else {
rownames.ret <- ret
}

ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2)))
rownames(ci) <- rownames.ret

class(ci) <- c("ci.coords", "matrix")
attr(ci, "conf.level") <- conf.level
attr(ci, "boot.n") <- boot.n
attr(ci, "boot.stratified") <- boot.stratified
attr(ci, "roc") <- roc
return(ci)
}
6 changes: 1 addition & 5 deletions R/coords.R
Expand Up @@ -133,11 +133,7 @@ coords.roc <- function(roc, x, input=c("threshold", "specificity", "sensitivity"
# match input
input <- match.arg(input)
# match return
valid.ret.args <- c("threshold", "specificity", "sensitivity", "accuracy", "tn", "tp", "fn", "fp", "npv", "ppv", "1-specificity", "1-sensitivity", "1-accuracy", "1-npv", "1-ppv")
ret <- replace(ret, ret=="t", "threshold")
ret <- replace(ret, ret=="npe", "1-npv")
ret <- replace(ret, ret=="ppe", "1-ppv")
ret <- match.arg(ret, valid.ret.args, several.ok=TRUE)
ret <- roc.utils.match.coords.ret.args(ret)
# make sure the sort of roc is correct
roc <- sort(roc)

Expand Down
4 changes: 3 additions & 1 deletion R/onLoad.R
Expand Up @@ -21,11 +21,13 @@
# Generate progressbar option with smart default values
if (is.null(getOption("pROCProgress"))) {
if (interactive()) {
# Check the presence of tcltk
tcltk.present <- length(find.package("tcltk", quiet=TRUE)) > 0
if (!is.null(getOption("STERM")) && getOption("STERM") == "iESS")
options("pROCProgress" = list(name = "text", width = NA, char = "=", style = 1))
else if (.Platform$OS.type == "windows")
options("pROCProgress" = list(name = "win", width = 300))
else if ("tcltk" %in% installed.packages()[,1] && Sys.getenv("DISPLAY") != "")
else if (tcltk.present && Sys.getenv("DISPLAY") != "")
options("pROCProgress" = list(name = "tk", width = 300))
else
options("pROCProgress" = list(name = "text", width = NA, char = "=", style = 3))
Expand Down
17 changes: 17 additions & 0 deletions R/print.R
Expand Up @@ -168,6 +168,23 @@ print.ci.se <- function(x, digits=max(3, getOption("digits") - 3), ...) {
invisible(x)
}

print.ci.coords <- function(x, digits=max(3, getOption("digits") - 3), ...) {
cat(attr(x, "conf.level")*100, "% CI", sep="")
cat(" (", attr(x, "boot.n"), " ", ifelse(attr(x, "boot.stratified"), "stratified", "non-stratified"), " bootstrap replicates):\n", sep="")

# Back to the standard object we'll print
unattr.coords <- x
class(unattr.coords) <- "matrix"
attr(unattr.coords, "conf.level") <- NULL
attr(unattr.coords, "boot.n") <- NULL
attr(unattr.coords, "boot.stratified") <- NULL
attr(unattr.coords, "roc") <- NULL

class(unattr.coords) <- "matrix"
print(unattr.coords, digits=digits)
invisible(x)
}

print.dataline <- function(x) {
# Case / Controls call
if ("cases" %in% names(x$call) && "controls" %in% names(x$call)) {
Expand Down
20 changes: 18 additions & 2 deletions R/roc.utils.R
Expand Up @@ -112,8 +112,14 @@ roc.utils.max.thresholds.idx <- function(thresholds, sp, se) {
# Define which progress bar to use
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, ...) {
if (name == "tk") { # load tcltk if possible
if (!require(tcltk))
stop("Package tcltk not available, required with progress='tk'")
if (!require(tcltk)) {
# If tcltk cannot be loaded fall back to default text progress bar
name <- "text"
style <- 3
char <- "="
width <- NA
warning("Package tcltk required with progress='tk' but could not be loaded. Falling back to text progress bar.")
}
}
if (name == "none")
progress_none()
Expand Down Expand Up @@ -150,3 +156,13 @@ sort.smooth.roc <- function(roc) {
}
return(roc)
}

# Arguments which can be returned by coords
roc.utils.match.coords.ret.args <- function(x) {
valid.ret.args <- c("threshold", "specificity", "sensitivity", "accuracy", "tn", "tp", "fn", "fp", "npv", "ppv", "1-specificity", "1-sensitivity", "1-accuracy", "1-npv", "1-ppv")
x <- replace(x, x=="t", "threshold")
x <- replace(x, x=="npe", "1-npv")
x <- replace(x, x=="ppe", "1-ppv")
match.arg(x, valid.ret.args, several.ok=TRUE)
}

0 comments on commit deda1b3

Please sign in to comment.