knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=8, fig.height=5, fig.path="figs-introduction/" )
This tutorial provides a worked example of outbreak reconstruction using outbreaker2. For installation guidelines, a general overview of the package's functionalities as well as other resources, see the 'overview' vignette:
vignette("Overview", package = "outbreaker2")
We will be analysing a small simulated outbreak distributed with the package,
fake_outbreak
. This dataset contains simulated dates of onsets, partial
contact tracing data and pathogen genome sequences for 30 cases:
library(ape) library(outbreaker2) col <- "#6666cc" fake_outbreak
Here, we will use the dates of case isolation $sample
, DNA sequences $dna
, contact tracing data $ctd
and the empirical distribution of the generation time $w
, which can be visualised as:
plot(fake_outbreak$w, type = "h", xlim = c(0, 5), lwd = 30, col = col, lend = 2, xlab = "Days after infection", ylab = "p(new case)", main = "Generation time distribution")
By default, outbreaker2 uses the settings defined by create_config()
; see
the documentation of this function for details. Note that the main function of
outbreaker2 is called outbreaker
(without number). The function's arguments are:
args(outbreaker)
The only mandatory input really is the data. For most cases, customising the
method will be done through config
and the function create_config()
, which
creates default and alters settings such as prior parameters, length and rate of
sampling from the MCMC, and definition of which parameters should be estimated
('moved'). The last arguments of outbreaker
are used to specify custom prior,
likelihood, and movement functions, and are detailed in the 'Customisation'
vignette.
Let us run the analysis with default settings:
dna <- fake_outbreak$dna dates <- fake_outbreak$sample ctd <- fake_outbreak$ctd w <- fake_outbreak$w data <- outbreaker_data(dna = dna, dates = dates, ctd = ctd, w_dens = w) ## we set the seed to ensure results won't change set.seed(1) res <- outbreaker(data = data)
This analysis will take around 40 seconds on a modern computer. Note that outbreaker2 is slower than outbreaker for the same number of iterations, but the two implementations are actually different. In particular, outbreaker2 performs many more moves than the original package for each iteration of the MCMC, resulting in more efficient mixing. In short: outbreaker2 is slower, but it requires far less iterations.
Results are stored in a data.frame
with the special class outbreaker_chains
:
class(res) dim(res) res
Each row of res
contains a sample from the MCMC. For each, informations about
the step (iteration of the MCMC), log-values of posterior, likelihood and
priors, and all parameters and augmented data are returned. Ancestries
(i.e. indices of the most recent ancestral case for a given case), are indicated
by alpha_[index of the case]
, dates of infections by t_inf_[index of the
case]
, and number of generations between cases and their infector / ancestor by
kappa_[index of the case]
:
names(res)
Results can be visualised using plot
, which has several options and can be
used to derive various kinds of graphics (see ?plot.outbreaker_chains
). The
basic plot shows the trace of the log-posterior values, which is useful to
assess mixing:
plot(res)
The second argument of plot
can be used to visualise traces of any
other column in res
:
plot(res, "prior") plot(res, "mu") plot(res, "t_inf_15")
burnin
can be used to discard the first iterations prior to mixing:
## compare this to plot(res) plot(res, burnin = 2000)
type
indicates the type of graphic to plot; roughly:
trace
for traces of the MCMC (default)
hist
, density
to assess distributions of quantitative values
alpha
, network
to visualise ancestries / transmission tree; note that
network
opens up an interactive plot and requires a web browser with
Javascript enabled; the argument min_support
is useful to select only the
most supported ancestries and avoid displaying too many links
kappa
to visualise the distributions generations between cases and their
ancestor / infector
Here are a few examples:
plot(res, "mu", "hist", burnin = 2000) plot(res, "mu", "density", burnin = 2000) plot(res, type = "alpha", burnin = 2000) plot(res, type = "t_inf", burnin = 2000) plot(res, type = "kappa", burnin = 2000) plot(res, type = "network", burnin = 2000, min_support = 0.01)
summary
The summary of results derives various distributional statistics for posterior, likelihood and prior densities, as well as for the quantitative parameters. It also builds a consensus tree, by finding for each case the most frequent infector / ancestor in the posterior samples. The corresponding frequencies are reported as 'support'. The most frequent value of kappa is also reported as 'generations':
summary(res)
As said before, most customisation can be achieved via create_config
.
In the following, we make the following changes to the defaults:
increase the number of iterations to 30,000
set the sampling rate to 20
use a star-like initial tree
disable to movement of kappa
, so that we assume that all cases have
observed
set a lower rate for the exponential prior of mu
(10 instead of 1000)
config2 <- create_config(n_iter = 3e4, sample_every = 20, init_tree ="star", move_kappa = FALSE, prior_mu = 10) set.seed(1) res2 <- outbreaker(data, config2) plot(res2) plot(res2, burnin = 2000)
We can see that the burnin is around 2,500 iterations (i.e. after the initial
step corresponding to a local optimum). We get the consensus tree from the new
results, and compare the inferred tree to the actual ancestries stored in the
dataset (fake_outbreak$ances
):
summary(res2, burnin = 3000) tree2 <- summary(res2, burnin = 3000)$tree comparison <- data.frame(case = 1:30, inferred = paste(tree2$from), true = paste(fake_outbreak$ances), stringsAsFactors = FALSE) comparison$correct <- comparison$inferred == comparison$true comparison mean(comparison$correct)
Let's visualise the posterior trees:
plot(res2, type = "network", burnin = 3000, min_support = 0.01)
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.