R/linearQuadraticModel.R

In RadioGx: Analysis of Large-Scale Radio-Genomic Data

#' Fit linear-quadratic curves to dose-response data
#'
#' This function fits a linear-quadratic curve to dose-response data.
#'
#' @examples
#' linearQuadraticModel(c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
#'  c(1.1, 0.8, 0.7, 0.45, 0.15, -0.1, -0.1, -0.4, -0.65, -0.75, -1.1))
#'
#' @param D vector of radiation doses
#' @param SF vector of survival fractions corresponding to the doses
#' @param lower_bounds vector of length 2 containing minimum allowed values of
#'   fitted alpha and beta, respectively
#' @param upper_bounds vector of length 2 containing maximum allowed values of
#'   fitted alpha and beta, respectively
#' @param scale parameter of the assumed error distribution of the data; see
#'   sdetails
#' @param family family of distributions of the error terms in the data;
#'   currently supported options are "normal" and "cauchy"
#' @param median_n see details
#' @param trunc should survival fractions be truncated downward to 1? Defaults
#'   to FALSE.
#' @param verbose 'verbose' outputs warnings that are otherwised suppressed when
#'   the function sanity-checks user inputs. 'median_n' denotes the number of
#'   distributions from family 'family' that are medianned. (Note that setting
#'   n = 1 (the default) is equivalent to using a simple normal or cauchy
#'   distribution without taking any medians.)
#'
#' @return \code{numeric} The estimated alpha and beta values
#' @export
linearQuadraticModel <- function (D,
                                  SF,
                                  lower_bounds = c(0, 0),
                                  upper_bounds = c(1, 1),
                                  scale = 5,
                                  family = c("normal", "Cauchy"),
                                  median_n = 1,
                                  trunc = FALSE,
                                  verbose = FALSE) {
  match.arg(family)

  CoreGx::.sanitizeInput(x = D,
                         y = SF,
                         x_as_log = FALSE,
                         y_as_log = FALSE,
                         y_as_pct = FALSE,
                         trunc = trunc,
                         verbose = verbose)

  DSF <- CoreGx::.reformatData(x = D,
                               y = SF,
                               x_to_log = FALSE,
                               y_to_log = TRUE,
                               y_to_frac = FALSE,
                               trunc = trunc)
  D <- DSF[["x"]]
  SF <- DSF[["y"]]

  if (!(all(lower_bounds < upper_bounds))) {
    if (verbose == 2) {
      message("lower_bounds:")
      message(lower_bounds)
      message("upper_bounds:")
      message(upper_bounds)
    }
    stop ("All lower bounds must be less than the corresponding upper_bounds.")
  }

  if(!(0 %in% D) || SF[D==0] != 0) {
    D <- c(0,D)
    SF <- c(0,SF)
  }

  gritty_guess <- .makeGrittyGuess(lower_bounds = lower_bounds,
                                   upper_bounds = upper_bounds,
                                   D = D,
                                   SF = SF)

  guess <- CoreGx::.fitCurve(x = D,
                              y = SF,
                              f = .linearQuadratic,
                              density = c(100, 100),
                              step = c(0.005, 0.005),
                              precision = 0.005,
                              lower_bounds = lower_bounds,
                              upper_bounds = upper_bounds,
                              scale = scale,
                              family = family,
                              median_n = median_n,
                              trunc = FALSE,
                              verbose = verbose,
                              gritty_guess = gritty_guess,
                              span = 0.1)

  names(guess) <- c("alpha", "beta")

  return(guess)
}