R/sigmoid.R

In hemoshear/assayr2: Tools for analyzing assay data

#' Make a robust DRC estimate using a method based on Nguyen et al. 2014
#'
#' @param object An object of class drc fit with \code{drm()} and the argument \code{fct = LL.4()}.
#' @param form A formula associated with the object.
#' @param conv Numeric value representing the convergence criterion (manhattan distance between weight vectors).
#' @param maxits Numeric value representing maximum number of iterations.
#' @param verbose Logical. If TRUE reports the manhattan distance for each iteration.
#' @param drm_error Logical. Determines if drc() convergence failure results in a error (TRUE) or warning (FALSE)
#' @examples
#' mtcars2 <- mtcars
#' mtcars2$wt[c(10:12)] <- mtcars2$wt[c(10:12)] * 2
#' object <- drc::drm(disp~wt, data=mtcars2, fct=drc::LL.4())
#' robustify_drc(object, disp ~ wt, verbose = TRUE)
#' @importFrom stats formula resid weights
#' @export
robustify_drc <- function(object, form, conv=0.01, maxits=100, verbose=FALSE,
                         drm_error=FALSE) {
  environment(form) <- environment()
  d <- object$data
  ow <- object$weights
  d$weights <- .biweightMean(resid(object))
  mdist <- sum(abs(ow - d$weights))
  if (verbose) {
    print(paste0("Manhattan distance = ", mdist))
  }
  if (mdist > conv & maxits > 0) {
    maxits <- maxits - 1
    # fit <- drc::drm(formula, weights=weights, data=d, fct=drc::LL.4(),
    #                 control = drc::drmc(errorm=drm_error, useD=deriv))
    object <- .tryFit(form, data=d, w=d$weights, drm_error=drm_error)
    robustify_drc(object, form, conv=conv, maxits=maxits, verbose = verbose)
  }
  else { return(object) }
}

.tryFit <- function(form, data, w=NULL, drm_error=FALSE) {
  if (is.null(w)) {
    w <- rep(1, nrow(data))
  }
  environment(form) <- environment()
  object <- try(drc::drm(form, data = data, weights = w, fct = drc::LL.4(),
                      control = drc::drmc(errorm = drm_error)), silent = TRUE)
  if (class(object) == "list" | class(object) == "try-error") {
    object <- drc::drm(form, data = data, weights = w, fct = drc::LL.4(),
                    control = drc::drmc(errorm = drm_error, useD = TRUE))
  }
  return(object)
}
# attempt to remove the 'form' arguement and pull it directly from 'object'
drc_robust <- function(object,
                       # form, 
                       conv = 0.000001, maxits = 100, verbose = FALSE,
                         drm_error = FALSE, useD = FALSE) {
  
  form <- formula(object) # trying to figure a way to remove the 'form' argument
  environment(form) <- environment() # bc formula objects point to their original calling environment

  d <- object$data
  d$weights <- .biweightMean(resid(object)) # down weight outliers via Tukey's bi-weight mean
  avg_weight_delta <- sum(abs(object$weights - d$weights)) / nrow(d) # average difference between weights
  print(avg_weight_delta)
  if (verbose) {
    print(paste0("Average manhattan distance between weights = ", avg_weight_delta))
  }
  if ( (avg_weight_delta > conv) & (maxits > 0) ) {
    print("enter")
    maxits <- maxits - 1
    # object <- .tryFit(form, data = d, w = d$weights, drm_error = drm_error)
    object <- try(
      drc::drm(form, data = d, weights = weights, fct = drc::LL.4(),
             control = drc::drmc(useD = useD)
             ),
      silent = TRUE
    )
    if (class(object) == "list" | class(object) == "try-error") { # come back to class(object) ==  "list"; not sure why that is included
      stop(glue::glue("Error: drc::drm() failed to converge, try setting useD = {!useD}"))
    }
    drc_robust(object, form, conv = conv, maxits = maxits, verbose = verbose, useD = useD)
  }
  else {
    return(object)
  }
}

