R/power.roc.test.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) 2011-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/>.

power.roc.test <- function(...) {
  UseMethod("power.roc.test")
}

power.roc.test.roc <- function(roc1, roc2, sig.level = 0.05, power = NULL, kappa = NULL,
                               alternative = c("two.sided", "one.sided"),
                               reuse.auc = TRUE, method = c("delong", "bootstrap", "obuchowski"),
                               ...) {
  # Basic sanity checks
  if (!is.null(power) && (power < 0 || power > 1)) {
    stop("'power' must range from 0 to 1")
  }
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1)) {
    stop("'sig.level' must range from 0 to 1")
  }

  # check that the AUC of roc1 was computed, or do it now
  if (is.null(roc1$auc) | !reuse.auc) {
    roc1$auc <- auc(roc1, ...)
  }
  if (!is.null(attr(roc1$auc, "partial.auc.correct")) && attr(roc1$auc, "partial.auc.correct")) {
    stop("Cannot compute power with corrected partial AUCs")
  }
  roc1 <- roc_utils_unpercent(roc1)

  if (!missing(roc2) && !is.null(roc2)) {
    alternative <- match.arg(alternative)
    if (!is.null(sig.level) && alternative == "two.sided") {
      sig.level <- sig.level / 2
    }

    if ("roc" %in% class(roc2)) {
      # check that the AUC of roc2 was computed, or do it now
      if (is.null(roc2$auc) | !reuse.auc) {
        roc2$auc <- auc(roc2, ...)
      }
      if (!is.null(attr(roc2$auc, "partial.auc.correct")) && attr(roc2$auc, "partial.auc.correct")) {
        stop("Cannot compute power with corrected partial AUCs")
      }
      roc2 <- roc_utils_unpercent(roc2)

      # Make sure the ROC curves are paired
      rocs.are.paired <- are.paired(roc1, roc2)
      if (!rocs.are.paired) {
        stop("The sample size for a difference in AUC cannot be applied to unpaired ROC curves yet.")
      }
      # Make sure the AUC specifications are identical
      attr1 <- attributes(roc1$auc)
      attr1$roc <- NULL
      attr2 <- attributes(roc2$auc)
      attr2$roc <- NULL
      if (!identical(attr1, attr2)) {
        stop("Different AUC specifications in the ROC curves.")
      }

      # check that the same region was requested in auc. Otherwise, issue a warning
      if (!identical(attributes(roc1$auc)[names(attributes(roc1$auc)) != "roc"], attributes(roc2$auc)[names(attributes(roc2$auc)) != "roc"])) {
        warning("Different AUC specifications in the ROC curves. Enforcing the inconsistency, but unexpected results may be produced.")
      }

      ncontrols <- length(roc1$controls)
      ncases <- length(roc1$cases)
      if (is.null(kappa)) {
        kappa <- ncontrols / ncases
      }

      # Power test
      if (is.null(power)) {
        if (is.null(sig.level)) {
          stop("'sig.level' or 'power' must be provided.")
        }
        zalpha <- qnorm(1 - sig.level)
        zbeta <- zbeta.obuchowski(roc1, roc2, zalpha, method = method, ...)
        power <- pnorm(zbeta)
      }
      # sig.level
      else if (is.null(sig.level)) {
        zbeta <- qnorm(power)
        zalpha <- zalpha.obuchowski(roc1, roc2, zbeta, method = method, ...)
        sig.level <- 1 - pnorm(zalpha)
      }
      # Sample size
      else {
        zalpha <- qnorm(1 - sig.level)
        zbeta <- qnorm(power)
        ncases <- ncases.obuchowski(roc1, roc2, zalpha, zbeta, method = method, ...)
        ncontrols <- kappa * ncases
      }

      # Restore sig.level if two.sided
      if (alternative == "two.sided") {
        sig.level <- sig.level * 2
      }
      return(structure(list(ncases = ncases, ncontrols = ncontrols, auc1 = roc1$auc, auc2 = roc2$auc, sig.level = sig.level, power = power, alternative = alternative, method = "Two ROC curves power calculation"), class = "power.htest"))
    } else {
      stop("'roc2' must be an object of class 'roc'.")
    }
  } else {
    ncontrols <- length(roc1$controls)
    ncases <- length(roc1$cases)
    if (!is.null(sig.level) && !is.null(power)) {
      if (is.null(kappa)) {
        kappa <- ncontrols / ncases
      }
      ncontrols <- ncases <- NULL
    }
    auc <- auc(roc1)
    # TODO: implement this with var() and cov() for the given ROC curve
    return(power.roc.test.numeric(ncontrols = ncontrols, ncases = ncases, auc = auc, sig.level = sig.level, power = power, alternative = alternative, kappa = kappa, ...))
  }
}

