This vignette shows a simple and typical use of
phylodyn and its main functions BNPR and BNPR_PS.
We start by loading the
We need to set the true effective population size trajectory function. For this example, we choose a seasonal growth/collapse trajectory (already implemented in
traj = logistic_traj
We set the end of the sampling interval and the number of samples we should expect.
samp_end = 48 nsamps = 500
We calculate the proportionality constant necessary to produce the correct expected number of samples. We then produce a set of samples by drawing from an inhomogeneous Poisson process with intensity proportional to the effective population size trajectory.
Cprop = nsamps/integrate(traj_beta, 0, samp_end, traj=traj, beta=1)$value samp_times = pref_sample(traj, c=Cprop, lim=c(0,samp_end), beta=1)
We simulate a genealogy based on our sample using the coalescent.
gene = coalsim(samp_times = samp_times, n_sampled = rep(1, length(samp_times)), traj = traj, lower_bound = 10)
We use BNPR and BNPR_PS to calculate approximate marginals (without and with a sampling model).
res_BNPR = BNPR(data = gene, lengthout = 100) res_BNPR_PS = BNPR_PS(data = gene, lengthout = 100)
We plot the results.
par(mfrow=c(1,2)) plot_BNPR(res_BNPR, traj = traj, main="BNPR", ylim = c(1, 700)) plot_BNPR(res_BNPR_PS, traj = traj, main="BNPR-PS", ylim = c(1, 700))
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.
M. S. Gill, P. Lemey, N. R. Faria, A. Rambaut, B. Shapiro, and M. A. Suchard. Improving Bayesian population dynamics inference: a coalescent-based model for multiple loci. Molecular biology and evolution, 30(3):713–724, 2013.
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.