R/getFit.R

Defines functions getFit

Documented in getFit

#' 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)
}
ornelles/virustiter documentation built on March 15, 2024, 9:28 a.m.