Many anthropologists focus on inferring the historical demography of human populations. One aspect of population history that has been of particular interest, especially with regard to human history, is changes in population size. The coalescent, as originally conceived, assumes a population does not change in size. We know that for many populations and species this is not true, and could not be. There have been many extensions of the coalescent to include populations that grow or shrink with deterministic demographic functions (Griffiths and Tavare 1994; Donnelly and Tavare 1995; Wilson and Balding 1998; Beaumont 1999), or suddenly shift between demographic functions (e.g., Shapiro et al. 2004). These extensions of the basic coalescent models may be used to reveal the past dynamics of a population history.

Perhaps the most popular demographic model in use today is the Bayesian skyline plot (BSP), which allows the effective population size to change in a piecewise fashion at coalescent events (Ho and Shapiro 2011; Drummond et al. 2012). The predecessors to the BSP, the classic (Pybus et al. 2000) and generalized skyline plots (Strimmer and Pybus 2001), estimate changes in effective population size over time based on a single genealogy in a likelihood framework. The BSP improves these models by estimating effective population size from a distribution of genealogies generated using MCMC algorithms. This both integrates over the uncertainty in genealogies as well as allowing for the calculation of credibility intervals for effective population size backward in time (Heled and Drummond 2008). Another advantage of a BSP, perhaps one that is just as pivotal for its popularity among authors, is the ability to render its output in a visually pleasing graphical format. This is the familiar plot found at the center of many studies of historical human population dynamics (e.g., Kitchen et al. 2008; Mulligan et al. 2008; Atkinson et al. 2009).

Under a BSP demographic model, effective population size is allowed to change an arbitrary number of times. The BSP model assumes that effective population size remains constant between change points, but can instantaneously change at coalescent events (Drummond et al. 2005). The change points are determined by grouping neighboring coalescent events such that each group is associated with a single constant population size that persists across all coalescent events, with changes occurring at the transition from one group to another. The minimum number of groups is 1, which reduces the BSP to a constant population demographic function, and the maximum is n-1 for a sample of n individuals, which provides for as many changes in N_e as there are coalescent events. The number of groups is fixed a priori, and though for most datasets an intermediate number of groups (such as 5 or 10) does not affect the analysis, an excessive number of groups will inhibit efficient MCMC mixing, and too few groups will not capture complex population histories. In these extreme cases, it is necessary to evaluate manually how the number of groups affects the fit of the model, which is often a slow and lengthy process.

It should be noted that a BSP is not the only Bayesian demographic model prior available in BEAST. More stringent model priors exist, such as constant population size, exponential growth, logistic growth, or expansion growth. These models can even be combined manually, and tested against each other (Pybus and Rambaut 2002). These simpler models often fit the data better than models allowing for “free” population size change, such as a BSP. However, testing each and every simple model can be extremely time consuming, with no guarantee that an adequately fitting model will be found (Drummond et al. 2005). The BSP provides an easy alternative, in which user input is reduced to deciding on the number of groups. It is noteworthy that, as a BSP is in itself a model prior, it is recommended that the fit of the empirical data to the model be tested by performing a Bayes factor model comparison with at least a few other demographic priors, such as the constant population model and the exponential model. This comparative step is important to ensure that biological conclusions are driven by the data, and not by prior model selection.

The nature of the basic BSP, in which the number of transitions in effective population size is determined a priori with no governing principle, can be problematic. Firstly, a poor choice in number of groups may lead to large credibility intervals or even prevent convergence of the analysis on an accurate estimate (Heled and Drummond 2008). Secondly, a few large transitions in effective population size between groups (in the fashion of steps) are an artificial representation of the historical reality of a natural population, in which transitions are expected to be gradual. To address this problem, extensions on the BSP have been made to specifically remove the necessity of a strong prior determining the number of transitions over evolutionary time. The Bayesian Skyride method (Minin et al. 2008) accomplished this by modifying the classic skyline plots (from Pybus et al. 2000) by imposition of a Gaussian Markov random field (GMRF) smoothing prior on the algorithm which estimates the population size trajectory, and “smooths” the abrupt transitions between intervals. An extended Bayesian skyline plot (EBSP, Heled and Drummond 2008) is a modification to the BSP supported in BEAST, in which the genetic data are referenced at each coalescent event to estimate a new effective population size by means of Bayesian stochastic variable selection. In an EBSP, transitions in effective population size are not reported for each coalescent event, but instead, there is a given probability calculated from the likelihood function at each transition, for which the effective population will either remain the same as the previous interval, or it will change to reflect the newly calculated parameter.

By their nature, Bayesian analyses revolve around model improvement as means to approximate natural processes, even when data are uninformative or incomplete. Under these circumstances, a poorly resolved likelihood function will still yield parameter estimation, which will be driven by the prior distribution rather than the data, resulting in incorrect biological inference (Konigsberg and Frankenberg 2013). For this reason, it is important to test the statistical power of any analysis. A commonly used metric of statistical power is the Effective Sample Size (ESS), which estimates the average number of independent (non-correlated) data points in the posterior distribution of sampled genealogies, ensuring that the MCMC chain has sampled a diverse mix of genealogies. While there is no hard-limit on how large an appropriate ESS should be, values under 200 are not recommended (Kuhner 2009). Software packages such as TRACER (Drummond and Rambaut 2007) provide the tools for ESS estimation, and other simple qualitative estimations of analysis appropriateness, including comparisons of the outputs of multiple independent replicate analyses, which provides a visual check of the convergence of posterior distributions on similar values across the runs.

Bayesian Skyline plots in BEAST can be extremely powerful tools for demographic analysis based on genetic data:

• BSPs can condense complex data into simple-to-interpret displays which are easily accessible for anyone familiar with basic plots.
• Thoughtful prior selection is fundamental. Inadequate priors may drive the analysis into converging on intractable likelihoods, in which case the BSP will still provide parameter estimations, which reflect prior selection instead of the information contained in the data.
• Post hoc scrutiny of prior appropriateness is paramount, in particular estimating (and reporting) the ESS metric. Replication and contrast of independent replicate analyses is also recommended.