power.roc.test.numeric <- function(auc = NULL, ncontrols = NULL, ncases = NULL, sig.level = 0.05, power = NULL, kappa = 1, alternative = c("two.sided", "one.sided"), ...) {
  # basic sanity checks
  if (!is.null(ncases) && ncases < 0) {
    stop("'ncases' must be positive")
  }
  if (!is.null(ncontrols) && ncontrols < 0) {
    stop("'ncontrols' must be positive")
  }
  if (!is.null(kappa) && kappa < 0) {
    stop("'kappa' must be positive")
  }
  if (!is.null(power) && (power < 0 || power > 1)) {
    stop("'power' must range from 0 to 1")
  }
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1)) {
    stop("'sig.level' must range from 0 to 1")
  }

  # Complete ncontrols and ncases with kappa
  if (is.null(ncontrols) && !is.null(ncases) && !is.null(kappa)) {
    ncontrols <- kappa * ncases
  } else if (is.null(ncases) && !is.null(ncontrols) && !is.null(kappa)) {
    ncases <- ncontrols / kappa
  }

  alternative <- match.arg(alternative)
  if (alternative == "two.sided" && !is.null(sig.level)) {
    sig.level <- sig.level / 2
  }

  # determine AUC
  if (is.null(auc)) {
    if (is.null(ncontrols) || is.null(ncases)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(power)) {
      stop("'power' or 'auc' must be provided.")
    } else if (is.null(sig.level)) {
      stop("'sig.level' or 'auc' must be provided.")
    }
    kappa <- ncontrols / ncases
    zalpha <- qnorm(1 - sig.level)
    zbeta <- qnorm(power)

    tryCatch(
      root <- uniroot(power_roc_test_optimize_auc_function, interval = c(0.5, 1 - 1e-16), ncontrols = ncontrols, ncases = ncases, zalpha = zalpha, zbeta = zbeta),
      error = function(e) {
        stop(sprintf("AUC could not be solved:\n%s", e))
      }
    )
    auc <- root$root
  }

  # Determine number of patients (sample size)
  else if (is.null(ncases) && is.null(ncontrols)) {
    if (is.null(power)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(kappa)) {
      stop("'kappa' must be provided.")
    } else if (is.null(auc)) {
      stop("'auc' or 'ncases' and 'ncontrols' must be provided.")
    } else if (is.null(sig.level)) {
      stop("'sig.level' or 'ncases' and 'ncontrols' must be provided.")
    }

    theta <- as.numeric(auc)
    Vtheta <- var_theta_obuchowski(theta, kappa)

    ncases <- solve.nd(
      zalpha = qnorm(1 - sig.level),
      zbeta = qnorm(power),
      v0 = 0.0792 * (1 + 1 / kappa),
      va = Vtheta,
      delta = theta - 0.5
    )
    ncontrols <- kappa * ncases
  }

  # Determine power
  else if (is.null(power)) {
    if (is.null(ncontrols) || is.null(ncases)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(auc)) {
      stop("'auc' or 'power' must be provided.")
    } else if (is.null(sig.level)) {
      stop("'sig.level' or 'power' must be provided.")
    }

    kappa <- ncontrols / ncases
    theta <- as.numeric(auc)
    Vtheta <- var_theta_obuchowski(theta, kappa)

    zbeta <- solve.zbeta(
      nd = ncases,
      zalpha = qnorm(1 - sig.level),
      v0 = 0.0792 * (1 + 1 / kappa),
      va = Vtheta,
      delta = theta - 0.5
    )
    power <- pnorm(zbeta)
  }

  # Determine sig.level
  else if (is.null(sig.level)) {
    if (is.null(ncontrols) || is.null(ncases)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(auc)) {
      stop("'auc' or 'sig.level' must be provided.")
    } else if (is.null(power)) {
      stop("'power' or 'sig.level' must be provided.")
    }

    kappa <- ncontrols / ncases
    theta <- as.numeric(auc)
    Vtheta <- var_theta_obuchowski(theta, kappa)

    zalpha <- solve.zalpha(
      nd = ncases,
      zbeta = qnorm(power),
      v0 = 0.0792 * (1 + 1 / kappa),
      va = Vtheta,
      delta = theta - 0.5
    )
    sig.level <- 1 - pnorm(zalpha)
  } else {
    stop("One of 'power', 'sig.level', 'auc', or both 'ncases' and 'ncontrols' must be NULL.")
  }
  # Restore sig.level if two.sided
  if (alternative == "two.sided") {
    sig.level <- sig.level * 2
  }
  return(structure(list(ncases = ncases, ncontrols = ncontrols, auc = auc, sig.level = sig.level, power = power, method = "One ROC curve power calculation"), class = "power.htest"))
}

