R/smooth.R

In pROC: Display and Analyze ROC Curves

# pROC: Tools Receiver operating characteristic (ROC curves) with
# (partial) area under the curve, confidence intervals and comparison.
# Copyright (C) 2010-2014 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/>.

smooth <- function(...) {
  UseMethod("smooth")
}

smooth.default <- function(...) {
  stats::smooth(...)
}

smooth.smooth.roc <- function(smooth.roc, ...) {
  roc <- attr(smooth.roc, "roc")
  if (is.null(roc)) {
    stop("Cannot smooth a ROC curve generated directly with numeric 'density.controls' and 'density.cases'.")
  }
  smooth.roc(roc, ...)
}

smooth.roc <- function(roc, method = c("binormal", "density", "fitdistr", "logcondens", "logcondens.smooth"), n = 512, bw = "nrd0",
                       density = NULL, density.controls = density, density.cases = density,
                       start = NULL, start.controls = start, start.cases = start,
                       reuse.auc = TRUE, reuse.ci = FALSE,
                       ...) {
  method <- match.arg(method)

  if (is.ordered(roc$original.predictor) && (method == "density" || method == "fitidstr")) {
    stop("ROC curves of ordered predictors can be smoothed only with binormal smoothing.")
  }

  if (method == "binormal") {
    sesp <- smooth_roc_binormal(roc, n)
  } else if (method == "fitdistr") {
    sesp <- smooth_roc_fitdistr(roc, n, density.controls, density.cases, start.controls, start.cases, ...)
  } else if (method == "density") {
    sesp <- smooth_roc_density(roc, n, density.controls, density.cases, bw, ...)
  } else if (method == "logcondens") {
    sesp <- smooth_roc_logcondens(roc, n)
  } else if (method == "logcondens.smooth") {
    sesp <- smooth_roc_logcondens_smooth(roc, n)
  } else {
    stop(sprintf(
      "Impossible smooth method value '%s'. Please report this bug to <%s>.",
      method, utils::packageDescription("pROC")$BugReports
    ))
  }

  class(sesp) <- "smooth.roc"
  sesp <- sort_smooth_roc(sesp) # sort se and sp
  # anchor SE/SP at 0/100
  sesp$specificities <- c(0, as.vector(sesp$specificities), ifelse(roc$percent, 100, 1))
  sesp$sensitivities <- c(ifelse(roc$percent, 100, 1), as.vector(sesp$sensitivities), 0)
  attr(sesp, "roc") <- roc # keep the original roc. May be useful in CI.
  sesp$percent <- roc$percent # keep some basic roc specifications
  sesp$direction <- roc$direction
  sesp$call <- match.call()
  # keep smoothing arguments (for print and bootstrap)
  sesp$smoothing.args <- list(...)
  sesp$smoothing.args$method <- method
  sesp$smoothing.args$n <- n
  sesp$smoothing.args$start.controls <- start.controls
  sesp$smoothing.args$start.cases <- start.cases
  sesp$smoothing.args$density.controls <- density.controls
  sesp$smoothing.args$density.cases <- density.cases
  sesp$smoothing.args$bw <- bw
  # complete fit.controls/cases if a function was passed as densfun
  if (method == "fitdistr") {
    if (is.null(sesp$fit.controls$densfun)) {
      if (missing(density.controls)) {
        sesp$fit.controls$densfun <- deparse(substitute(density))
      } else {
        sesp$fit.controls$densfun <- deparse(substitute(density.controls))
      }
    }
    if (is.null(sesp$fit.cases$densfun)) {
      if (missing(density.cases)) {
        sesp$fit.cases$densfun <- deparse(substitute(density))
      } else {
        sesp$fit.cases$densfun <- deparse(substitute(density.cases))
      }
    }
  }

  # if there was an auc and a ci, re-do them
  if (!is.null(roc$auc) && reuse.auc) {
    args <- attributes(roc$auc)
    args$roc <- NULL
    args$smooth.roc <- sesp
    sesp$auc <- do.call("auc.smooth.roc", args)
  }
  if (!is.null(roc$ci) && reuse.ci) {
    args <- attributes(roc$ci)
    args$roc <- NULL
    args$smooth.roc <- sesp
    sesp$ci <- do.call(paste(class(roc$ci), "smooth.roc", sep = "."), args)
  }

  return(sesp)
}

