rsu.pfree.rs.Rd
\name{rsu.pfree.rs}
\alias{rsu.pfree.rs}
\title{
Calculate the probability of freedom for given population sensitivity and probability of introduction
}
\description{
Calculates the posterior probability (confidence) of disease freedom (negative predictive value) for one or more population sensitivity (se.p) estimates, over one or more time periods.
}
\usage{
rsu.pfree.rs(se.p, p.intro = 0, prior = 0.5, by.time = TRUE)
}
\arguments{
\item{se.p}{scalar, vector or matrix representing the population sensitivity estimates. \code{se.p} will be scalar if you're calculating the posterior probability of disease freedom for a single time period. If \code{se.p} is a vector set \code{by.time = TRUE} if the \code{se.p} estimates are for separate time periods. Set \code{by.time = FALSE} if the \code{se.p} estimates are variations (iterations) within a single time period. If \code{se.p} is a matrix, columns represent consecutive time periods and rows represent multiple \code{se.p} estimates per time period.}
\item{p.intro}{scalar, vector or matrix representing the probability of disease introduction per time period. If \code{p.intro} is scalar this value is applied across all \code{se.p} values and time periods. If \code{p.intro} is a vector set \code{by.time = TRUE} if the \code{p.intro} estimates are for separate time periods. Set \code{by.time = FALSE} if the \code{p.intro} estimates are variations (iterations) within a single time period. If \code{p.intro} is a matrix it should have the same dimensions as \code{se.p} with columns representing time periods and rows representing multiple \code{p.intro} estimates per time period.}
\item{prior}{scalar or vector of the same length as the number of rows of \code{se.p} representing the prior probability of disease freedom before surveillance.}
\item{by.time}{logical, representing the type of analysis. See details, below.}
}
\details{
The \code{by.time} argument is used for two specific circumstances.
Use \code{by.time = TRUE} if the \code{se.p} estimates are a vector of values for consecutive time periods. Use \code{by.time = FALSE} if the \code{se.p} estimates are a vector of multiple values (iterations) for a single time period.
Use \code{by.times = TRUE} if \code{se.p} is a symmetrical matrix and \code{p.intro} is a vector of values representing the probability of disease introduction over consecutive time periods. Use \code{by.time = FALSE} if \code{se.p} is a symmetrical matrix (with columns for time periods and rows representing estimates of \code{se.p} within each time period) and \code{p.intro} is a vector of values corresponding to multiple values for a single time period that are the same across all periods.
}
\value{
A list comprised of six elements:
\item{PFree}{The posterior probability of disease freedom.}
\item{SeP}{The population sensitivity.}
\item{PIntro}{The probability of disease introduction (as entered by the user).}
\item{Discounted prior}{The discounted prior confidence of disease freedom.}
\item{Equilibrium PFree}{The equilibrium probability of disease freedom.}
\item{Equilibrium prior}{The equilibrium discounted prior probability of disease freedom.}
}
\references{
Martin P, Cameron A, Greiner M (2007). Demonstrating freedom from disease using multiple complex data sources 1: A new methodology based on scenario trees. Preventive Veterinary Medicine 79: 71 - 97.
Martin P, Cameron A, Barfod K, Sergeant E, Greiner M (2007). Demonstrating freedom from disease using multiple complex data sources 2: Case study - classical swine fever in Denmark. Preventive Veterinary Medicine 79: 98 - 115.
}
\examples{
## EXAMPLE 1:
## You have estimated herd-sensitivity for 20 herds for a disease of concern,
## all returned negative results. What is the confidence of disease freedom
## for these herds, assuming that based on other data, 20\% of herds in the
## region are estimated to be disease positive?
## Generate 20 herd sensitivity estimates, using random values between 70\%
## and 95\%:
herd.sens <- runif(n = 20, min = 0.70, max = 0.95)
## The background herd prevalence is 0.20, so the prior confidence of freedom
## is 1 - 0.2 = 0.8. For this example we assume the prior is applicable at
## the time of sampling so p.intro = 0 (the default) and we are carrying out
## an analysis using multiple estimates of population sensitivities for a
## single time period so we set by.time = FALSE.
rval.df <- rsu.pfree.rs(se.p = herd.sens, p.intro = 0, prior = 0.80,
by.time = FALSE)
rval.df <- data.frame(SeP = rval.df$SeP, PFree = rval.df$PFree)
range(rval.df$SeP)
## The herd-level probability of disease freedom ranges from about 0.93 to
## 0.99 depending on individual herd level sensitivity values.
## EXAMPLE 2:
## You have analysed 12 months of surveillance data for disease X, to provide
## 12 monthly estimates of population sensitivity. In addition, based on
## previous data, the monthly probability of the introduction of disease is
## estimated to be in the range of 0.005 (0.5\%) to 0.02 (2\%). The prior
## confidence of disease freedom is assumed to be 0.5 (i.e., uninformed).
## What is your level of confidence of disease freedom at the end of the 12
## month surveillance period?
## Generate 12, monthly estimates of se.p and p.intro:
pop.sens <- runif(n = 12, min = 0.40, max = 0.70)
pintro <- runif(n = 12, min = 0.005, max = 0.020)
## For this example we're analysing a single population over multiple time
## periods, so we set by.time = TRUE:
rval.df <- rsu.pfree.rs(se.p = pop.sens, p.intro = pintro, prior = 0.50,
by.time = TRUE)
rval.df <- data.frame(mnum = 1:12, mchar = seq(as.Date("2020/1/1"),
by = "month", length.out = 12), SeP = t(rval.df$SeP),
PFree = t(rval.df$PFree))
## Plot the probability of disease freedom as a function of time:
plot(x = rval.df$mnum, y = rval.df$PFree, xlim = c(1,12), ylim = c(0,1),
xlab = "Month", ylab = "Probability of disease freedom",
pch = 16, type = "b", xaxt = "n")
axis(side = 1, at = rval.df$mnum,
labels = format(rval.df$mchar, format = "\%b"))
abline(h = 0.95, lty = 2)
\dontrun{
library(ggplot2); library(scales)
ggplot(data = rval.df, aes(x = mchar, y =PFree)) +
geom_line(col = "black") +
scale_x_date(breaks = date_breaks("1 month"), labels = date_format("\%b"),
name = "Month") +
scale_y_continuous(limits = c(0,1), name = "Probability of disease freedom") +
geom_hline(yintercept = 0.95, linetype = "dashed") +
theme_bw()
}
## The estimated probability of disease freedom (Pfree) increases over time
## from about 0.70 (or less) to >0.99, depending on the actual se.p values
## generated by simulation.
## EXAMPLE 3:
## Extending the above example, instead of a simple deterministic estimate,
## you decide to use simulation to account for uncertainty in the monthly
## se.p and p.intro estimates.
## For simplicity, we generate 1200 random estimates of se.p and coerce them
## into a matrix with 12 columns and 100 rows:
pop.sens <- matrix(runif(n = 1200, min = 0.40, max = 0.70), nrow = 100)
## For p.intro we generate a vector of 100 random values, which will then be
## used across all time periods:
pintro <- runif(n = 100, min = 0.005, max = 0.020)
## For this example, because se.p is a matrix and p.intro is a vector matching
## one of the dimensions of se.p, by.time is ignored:
rval.df <- rsu.pfree.rs(se.p = pop.sens, p.intro = pintro, prior = 0.5,
by.time = TRUE)
## Calculate 95\% confidence intervals for the probability of disease freedom:
rval.df <- apply(rval.df$PFree, FUN = quantile, MARGIN = 2,
probs = c(0.025,0.5,0.975))
rval.df <- data.frame(mnum = 1:12, mchar = seq(as.Date("2020/1/1"),
by = "month", length.out = 12), t(rval.df))
## Plot the probability of disease freedom as a function of time. Dashed lines
## show the lower and upper bound of the confidence interval around the
## probability of disease freedom estimates:
plot(x = rval.df$mnum, y = rval.df$X50., xlim = c(1,12), ylim = c(0,1),
xlab = "Month", ylab = "Probability of disease freedom",
type = "l", lwd = 2, xaxt = "n")
axis(side = 1, at = rval.df$mnum, labels = format(rval.df$mchar, format = "\%b"))
lines(x = rval.df$mnum, y = rval.df$X2.5., type = "l", lty = 2)
lines(x = rval.df$mnum, y = rval.df$X97.5., type = "l", lty = 2)
\dontrun{
library(ggplot2); library(scales)
ggplot(data = rval.df, aes(x = mchar, y = X50.)) +
geom_line(col = "black") +
geom_ribbon(aes(ymin = X2.5., ymax = X97.5.), alpha = 0.25) +
scale_x_date(breaks = date_breaks("1 month"), labels = date_format("\%b"),
name = "Month") +
scale_y_continuous(limits = c(0,1), name = "Probability of disease freedom") +
theme_bw()
}
## The median probability of disease freedom increases over time from about
## 0.7 (or less) to >0.99, depending on the actual se.p values generated by
## simulation.
}
\keyword{methods}