#' 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 places
#' 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{getTiter}} to derive the titer and includes
#' 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"), axes = TRUE,
xlab = NULL, ylab = NULL, frame.plot = axes, digits = 2, ...)
{
# 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) && !is(fm, "glm"))
main <- names(fm)
else 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 <- 1/getTiter(mod, level = NULL)
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/" * .(units))
else
xlab <- "Titer = NaN"
}
else
xlab <- sprintf(sprintf("One IU = %%0.%dg", digits), xpp)
}
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, digits, ...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.