R/predictIC50.R

In lujunyan1118/DrugScreenExplorer_dev: A pipeline for processing drug sensitivity screen data

#' Predict responses based on input dose information and fitted log-logistic model
#'
#' This function predicts the response based on input fitted log-logistic model and concentration values.
#' @param object an object of class 'fitIC50', generated by fitIC50 function
#' @param newdata a numeric vector or data fram with one column that contains the new concentration values
#' If omitted, the concentrations used to fit the model are used.
#' @param se.fit not used, mainly for geom_smooth purpose
#' @param level not used, mainly for geom_smooth purpose
#' @param interval not used, mainly for geom_smooth purpose
#' @param ... further arguments to be passed to ggplot when using geom_smooth().
#' @export
#' @return a vector of predicted response based on a fitted model and input concentration values.
#' @examples
#' # create an example dose-response table
#' doseTab <- data.frame(viability = c(1.0, 0.95, 0.5, 0.2, 0.1),
#'                       concentration = c(0.001, 0.01, 0.1, 1, 10))
#'
#' # fit a four-parameter logistic model using dr4pl package
#' mm <- fitIC50(viability ~ concentration, data = doseTab)
#'
#' # Predict response using the input concentrations
#' predict(mm)
#'
#' # Predict response using new concentrations
#' predict(mm, newdata = c(0.25, 2.5, 5, 15))

predict.fitIC50 <- function(object, newdata = NULL, se.fit = FALSE, level = 0.95,
                            interval = c("none", "confidence", "prediction"), ...) {

  if (is.null(newdata))
    newdata <- object$model else newdata <- newdata

    params <- object$parm_fit$parameters
    logDose <- object$logDose

    a <- min(params[1], params[4])
    d <- max(params[4], params[1])
    c <- params[2]
    b <- params[3]

    predY <- function(x) {
      y = d + (a - d)/(1 + (x/c)^b)
      names(y) <- NULL
      return(y)
    }

    if (is.vector(newdata)) {
      conc <- newdata
    } else if (is.data.frame(newdata)) {
      if (ncol(newdata) == 1) {
        conc <- newdata[, 1]
      } else {
        concName <- as.character(object$formula)[3]
        conc <- newdata[, concName]
      }
    }

    if (!is.null(logDose))
      conc <- logDose^conc
    res <- vapply(conc, predY, numeric(1))

    return(res)
}