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') } }
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.
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.
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.