R/plotFit.R

Defines functions plotFit

Documented in plotFit

#' Plot Titer Fit from GLM Model
#' 
#' Prepare an informative plot of viral titer from the GLM fit.
#' 
#' @param fm fitted model or list of fitted models from \code{\link{getFit}}.
#' @param main plot title as a character, vector of characters or list of
#'  characters to be used as the plot title(s). Values will be replicated as
#'  necessary. If \code{NULL}, the plot date will be used.
#' @param pch,col,col.pch plot character and color for data points. The default
#'  color for the points and the fitted line is \code{col} but these can be
#'  changed with the appropriate argument. \code{col} will also be used as the
#'  default color for other plot elements.
#' @param lty.fit,col.fit line type and color for the GLM best-fit line.  Set
#'  \code{lty.fit} to \code{NA} to exclude the best-fit line.
#' @param lty.ref,col.ref line type and color for the value on the x-axis
#'  that intersect the 63% value on the best-fit line. Set \code{lty.ref} to
#'  \code{NA} to exclude the reference line.
#' @param xlim,ylim optional parameters to override the default parameters of
#'  the full range of x values for \code{xlim} and \code{c(0, 1)} for
#'  \code{ylim}.
#' @param ann a logical value indicating whether the default annotation (title
#'  and x and y axis labels) should appear on the plot.
#' @param xlab a label for the x axis. If the \code{unit} attribute associated
#'  with \code{fm} is a volume (ml, ul, nl, pl or fl), the default will be the
#'  titer as IU per ml. Otherwise, the default will be the x-value
#'  corresponding to an MOI of 1 (63\% point). 
#' @param ylab a label for the y axis, defaults to "Infected fraction".
#' @param axes a logical value indicating whether both axes should be drawn on
#'  the plot. 
#' @param frame.plot a logical value indicating whether a box should be drawn
#'  around the plot.
#' @param digits integer indicating the number of significant decimal palces 
#'  to use for the calculated x-axis label
#' @param ... arguments to be passed to \code{\link[graphics]{plot.default}}
#' 
#' @details
#' 
#' Base graphics are used to prepare a plot with \code{log = "x"} from a single 
#' fitted model or a list of fitted models generated by \code{\link{getFit}}. 
#' The function calls \code{\link{getEC63}} to derive the titer and include 
#' this information on the plot by default. The graphic parameters can be 
#' used to change or exclude many elements of the plot. See the help 
#' information in \code{\link[graphics]{plot.default}} for details.
#' 
#' If \code{fm} is a list of fitted models, each will be plotted. Use 
#' \code{par(ask = TRUE)} to see each in turn or use 
#' \code{par(mfrow = c(nr, nc))} to place \code{nr x nc} plots on one device.
#' 
#' Additional curves can be added with the \code{\link{addOneFit}} as shown in 
#' the examples.
#' 
#' @return
#' 
#' No value is returned. This function is called for the side-effect of 
#' producing a plot (or plots).
#' 
#' @examples
#' # sample result from tally() with two sets of data
#'   plotFit_data <- read.csv(system.file("extdata", "plotFit_data.csv", 
#'     package = "virustiter"))
#'   fm <- getFit(plotFit_data, by = "row")
#'
#' # Default format and example with minimal annotation
#'   plotFit(fm$A)
#'   plotFit(fm$A, ann = FALSE, axes = FALSE, lty.ref = NA)
#'
#' # Example of two plots
#'   plotFit(fm$A, "Two fits", xlab = "Multiplicity (ml/cell)")
#'   addOneFit(fm$B)
#'   legend("topleft", legend = c("A", "B"), lty = 1, col = c(1, 2), bty = "n")
#'
#' # With two panel figure 
#'   par(mfrow = c(2, 1))
#'   plotFit(fm, main = names(fm))
#'
#' @export
#'
plotFit <- function(fm, main = NULL, pch = 1, col = 1, col.pch = col,
	lty.fit = 1, col.fit = col, lty.ref = 2, col.ref = "gray",
	xlim = NULL, ylim = NULL, ann = par("ann"), xlab = NULL, ylab = NULL,
	axes = TRUE, frame.plot = axes, digits = 3, ...) 
{
  # argument check
    if (!(is(fm, "glm") | all(sapply(fm, is, "glm"))))
      stop("requires a glm model generated by getFit()")

  # ensure main is appropriate and pad with NULL if need be
    if (is.null(main)) main <- list(NULL)
    main <- as.list(main)
    main <- c(main, rep(list(NULL), max(0, length(fm) - length(main))))
    main <- main[seq_along(fm)]

  # working function to dispatch on each glm model
    .plotFit <- function(mod, main = main, pch = pch, col = col,
    	col.pch = col.pch, lty.fit = lty.fit, col.fit = col.fit, lty.ref = lty.ref,
    	col.ref = col.ref, xlim = xlim, ylim = ylim, ann = ann, xlab = xlab,
    	ylab = ylab, axes = axes, frame.plot = frame.plot, digits = digits, ...)
    {
    # generate plot title
      if (is.null(main[[1]])) main <- paste("Plot generated", Sys.Date())

    # local data and helper functions
      x <- exp(mod$model[[2]])
      y <- prop.table(mod$model[[1]], 1)[,1]
      units <- attr(mod, "unit")[1]
      localTitle <- function(..., col, bg, pch, cex, lty, lwd) title(...)
      localAxis <- function(..., col, bg, pch, cex, lty, lwd) Axis(...)
      localBox <- function(..., col, bg, pch, cex, lty, lwd) box(...)

    # glm fit line
      xlo <- min(x[x > 0]); xhi <- max(x)
      xp <- exp(seq(log(xlo), log(xhi), length = 101))
      yp <- predict(mod, data.frame(x = xp), type = "response")

    # reference lines
      xpp <- getEC63(mod)[1]
      ypp <- 1 - exp(-1)

    # annotation
      if (is.null(xlab)) {
        if (grepl("[munpf]l", units, ignore.case = TRUE)) { # units are volume
          val <- 1/xpp
          expt <- floor(log10(abs(val)))
          mant <- round(val/10^expt, digits)
          if (is.finite(xpp))
            xlab <- bquote("Titer = " * .(mant) %*% 10^.(expt) * " IU/ml")
          else
            xlab <- "Titer = NaN"
        }
        else
          xlab <- sprintf("One IU = %0.3g %s", xpp, units)
      }
      if (is.null(ylab)) ylab <- "Infected fraction"
      if (is.null(xlim)) xlim <- range(x[is.finite(x)])
      if (is.null(ylim)) ylim <- c(0, 1)
      
    # plot with base graphics
      plot(y ~ x, subset = x > 0, log = "x", xlim = xlim, ylim = ylim,
        xlab = "", ylab = "", axes = FALSE, ann = FALSE, frame.plot = FALSE,
        pch = pch, col = col.pch,  ...)
      if (!is.na(lty.fit)) lines(xp, yp, lty = lty.fit, col = col.fit)
      if (!is.na(lty.ref)) lines(c(xlo, xpp, xpp), c(ypp, ypp, ylim[1]),
        lty = lty.ref, col = col.ref)
      if (axes) localAxis(x, side = 1, ...)
      if (axes) localAxis(y, side = 2, las = 1, ...)
      if (frame.plot) localBox(...)
      if (ann) localTitle(main = main, xlab = xlab, ylab = ylab, ...)
    }
  # dispatch function
    if(is(fm, "list")) {
      ll <- function(x) as.list(list(x)) # prevents zero-length arguments
      junk <- Map(.plotFit, fm, main = main, pch = ll(pch), col = ll(col),
      	col.pch = ll(col.pch), lty.fit = ll(lty.fit), col.fit = ll(col.fit),
      	lty.ref = ll(lty.ref), col.ref = ll(col.ref),
        xlim = ll(xlim), ylim = ll(ylim), ann = ll(ann), xlab = ll(xlab),
        ylab = ll(ylab), axes = ll(axes), frame.plot = ll(frame.plot),
        digits = ll(digits), ...)
    }
    else 
      junk <- .plotFit(fm, main, pch, col, col.pch, lty.fit, col.fit, lty.ref,
      	col.ref, xlim, ylim, ann, xlab, ylab, axes, frame.plot, ...)
}
ornelles/virustiter documentation built on March 15, 2024, 9:28 a.m.