---
title: "Using custom priors, likelihood, or movements in outbreaker2"
author: "Thibaut Jombart"
date: "`r Sys.Date()`"
output:
   rmarkdown::html_vignette:
     toc: true
     toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Using custom priors, likelihood, or movements in 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-customisation/"
)
```


In this vignette, we show how custom functions for priors, likelihood, or
movement of parameters and augmented data can be used in *outbreaker2*. In all
these functions, the process will be similar:

1. write your own function with the right arguments
2. pass this function as an argument to a `custom...` function
3. pass the result to *outbreaker2*

Note that 2-3 can be a single step if passing the function to the arguments of
*outbreaker2* directly. Also note that **all priors and likelihoods are expected
on a log scale**. Finally, also note that while the various `custom...`
functions will try to some extent to check that the provided functions are
valid, such tests are very difficult to implement. In short: you are using these
custom features at your own risks - make sure these functions work before
passing them to *outbreaker2*.




<br>

# Customising priors

Priors of *outbreaker2* must be a function of an `outbreaker_param` list (see
`?outbreaker_param`). Here, we decide to use a step function rather than the
default Beta function as a prior for *pi*, the reporting probability, and a flat
prior between 0 and 1 for the mutation rate (which is technically a probability
in the basic genetic model used in *outbreaker2*).


We start by defining two functions: an auxiliary function `f` which returns
values on the natural scale, and which we can use for plotting the prior
distribution, and then a function `f_pi` which will be used for the
customisation.

```{r, f_pi}
f <- function(pi) {
    ifelse(pi < 0.8, 0, 5)
}

f_pi <- function(param) { 
    log(f(param$pi))
}

plot(f, type = "s", col = "blue", 
     xlab = expression(pi), ylab = expression(p(pi)), 
     main = expression(paste("New prior for ", pi)))

```


While `f` is a useful function to visualise the prior, `f_pi` is the function
which will be passed to `outbreaker`. To do so, we pass it to `custom_priors`:

```{r, custom_prior}
library(outbreaker2)

f_mu <- function(param) {
  if (param$mu < 0 || param$mu > 1) {
    return(-Inf)
  } else {
    return(0.0)
  }
  
}

priors <- custom_priors(pi = f_pi, mu = f_mu)
priors

```

Note that `custom_priors` does more than just adding the custom function to a
list. For instance, the following customisations are all wrong, and rightfully
rejected:

```{r, wrong_prior, error = TRUE}

## wrong: not a function
## should be pi = function(x){0.0}
custom_priors(pi = 0.0)

## wrong: two arguments
custom_priors(pi = function(x, y){0.0})

```

We can now use the new priors to run `outbreaker` on the `fake_outbreak` data
(see [*introduction vignette*](introduction.html)):

```{r, run_custom_priors}

dna <- fake_outbreak$dna
dates <- fake_outbreak$sample
w <- fake_outbreak$w
data <- outbreaker_data(dna = dna, dates = dates, w_dens = w)

## we set the seed to ensure results won't change
set.seed(1)


res <- outbreaker(data = data, priors = priors)

```

We can check the results first by looking at the traces, and then by plotting
the posterior distributions of `pi` and `mu`, respectively:

```{r, traces_custom_priors}

plot(res)
plot(res, "pi", burnin = 500)
plot(res, "mu", burnin = 500)
plot(res, "pi", type = "density", burnin = 500)
plot(res, "mu", type = "hist", burnin = 500)

```

Note that we are using density and histograms here for illustrative purposes,
but there is no reason to prefer one or the other for a specific
parameter.


Interestingly, the trace of `pi` suggests that the MCMC oscillates between two
different states, on either bound of the interval on which the prior is positive
(it is `-Inf` outside (0.8; 1)). This may be a consequence of the step function,
which causes sharp 'cliffs' in the posterior landscape. What shall one do to
derive good samples from the posterior distribution in this kind of situation?
There are several options, which in fact apply to typical cases of multi-modal
posterior distributions:

- Avoid 'cliffs', i.e. sharp drops in the posterior landscape, typically created
  by using step-functions in likelihoods and in priors.

- Use larger samples, i.e. run more MCMC iterations.

- Use a different sampler, better than Metropolis-Hasting at deriving samples
  from multi-modal distributions.


Because we know what the real transmission tree is for this dataset, we can
assess how the new priors impacted the inference of the transmission tree.

```{r, tree_custom_priors}

