## 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)
if (FALSE) {
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)
if (FALSE) {
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.
Run the code above in your browser using DataLab