knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=7, fig.height=5 )
incidence implements functions and classes to compute, handle, visualise and model incidences from
dates data. This vignette provides an overview of current features. It largely reproduces the
content of REAME.md
.
To install the current stable, CRAN version of the package, type:
install.packages("incidence")
To benefit from the latest features and bug fixes, install the development, github version of the package using:
devtools::install_github("reconhub/incidence")
Note that this requires the package devtools installed.
The main functions of the package include:
incidence
: compute incidence from dates in various formats; any fixed
time interval can be used; the returned object is an instance of the (S3)
class incidence.
plot
: this method (see ?plot.incidence
for details) plots incidence
objects, and can also add predictions of the model(s) contained in an
incidence_fit object (or a list of such objects).
fit
: fit one or two exponential models (i.e. linear regression on
log-incidence) to an incidence object; two models are calibrated only if a
date is provided to split the time series in two (argument split
); this is
typically useful to model the two phases of exponential growth, and decrease
of an outbreak; each model returned is an instance of the (S3) class
incidence_fit, each of which contains various useful information
(e.g. growth rate r, doubling/halving time, predictions and confidence
intervals); results can be plotted using plot
, or added to an existing
uncudence
plot using the piping-friendly function add_incidence_fit
.
fit_optim_split
: finds the optimal date to split the time series in two,
typically around the peak of the epidemic.
[
: lower-level subsetting of incidence objects, permitting to specify
which dates and groups to retain; uses a syntax similar to matrices,
i.e. x[i, j]
, where x
is the incidence object, i
a subset of dates,
and j
a subset of groups.
subset
: subset an incidence object by specifying a time window.
pool
: pool incidence from different groups into one global incidence
time series.
cumulate
: computes cumulative incidence over time from and incidence
object.
as.data.frame
: converts an incidence object into a data.frame
containing dates and incidence values.
bootstrap
: generates a bootstrapped incidence object by re-sampling,
with replacement, the original dates of events.
find_peak
: locates the peak time of the epicurve.
estimate_peak
: uses bootstrap to estimate the peak time (and related
confidence interval) of a partially observed outbreak.
This example uses the simulated Ebola Virus Disease (EVD) outbreak from the package outbreaks. We will compute incidence for various time steps, calibrate two exponential models around the peak of the epidemic, and analyse the results.
First, we load the data:
library(outbreaks) library(ggplot2) library(incidence) dat <- ebola_sim$linelist$date_of_onset class(dat) head(dat)
We compute the daily incidence:
i <- incidence(dat) i plot(i)
The daily incidence is quite noisy, but we can easily compute other incidence using larger time intervals:
# weekly, starting on Monday (ISO week, default) i.7 <- incidence(dat, interval = "1 week") plot(i.7) # semi-weekly, starting on Saturday i.14 <- incidence(dat, interval = "2 saturday weeks") plot(i.14, border = "white") ## monthly i.month <- incidence(dat, interval = "1 month") plot(i.month, border = "white")
incidence
can also compute incidence by specified groups using the groups
argument. For instance, we can compute incidence by gender:
i.7.sex <- incidence(dat, interval = "1 week", groups = ebola_sim$linelist$gender) i.7.sex plot(i.7.sex, stack = TRUE, border = "grey")
We can do the same for hospitals, using the 'clean' version of the dataset, with some customization of the legend:
i.7.hosp <- with(ebola_sim_clean$linelist, incidence(date_of_onset, interval = "week", groups = hospital)) i.7.hosp head(get_counts(i.7.hosp)) plot(i.7.hosp, stack=TRUE) + theme(legend.position= "top") + labs(fill="")
incidence
objectsincidence
objects can be manipulated easily. The [
operator implements
subsetting of dates (first argument) and groups (second argument). For
instance, to keep only the peak of the distribution:
i[100:250] plot(i[100:250])
Or to keep every other week:
i.7[c(TRUE,FALSE)] plot(i.7[c(TRUE,FALSE)])
Some temporal subsetting can be even simpler using subset
, which permits to
retain data within a specified time window:
i.tail <- subset(i, from=as.Date("2015-01-01")) i.tail plot(i.tail, border="white")
Subsetting groups can also matter. For instance, let's try and visualise the incidence based on onset of symptoms by outcome:
i.7.outcome <- incidence(dat, 7, groups=ebola_sim$linelist$outcome) i.7.outcome plot(i.7.outcome, stack = TRUE, border = "grey")
By default, incidence
treats missing data (NA) as a separate group (see
argument na_as_group
). We could disable this to retain only known outcomes,
but alternatively we can simply subset the object to exclude the last (3rd)
group:
i.7.outcome[,1:2] plot(i.7.outcome[,1:2], stack = TRUE, border = "grey")
Groups can also be collapsed into a single time series using pool
:
i.pooled <- pool(i.7.outcome) i.pooled identical(i.7$counts, i.pooled$counts)
Incidence data, excluding zeros, can be modelled using log-linear regression of the form: log(y) = r x t + b
where y is the incidence, r is the growth rate, t is the number of days since a specific point in time (typically the start of the outbreak), and b is the intercept.
Such model can be fitted to any incidence object using fit
. Of course, a
single log-linear model is not sufficient for modelling our epidemic curve, as
there is clearly an growing and a decreasing phase. As a start, we can
calibrate a model on the first 20 weeks of the epidemic:
plot(i.7[1:20]) early.fit <- fit(i.7[1:20]) early.fit
The resulting objects (known as incidence_fit
objects) can be plotted, in
which case the prediction and its confidence interval is displayed:
plot(early.fit)
However, a better way to display these predictions is adding them to the
incidence plot using the argument fit
:
plot(i.7[1:20], fit = early.fit)
In this case, we would ideally like to fit two models, before and after the peak of the epidemic. This is possible using the following approach, if you know what date to use to split the data in two phases:
fit.both <- fit(i.7, split=as.Date("2014-10-15")) fit.both plot(i.7, fit=fit.both)
This is much better, but the splitting date is not completely optimal. To look for the best possible splitting date (i.e. the one maximizing the average fit of both models), we use:
best.fit <- fit_optim_split(i.7) best.fit plot(i.7, fit=best.fit$fit)
These models are very good approximation of these data, showing a doubling time
of r round(get_info(best.fit$fit, "doubling"), 1)
days during the first
phase, and a halving time of r round(get_info(best.fit$fit, "halving"), 1)
days during the second.
To access these parameters, you can use the get_info()
function.
The possible parameters are:
For "r", "doubling", and "halving", you can also add ".conf" to get the confidence intervals. Here's how you can get the doubling and halving times of the above epi curve:
get_info(best.fit$fit, "doubling") # doubling time get_info(best.fit$fit, "doubling.conf") # confidence interval get_info(best.fit$fit, "halving") get_info(best.fit$fit, "halving.conf")
Note that fit
will also take groups into account if incidence has been
computed for several groups:
best.fit2 <- fit_optim_split(i.7.sex) best.fit2 plot(i.7.sex, fit=best.fit2$fit)
Using get_info()
on this fit object will return all groups together:
get_info(best.fit2$fit, "doubling") # doubling time get_info(best.fit2$fit, "doubling.conf") # confidence interval get_info(best.fit2$fit, "halving") get_info(best.fit2$fit, "halving.conf")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.