summary(res, burnin = 500)
tree <- summary(res, burnin = 500)$tree

comparison <- data.frame(case = 1:30,
                       	 inferred = paste(tree$from),
			 true = paste(fake_outbreak$ances),
			 stringsAsFactors = FALSE)
			 
comparison$correct <- comparison$inferred == comparison$true
comparison
mean(comparison$correct)

```





<br>

# Customizing likelihood

Likelihood functions customisation works identically to prior functions. The
only difference is that custom functions will take two arguments (`data` and
`param`) instead of one in the prior functions. The function used to specify
custom likelihood is `custom_likelihoods`. Each custom function will correspond
to a specific likelihood component:

```{r, likelihood_components}

custom_likelihoods()

```

see `?custom_likelihoods` for details of these components, and see the section
'Extending the model' for new, other components. As for `custom_priors`, a few
checks are performed by `custom_likelihoods`:

```{r, wrong_likelihood, error = TRUE}

## wrong: not a function
custom_likelihoods(genetic = "fubar")

## wrong: only one argument
custom_likelihoods(genetic = function(x){ 0.0 })

```

A trivial customisation is to disable some or all of the likelihood components
of the model by returning a finite constant. Here, we apply this to two cases:
first, we will disable all likelihood components as a sanity check, making sure
that the transmission tree landscape is explored freely by the MCMC. Second, we
will recreate the [Wallinga & Teunis (1994)](http://dx.doi.org/10.1093/aje/kwh255)
model, by disabling specific components.



## A null model

```{r, null_model}

f_null <- function(data, param) {
   return(0.0)
}

null_model <- custom_likelihoods(genetic = f_null,
                                 timing_sampling = f_null,
                                 timing_infections = f_null,
                                 reporting = f_null,
                                 contact = f_null)

null_model

```


We also specify settings via the `config` argument to avoid detecting imported
cases, reduce the number of iterations and sampling each of them:

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

null_config <- list(find_import = FALSE,
n_iter = 500,
sample_every = 1)

set.seed(1)

res_null <- outbreaker(data = data,
config = null_config,
likelihoods = null_model)

```


```{r, res_null_model}

plot(res_null)
plot(res_null, "pi")
plot(res_null, "mu")

```

By typical MCMC standards, these traces look appaling, as they haven't reach
stationarity (i.e. same mean and variance over time), and are grossly
autocorrelated in parts. Fair enough, as these are only the first 500 iterations
of the MCMC, so that autocorrelation is expected. In fact, what we observe here
literally is the random walk across the posterior landscape, which in this case
is only impacted by the priors.


We can check that transmission trees are indeed freely explored:

```{r, null_trees}

plot(res_null, type = "alpha")

```

Do **not** try to render the corresponding network using `plot(..., type =
"network")` as the force-direction algorithm will go insane. However, this
network can be visualised using *igraph*, extracting the edges and nodes from
the plot (without displaying it):

```{r, null_net}

## extract nodes and edges from the visNetwork object
temp <- plot(res_null, type = "network", min_support = 0)
class(temp)
head(temp$x$edges)
head(temp$x$nodes)

## make an igraph object
library(igraph)

net_null <- graph.data.frame(temp$x$edges,
                             vertices = temp$x$nodes[1:4])

plot(net_null, layout = layout.circle,
     main = "Null model, posterior trees")

```

We can derive similar diagnostics for the number of generations betweens cases
(`kappa`), only constrained by default settings to be between 1 and 5, and for
the infection dates (*t_inf*):

```{r, res_null_diag}

plot(res_null, type = "kappa")
plot(res_null, type = "t_inf")

```

Finally, we can verify that the distributions of `mu` and `pi` match their
priors, respectively an exponential distribution with rate 1000 and a beta with
parameters 10 and 1. Here, we get a qualitative assessment by comparing the
observed distribution (histograms) to the densities of similar sized random
samples from the priors:

```{r, res_null_priors}

par(xpd=TRUE)
hist(res_null$mu, prob = TRUE, col = "grey",
     border = "white",
     main = "Distribution of mu")

invisible(replicate(30,
     points(density(rexp(500, 1000)), type = "l", col = "blue")))


hist(res_null$pi, prob = TRUE, col = "grey",
     border = "white", main = "Distribution of pi")

invisible(replicate(30,
     points(density(rbeta(500, 10, 1)), type = "l", col = "blue")))

```