# delete contigent on buidling a test taht breaks the function abo
# .tryFit <- function(form, data, w = NULL, drm_error = FALSE) {
#   if (is.null(w)) {
#     w <- rep(1, nrow(data))
#   }
#   environment(form) <- environment()
#   fit <- try(
#     drc::drm(form, data = data, weights = w,
#              fct = drc::LL.4(),
#              control = drc::drmc(errorm = drm_error)),
#     silent = TRUE)
#   if (class(fit) == "list" | class(fit) == "try-error") {
#     fit <- drc::drm(form,
#       data = data, weights = w, fct = drc::LL.4(),
#       control = drc::drmc(errorm = drm_error, useD = TRUE)
#     )
#     warning("just kidding, it worked")
#   }
#   return(fit)
# }

# export the reference for the paper in package documentation

.biweightMean <- function(r) { # weighting scheme from Nguyen et al. 2014
  mr <- mean(abs(r))
  lim <- 6 * mr
  w <- sapply(r, function(x) {
    ifelse(abs(x) < lim, (1 - (x / lim)^2)^2, 0)
  })
  return(w)
}

#' Get the y value of the upper asymptote in a 4-parmeter log-logistic function
#'
#' @param object An object of class "drc" generated from \code{drm()} and the argument \code{fct = LL.4()}.
#' @param CI95 Logical. Report 95\% confidence interval.
#' @return The y intercept of the upper asymptote.
#' @examples
#' fit <- drc::drm(disp~wt, data=mtcars, fct=drc::LL.4())
#' getUpperAsym(fit)
#' @export
getUpperAsym <- function(object, CI95 = FALSE) {
  res <- cbind(confint(object), object$coefficients)
  colnames(res)[3] <- "estimate"
  res <- res[c(
    "c:(Intercept)",
    "d:(Intercept)"
  ), ]
  ind <- which.max(res[, 3])
  res <- res[ind, ]
  if (!CI95) {
    res <- res["estimate"]
  }
  return(res)
}

#' Get the y value of the lower asymptote in a 4-parmeter log-logistic function
#'
#' @param object An object of class "drc" generated from \code{drm()} and the argument \code{fct = LL.4()}.
#' @param CI95 Logical. Report 95\% confidence interval.
#' @return The y intercept of the lower asymptote.
#' @examples
#' fit <- drc::drm(disp~wt, data=mtcars, fct=drc::LL.4())
#' getLowerAsym(fit)
#' @export
getLowerAsym <- function(object, CI95 = FALSE) {
  res <- cbind(confint(object), object$coefficients)
  colnames(res)[3] <- "estimate"
  res <- res[c(
    "c:(Intercept)",
    "d:(Intercept)"
  ), ]
  ind <- which.min(res[, 3])
  res <- res[ind, ]
  if (!CI95) {
    res <- res["estimate"]
  }
  return(res)
}

#' Get the EC50 value of a 4-parmeter log-logistic function
#'
#' @param object An object of class "drc" generated from \code{drm()} and the argument \code{fct = LL.4()}.
#' @param CI95 Logical. Report 95\% confidence interval.
#' @return The x value corresponding to point at which the function is rotationally symmetric.
#' @examples
#' fit <- drc::drm(disp~wt, data=mtcars, fct=drc::LL.4())
#' getEC50(fit)
#' @export
getEC50 <- function(object, CI95 = FALSE) {
  res <- cbind(confint(object), object$coefficients)
  colnames(res)[3] <- "estimate"
  res <- res["e:(Intercept)", ]
  if (!CI95) {
    res <- res["estimate"]
  }
  return(res)
}

#' Get the Hill coefficient of a 4-parmeter log-logistic function
#'
#' @param object An object of class "drc" generated from \code{drm()} and the argument \code{fct = LL.4()}.
#' @param CI95 Logical. Report 95\% confidence interval.
#' @return The Hill coefficient:
#' \eqn{y=c+\frac{d-c}{1+10^{(ln(\tilde{e})-x) \times Hill}}}.
#' Where c = upper asymptote, d = lower asymptote, and \eqn{\tilde{e}} = the EC50 value.
#' @examples
#' fit <- drc::drm(disp~wt, data=mtcars, fct=drc::LL.4())
#' getEC50(fit)
#' @export
getHillSlope <- function(object, CI95 = FALSE) {
  res <- cbind(confint(object), object$coefficients)
  colnames(res)[3] <- "estimate"
  b <- res["b:(Intercept)", ]
  if (res[3, 3] > res[2, 3]) {
    hill <- (-1 * b) / log(10)
  } else {
    hill <- b / log(10)
  }
  if (!CI95) {
    hill <- hill["estimate"]
  }
  return(hill)
}

