Bayesian <scp>StairwayPlot</scp> for Inferring Single Population Demographic Histories From Site Frequency Spectra

The StairwayPlot approach provides an elegant, flexible and powerful method to estimate complex demographic histories of single populations from site frequency spectrum data. It uses expected coalescent times to compute the expected site frequency spectrum within a multinomial likelihood function. Population sizes are allowed to vary freely between coalescent events but are constant within each interval. Here, we implement the StairwayPlot approach in the Bayesian software package RevBayes. We use approaches developed for Bayesian Skyline Plots, which include independent and identically distributed (i.i.d.) population sizes, Gaussian Markov random fields and Horseshoe Markov random fields as prior distributions on population sizes. Furthermore, we implement a recently developed approach for computing the leave-one-out cross-validation probability for efficient model selection. We compare inference from our Bayesian implementation to the original Maximum Likelihood implementation, StairwayPlot2. Our results show that our Bayesian implementation in RevBayes performs comparable to StairwayPlot2 in terms of parameter accuracy, which is expected given that both use the same underlying likelihood function. From our set of prior models, the Gaussian Markov random field prior performed best for smoothly varying demographic histories, while the Horseshoe Markov random field performs best for abruptly changing demographic histories. We conclude the study by exploring several choices often faced in empirical studies, including the estimate of the total sequence length, the assumed mutation rate, as well as biases through mis-calling ancestral alleles. We show using our empirical example that as few as 10 diploid individuals are sufficient to infer complex demographic histories, but at least 500 k single nucleotide polymorphisms (SNPs) are required.