From 570a08f316c4163dc2c14c6194da2635ebfb5aee Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Wed, 14 Dec 2011 14:39:18 +0100 Subject: [PATCH 01/16] Initial work on ci.coords --- R/ci.coords.R | 120 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 R/ci.coords.R diff --git a/R/ci.coords.R b/R/ci.coords.R new file mode 100644 index 0000000..91776ec --- /dev/null +++ b/R/ci.coords.R @@ -0,0 +1,120 @@ +# 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 . + +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].") + + # 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="SE confidence interval", label="Bootstrap in progress...", ...) + + if (boot.stratified) { + perfs <- rdply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call), .progress=progress)[-1] + } + else { + perfs <- rdply(boot.n, nonstratified.ci.smooth.se(roc, x, input, ret, best.method, best.weights,smooth.roc.call), .progress=progress)[-1] + } + + 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))),] + } + + ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2))) + rownames(ci) <- paste(specificities, ifelse(roc$percent, "%", ""), sep="") + + 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), + specificities = seq(0, 1, .1) * ifelse(roc$percent, 100, 1), + 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].") + + if(class(progress) != "list") + progress <- roc.utils.get.progress.bar(progress, title="SE confidence interval", label="Bootstrap in progress...", ...) + + if (boot.stratified) { + perfs <- rdply(boot.n, stratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress)[-1] + } + else { + perfs <- rdply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress)[-1] + } + + 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))),] + } + + ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2))) + rownames(ci) <- paste(specificities, ifelse(roc$percent, "%", ""), sep="") + + 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") <- roc + return(ci) +} From 83df85d139b9e9b21effc8a8b182f3155803eeb3 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Wed, 21 Dec 2011 15:52:31 +0100 Subject: [PATCH 02/16] pROC now also depends on methods package ('is') --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index abf443a..e2d32d1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -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 From cf4803146dcc4ca61fa9157deddc6f250a932600 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 22 Dec 2011 13:26:05 +0100 Subject: [PATCH 03/16] Plyr not used directly in roc function --- man/roc.Rd | 6 ------ 1 file changed, 6 deletions(-) diff --git a/man/roc.Rd b/man/roc.Rd index a112f55..4420b8e 100644 --- a/man/roc.Rd +++ b/man/roc.Rd @@ -230,16 +230,10 @@ plot=FALSE, smooth.method="binormal", ci.method=NULL, density=NULL, ...) (2011) ``pROC: an open-source package for R and S+ to analyze and compare ROC curves''. \emph{BMC Bioinformatics}, \bold{7}, 77. DOI: \href{http://dx.doi.org/10.1186/1471-2105-12-77}{10.1186/1471-2105-12-77}. - - Hadley Wickham (2011) ``The Split-Apply-Combine Strategy for Data Analysis''. \emph{Journal of Statistical Software}, \bold{40}, 1--29. - URL: \href{http://www.jstatsoft.org/v40/i01}{www.jstatsoft.org/v40/i01}. - } \seealso{ \code{\link{auc}}, \code{\link{ci}}, \code{\link{plot.roc}}, \code{\link{print.roc}}, \code{\link{roc.test}} - - CRAN package \pkg{plyr}, employed in this function. } \examples{ From 7be3e8e7625f329d90dabcfb507e1fb78f9ae301 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 22 Dec 2011 13:35:41 +0100 Subject: [PATCH 04/16] Typo --- man/smooth.roc.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/smooth.roc.Rd b/man/smooth.roc.Rd index 5d7928f..9bc641f 100644 --- a/man/smooth.roc.Rd +++ b/man/smooth.roc.Rd @@ -232,7 +232,7 @@ reuse.auc=TRUE, reuse.ci=FALSE, ...) require additional packages. If not available, the following message will be displayed with the required command to install the package: \dQuote{Package ? not available, required with method='?'. - Please install it with 'install.packages("?")'." + Please install it with 'install.packages("?")'. } } From d923552426b5f6d31b7767105b28b09b9d691458 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 22 Dec 2011 13:46:29 +0100 Subject: [PATCH 05/16] Properly citing MASS --- man/smooth.roc.Rd | 3 +++ 1 file changed, 3 insertions(+) diff --git a/man/smooth.roc.Rd b/man/smooth.roc.Rd index 9bc641f..07f2159 100644 --- a/man/smooth.roc.Rd +++ b/man/smooth.roc.Rd @@ -248,6 +248,9 @@ reuse.auc=TRUE, reuse.ci=FALSE, ...) (2011) ``pROC: an open-source package for R and S+ to analyze and compare ROC curves''. \emph{BMC Bioinformatics}, \bold{7}, 77. DOI: \href{http://dx.doi.org/10.1186/1471-2105-12-77}{10.1186/1471-2105-12-77}. + + William N. Venables, Brian D. Ripley (2002). ``Modern Applied Statistics with S''. New York, Springer. + \href{http://books.google.ch/books?id=974c4vKurNkC}{Google books}. Kelly H. Zou, W. J. Hall and David E. Shapiro (1997) ``Smooth non-parametric receiver operating characteristic (ROC) curves for From 0dc7c858ed9e8b82ea92e4e21ddf3e59c35bebbc Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Fri, 20 Jan 2012 18:01:11 +0100 Subject: [PATCH 06/16] More on ci.coords --- R/bootstrap.R | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/R/bootstrap.R b/R/bootstrap.R index c4cffb3..4075253 100644 --- a/R/bootstrap.R +++ b/R/bootstrap.R @@ -591,3 +591,34 @@ nonstratified.ci.thresholds <- function(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) + roc$sensitivities <- perfs[2,] + roc$specificities <- perfs[1,] + + as.numeric(coords.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"))) +} + +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) + roc$sensitivities <- perfs[2,] + roc$specificities <- perfs[1,] + + as.numeric(auc.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"))) +} + From fb3be140530432903fe43282db436e979e8d470b Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 8 Mar 2012 13:03:10 +0100 Subject: [PATCH 07/16] Removing installed.packages() call in .onLoad --- R/onLoad.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/onLoad.R b/R/onLoad.R index dd7d51d..10b0ed9 100644 --- a/R/onLoad.R +++ b/R/onLoad.R @@ -25,7 +25,7 @@ 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 (suppressPackageStartupMessages(suppressWarnings(require(tcltk))) && Sys.getenv("DISPLAY") != "") options("pROCProgress" = list(name = "tk", width = 300)) else options("pROCProgress" = list(name = "text", width = NA, char = "=", style = 3)) From acc90787975896c1554cb9deae4e19c8ccc9c712 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Fri, 9 Mar 2012 11:28:32 +0100 Subject: [PATCH 08/16] Better fix for slow startup. Use find.packages and fallback to text progress bar if tcltk is installed but cannot be loaded. --- NEWS | 3 +++ R/onLoad.R | 4 +++- R/roc.utils.R | 10 ++++++++-- 3 files changed, 14 insertions(+), 3 deletions(-) diff --git a/NEWS b/NEWS index de1c097..9932f38 100644 --- a/NEWS +++ b/NEWS @@ -3,6 +3,9 @@ * New 'cov' and 'var' functions supports new "obuchowski" method * '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" diff --git a/R/onLoad.R b/R/onLoad.R index 10b0ed9..e7bb2aa 100644 --- a/R/onLoad.R +++ b/R/onLoad.R @@ -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 (suppressPackageStartupMessages(suppressWarnings(require(tcltk))) && 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)) diff --git a/R/roc.utils.R b/R/roc.utils.R index cee495c..ae97b0e 100644 --- a/R/roc.utils.R +++ b/R/roc.utils.R @@ -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() From 3ebc4d3b295013c23d6547304cb4459c82873ce2 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Tue, 13 Mar 2012 16:23:50 +0100 Subject: [PATCH 09/16] Documenting the behavior when the value of predictor == threshold. --- man/roc.Rd | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/man/roc.Rd b/man/roc.Rd index 4420b8e..9881a58 100644 --- a/man/roc.Rd +++ b/man/roc.Rd @@ -65,9 +65,9 @@ plot=FALSE, smooth.method="binormal", ci.method=NULL, density=NULL, ...) \dQuote{auto} (default): automatically define in which group the median is higher and take the direction accordingly. \dQuote{>}: if the predictor values for the control group are - higher than the values of the case group. \dQuote{<}: if the predictor - values for the control group are higher or equal than the values of - the case group. + higher than the values of the case group (controls > t >= cases). + \dQuote{<}: if the predictor values for the control group are higher + or equal than the values of the case group (controls < t <= cases). } \item{smooth}{if TRUE, the ROC curve is passed to \code{\link{smooth}} to be smoothed. From 34cbcd90410843ea0477024fab61c78f9b13751b Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Wed, 14 Mar 2012 09:55:11 +0100 Subject: [PATCH 10/16] methods and utils packages need to be imported under some specific conditions --- NAMESPACE | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 9568454..2b2ef08 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -20,3 +20,5 @@ export( smooth, smooth.roc, smooth.default, smooth.smooth.roc, var, var.default, var.auc, var.smooth.roc, var.roc ) + +import(methods, utils) From ef21bccc9e9368d11ad67a54cc6d2fb0f7d8aee5 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 15 Mar 2012 14:33:00 +0100 Subject: [PATCH 11/16] ci.coords basically working now --- R/bootstrap.R | 16 ++++++++++++++-- R/ci.coords.R | 28 ++++++++++++++++++++++------ R/coords.R | 6 +----- R/print.R | 17 +++++++++++++++++ R/roc.utils.R | 10 ++++++++++ 5 files changed, 64 insertions(+), 13 deletions(-) diff --git a/R/bootstrap.R b/R/bootstrap.R index 4075253..fbabe8b 100644 --- a/R/bootstrap.R +++ b/R/bootstrap.R @@ -600,10 +600,16 @@ stratified.ci.coords <- function(roc, x, input, ret, best.method, best.weights) 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, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"))) + 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) { @@ -616,9 +622,15 @@ nonstratified.ci.coords <- function(roc, x, input, ret, best.method, best.weight 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(auc.roc(roc, partial.auc=attr(roc$auc, "partial.auc"), partial.auc.focus=attr(roc$auc, "partial.auc.focus"), partial.auc.correct=attr(roc$auc, "partial.auc.correct"))) + as.numeric(coords.roc(roc, x=x, input=input, ret=ret, best.method=best.method, best.weights=best.weights)) } diff --git a/R/ci.coords.R b/R/ci.coords.R index 91776ec..517db40 100644 --- a/R/ci.coords.R +++ b/R/ci.coords.R @@ -60,7 +60,7 @@ ci.coords.smooth.roc <- function(smooth.roc, perfs <- rdply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call), .progress=progress)[-1] } else { - perfs <- rdply(boot.n, nonstratified.ci.smooth.se(roc, x, input, ret, best.method, best.weights,smooth.roc.call), .progress=progress)[-1] + perfs <- rdply(boot.n, nonstratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights,smooth.roc.call), .progress=progress)[-1] } if (any(is.na(perfs))) { @@ -69,7 +69,7 @@ ci.coords.smooth.roc <- function(smooth.roc, } ci <- t(apply(perfs, 2, quantile, probs=c(0+(1-conf.level)/2, .5, 1-(1-conf.level)/2))) - rownames(ci) <- paste(specificities, ifelse(roc$percent, "%", ""), sep="") + rownames(ci) <- ret class(ci) <- "ci.coords" attr(ci, "conf.level") <- conf.level @@ -92,15 +92,23 @@ ci.coords.roc <- function(roc, ) { 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="SE confidence interval", label="Bootstrap in progress...", ...) if (boot.stratified) { - perfs <- rdply(boot.n, stratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress)[-1] + perfs <- raply(boot.n, stratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress, .drop=FALSE) } else { - perfs <- rdply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress)[-1] + perfs <- raply(boot.n, nonstratified.ci.coords(roc, x, input, ret, best.method, best.weights), .progress=progress, .drop=FALSE) } if (any(is.na(perfs))) { @@ -108,10 +116,18 @@ ci.coords.roc <- function(roc, 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) <- paste(specificities, ifelse(roc$percent, "%", ""), sep="") + rownames(ci) <- rownames.ret - class(ci) <- "ci.coords" + class(ci) <- c("ci.coords", "matrix") attr(ci, "conf.level") <- conf.level attr(ci, "boot.n") <- boot.n attr(ci, "boot.stratified") <- boot.stratified diff --git a/R/coords.R b/R/coords.R index 4bc0e93..0c0a515 100644 --- a/R/coords.R +++ b/R/coords.R @@ -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) diff --git a/R/print.R b/R/print.R index 7452a6a..a97fe4e 100644 --- a/R/print.R +++ b/R/print.R @@ -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)) { diff --git a/R/roc.utils.R b/R/roc.utils.R index cee495c..d95a2a4 100644 --- a/R/roc.utils.R +++ b/R/roc.utils.R @@ -150,3 +150,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) +} + From b6606efe85d30d4b9f7dc536ede3e32f1270a6c2 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 15 Mar 2012 14:52:26 +0100 Subject: [PATCH 12/16] ci.coords for smooth.roc curves --- R/bootstrap.R | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++ R/ci.coords.R | 23 ++++++++++++++++++---- 2 files changed, 72 insertions(+), 4 deletions(-) diff --git a/R/bootstrap.R b/R/bootstrap.R index fbabe8b..845fbca 100644 --- a/R/bootstrap.R +++ b/R/bootstrap.R @@ -634,3 +634,56 @@ nonstratified.ci.coords <- function(roc, x, input, ret, best.method, best.weight 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)) +} diff --git a/R/ci.coords.R b/R/ci.coords.R index 517db40..af68c8a 100644 --- a/R/ci.coords.R +++ b/R/ci.coords.R @@ -41,6 +41,14 @@ ci.coords.smooth.roc <- function(smooth.roc, ) { 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)) @@ -57,10 +65,10 @@ ci.coords.smooth.roc <- function(smooth.roc, progress <- roc.utils.get.progress.bar(progress, title="SE confidence interval", label="Bootstrap in progress...", ...) if (boot.stratified) { - perfs <- rdply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call), .progress=progress)[-1] + perfs <- raply(boot.n, stratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights, smooth.roc.call), .progress=progress) } else { - perfs <- rdply(boot.n, nonstratified.ci.smooth.coords(roc, x, input, ret, best.method, best.weights,smooth.roc.call), .progress=progress)[-1] + 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))) { @@ -68,8 +76,16 @@ ci.coords.smooth.roc <- function(smooth.roc, 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) <- ret + rownames(ci) <- rownames.ret class(ci) <- "ci.coords" attr(ci, "conf.level") <- conf.level @@ -83,7 +99,6 @@ 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), - specificities = seq(0, 1, .1) * ifelse(roc$percent, 100, 1), conf.level = 0.95, boot.n = 2000, boot.stratified = TRUE, From 4452abc2393ae0fc094bea830870a04a526cfb65 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 15 Mar 2012 15:05:26 +0100 Subject: [PATCH 13/16] Call ci.coords directly with ci. --- R/ci.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/ci.R b/R/ci.R index 64ae158..cad9dd9 100644 --- a/R/ci.R +++ b/R/ci.R @@ -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") @@ -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") From 7e459520789bf34b7ff69d71dc3c5de033ab5085 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 15 Mar 2012 15:10:00 +0100 Subject: [PATCH 14/16] Correct bootstrap message for ci.coords --- R/ci.coords.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/ci.coords.R b/R/ci.coords.R index af68c8a..748b943 100644 --- a/R/ci.coords.R +++ b/R/ci.coords.R @@ -62,7 +62,7 @@ ci.coords.smooth.roc <- function(smooth.roc, 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="SE confidence interval", label="Bootstrap in progress...", ...) + 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) @@ -117,7 +117,7 @@ ci.coords.roc <- function(roc, } if(class(progress) != "list") - progress <- roc.utils.get.progress.bar(progress, title="SE confidence interval", label="Bootstrap in progress...", ...) + 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) From 7c38eec481b4cecd89503aeffa49553511040c59 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 15 Mar 2012 15:11:48 +0100 Subject: [PATCH 15/16] Documenting ci.coords --- man/ci.Rd | 13 ++++++------- man/coords.Rd | 2 +- man/pROC-package.Rd | 3 ++- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/man/ci.Rd b/man/ci.Rd index 7d96a38..31c75ed 100644 --- a/man/ci.Rd +++ b/man/ci.Rd @@ -17,8 +17,8 @@ } \usage{ ci(...) -\S3method{ci}{roc}(roc, of = c("auc", "thresholds", "sp", "se"), ...) -\S3method{ci}{smooth.roc}(smooth.roc, of = c("auc", "sp", "se"), ...) +\S3method{ci}{roc}(roc, of = c("auc", "thresholds", "sp", "se", "coords"), ...) +\S3method{ci}{smooth.roc}(smooth.roc, of = c("auc", "sp", "se", "coords"), ...) \S3method{ci}{formula}(formula, data, ...) \S3method{ci}{default}(response, predictor, ...) } @@ -33,7 +33,7 @@ ci(...) response~predictor for the \code{\link{roc}} function. } \item{of}{The type of confidence interval. One of \dQuote{auc}, - \dQuote{thresholds}, \dQuote{sp} or \dQuote{se}. Note that + \dQuote{thresholds}, \dQuote{sp}, \dQuote{se} or \dQuote{coords}. Note that confidence interval on \dQuote{thresholds} are not available for smoothed ROC curves. } @@ -54,13 +54,12 @@ ci(...) This function is typically called from \code{\link{roc}} when \code{ci=TRUE} (not by default). Depending on the \code{of} argument, the specific \code{ci} functions \code{\link{ci.auc}}, \code{\link{ci.thresholds}}, - \code{\link{ci.sp}} or \code{\link{ci.se}} are called. + \code{\link{ci.sp}}, \code{\link{ci.se}} or \code{\link{ci.coords}} are called. } \value{ The return value of the specific \code{ci} functions -\code{\link{ci.auc}}, \code{\link{ci.thresholds}}, \code{\link{ci.sp}} -or \code{\link{ci.se}}, depending on the +\code{\link{ci.auc}}, \code{\link{ci.thresholds}}, \code{\link{ci.sp}}, \code{\link{ci.se}} or \code{\link{ci.coords}}, depending on the \code{of} argument. } @@ -73,7 +72,7 @@ or \code{\link{ci.se}}, depending on the \seealso{ \code{\link{roc}}, \code{\link{auc}}, \code{\link{ci.auc}}, - \code{\link{ci.thresholds}}, \code{\link{ci.sp}}, \code{\link{ci.se}} + \code{\link{ci.thresholds}}, \code{\link{ci.sp}}, \code{\link{ci.se}}, \code{\link{ci.coords}} } \examples{ diff --git a/man/coords.Rd b/man/coords.Rd index d916077..c5d8621 100644 --- a/man/coords.Rd +++ b/man/coords.Rd @@ -207,7 +207,7 @@ best.weights=c(1, 0.5), ...) } \seealso{ -\code{\link{roc}} +\code{\link{roc}}, \code{\link{ci.coords}} } \examples{ data(aSAH) diff --git a/man/pROC-package.Rd b/man/pROC-package.Rd index 20eb724..ce681d4 100644 --- a/man/pROC-package.Rd +++ b/man/pROC-package.Rd @@ -70,10 +70,11 @@ \code{\link{auc}} \tab Compute the area under the ROC curve \cr \code{\link{ci}} \tab Compute confidence intervals of a ROC curve \cr \code{\link{ci.auc}} \tab Compute the CI of the AUC \cr + \code{\link{ci.coords}} \tab Compute the CI of arbitrary coordinates \cr \code{\link{ci.se}} \tab Compute the CI of sensitivities at given specificities \cr \code{\link{ci.sp}} \tab Compute the CI of specificities at given sensitivities \cr \code{\link{ci.thresholds}} \tab Compute the CI of specificity and sensitivity of thresholds \cr - \code{\link{coords}} \tab Coordinates of a ROC curve \cr + \code{\link{coords}} \tab Coordinates of the ROC curve \cr \code{\link[=cov.roc]{cov}} \tab Covariance between two AUCs\cr \code{\link{has.partial.auc}} \tab Determine if the ROC curve have a partial AUC\cr \code{\link{lines.roc}} \tab Add a ROC line to a ROC plot \cr From 89e8ca8f60e054e9a0aa4cb262812172bc54b4b4 Mon Sep 17 00:00:00 2001 From: Xavier Robin Date: Thu, 15 Mar 2012 15:15:57 +0100 Subject: [PATCH 16/16] News and namespace for ci.coords --- NAMESPACE | 1 + NEWS | 1 + 2 files changed, 2 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 2b2ef08..f5f5716 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -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, diff --git a/NEWS b/NEWS index 9932f38..aad40b0 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,7 @@ 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)