# 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.

## Try the phylodyn package in your browser

Any scripts or data that you put into this service are public.

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