diff options
| -rw-r--r-- | R/WP.R | 71 |
1 files changed, 30 insertions, 41 deletions
@@ -9,33 +9,29 @@ source("WP_unknown.R") #' #' This method is based on a Poisson transmission model, and hence may be most most valid at the beginning #' of an epidemic. In their model, the serial distribution is assumed to be discrete with a finite number -#' of posible values. In this implementation, if the serial distribution is assumed known, it is taken to +#' of posible values. In this implementation, if \code{mu} is not {NA}, the serial distribution is taken to #' be a discretized version of a gamma distribution with mean \code{mu}, shape parameter one, and largest -#' possible value based on parameter \code{tol}. When the serial distribution is unknown, the function -#' implements a grid search algorithm to find the maximum likelihood estimator over all possible gamma -#' distributions with unknown mean and variance, restricting these to a prespecified grid (see -#' \code{search} parameter). +#' possible value based on parameter \code{tol}. When \code{mu} is \code{NA}, the function implements a +#' grid search algorithm to find the maximum likelihood estimator over all possible gamma distributions +#' with unknown mean and variance, restricting these to a prespecified grid (see \code{search} parameter). #' -#' When the serial distribution is taken to be \code{known}, sensitivity testing of the parameter \code{mu} -#' is strongly recommended. If the serial distribution is \code{unknown}, the likelihood function can be -#' flat near the maximum, resulting in numerical instability of the optimizer. When the serial distribution -#' is \code{unkown} the implementation takes considerably longer to run. Users should be careful about units -#' of time (e.g. are counts observed daily or weekly?) when implementing. +#' When the serial distribution is known (i.e., \code{mu} is not \code{NA}), sensitivity testing of \code{mu} +#' is strongly recommended. If the serial distribution is unknown (i.e., \code{mu} is \code{NA}), the +#' likelihood function can be flat near the maximum, resulting in numerical instability of the optimizer. +#' When \code{mu} is \code{NA}, the implementation takes considerably longer to run. Users should be careful +#' about units of time (e.g. are counts observed daily or weekly?) when implementing. #' #' The model developed in White and Pagano (2008) is discrete, and hence the serial distribution is finite #' discrete. In our implementation, the input value \code{mu} is that of a continuous distribution. The -#' algorithm when \code{method="known"} disretizes this input, and hence the mean of the serial distribution -#' returned in the list \code{SD} will differ from \code{mu} somewhat. That is to say, if the user notices that -#' the input \code{mu} and out put mean of \code{SD} are different, this is to be expected, and is caused by -#' the discretization. +#' algorithm discretizes this input when \code{mu} is not \code{NA}, and hence the mean of the serial +#' distribution returned in the list \code{SD} will differ from \code{mu} somewhat. That is to say, if the +#' user notices that the input \code{mu} and output mean of \code{SD} are different, this is to be expected, +#' and is caused by the discretization. #' #' @param NT Vector of case counts #' @param mu Mean of the serial distribution (needs to match case counts in time units; for example, if case #' counts are weekly and the serial distribution has a mean of seven days, then \code{mu} should be #' set to one). The default value of \code{mu} is set to \code{NA}. -#' @param method Variable taking one of two possible values: \code{known} or \code{unknown}. If "known", the -#' serial distribution is assumed to be gamma with rate 1/\code{mu} and shape equal to one, if -#' "unknown" then the serial distribution is gamma with unknown parameters. Defaults to "unknown" #' @param search List of default values for the grid search algorithm; the list includes three elements: the #' first is \code{B} which is the length of the grid in one dimension, the second is #' \code{scale.max} which is the largest possible value of the scale parameter, and the third is @@ -47,8 +43,8 @@ source("WP_unknown.R") #' original gamma distribution has cumulative probability of no more than 0.999 at this maximum). #' #' @return WP returns a list containing the following components: \code{Rhat} is the estimate of R0, \code{SD} -#' is either the discretized serial distribution (if \code{method="known"}) or the estimated -#' discretized serial distribution (if \code{method="unknown"}), and \code{inputs} is a list of the +#' is either the discretized serial distribution (if \code{mu} is not \code{NA}) or the estimated +#' discretized serial distribution (if \code{mu} is \code{NA}), and \code{inputs} is a list of the #' original input variables \code{NT, mu, method, search, tol}. The list also returns the variable #' \code{check}, which is equal to the number of non-unique maximum likelihood estimators. The serial #' distribution \code{SD} is returned as a list made up of \code{supp} the support of the distribution @@ -59,16 +55,16 @@ source("WP_unknown.R") #' ## Illustrate on weekly data ## #' ## ===================================================== ## #' -#' NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4) +#' NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4) #' ## obtain Rhat when serial distribution has mean of five days -#' res1 <- WP(NT=NT, mu=5/7, method="known") +#' res1 <- WP(NT=NT, mu=5/7) #' res1$Rhat #' ## obtain Rhat when serial distribution has mean of three days -#' res2 <- WP(NT=NT, mu=3/7, method="known") +#' res2 <- WP(NT=NT, mu=3/7) #' res2$Rhat #' ## obtain Rhat when serial distribution is unknown #' ## NOTE: this implementation will take longer to run -#' res3 <- WP(NT=NT) +#' res3 <- WP(NT=NT) #' res3$Rhat #' ## find mean of estimated serial distribution #' serial <- res3$SD @@ -78,33 +74,26 @@ source("WP_unknown.R") #' ## Compute Rhat using only the first five weeks of data ## #' ## ========================================================= ## #' -#' res4 <- WP(NT=NT[1:5], mu=5/7, method="known") # serial distribution has mean of five days +#' res4 <- WP(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days #' res4$Rhat #' #' @export -WP <- function(NT, mu="NA", method="unknown", search=list(B=100, shape.max=10, scale.max=10), tol=0.999) { - if (method == "unknown") { +WP <- function(NT, mu=NA, search=list(B=100, shape.max=10, scale.max=10), tol=0.999) { + if (is.na(mu)) { print("You have assumed that the serial distribution is unknown.") res <- WP_unknown(NT=NT, B=search$B, shape.max=search$shape.max, scale.max=search$scale.max, tol=tol) Rhat <- res$Rhat p <- res$p range.max <- res$range.max JJ <- res$JJ - } - - if (method == "known") { - if (mu=="NA") { - res <- "NA" - print("For method=known, the mean of the serial distribution must be specified.") - } else { - print("You have assumed that the serial distribution is known.") - range.max <- ceiling(qexp(tol, rate=1/mu)) - p <- diff(pexp(0:range.max, 1/mu)) - p <- p / sum(p) - res <- WP_known(NT=NT, p=p) - Rhat <- res - JJ <- NA - } + } else { + print("You have assumed that the serial distribution is known.") + range.max <- ceiling(qexp(tol, rate=1/mu)) + p <- diff(pexp(0:range.max, 1/mu)) + p <- p / sum(p) + res <- WP_known(NT=NT, p=p) + Rhat <- res + JJ <- NA } return(list(Rhat=Rhat, check=length(JJ), SD=list(supp=1:range.max, pmf=p))) |