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:
custom...
functionNote 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.
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.
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
:
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:
## 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):
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:
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.
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)
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:
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
:
## 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) model, by disabling specific components.
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:
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)
plot(res_null) plot(res_null, "pi") plot(res_null, "mu")
By typical MCMC standards, these traces look appalling, 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:
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):
## 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 between cases
(kappa
), only constrained by default settings to be between 1 and 5, and for
the infection dates (t_inf):
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:
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")))
We can use data and likelihood customisation to change the default outbreaker2 model into a Wallinga & Teunis (1994) 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:
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:
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:
set.seed(1) res_wt <- outbreaker(data = onset_data, config = wt_config, likelihoods = wt_model)
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:
## 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")
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. 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.
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:
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 above. 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.
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.
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:
## 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:
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.
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
).
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.
N_ACCEPT <- 0 N_REJECT <- 0 set.seed(1) res_new_move <- outbreaker(data, list(move_kappa = FALSE), moves = moves) N_ACCEPT N_REJECT
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
):
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.