R/predict.nls.R

In Novartis/xgxr: Exploratory Graphics for Pharmacometrics

#' predict.nls
#'
#' @param object Object of class inheriting from "nls"
#' @param newdata An optional data frame in which to look for variables with which to predict.
#' If omitted, the fitted values are used.
#' @param se.fit A switch indicating if standard errors are required.
#' @param interval Type of interval calculation, "none" or "confidence"
#' @param level Level of confidence interval to use
#' @param ... additional arguments affecting the predictions produced.
#'
#' @return \code{predict.nls} produces a vector of predictions or a matrix of predictions and
#' bounds with column names \code{fit}, \code{lwr}, and \code{upr} if interval is set.
#'
#' If \code{se.fit} is \code{TRUE}, a list with the following components is returned:
#'
#' \item{fit}{vector or matrix as above}
#'
#' \item{se.fit}{standard error of predicted means}
#'
#' \item{residual.scale}{residual standard deviations}
#'
#' \item{df}{degrees of freedom for residual}
#'
#' @examples
#'
#' set.seed(12345)
#' data_to_plot <- data.frame(x1 = rep(c(0, 25, 50, 100, 200, 400, 600), 10)) %>%
#'   dplyr::mutate(AUC = x1*rlnorm(length(x1), 0, 0.3),
#'          x2 = x1*stats::rlnorm(length(x1), 0, 0.3),
#'          Response = (15 + 50*x2/(20+x2))*stats::rlnorm(length(x2), 0, 0.3))
#'
#'
#' gg <- ggplot2::ggplot(data = data_to_plot, ggplot2::aes(x = AUC, y = Response)) +
#'   ggplot2::geom_point() +
#'   xgx_geom_smooth(method = "nls",
#'                   method.args = list(formula = y ~ E0 + Emax* x / (EC50 + x),
#'                                      start = list(E0 = 15, Emax = 50, EC50 = 20) ),
#'                   color = "black", size = 0.5, alpha = 0.25)
#' gg
#'
#' mod <- stats::nls(formula = Response ~ E0 + Emax * AUC / (EC50 + AUC),
#' data = data_to_plot,
#' start = list(E0 = 15, Emax = 50, EC50 = 20))
#'
#' predict(mod)
#'
#' predict(mod, se.fit = TRUE)
#'
#' predict(mod,
#'             newdata = data.frame(AUC = c(0, 25, 50, 100, 200, 400, 600)),
#'             se.fit = TRUE)
#'
#' predict(mod,
#'             newdata = data.frame(AUC = c(0, 25, 50, 100, 200, 400, 600)),
#'             se.fit = TRUE, interval = "confidence", level = 0.95)
#'
#' predict(mod,
#'             newdata = data.frame(AUC = c(0, 25, 50, 100, 200, 400, 600)),
#'             se.fit = TRUE, interval = "confidence", level = 0.95)
#'
#' @importFrom Deriv Deriv
#' @importFrom stats nls
predict.nls <- function(object, newdata = NULL, se.fit = FALSE, interval = "none", level = 0.95, ...){

  pred <- list()

  # function to calculate gradient wrt model parameters
  # value is the function value
  # grad is the gradient
  fun_grad <- function(form, newdata, pars, se.fit){

    # extract the model parameters to the local environment
    list2env(pars %>% as.list(), envir = environment())

    ret <- list()

    if(se.fit){
      ret$grad <- list()
    }


    for(i in 1:length(newdata[,1])){

      if(length(newdata[1,]) > 1){
        for(j in names(newdata[i,])){
          assign(j, newdata[i,j])
        }
      }else{
        j = names(newdata[1])
        assign(j, newdata[i,j])
      }

      ret$value[i] <- eval(form[[3L]]) # this is the value of the formula

      if(se.fit){
        ret$grad[[i]] <- eval(Deriv::Deriv(form, names(pars), cache.exp = FALSE)) %>% as.list()

        if(is.null(names(ret$grad[[i]]))){
          names(ret$grad[[i]]) <- names(pars)
        }
      }
    }


    if(se.fit){
      ret$grad <- dplyr::bind_rows(ret$grad) %>% as.matrix
    }

    return(ret)
  }

  if(is.null(newdata)){
    fg <- list()
    fg$value <- as.numeric(object$m$fitted())
    fg$grad <- object$m$gradient()

  }else{
    fg <- fun_grad(form = object$m$formula(), newdata, pars = object$m$getPars(), se.fit)

  }


  f.new <- fg$value # value of function

  if(se.fit){
    pred$fit <- f.new

    grad.new <- fg$grad # value of gradient

    vcov <- vcov(object)
    GS = rowSums((grad.new%*%vcov)*grad.new)

    if(interval == "confidence"){

      alpha = 1 - level
      deltaf <- sqrt(GS)*qt(1 - alpha/2, df = summary(object)$df[2])

      pred$fit <- data.frame(fit = pred$fit)
      pred$fit$lwr <- f.new - deltaf
      pred$fit$upr <- f.new + deltaf
    }

    pred$se.fit <- sqrt(GS)
    pred$df <- summary(object)$df[2]

  }else{
    pred <- f.new
  }

  return(pred)
}