power.roc.test.list <- function(parslist, ncontrols = NULL, ncases = NULL, sig.level = 0.05, power = NULL, kappa = 1, alternative = c("two.sided", "one.sided"), ...) {
  # basic sanity checks
  if (!is.null(ncases) && ncases < 0) {
    stop("'ncases' must be positive")
  }
  if (!is.null(ncontrols) && ncontrols < 0) {
    stop("'ncontrols' must be positive")
  }
  if (!is.null(kappa) && kappa < 0) {
    stop("'kappa' must be positive")
  }
  if (!is.null(power) && (power < 0 || power > 1)) {
    stop("'power' must range from 0 to 1")
  }
  if (!is.null(sig.level) && (sig.level < 0 || sig.level > 1)) {
    stop("'sig.level' must range from 0 to 1")
  }


  # Complete ncontrols and ncases with kappa
  if (is.null(ncontrols) && !is.null(ncases) && !is.null(kappa)) {
    ncontrols <- kappa * ncases
  } else if (is.null(ncases) && !is.null(ncontrols) && !is.null(kappa)) {
    ncases <- ncontrols / kappa
  }

  # Warn if anything is passed with ...
  if (length(list(...)) > 0) {
    warning(paste("The following arguments were ignored:", paste(names(list(...)), collapse = ", ")))
  }

  alternative <- match.arg(alternative)
  if (alternative == "two.sided" && !is.null(sig.level)) {
    sig.level <- sig.level / 2
  }

  # Check required elements of parslist
  required <- c("A1", "B1", "A2", "B2", "rn", "ra", "delta")
  if (any(!required %in% names(parslist))) {
    stop(paste("Missing parameter(s):", paste(required[!required %in% names(parslist)], collapse = ", ")))
  }

  # Determine number of patients (sample size)
  if (is.null(ncases) && is.null(ncontrols)) {
    if (is.null(power)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(kappa)) {
      stop("'kappa' must be provided.")
    } else if (is.null(sig.level)) {
      stop("'sig.level' or 'ncases' and 'ncontrols' must be provided.")
    }

    zalpha <- qnorm(1 - sig.level)
    zbeta <- qnorm(power)
    ncases <- ncases.obuchowski.params(parslist, zalpha, zbeta, kappa)
    ncontrols <- kappa * ncases
  }

  # Determine power
  else if (is.null(power)) {
    if (is.null(ncontrols) || is.null(ncases)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(sig.level)) {
      stop("'sig.level' or 'power' must be provided.")
    }
    kappa <- ncontrols / ncases

    zalpha <- qnorm(1 - sig.level)
    zbeta <- zbeta.obuchowski.params(parslist, zalpha, ncases, kappa)
    power <- pnorm(zbeta)
  }

  # Determine sig.level
  else if (is.null(sig.level)) {
    if (is.null(ncontrols) || is.null(ncases)) {
      stop("'ncontrols' and 'ncases' (or one of these with 'kappa') or 'auc' must be provided.")
    } else if (is.null(power)) {
      stop("'power' or 'sig.level' must be provided.")
    }
    kappa <- ncontrols / ncases

    zbeta <- qnorm(power)
    zalpha <- zalpha.obuchowski.params(parslist, zbeta, ncases, kappa)
    sig.level <- 1 - pnorm(zalpha)
  } else {
    stop("One of 'power', 'sig.level', 'auc', or both 'ncases' and 'ncontrols' must be NULL.")
  }
  # Restore sig.level if two.sided
  if (alternative == "two.sided") {
    sig.level <- sig.level * 2
  }
  return(structure(list(ncases = ncases, ncontrols = ncontrols, sig.level = sig.level, power = power, method = "Two ROC curves power calculation"), class = "power.htest"))
}


