#' Get Titer from Fitted Model
#'
#' Extract the viral titer \emph{or} the volume required for one
#' infectious unit with confidence intervals from a fit object or
#' a list of fit objects generated by \code{\link{getFit}}.
#'
#' @param fm fitted model \emph{or} list of fitted models generated by
#' \code{\link{getFit}}.
#' @param level this specifies the confidence level (default of 0.95).
#' Use \code{NULL} to return only the titer (or EC63 value).
#' @param digits \code{Integer} number of significant digits to report
#' (default = 3)
#'
#' @return
#'
#' For \code{getTiter}, a named numeric vector with the estimated titer in
#' IU per ml.
#'
#' For \code{getEC63}, a named numeric vector with the volume corresponding to
#' one infectious unit.
#'
#' If \code{level} was provided, the return value is a matrix with additional
#' columns providing the desired +/- confidence intervals.
#'
#' @export
#'
getTiter <- function(fm, level = 0.95, digits = 3)
{
# option checks
digits <- as.integer(digits[1])
if (!is.null(level)) {
level <- as.numeric(level)[1]
if (level <= 0 | level >= 1)
stop("'level' must be NULL or a single value between 0 and 1")
}
# working function
warnUnits <- FALSE
.fun <- function(fm, digits, level) {
units <- attr(fm, "unit")
if (!units %in% c("ml", "ul", "nl", "pl")) warnUnits <<- TRUE
mult <- switch(units, ml = 1, ul = 1e3, nl = 1e6, pl = 1e9, 1)
return(signif(mult/getEC63(fm, level)[c(1,3,2)], digits))
}
# apply according to argument 'fm'
if (!is.null(level)) {
if (is(fm, "list")) ret <- t(sapply(fm, .fun, digits, level))
else ret <- .fun(fm, digits, level)
}
else {
if (is(fm, "list")) ret <- t(sapply(fm, .fun, digits, level))[,1]
else ret <- .fun(fm, digits, level)[1]
}
# indicate units warning and return named vector or array
if (warnUnits)
warning("no volume attribute found in '", deparse(substitute(fm)), "'")
return(ret)
}
#' @name getEC63
#' @rdname getTiter
#'
#' @export
#'
getEC63 <- function(fm, level = 0.95)
{
if (!is.null(level)) {
level <- as.numeric(level)[1]
if (level <= 0 | level >= 1)
stop("'level' must be NULL or a single value between 0 and 1")
}
# working function
.coef <- function(fm, level) {
cf <- numeric(3)
if (is.null(fm))
cf <- c(NA, NA, NA)
else {
cf[1] <- exp(-coef(fm))
if (!is.null(level)) {
cfVal <- try(confint(fm, level = level))
if(class(cfVal)[1] != "try-error")
cf[c(3,2)] <- exp(-cfVal)
names(cf) <- c("est", "lo.CI", "hi.CI")
}
}
if (is.null(level)) return(cf[1]) else return(cf)
}
if ("list" %in% class(fm))
ret <- t(sapply(fm, .coef, level))
else
ret <- .coef(fm, level)
return(ret)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.