<br>

## A model using symptom onset only

We can use data and likelihood customisation to change the default *outbreaker2*
model into a [Wallinga & Teunis (1994)](http://dx.doi.org/10.1093/aje/kwh255)
model. To do so, we need to:

- Remove the DNA sequences from the data; alternatively we could also specify a
  'null' function (i.e. returning a finite constant, as above) for the genetic
  likelihood.

- Disable all likelihood components other than `timing_infections` using
  `custom_likelihoods`. This means that the dates provided will be treated as
  dates of symptom onset, and the timing distribution `w` will be taken as the
  serial interval.

- Disable the detection of imported cases, and forcing all `kappa` values to be
  1.



While these are fairly major changes, they are straightforward to implement. We
first create the dataset and custom likelihood functions:

```{r, wt_custom}

onset_data <- outbreaker_data(dates = fake_outbreak$onset,
	      	              w_dens = fake_outbreak$w)

wt_model <- custom_likelihoods(timing_sampling = f_null,
                               reporting = f_null)

```

To fix parameters or augmented data (here, fix all `kappa` values to 1), we set
the initial states to the desired values and disable the corresponding moves:

```{r, wt_config}

wt_config <- create_config(init_kappa = 1, move_kappa = FALSE,
                           init_pi = 1, move_pi = FALSE,
                           move_mu = FALSE)

```

We can now run the analyses for this new model:
```{r, run_wt, cache = TRUE}

set.seed(1)
res_wt <- outbreaker(data = onset_data,
                     config = wt_config,
                     likelihoods = wt_model)
		       
```


```{r, res_wt}

plot(res_wt)
plot(res_wt, burnin = 500)
plot(res_wt, burnin = 500, type = "alpha")
summary(res_wt)

```

As before for the 'null' model, the transmission tree is very poorly resolved in
this case. We use the same approach to visualise it: extract nodes and edges
from the `visNetork` object, use this information to create an `igraph` object,
and visualise the result using a circular layout:

```{r, wt_net}

## extract nodes and edges from the visNetwork object
temp <- plot(res_wt, type = "network", min_support = 0.05)
class(temp)
head(temp$x$edges)
head(temp$x$nodes)

## make an igraph object

net_wt <- graph.data.frame(temp$x$edges,
                             vertices = temp$x$nodes[1:4])
			     
plot(net_wt, layout = layout.circle,
     main = "WT model, posterior trees")

```






<br>

# Customising movements

Customising movements works in similar ways to priors and likelihoods. In
practice, this type of customisation is more complex as it most likely will
require evaluation of likelihoods and priors, and therefore require the user to
know which functions to all, and how. These are documented in the [API
vignette](Rcpp_API.html). In the following, we provide two examples:

- a (fake) Gibbs sampler for the movement of the mutation rate `mu`

- a new 'naive' movement of ancestries in which infectors are picked at random
  from all cases

But before getting into these, we need to explicit how movements are happening
in *outbreaker2*.



## Movements in *outbreaker2*

At the core of the `outbreaker` function, movements are implemented as a list of
functions, which are all evaluated in turn during every iteration of the
MCMC. All movement functions must obey two rules:

- The first argument must be an `outbreaker_param` object (typically called
  `param` in the original code); see `?create_param` for details.

- All movement functions must return a valid, `outbreaker_param` object.


However, a new difficulty compared to prior or likelihood customisation is that
different movements may require different components of the model, and a
different set of arguments after `param`. In fact, this can be seen by examining
the arguments of all the default movement functions:

```{r, move_defaults}

lapply(custom_moves(),  args)

```

To handle this difficulty, *outbreaker2* transforms every movement function
before running the MCMC into a new function with a single parameter `param`,
attaching all other required argument to the function's environment. The
function achieving this transformation is called `bind_moves`. This function
'knows' what these components are for known moves listed aboves. For new,
unknown moves, it attaches the following components, provided they are used as
arguments of the new function:

- `data`: the processed data; see `?outbreaker_data`

- `config`: the configuration list; see `create_config`

- `likelihoods`: a list of custom likelihood functions; see
  `?custom_likelihoods`

- `priors`: a list of custom prior functions; see `?custom_priors`


See examples in `?bind_moves` for details of how this process works.





## A (fake) Gibbs sampler for `mu`

A Gibbs sampler supposes that the conditional distribution of a parameter is
known and can directly be sampled from. Here, we use this principle to provide a
toy example of custom movement for `mu`, assuming that this conditional
distribution is always an Exponential distribution with a rate of 1000. This is
easy to implement; to make sure that the function is actually used, we set a
global variable changed when the function is called.


```{r, custom_move_mu}

move_mu <- function(param, config) {

  NEW_MOVE_HAS_BEEN_USED <<- TRUE
  param$mu <- rexp(1, 1000)  
  return(param)
  
}

moves <- custom_moves(mu = move_mu)

quick_config <- list(n_iter = 500, sample_every = 1, find_import = FALSE)

```

Note that the new movement function `move_mu` has two arguments, and that we do
not specify `config`. Internally, what happens is:

```{r, bind_moves}

## bind quick_config to function
move_mu_intern <- bind_to_function(move_mu, config = quick_config)

## new function has just one argument
move_mu_intern

## 'config' is in the function's environment
names(environment(move_mu_intern))

## 'config' is actually 'quick_config'
identical(environment(move_mu_intern)$config, quick_config)

```


We perform a quick run using this new movement:

```{r, run_custom_move_mu}

NEW_MOVE_HAS_BEEN_USED <- FALSE

set.seed(1)
res_move_mu <- outbreaker(data, quick_config, moves = moves)
NEW_MOVE_HAS_BEEN_USED

plot(res_move_mu)
plot(res_move_mu, "pi")
plot(res_move_mu, "mu")

```

This short, full trace, clearly hasn't mixed well (as is to be expected). But
while we see the effect of accept/reject sampling (Metropolis algorithm) for
`pi` with a lot of autocorrelation, the trace of `mu` shows complete
independence between successive values. While the Gibbs sampler used here is not
correct, this result is: a Gibbs sampler will be more efficient than the
classical Metropolis(-Hasting) algorithm for a given number a iterations.




<br>

## A new movement of ancestries

Moves of ancestries are done in two ways in outbreaker: by picking ancestors at
random from any prior case, and by swapping cases from a transmission
link. Here, we implement a new move, which will propose infectors which have
been infected on the same day of the current infector. As before, we will use
global variables to keep track of the resulting movements (see `N_ACCEPT` and
`N_REJECT`).


```{r, new_move_ances}

api <- get_cpp_api()

new_move_ances <- function(param, data, custom_likelihoods = NULL) {

for (i in 1:data$N) {
    current_ances <- param$alpha[i]
    if (!is.na(current_ances)) {
      ## find cases infected on the same days
      current_t_inf <- param$t_inf[current_ances]
      pool <- which(param$t_inf == current_t_inf)
      pool <- setdiff(pool, i)
      
      if (length(pool) > 0) {
        ## propose new ancestor
        current_ll <- api$cpp_ll_all(data, param, i = i, custom_likelihoods)
        param$alpha[i] <- sample(pool, 1)
        new_ll <- api$cpp_ll_all(data, param, i = i, custom_likelihoods)

        ## likelihood ratio - no correction, move is symmetric
        ratio <- exp(new_ll - current_ll)

        ## accept / reject
        if (runif(1) <= ratio) { # accept
          N_ACCEPT <<- N_ACCEPT + 1
        } else { # reject
          N_REJECT <<- N_REJECT + 1
          param$alpha[i] <- current_ances
        }
      }
    }
  }
  return(param)
}

moves <- custom_moves(new_move = new_move_ances)

```

We can now use this new move in our transmission tree reconstruction. We will
use a shorter chain than the defaults as this new move is likely to be computer
intensive.



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

N_ACCEPT <- 0
N_REJECT <- 0

set.seed(1)
res_new_move <- outbreaker(data, list(move_kappa = FALSE), moves = moves)

N_ACCEPT
N_REJECT

```


```{r, res_new_move}

plot(res_new_move)
plot(res_new_move, type = "alpha")
summary(res_new_move)

```

Results show a switch to a new mode at about 5000 iterations. Let us compare the
consensus tree to the actual one (store in `fake_outbreak$ances`):

```{r, check_new_move}

summary(res_new_move, burnin = 5000)
tree2 <- summary(res_new_move, burnin = 5000)$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)

```