R/deprecated/tcplFit2.R

In USEPA/CompTox-ToxCast-tcplFit2: A Concentration-Response Modeling Utility

#' Concentration-response curve fitting
#'
#' Concentration response curve fitting using the methods from BMDExpress
#'
#' All models are equal to 0 at 0 concentration (zero background).
#' To add more models in the future, write a fit____ function, and add the model
#' name to the fitmodels and modelnames vectors.
#'
#' @param conc Vector of concentrations (NOT in log units).
#' @param resp Vector of responses.
#' @param cutoff Desired cutoff. If no absolute responses > cutoff and
#'   force.fit = FALSE, will only fit constant model.
#' @param force.fit If force.fit = TRUE, will fit all models regardless of cutoff.
#' @param bidirectional If bidirectional = FALSE, will only give positive fits.
#' @param verbose If verbose = TRUE, will print optimization details and aics.
#' @param do.plot If do.plot = TRUE, will generate a plot comparing model curves.
#' @param fitmodels Vector of model names to try fitting. Missing models still
#'   return a skeleton output filled with NAs.
#' @param ... Other fitting parameters (deprecated).
#'
#' @import RColorBrewer
#' @import graphics
#' @importFrom stats median
#'
#' @return List of 11 elements. First 10 elements are the output generated by each
#'   fit function with their given model names. Last element is "modelnames": a
#'   vector of model names so other functions can easily cycle through the output.
#' @export
#'
#' @examples
#' conc = c(.03,.1,.3,1,3,10,30,100)
#' resp = c(0,.1,0,.2,.6,.9,1.1,1)
#' output = httrFit(conc,resp, .8, fitmodels = c("cnst", "hill"),verbose = TRUE,
#'   do.plot = TRUE)
tcplFit2 <- function(conc, resp, cutoff, force.fit = FALSE, bidirectional = TRUE, verbose = FALSE, do.plot = FALSE,
                    fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"),
                    ...) {

  logc = log10(conc)
  rmds <- tapply(resp, logc, median)

  #first decide which of possible models will be fit
  modelnames = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5")
  #decide whether to run each model, then use generic functions to run model by name
  for(model in modelnames){
    #only fit when four or more concentrations, the model is in fitmodels, and
    #( either one response is above cutoff OR force.fit == T OR it's the constant model.)
    to.fit = (length(rmds) >= 4 && model %in% fitmodels && (length(which(abs(rmds) >= cutoff)) > 0 || force.fit ||
                                                                     model == "cnst") )
    fname = paste0("fit",model) #requires each model function have name "fit____" where ____ is the model name
    #use do.call to call fit function; cnst has different inputs than others.
    if(model == "cnst") assign(model, do.call(fname, list(conc = conc, resp = resp, nofit = !to.fit))) else{
      assign(model, do.call(fname, list(conc = conc, resp = resp, bidirectional = bidirectional, verbose = verbose,
                                        nofit = !to.fit)))
    }
  }
  #optionally print out AICs
  if(verbose) {
    print("aic values:")
    aics = sapply(modelnames, function(x){get(x)[["aic"]]})
    names(aics) = modelnames
    print(aics)
    cat("Winner: ", modelnames[which.min(aics)])
  }

  #optionallky plot all models if there's at least one model to plot
  shortnames = modelnames[modelnames != "cnst"]
  successes = sapply(shortnames, function(x){get(x)[["success"]]})
  if(do.plot && sum(successes, na.rm = T) == length(shortnames)){
    resp = resp[order(logc)]
    #par(xpd = T)
    cols = c("black",brewer.pal(9,"Set1"))
    n = length(logc)
    allresp =  c(resp ,sapply(shortnames, function(x){get(x)[["modl"]][order(logc)]}))
    logc = logc[order(logc)]
    plot(rep(logc,length.out =length(allresp)), allresp, col = rep(cols,each = n), pch = 16)

    for(i in 1:length(allresp)){
      points(logc,allresp[((i-1)*n + 1):(i*n)],col = cols[i], type = "l")
    }

    legend("top", legend = c("resp", shortnames), col = cols, pch = 16, ncol = 10, inset = c(0,-.1))
  }

  #put all the model outputs into one list and return
  out <- c(mget(modelnames),
           list(modelnames  = modelnames, ...))

  return(out)
}