#' Fit Infectivity Results with Binomial GLM Model
#'
#' Fit the results in a \code{data.frame} produced by either \code{\link{tally}} or
#' \code{\link{score}}.
#'
#' @param obj A \code{data.frame} produced by \code{\link{tally}} or by \code{\link{score}}.
#' @param by A character string indicating the organization of the results
#' (see details.)
#'
#' @details
#' This function will fit a GLM binomial model using a two-column response
#' where the fit is weighted by the total numbers of cases. Data are fit by
#' the \code{\link{glm}} function using a complementary log-log model. The data
#' can be organized as a single, sequential series of values \emph{or} by
#' rows \emph{or} columns. If the argument \code{by} is \code{"sequential"},
#' all the data are assumed to represent a single assay. If \code{by} is
#' either \code{"row"} or \code{"column"}, fits will be generated for data
#' partitioned by row or column, respectively.
#'
#' The argument \code{obj} must be produced either by \code{\link{tally}} or by
#' \code{\link{score}}. If a variable named \code{"positive"} is present, the data
#' are assumed to have been produced by \code{\link{score}} and the function
#' \code{\link{tally}} will be called before proceeding. If variables named
#' \code{"pos"} and \code{"neg"} are present, the data are assumed to be the
#' product of \code{\link{tally}}. In the case of a data frame produced by
#' \code{\link{score}}, the data must have a variable named \code{"moi"} or \code{"x"}
#' representing the multiplicity of infection. Data produced by \code{\link{tally}}
#' will always have a variable named \code{"x"} representing the multiplicity.
#'
#' @return
#'
#' A fitted \code{glm} model or list of models. Each fitted model has the unit
#' of measure attached as an attribute with the text value from \code{unit}.
#'
#' @export
#'
getFit <- function(obj, by = c("sequential", "column", "row"))
{
if (missing(obj))
stop("getFit requires a data.frame from tally() or score()")
if (all(c("pos", "neg") %in% names(obj))) # data.frame produced by tally()
res <- obj
else if ("positive" %in% names(obj)) # data.frame from score()
res <- tally(obj)
else
stop(deparse(substitute(obj)), " has not been generated by tally() or score()")
# confirm availability of moi variable (change name for fitting purposes...)
if (!"x" %in% names(res))
stop(deparse(substitute(obj)), " requires either 'moi' or 'x' as the multiplicity")
if (!is.numeric(res$x))
stop("unable to use non-numeric values for 'x'")
# check for unit string
if ("unit" %in% names(res))
unit <- as.character(res$unit[1])
else
unit <- "none"
# determine how proceed from 'by'
by <- match.arg(by)
fmList <- list()
if (by == "sequential") {
form <- as.formula(cbind(pos, neg) ~ offset(log(x)))
fm <- try(glm(form, data = res, subset = res$x > 0, family = binomial("cloglog")))
if (class(fm)[1] == "try-error")
fmList[[1]] <- NULL
else
fmList[[1]] <- fm
}
else if (by == "column") {
stopifnot("column" %in% names(res))
for (i in levels(as.factor(res[[by]]))) {
form <- as.formula(cbind(pos, neg) ~ offset(log(x)))
fm <- try(glm(form, family = binomial("cloglog"), data = res,
subset = res[[by]] == i & res$x > 0), silent = TRUE)
if (class(fm)[1] == "try-error")
fmList[[i]] <- NULL
else
fmList[[i]] <- fm
}
}
else { # by == "row"
stopifnot("row" %in% names(res))
for (i in levels(as.factor(res[[by]]))) {
form <- as.formula(cbind(pos, neg) ~ offset(log(x)))
fm <- try(glm(form, family = binomial("cloglog"), data = res,
subset = res[[by]] == i & res$x > 0), silent = TRUE)
if (class(fm)[1] == "try-error")
fmList[[i]] <- NULL
else
fmList[[i]] <- fm
}
}
for (i in seq_along(fmList))
attr(fmList[[i]], "unit") <- unit
# return fit or list of fits
if (length(fmList) == 1)
return(fmList[[1]])
else
return(fmList)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.