smooth_roc_density <- function(roc, n, density.controls, density.cases, bw,
                               # catch args for density
                               cut = 3, adjust = 1, kernel = window, window = "gaussian",
                               percent = roc$percent, direction = roc$direction,
                               ...) {
  if (!is.numeric(density.controls) || !is.numeric(density.cases)) {
    predictor <- c(roc$controls, roc$cases)
    if (is.character(bw)) {
      bw <- match.fun(paste("bw", bw, sep = "."))(predictor)
    }
    bw <- bw * adjust
    from <- min(predictor) - (cut * bw)
    to <- max(predictor) + (cut * bw)
  }
  if (mode(density.controls) == "function") {
    density.controls <- density.controls(roc$controls, n = n, from = from, to = to, bw = bw, kernel = kernel, ...)
    if (!is.numeric(density.controls)) {
      if (is.list(density.controls) && !is.null(density.controls$y) && is.numeric(density.controls$y)) {
        density.controls <- density.controls$y
      } else {
        stop("The 'density' function must return a numeric vector or a list with a 'y' item.")
      }
    }
  } else if (is.null(density.controls)) {
    density.controls <- suppressWarnings(density(roc$controls, n = n, from = from, to = to, bw = bw, kernel = kernel, ...))$y
  } else if (!is.numeric(density.controls)) {
    stop("'density.controls' must be either NULL, a function or numeric values of density (over the y axis).")
  }
  if (mode(density.cases) == "function") {
    density.cases <- density.cases(roc$cases, n = n, from = from, to = to, bw = bw, kernel = kernel, ...)
    if (!is.numeric(density.cases)) {
      if (is.list(density.cases) && !is.null(density.cases$y) && is.numeric(density.cases$y)) {
        density.cases <- density.cases$y
      } else {
        stop("The 'density' function must return a numeric vector or a list with a 'y' item.")
      }
    }
  } else if (is.null(density.cases)) {
    density.cases <- suppressWarnings(density(roc$cases, n = n, from = from, to = to, bw = bw, kernel = kernel, ...))$y
  } else if (!is.numeric(density.cases)) {
    stop("'density.cases' must be either NULL, a function or numeric values of density (over the y axis).")
  }
  if (length(density.controls) != length(density.cases)) {
    stop("Length of 'density.controls' and 'density.cases' differ.")
  }

  perfs <- sapply((1:length(density.controls)) + .5, roc_utils_perfs_dens, x = (1:length(density.controls)) + .5, dens.controls = density.controls, dens.cases = density.cases, direction = direction)

  return(list(
    sensitivities = perfs[2, ] * ifelse(percent, 100, 1),
    specificities = perfs[1, ] * ifelse(percent, 100, 1)
  ))
}

smooth_roc_binormal <- function(roc, n) {
  df <- data.frame(sp = qnorm(roc$sp * ifelse(roc$percent, 1 / 100, 1)), se = qnorm(roc$se * ifelse(roc$percent, 1 / 100, 1)))
  df <- df[apply(df, 1, function(x) all(is.finite(x))), ]
  if (dim(df)[1] <= 1) { # ROC curve or with only 1 point
    stop("ROC curve not smoothable (not enough points).")
  }
  model <- lm(sp ~ se, df)
  if (any(is.na(model$coefficients[2]))) {
    stop("ROC curve not smoothable (not enough points).")
  }
  se <- qnorm(seq(0, 1, 1 / (n - 1)))
  sp <- predict(model, data.frame(se))

  return(list(
    sensitivities = pnorm(se) * ifelse(roc$percent, 100, 1),
    specificities = pnorm(sp) * ifelse(roc$percent, 100, 1),
    model = model
  ))
}