#' Generate a data frame with y values in increments of a given x range.
#'
#' Useful for plotting fitted curves.
#'
#' @param object An model object with a \code{predict()} method generated from \code{drm()} and the argument \code{fct = LL.4()}.
#' @param x_range An optional 2-element numeric vector giving the x-range. If `NULL`, will
#' calculate range from `object`.
#' @param logrange Logical. If TRUE increments are in log space.
#' @param npoint Numeric. The number of increments in the x range.
#' @return A data frame with predicted x and y values pulled from teh model call.
#' @examples
#' fit_drc <- drc::drm(disp~wt, data=mtcars, fct=drc::LL.4())
#' curve_drc <- get_curves(fit_drc, range(mtcars$wt), logrange = FALSE)
#' plot(wt~disp, data=curve_drc) # names come from model call
#'
#' fit_lm <- lm(disp~wt, data = mtcars)
#' curve_lm <- get_curves(fit_lm, range(mtcars$wt), logrange = FALSE)
#' plot(wt~disp, data=curve_lm)
#' @importFrom stats predict
#' @export
get_curves <- function(object, x_range = NULL, logrange = TRUE, npoint = 100) {
  UseMethod("get_curves")
}

#' @export
get_curves.drc <- function(object, x_range = NULL, logrange = TRUE, npoint = 100) {
  # this is the only  difference between lm and drc, where the data is stored
  # get column names
  x_name <- names(object$data)[1]
  y_name <- names(object$data)[2]
  if (is.null(x_range)) {
    x_range <- range(object$data[x_name])
  }

  # core is the same
  x_range %<>% sort()

  if (logrange) {
    x_range <- log(x_range)
  }
  xs <- seq(x_range[1], x_range[2], length.out = npoint)

  if (logrange) {
    xs <- exp(xs)
  }

  newdata <- as.data.frame(xs)
  ys <- predict(object, newdata)

  data.frame(x = xs, y = ys) %>% setNames(c(eval(x_name), eval(y_name)))
}
#' @export
get_curves.lm <- function(object, x_range = NULL, logrange = TRUE, npoint = 100) {
  # this is the only  difference between lm and drc, where the data is stored
  # get column names
  x_name <- names(object$model)[2]
  y_name <- names(object$model)[1]
  # get range data form model object
  if (is.null(x_range)) {
    x_range <- range(object$model[x_name])
  }

  # core is the same
  x_range %<>% sort()

  if (logrange) {
    x_range <- log(x_range)
  }
  xs <- seq(x_range[1], x_range[2], length.out = npoint)

  if (logrange) {
    xs <- exp(xs)
  }

  newdata <- as.data.frame(xs) %>%
    setNames(eval(x_name))
  ys <- predict(object, newdata)

  data.frame(x = xs, y = ys) %>% setNames(c(eval(x_name), eval(y_name)))
}

#' Get the bottom and top of a log-logistic curve in a given x range.
#'
#' @param object An object of class "drc" generated from \code{drm()} and the argument \code{fct = LL.4()}.
#' @param high_x Numeric. The upper bound of the x range.
#' @param low_x Numeric. The lower bound of the x range.
#' @return A named vector reporting the top and bottom of the curve (toc, boc)
#' @examples
#' fit <- drc::drm(disp~wt, data=mtcars, fct=drc::LL.4())
#' getYBounds(fit, 2, 5)
#' @importFrom stats predict
#' @export
getYBounds <- function(object, high_x, low_x) {
  y1 <- predict(object, data.frame(high_x))
  y2 <- predict(object, data.frame(low_x))
  ys <- sort(c(y1, y2))
  res <- c(boc = ys[1], toc = ys[2])
  return(res)
}