#### HIDDEN FUNCTIONS ####

# A function to 'optimize' auc
power_roc_test_optimize_auc_function <- function(x, ncontrols, ncases, zalpha, zbeta) {
  kappa <- ncontrols / ncases
  Vtheta <- var_theta_obuchowski(x, kappa)
  (zalpha * sqrt(0.0792 * (1 + 1 / kappa)) + zbeta * sqrt(Vtheta))^2 / (x - 0.5)^2 - ncases
}

# Compute variance of a delta from a 'covvar' list (see 'covvar' below)
var_delta_covvar <- function(covvar) {
  covvar$var1 + covvar$var2 - 2 * covvar$cov12
}

# Compute variance of a delta from a 'covvar' list (see 'covvar' below)
# under the null hypothesis
# roc1 taken as reference.
var0.delta.covvar <- function(covvar) {
  2 * covvar$var1 - 2 * covvar$cov12
}

# Compute the number of cases with Obuchowski formula and var(... method=method)
ncases.obuchowski <- function(roc1, roc2, zalpha, zbeta, method, ...) {
  delta <- roc1$auc - roc2$auc
  covvar <- covvar(roc1, roc2, method, ...)
  v0 <- var0.delta.covvar(covvar)
  va <- var_delta_covvar(covvar)
  nd <- solve.nd(
    zalpha = zalpha,
    zbeta = zbeta,
    v0 = v0, va = va,
    delta = delta
  )
  return(nd)
}

# Compute the number of cases with Obuchowski formula from params
ncases.obuchowski.params <- function(parslist, zalpha, zbeta, kappa) {
  covvar <- list(
    var1 = var_params_obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
    var2 = var_params_obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
    cov12 = cov_params_obuchowski(parslist$A1, parslist$B1, parslist$A2, parslist$B2, parslist$rn, parslist$ra, kappa, parslist$FPR11, parslist$FPR12, parslist$FPR21, parslist$FPR22)
  )
  v0 <- var0.delta.covvar(covvar)
  va <- var_delta_covvar(covvar)
  nd <- solve.nd(
    zalpha = zalpha,
    zbeta = zbeta,
    v0 = v0, va = va,
    delta = parslist$delta
  )
  return(nd)
}

# Compute the z alpha with Obuchowski formula and var(... method=method)
zalpha.obuchowski <- function(roc1, roc2, zbeta, method, ...) {
  delta <- roc1$auc - roc2$auc
  ncases <- length(roc1$cases)
  covvar <- covvar(roc1, roc2, method, ...)
  v0 <- var0.delta.covvar(covvar)
  va <- var_delta_covvar(covvar)
  zalpha <- solve.zalpha(
    nd = ncases,
    zbeta = zbeta,
    v0 = v0, va = va,
    delta = delta
  )
  return(zalpha)
}

# Compute the z alpha with Obuchowski formula from params
zalpha.obuchowski.params <- function(parslist, zbeta, ncases, kappa) {
  covvar <- list(
    var1 = var_params_obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
    var2 = var_params_obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
    cov12 = cov_params_obuchowski(parslist$A1, parslist$B1, parslist$A2, parslist$B2, parslist$rn, parslist$ra, kappa, parslist$FPR11, parslist$FPR12, parslist$FPR21, parslist$FPR22)
  )
  v0 <- var0.delta.covvar(covvar)
  va <- var_delta_covvar(covvar)
  zalpha <- solve.zalpha(
    nd = ncases,
    zbeta = zbeta,
    v0 = v0, va = va,
    delta = parslist$delta
  )
  return(zalpha)
}

# Compute the z beta with Obuchowski formula and var(... method=method)
zbeta.obuchowski <- function(roc1, roc2, zalpha, method, ...) {
  delta <- roc1$auc - roc2$auc
  ncases <- length(roc1$cases)
  covvar <- covvar(roc1, roc2, method, ...)
  v0 <- var0.delta.covvar(covvar)
  va <- var_delta_covvar(covvar)
  zbeta <- solve.zbeta(
    nd = ncases,
    zalpha = zalpha,
    v0 = v0, va = va,
    delta = delta
  )
  return(zbeta)
}