smooth_roc_fitdistr <- function(roc, n, densfun.controls, densfun.cases, start.controls, start.cases, ...) {
  load.suggested.package("MASS")

  densfuns.list <- list(
    beta = "dbeta", cauchy = "dcauchy", "chi-squared" = "dchisq", exponential = "dexp", f = "df",
    gamma = "dgamma", geometric = "dgeom", "log-normal" = "dlnorm", lognormal = "dlnorm",
    logistic = "dlogis", "negative binomial" = "dnbinom", normal = "dnorm", poisson = "dpois",
    t = "dt", weibull = "dweibull"
  )

  if (is.null(densfun.controls)) {
    densfun.controls <- "normal"
  } else if (is.character(densfun.controls)) {
    densfun.controls <- match.arg(densfun.controls, names(densfuns.list))
  }

  if (is.null(densfun.cases)) {
    densfun.cases <- "normal"
  } else if (is.character(densfun.cases)) {
    densfun.cases <- match.arg(densfun.cases, names(densfuns.list))
  }

  fit.controls <- MASS::fitdistr(roc$controls, densfun.controls, start.controls, ...)
  fit.cases <- MASS::fitdistr(roc$cases, densfun.cases, start.cases, ...)

  # store function name in fitting results
  if (mode(densfun.controls) != "function") {
    fit.controls$densfun <- densfun.controls
  }
  if (mode(densfun.cases) != "function") {
    fit.cases$densfun <- densfun.cases
  }

  x <- seq(min(c(roc$controls, roc$cases)), max(c(roc$controls, roc$cases)), length.out = n)

  # get the actual function name for do.call
  if (is.character(densfun.controls)) {
    densfun.controls <- match.fun(densfuns.list[[densfun.controls]])
  }
  if (is.character(densfun.cases)) {
    densfun.cases <- match.fun(densfuns.list[[densfun.cases]])
  }

  # ... that should be passed to densfun
  if (length(list(...)) > 0 && any(names(list(...)) %in% names(formals(densfun.controls)))) {
    dots.controls <- list(...)[names(formals(densfun.controls))[match(names(list(...)), names(formals(densfun.controls)))]]
  } else {
    dots.controls <- list()
  }
  if (length(list(...)) > 0 && any(names(list(...)) %in% names(formals(densfun.cases)))) {
    dots.cases <- list(...)[names(formals(densfun.cases))[match(names(list(...)), names(formals(densfun.cases)))]]
  } else {
    dots.cases <- list()
  }

  density.controls <- do.call(densfun.controls, c(list(x = x), fit.controls$estimate, dots.controls))
  density.cases <- do.call(densfun.cases, c(list(x = x), fit.cases$estimate, dots.cases))

  perfs <- sapply(x, roc_utils_perfs_dens, x = x, dens.controls = density.controls, dens.cases = density.cases, direction = roc$direction)

  return(list(
    sensitivities = perfs[2, ] * ifelse(roc$percent, 100, 1),
    specificities = perfs[1, ] * ifelse(roc$percent, 100, 1),
    fit.controls = fit.controls, fit.cases = fit.cases
  ))
}

smooth_roc_logcondens <- function(roc, n) {
  load.suggested.package("logcondens")

  sp <- seq(0, 1, 1 / (n - 1))
  logcondens <- logcondens::logConROC(roc$cases, roc$controls, sp)
  se <- logcondens$fROC

  return(list(
    sensitivities = se * ifelse(roc$percent, 100, 1),
    specificities = (1 - sp) * ifelse(roc$percent, 100, 1),
    logcondens = logcondens
  ))
}

smooth_roc_logcondens_smooth <- function(roc, n) {
  load.suggested.package("logcondens")

  sp <- seq(0, 1, 1 / (n - 1))
  logcondens <- logcondens::logConROC(roc$cases, roc$controls, sp)
  se <- logcondens$fROC.smooth

  return(list(
    sensitivities = se * ifelse(roc$percent, 100, 1),
    specificities = (1 - sp) * ifelse(roc$percent, 100, 1),
    logcondens = logcondens
  ))
}