diff options
| -rw-r--r-- | R/WP.R | 4 | ||||
| -rw-r--r-- | R/WP_internal.R | 8 | ||||
| -rw-r--r-- | R/seqB.R | 4 | ||||
| -rw-r--r-- | man/ID.Rd | 44 | ||||
| -rw-r--r-- | man/IDEA.Rd | 44 | ||||
| -rw-r--r-- | man/WP.Rd | 123 | ||||
| -rw-r--r-- | man/WP_known.Rd | 21 | ||||
| -rw-r--r-- | man/WP_unknown.Rd | 36 | ||||
| -rw-r--r-- | man/computeLL.Rd | 22 | ||||
| -rw-r--r-- | man/seqB.Rd | 117 |
10 files changed, 175 insertions, 248 deletions
@@ -43,8 +43,8 @@ #' \code{scale.max}, which is the largest possible value of the #' scale parameter; and the third is \code{shape.max}, which is #' the largest possible value of the shape parameter. Defaults to -#' \code{B=100, scale.max=10, shape.max=10}. For both shape and -#' scale, the smallest possible value is 1/\code{B}. +#' \code{B = 100, scale.max = 10, shape.max = 10}. For both shape +#' and scale, the smallest possible value is 1/\code{B}. #' @param tol Cutoff value for cumulative distribution function of the #' pre-discretization gamma serial distribution. Defaults to 0.999 #' (i.e. in the discretization, the maximum is chosen such that the diff --git a/R/WP_internal.R b/R/WP_internal.R index 420d0c0..dd10d29 100644 --- a/R/WP_internal.R +++ b/R/WP_internal.R @@ -9,7 +9,7 @@ #' #' @return The function returns the maximum likelihood estimator of R0. #' -#' @keywords internal +#' @noRd WP_known <- function(NT, p) { k <- length(p) TT <- length(NT) - 1 @@ -54,7 +54,7 @@ WP_known <- function(NT, p) { #' #' @importFrom stats pgamma qgamma #' -#' @keywords internal +#' @noRd WP_unknown <- function(NT, B = 100, shape.max = 10, scale.max = 10, tol = 0.999) { shape <- seq(0, shape.max, length.out = B + 1) @@ -91,14 +91,14 @@ WP_unknown <- function(NT, B = 100, shape.max = 10, scale.max = 10, #' This is a background/internal function called by \code{WP}. It computes the #' log-likelihood. #' -#' @param NT Vector of case counts. #' @param p Discretized version of the serial distribution. +#' @param NT Vector of case counts. #' @param R0 Basic reproductive ratio. #' #' @return This function returns the log-likelihood at the input variables and #' parameters. #' -#' @keywords internal +#' @noRd computeLL <- function(p, NT, R0) { k <- length(p) TT <- length(NT) - 1 @@ -11,7 +11,7 @@ #' posterior distribution. The prior distribution is the initial belief of the #' distribution of R0, which is the uninformative uniform distribution with #' values between zero and \code{kappa}. Users can change the value of -#' /code{kappa} only (i.e., the prior distribution cannot be changed from the +#' \code{kappa} only (i.e., the prior distribution cannot be changed from the #' uniform). As more case counts are observed, the influence of the prior #' distribution should lessen on the final estimate \code{Rhat}. #' @@ -41,7 +41,7 @@ #' \code{Rhat} is the estimate of R0 (the posterior mean), #' \code{posterior} is the posterior distribution of R0 from which #' alternate estimates can be obtained (see examples), and \code{group} -#' is an indicator variable (if \code{group=TRUE}, zero values of NT +#' is an indicator variable (if \code{group == TRUE}, zero values of NT #' were input and grouping was done to obtain \code{Rhat}). The variable #' \code{posterior} is returned as a list made up of \code{supp} (the #' support of the distribution) and \code{pmf} (the probability mass @@ -9,41 +9,39 @@ ID(NT, mu) \arguments{ \item{NT}{Vector of case counts.} -\item{mu}{Mean of the serial distribution. This 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 If case -counts are daily and the serial distribution has a mean of seven days, then \code{mu} should be set to seven.} +\item{mu}{Mean of the serial distribution. This 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. If case counts are daily and the serial distribution +has a mean of seven days, then \code{mu} should be set to seven.} } \value{ \code{ID} returns a single value, the estimate of R0. } \description{ -This function implements a least squares estimation method of R0 due to Fisman et al. (PloS One, 2013). -See details for implementation notes. +This function implements a least squares estimation method of R0 due to +Fisman et al. (PloS One, 2013). See details for implementation notes. } \details{ -The method is based on a straightforward incidence decay model. The estimate of R0 is the value which -minimizes the sum of squares between observed case counts and cases counts 'expected' under the model. +The method is based on a straightforward incidence decay model. The estimate +of R0 is the value which minimizes the sum of squares between observed case +counts and cases counts 'expected' under the model. -This method is based on an approximation of the SIR model, which is most valid at the beginning of an epidemic. -The method assumes that the mean of the serial distribution (sometimes called the serial interval) is known. -The final estimate can be quite sensitive to this value, so sensitivity testing is strongly recommended. -Users should be careful about units of time (e.g., are counts observed daily or weekly?) when implementing. +This method is based on an approximation of the SIR model, which is most +valid at the beginning of an epidemic. The method assumes that the mean of +the serial distribution (sometimes called the serial interval) is known. The +final estimate can be quite sensitive to this value, so sensitivity testing +is strongly recommended. Users should be careful about units of time (e.g., +are counts observed daily or weekly?) when implementing. } \examples{ -## ===================================================== ## -## Illustrate on weekly data ## -## ===================================================== ## - +# Weekly data: NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4) -## obtain Rhat when serial distribution has mean of five days -ID(NT=NT, mu=5/7) -## obtain Rhat when serial distribution has mean of three days -ID(NT=NT, mu=3/7) -## ========================================================= ## -## Compute Rhat using only the first five weeks of data ## -## ========================================================= ## +# Obtain R0 when the serial distribution has a mean of five days. +ID(NT, mu = 5 / 7) -ID(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days +# Obtain R0 when the serial distribution has a mean of three days. +ID(NT, mu = 3 / 7) } diff --git a/man/IDEA.Rd b/man/IDEA.Rd index 2dc8240..f4c7d14 100644 --- a/man/IDEA.Rd +++ b/man/IDEA.Rd @@ -9,42 +9,40 @@ IDEA(NT, mu) \arguments{ \item{NT}{Vector of case counts.} -\item{mu}{Mean of the serial distribution. This 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. If case -counts are daily and the serial distribution has a mean of seven days, then \code{mu} should be set to seven.} +\item{mu}{Mean of the serial distribution. This 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. If case counts are daily and the serial distribution has a +mean of seven days, then \code{mu} should be set to seven.} } \value{ \code{IDEA} returns a single value, the estimate of R0. } \description{ -This function implements a least squares estimation method of R0 due to Fisman et al. (PloS One, 2013). -See details for implementation notes. +This function implements a least squares estimation method of R0 due to +Fisman et al. (PloS One, 2013). See details for implementation notes. } \details{ -This method is closely related to that implemented in \code{ID}. The method is based on an incidence decay model. -The estimate of R0 is the value which minimizes the sum of squares between observed case counts and cases counts +This method is closely related to that implemented in \code{ID}. The method +is based on an incidence decay model. The estimate of R0 is the value which +minimizes the sum of squares between observed case counts and cases counts expected under the model. -This method is based on an approximation of the SIR model, which is most valid at the beginning of an epidemic. -The method assumes that the mean of the serial distribution (sometimes called the serial interval) is known. -The final estimate can be quite sensitive to this value, so sensitivity testing is strongly recommended. -Users should be careful about units of time (e.g., are counts observed daily or weekly?) when implementing. +This method is based on an approximation of the SIR model, which is most +valid at the beginning of an epidemic. The method assumes that the mean of +the serial distribution (sometimes called the serial interval) is known. The +final estimate can be quite sensitive to this value, so sensitivity testing +is strongly recommended. Users should be careful about units of time (e.g., +are counts observed daily or weekly?) when implementing. } \examples{ -## ===================================================== ## -## Illustrate on weekly data ## -## ===================================================== ## - +# Weekly data. NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4) -## obtain Rhat when serial distribution has mean of five days -IDEA(NT=NT, mu=5/7) -## obtain Rhat when serial distribution has mean of three days -IDEA(NT=NT, mu=3/7) -## ========================================================= ## -## Compute Rhat using only the first five weeks of data ## -## ========================================================= ## +# Obtain R0 when the serial distribution has a mean of five days. +IDEA(NT, mu = 5 / 7) -IDEA(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days +# Obtain R0 when the serial distribution has a mean of three days. +IDEA(NT, mu = 3 / 7) } @@ -14,81 +14,90 @@ WP( \arguments{ \item{NT}{Vector of case counts.} -\item{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}.} +\item{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}.} -\item{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 \code{shape.max}, which is the largest possible value of the shape parameter. Defaults to -\code{B=100, scale.max=10, shape.max=10}. For both shape and scale, the smallest possible -value is 1/\code{B}.} +\item{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 \code{shape.max}, which is +the largest possible value of the shape parameter. Defaults to +\code{B = 100, scale.max = 10, shape.max = 10}. For both shape +and scale, the smallest possible value is 1/\code{B}.} -\item{tol}{Cutoff value for cumulative distribution function of the pre-discretization gamma serial -distribution. Defaults to 0.999 (i.e. in the discretization, the maximum is chosen such that the -original gamma distribution has cumulative probability of no more than 0.999 at this maximum).} +\item{tol}{Cutoff value for cumulative distribution function of the +pre-discretization gamma serial distribution. Defaults to 0.999 +(i.e. in the discretization, the maximum is chosen such that the +original gamma distribution has cumulative probability of no more +than 0.999 at this maximum).} } \value{ -\code{WP} returns a list containing the following components: \code{Rhat} is the estimate of R0, - and \code{SD} 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}). 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) and \code{pmf} (the probability mass function). +\code{WP} returns a list containing the following components: + \code{Rhat} is the estimate of R0, and \code{SD} 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}). 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) and \code{pmf} (the + probability mass function). } \description{ -This function implements an R0 estimation due to White and Pagano (Statistics in Medicine, 2008). -The method is based on maximum likelihood estimation in a Poisson transmission model. -See details for important implementation notes. +This function implements an R0 estimation due to White and Pagano (Statistics +in Medicine, 2008). The method is based on maximum likelihood estimation in a +Poisson transmission model. See details for important implementation notes. } \details{ -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 \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 \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). +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 \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 \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 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. +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 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, +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 +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. } \examples{ -## ===================================================== ## -## Illustrate on weekly data ## -## ===================================================== ## - +# Weekly data. 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) + +# Obtain R0 when the serial distribution has a mean of five days. +res1 <- WP(NT, mu = 5 / 7) res1$Rhat -## obtain Rhat when serial distribution has mean of three days -res2 <- WP(NT=NT, mu=3/7) + +# Obtain R0 when the serial distribution has a mean of three days. +res2 <- WP(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) + +# Obtain R0 when the serial distribution is unknown. +# NOTE: This implementation will take longer to run. +res3 <- WP(NT) res3$Rhat -## find mean of estimated serial distribution + +# Find the mean of the estimated serial distribution. serial <- res3$SD sum(serial$supp * serial$pmf) -## ========================================================= ## -## Compute Rhat using only the first five weeks of data ## -## ========================================================= ## - -res4 <- WP(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days -res4$Rhat - } diff --git a/man/WP_known.Rd b/man/WP_known.Rd deleted file mode 100644 index aa6ab84..0000000 --- a/man/WP_known.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/WP_known.R -\name{WP_known} -\alias{WP_known} -\title{WP method background function WP_known} -\usage{ -WP_known(NT, p) -} -\arguments{ -\item{NT}{Vector of case counts.} - -\item{p}{Discretized version of the serial distribution.} -} -\value{ -The function returns the maximum likelihood estimator of R0. -} -\description{ -This is a background/internal function called by \code{WP}. It computes the maximum -likelihood estimator of R0 assuming that the serial distribution is known and finite discrete. -} -\keyword{internal} diff --git a/man/WP_unknown.Rd b/man/WP_unknown.Rd deleted file mode 100644 index 02a3280..0000000 --- a/man/WP_unknown.Rd +++ /dev/null @@ -1,36 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/WP_unknown.R -\name{WP_unknown} -\alias{WP_unknown} -\title{WP method background function WP_unknown} -\usage{ -WP_unknown(NT, B = 100, shape.max = 10, scale.max = 10, tol = 0.999) -} -\arguments{ -\item{NT}{Vector of case counts.} - -\item{B}{Length of grid for shape and scale (grid search parameter).} - -\item{shape.max}{Maximum shape value (grid \code{search} parameter).} - -\item{scale.max}{Maximum scale value (grid \code{search} parameter).} - -\item{tol}{cutoff value for cumulative distribution function of the serial distribution (defaults to 0.999).} -} -\value{ -The function returns \code{Rhat}, the maximum likelihood estimator of R0, as well as the maximum - likelihood estimator of the discretized serial distribution given by \code{p} (the probability mass - function) and \code{range.max} (the distribution has support on the integers one to \code{range.max}). - The function also returns \code{resLL} (all values of the log-likelihood) at \code{shape} (grid for - shape parameter) and at \code{scale} (grid for scale parameter), as well as \code{resR0} (the full - vector of maximum likelihood estimators), \code{JJ} (the locations for the likelihood for these), and - \code{J0} (the location for the maximum likelihood estimator \code{Rhat}). If \code{JJ} and \code{J0} - are not the same, this means that the maximum likelihood estimator is not unique. -} -\description{ -This is a background/internal function called by \code{WP}. It computes the maximum likelihood estimator -of R0 assuming that the serial distribution is unknown but comes from a discretized gamma distribution. -The function then implements a simple grid search algorithm to obtain the maximum likelihood estimator -of R0 as well as the gamma parameters. -} -\keyword{internal} diff --git a/man/computeLL.Rd b/man/computeLL.Rd deleted file mode 100644 index 7e9926a..0000000 --- a/man/computeLL.Rd +++ /dev/null @@ -1,22 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/computeLL.R -\name{computeLL} -\alias{computeLL} -\title{WP method background function computeLL} -\usage{ -computeLL(p, NT, R0) -} -\arguments{ -\item{p}{Discretized version of the serial distribution.} - -\item{NT}{Vector of case counts.} - -\item{R0}{Basic reproductive ratio.} -} -\value{ -This function returns the log-likelihood at the input variables and parameters. -} -\description{ -This is a background/internal function called by \code{WP}. It computes the log-likelihood. -} -\keyword{internal} diff --git a/man/seqB.Rd b/man/seqB.Rd index 184f476..0864294 100644 --- a/man/seqB.Rd +++ b/man/seqB.Rd @@ -9,82 +9,83 @@ seqB(NT, mu, kappa = 20) \arguments{ \item{NT}{Vector of case counts.} -\item{mu}{Mean of the serial distribution. This 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. If case -counts are daily and the serial distribution has a mean of seven days, then \code{mu} should be set to seven.} +\item{mu}{Mean of the serial distribution. This 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. If case counts are daily and the serial distribution has a +mean of seven days, then \code{mu} should be set to seven.} -\item{kappa}{Largest possible value of uniform prior (defaults to 20). This describes the prior belief on ranges of R0, -and should be set to a higher value if R0 is believed to be larger.} +\item{kappa}{Largest possible value of uniform prior (defaults to 20). This +describes the prior belief on ranges of R0, and should be set to +a higher value if R0 is believed to be larger.} } \value{ -\code{secB} returns a list containing the following components: \code{Rhat} is the estimate of R0 (the posterior mean), - \code{posterior} is the posterior distribution of R0 from which alternate estimates can be obtained (see examples), - and \code{group} is an indicator variable (if \code{group=TRUE}, zero values of NT were input and grouping was done - to obtain \code{Rhat}). The variable \code{posterior} is returned as a list made up of \code{supp} (the support of - the distribution) and \code{pmf} (the probability mass function). +\code{seqB} returns a list containing the following components: + \code{Rhat} is the estimate of R0 (the posterior mean), + \code{posterior} is the posterior distribution of R0 from which + alternate estimates can be obtained (see examples), and \code{group} + is an indicator variable (if \code{group == TRUE}, zero values of NT + were input and grouping was done to obtain \code{Rhat}). The variable + \code{posterior} is returned as a list made up of \code{supp} (the + support of the distribution) and \code{pmf} (the probability mass + function). } \description{ -This function implements a sequential Bayesian estimation method of R0 due to Bettencourt and Riberio (PloS One, 2008). -See details for important implementation notes. +This function implements a sequential Bayesian estimation method of R0 due to +Bettencourt and Riberio (PloS One, 2008). See details for important +implementation notes. } \details{ -The method sets a uniform prior distribution on R0 with possible values between zero and \code{kappa}, discretized to a fine grid. -The distribution of R0 is then updated sequentially, with one update for each new case count observation. -The final estimate of R0 is \code{Rhat}, the mean of the (last) posterior distribution. -The prior distribution is the initial belief of the distribution of R0; which in this implementation is the uninformative uniform -distribution with values between zero and \code{kappa}. Users can change the value of /code{kappa} only (i.e., the prior distribution -cannot be changed from the uniform). As more case counts are observed, the influence of the prior distribution should lessen on -the final estimate \code{Rhat}. +The method sets a uniform prior distribution on R0 with possible values +between zero and \code{kappa}, discretized to a fine grid. The distribution +of R0 is then updated sequentially, with one update for each new case count +observation. The final estimate of R0 is \code{Rhat}, the mean of the (last) +posterior distribution. The prior distribution is the initial belief of the +distribution of R0, which is the uninformative uniform distribution with +values between zero and \code{kappa}. Users can change the value of +\code{kappa} only (i.e., the prior distribution cannot be changed from the +uniform). As more case counts are observed, the influence of the prior +distribution should lessen on the final estimate \code{Rhat}. -This method is based on an approximation of the SIR model, which is most valid at the beginning of an epidemic. The method assumes -that the mean of the serial distribution (sometimes called the serial interval) is known. The final estimate can be quite sensitive -to this value, so sensitivity testing is strongly recommended. Users should be careful about units of time (e.g., are counts observed -daily or weekly?) when implementing. +This method is based on an approximation of the SIR model, which is most +valid at the beginning of an epidemic. The method assumes that the mean of +the serial distribution (sometimes called the serial interval) is known. The +final estimate can be quite sensitive to this value, so sensitivity testing +is strongly recommended. Users should be careful about units of time (e.g., +are counts observed daily or weekly?) when implementing. -Our code has been modified to provide an estimate even if case counts equal to zero are present in some time intervals. This is done -by grouping the counts over such periods of time. Without grouping, and in the presence of zero counts, no estimate can be provided. +Our code has been modified to provide an estimate even if case counts equal +to zero are present in some time intervals. This is done by grouping the +counts over such periods of time. Without grouping, and in the presence of +zero counts, no estimate can be provided. } \examples{ -## ===================================================== ## -## Illustrate on weekly data ## -## ===================================================== ## - +# Weekly data. NT <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4) -## obtain Rhat when serial distribution has mean of five days -res1 <- seqB(NT=NT, mu=5/7) + +## Obtain R0 when the serial distribution has a mean of five days. +res1 <- seqB(NT, mu = 5 / 7) res1$Rhat -## obtain Rhat when serial distribution has mean of three days -res2 <- seqB(NT=NT, mu=3/7) + +## Obtain R0 when the serial distribution has a mean of three days. +res2 <- seqB(NT, mu = 3 / 7) res2$Rhat -## ============================================================= ## -## Compute posterior mode instead of posterior mean and plot ## -## ============================================================= ## +# Compute posterior mode instead of posterior mean and plot. -Rpost <- res1$posterior +Rpost <- res1$posterior loc <- which(Rpost$pmf == max(Rpost$pmf)) -Rpost$supp[loc] # posterior mode -res1$Rhat # compare with posterior mean - -par(mfrow=c(2, 1), mar=c(2, 2, 1, 1)) -plot(Rpost$supp, Rpost$pmf, col="black", type="l", xlab="", ylab="") -abline(h=1/(20/0.01+1), col="red") -abline(v=res1$Rhat, col="blue") -abline(v=Rpost$supp[loc], col="purple") -legend("topright", legend=c("prior", "posterior", "posterior mean (Rhat)", "posterior mode"), - col=c("red", "black", "blue", "purple"), lty=1) -plot(Rpost$supp, Rpost$pmf, col="black", type="l", xlim=c(0.5, 1.5), xlab="", ylab="") -abline(h=1/(20/0.01+1), col="red") -abline(v=res1$Rhat, col="blue") -abline(v=Rpost$supp[loc], col="purple") -legend("topright", legend=c("prior", "posterior", "posterior mean (Rhat)", "posterior mode"), - col=c("red", "black", "blue", "purple"), lty=1) +Rpost$supp[loc] # Posterior mode. +res1$Rhat # Compare with the posterior mean. -## ========================================================= ## -## Compute Rhat using only the first five weeks of data ## -## ========================================================= ## +par(mfrow = c(2, 1), mar = c(2, 2, 1, 1)) -res3 <- seqB(NT=NT[1:5], mu=5/7) # serial distribution has mean of five days -res3$Rhat +plot(Rpost$supp, Rpost$pmf, col = "black", type = "l", xlab = "", ylab = "") +abline(h = 1 / (20 / 0.01 + 1), col = "red") +abline(v = res1$Rhat, col = "blue") +abline(v = Rpost$supp[loc], col = "purple") +legend("topright", + legend = c("Prior", "Posterior", "Posterior mean", "Posterior mode"), + col = c("red", "black", "blue", "purple"), lty = 1) } |