---

title: "Introduction to outbreaker2"

author: "Thibaut Jombart"

date: "`r Sys.Date()`"

output:

   rmarkdown::html_vignette:

     toc: true

     toc_depth: 2

vignette: >

  %\VignetteIndexEntry{Introduction to outbreaker2}

  %\VignetteEngine{knitr::rmarkdown}

  %\VignetteEncoding{UTF-8}

---

```{r, echo = FALSE}

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:

```{r, eval=FALSE}

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:

```{r, data}

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:

```{r, w}

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")

```

<br>

# Running the analysis with defaults

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:

```{r}

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:

```{r, first_run, cache = TRUE}

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`:

```{r}

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]`:

```{r}

names(res)

```

<br>

# Analysing the results

## Graphics 

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:

```{r, basic_trace}

plot(res)

```

The second argument of `plot` can be used to visualise traces of any

other column in `res`:

```{r, traces}

plot(res, "prior")

plot(res, "mu")

plot(res, "t_inf_15")

```

`burnin` can be used to discard the first iterations prior to mixing:

```{r, basic_trace_burn}

## 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:

```{r, many_plots}

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)

```

## Using `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':

```{r, summary}

summary(res)

```

<br>

# Customising settings and priors

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)

```{r, config2, cache = TRUE}

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`):

```{r, res2}

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:

```{r, net2}

plot(res2, type = "network", burnin = 3000, min_support = 0.01)

```