Producing results and figures using BNPR/MCMC

In phylodyn: Statistical Tools for Phylodynamics

This vignette shows a simple and typical use of phylodyn and its MCMC tools. If we increase the number of MCMC samples, we reproduce many of the results from Lan, et al. 2014.

We start by loading the phylodyn package.

set.seed(8675309)
library(phylodyn)

We need to set the true effective population size trajectory functions, as well as the lower bounds of the trajectories on the intervals we expect to use.

traj = list(logistic_traj,exp_traj,boombust_traj,bottleneck_traj)
lower_bounds = c(10, 0.01, 0.1, 0.1)

We simulate heterochronous sampling. We sample 10 individuals at t=0, and 40 more distributed uniformly between t=0 and t=8.

samp_times = c(0, sort(runif(40, 0, 8)))
n_sampled = c(10, rep(1, 40))

We simulate genealogies based on our sample using the coalescent and the four trajectories.

gene = list()
for (i in 1:4)
{
  gene[[i]] = coalsim(samp_times = samp_times, n_sampled = n_sampled, traj = traj[[i]], lower_bound = lower_bounds[i])
}

We first use BNPR to calculate approximate marginals. We use a helper function to generate the arguments for the BNPR function.

res_BNPR = list()
for (i in 1:4)
{
  res_BNPR[[i]] = BNPR(data = gene[[i]], lengthout = 100, prec_alpha = 0.01, prec_beta = 0.01)
}

We set the number of samples and burn-in parameters. For expediency in this vignette we set them to be small. To produce more meaningful results, increase both parameters. We also use all five implemented samplers.

nsamp = 500
nburnin = 100
alg = c("HMC", "splitHMC", "MALA", "aMALA", "ESS")

We package the data in a form we can pass to mcmc_sampling. We invoke the mcmc_sampling function to run MCMC on the data.

res_MCMC = list()
for (i in 1:4)
{
  res_MCMC[[i]] = list()
  data = list(samp_times=samp_times, n_sampled=n_sampled, coal_times=gene[[i]]$coal_times)
  for (j in 1:5)
  {
    res_MCMC[[i]][[j]] = mcmc_sampling(data = data, alg = alg[j],
                                       nsamp = nsamp, nburnin = nburnin)
  }
}

We plot and compare the results from BNPR versus the MCMC samplers.

par(mfrow=c(4,5))
traj_title = c("Logistic", "Exponential", "Boombust", "Bottleneck")
color = c("green", "red", "pink", "purple", "cyan")
for (i in 1:4)
{
  for (j in 1:5)
  {
    plot_BNPR(res_BNPR[[i]], traj[[i]])
    title(traj_title[i])

    lines(res_MCMC[[i]][[j]]$med_fun, pch="", col=color[j], lwd=2.5)
    lines(res_MCMC[[i]][[j]]$low_fun, pch="", col=color[j], lwd=2, lty=3)
    lines(res_MCMC[[i]][[j]]$up_fun,  pch="", col=color[j], lwd=2, lty=3)
    legend('bottomright',c('Truth','BNPR',alg[j]),col=c('black','blue',color[j]),lwd=c(1,2.5,2.5),bty='n')
  }
}

References

1. J. A. Palacios and V. N. Minin. Integrated nested Laplace approximation for Bayesian nonparametric phylodynamics. In Proceedings of the Twenty-Eighth International Conference on Uncertainty in Artificial Intelligence, pages 726–735, 2012.

2. Lan, S., Palacios, J. A., Karcher, M., Minin, V. N., & Shahbaba, B. (2014). An Efficient Bayesian Inference Framework for Coalescent-Based Nonparametric Phylodynamics. arXiv preprint arXiv:1412.0158.

phylodyn documentation built on May 29, 2017, 1:28 p.m.