# Compute the z beta with Obuchowski formula from params
zbeta.obuchowski.params <- function(parslist, zalpha, ncases, kappa) {
  covvar <- list(
    var1 = var_params_obuchowski(parslist$A1, parslist$B1, kappa, parslist$FPR11, parslist$FPR12),
    var2 = var_params_obuchowski(parslist$A2, parslist$B2, kappa, parslist$FPR21, parslist$FPR22),
    cov12 = cov_params_obuchowski(parslist$A1, parslist$B1, parslist$A2, parslist$B2, parslist$rn, parslist$ra, kappa, parslist$FPR11, parslist$FPR12, parslist$FPR21, parslist$FPR22)
  )
  v0 <- var0.delta.covvar(covvar)
  va <- var_delta_covvar(covvar)
  a <- va
  zbeta <- solve.zbeta(
    nd = ncases,
    zalpha = zalpha,
    v0 = v0, va = va,
    delta = parslist$delta
  )
  return(zbeta)
}

solve.zbeta <- function(nd, zalpha, v0, va, delta) {
  # Solve for z_\beta in Obuchowski formula:
  # See formula 2 in Obuchowsk & McClish 1997  (2 ROC curves)
  # or formula 2 in Obuchowski et al 2004 (1 ROC curve)
  # The formula is of the form:
  # nd = (z_alpha * sqrt(v0) - z_beta * sqrt(va)) / delta ^ 2
  # Re-organized:
  # z_beta = (sqrt(nd * delta ^ 2) - z_alpha * sqrt(v0)) / sqrt(va)
  # @param nd: number of diseased patients (or abornmal, N_A in Obuchowsk & McClish 1997)
  # @param zalpha: upper \alpha (sig.level) percentile of the standard normal distribution
  # @param v0 the null variance associated with z_alpha
  # @param va: the alternative variance associated with z_beta
  # @param delta: the difference in AUC
  return((sqrt(nd * delta^2) - zalpha * sqrt(v0)) / sqrt(va))
}

solve.nd <- function(zalpha, zbeta, v0, va, delta) {
  # Solve for number of diseased (abnormal) patients in Obuchowski formula:
  # See formula 2 in Obuchowsk & McClish 1997  (2 ROC curves)
  # or formula 2 in Obuchowski et al 2004 (1 ROC curve)
  # nd = (z_alpha * sqrt(v0) - z_beta * sqrt(va)) / delta ^ 2
  # @param zalpha: upper \alpha (sig.level) percentile of the standard normal distribution
  # @param zbeta: upper \beta (power) percentile of the standard normal distribution
  # @param v0 the null variance associated with z_alpha
  # @param va: the alternative variance associated with z_beta
  # @param delta: the difference in AUC
  return((zalpha * sqrt(v0) + zbeta * sqrt(va))^2 / delta^2)
}

solve.zalpha <- function(nd, zbeta, v0, va, delta) {
  # Solve for z_\alpha in Obuchowski formula:
  # See formula 2 in Obuchowsk & McClish 1997  (2 ROC curves)
  # or formula 2 in Obuchowski et al 2004 (1 ROC curve)
  # The formula is of the form:
  # nd = (z_alpha * sqrt(v0) - z_beta * sqrt(va)) / delta ^ 2
  # Re-organized:
  # z_alpha = (sqrt(nd * delta ^ 2) - z_beta * sqrt(va)) / sqrt(v0)
  # @param nd: number of diseased patients (or abornmal, N_A in Obuchowsk & McClish 1997)
  # @param zbeta: upper \beta (power) percentile of the standard normal distribution
  # @param v0 the null variance associated with z_alpha
  # @param va: the alternative variance associated with z_beta
  # @param delta: the difference in AUC
  return((sqrt(nd * delta^2) - zbeta * sqrt(va)) / sqrt(v0))
}

# Compute var and cov of two ROC curves by bootstrap in a single bootstrap run
covvar <- function(roc1, roc2, method, ...) {
  cov12 <- cov(roc1, roc2, boot.return = TRUE, method = method, ...)
  if (!is.null(attr(cov12, "resampled.values"))) {
    var1 <- var(attr(cov12, "resampled.values")[, 1])
    var2 <- var(attr(cov12, "resampled.values")[, 2])
    attr(cov12, "resampled.values") <- NULL
  } else {
    var1 <- var(roc1, method = method, ...)
    var2 <- var(roc2, method = method, ...)
  }
  ncases <- length(roc1$cases)
  return(list(var1 = var1 * ncases, var2 = var2 * ncases, cov12 = cov12 * ncases))
}