Stepping-stone Sampling Research Articles

Diversification Models Conflate Likelihood and Prior, and Cannot be Compared Using Conventional Model-Comparison Tools.

Time-calibrated phylogenetic trees are a tremendously powerful tool for studying evolutionary, ecological, and epidemiological phenomena. Such trees are predominantly inferred in a Bayesian framework, with the phylogeny itself treated as a parameter with a prior distribution (a "tree prior"). However, we show that the tree "parameter" consists, in part, of data, in the form of taxon samples. Treating the tree as a parameter fails to account for these data and compromises our ability to compare among models using standard techniques (e.g., marginal likelihoods estimated using path-sampling and stepping-stone sampling algorithms). Since accuracy of the inferred phylogeny strongly depends on how well the tree prior approximates the true diversification process that gave rise to the tree, the inability to accurately compare competing tree priors has broad implications for applications based on time-calibrated trees. We outline potential remedies to this problem, and provide guidance for researchers interested in assessing the fit of tree models. [Bayes factors; Bayesian model comparison; birth-death models; divergence-time estimation; lineage diversification].

Improved marginal likelihood estimation via power posteriors and importance sampling

Power posteriors have become popular in estimating the marginal likelihood of a Bayesian model. A power posterior is referred to as the posterior distribution that is proportional to the likelihood raised to a power b∈[0,1]. Important power-posterior-based algorithms include thermodynamic integration (TI) of Friel and Pettitt (2008) and steppingstone sampling (SS) of Xie et al. (2011). In this paper, it is shown that the Bernstein–von Mises (BvM) theorem holds for power posteriors under regularity conditions. Due to the BvM theorem, power posteriors, when adjusted by the square root of the auxiliary constant, have the same limit distribution as the original posterior distribution, facilitating the implementation of the modified TI and SS methods via importance sampling. Unlike the TI and SS methods that require repeated sampling from the power posteriors, the modified methods only need the original posterior output and hence, are computationally more efficient. Moreover, they completely avoid the coding efforts associated with sampling from the power posteriors. Primitive conditions, under which the TI and modified TI algorithms can produce consistent estimators of the marginal likelihood, are provided. The numerical efficiency of the proposed methods is illustrated using two models.

Reporting of New tick-borne encephalitis virus strains isolated in Eastern Siberia (Russia) in 1960-2011 and explaining them in an evolutionary context using Bayesian phylogenetic inference

Tick-borne encephalitis virus (TBEV) is one of the main tick-borne viral pathogens of humans. Infection may induce signs of meningitis, encephalitis, paralysis and high fever. TBEV is well studied by molecular phylogenetic methods. The present-day implementation of Bayesian phylogenetic models allows population dynamics to be tracked, providing changes in population size that were not directly observed. However, the description of the past population dynamics of TBEV is rare in the literature. In our investigation, we provide data on the dynamics of viral genetic diversity of TBEV in Zabaikalsky Krai (Eastern Siberia, Russia) revealed by the Bayesian coalescent inference in a BEAST program. As a data set, we used the envelope (E) protein partial gene sequences (1308 nt) of 38 TBEV strains (including six “886−84-like” or Baikalian subtype strains (TBEV-B)), isolated in Zabaikalsky Krai (Eastern Siberia, Russia) in 1960–1963 and 1995–2011. To increase estimations reliability, we compared 9 model combinations by Path sampling and Stepping-stone sampling methods. It has been shown that the genetic diversity decline in the population history of TBEV in the 1950s coincides with the date of the beginning of wide dichlorodiphenyltrichloroethane forest dusting in Siberia. We assumed that the TBEV population on the territory of Siberia went through a genetic bottleneck. Also, we provide data estimating the divergence time of TBEV-B strains and indicate the specific evolution rate of an ancestor lineage of the Baikalian subtype, illustrated on a phylogenetic tree, and reconstructed under a relaxed clock model.

A66 Tracing the evolutionary history of an emerging Salmonella 4,[5],12:i:- clone in the United States

Salmonellosis is one of the leading causes of foodborne disease worldwide, with an estimated one million cases a year in the United States. Salmonella 4,[5],12: i:-, a monophasic variant of Salmonella typhimurium, is an emerging serovar that has been associated with multiple foodborne outbreaks throughout the world, mostly attributed to pig and pig products. Recently, we have demonstrated that two distinct groups of Salmonella 4,[5],12:i:- circulate in the USA and Europe, with the majority of isolates recovered during recent years belonging to an emerging multidrug-resistant clade (Elnekave et al. 2018). We applied Bayesian phylodynamic reconstruction to uncover the evolutionary history of this clade. We used a dataset of whole-genome sequences of 1446 4,[5],12:i:- isolates from different sources (livestock, human, food products, and others) from the USA (n = 752) and Europe (n = 694), collected between 2008 and 2017 and belonging to the Multilocus Subtype 34, which was predominant in the emerging clade (Elnekave et al. 2018). A subset (n = 110) of Salmonella 4,[5],12:i: isolates was then randomly selected after stratifying by location and year of isolation in order to achieve balanced sampling. Evidence of temporal signal was confirmed by looking at root-to-tip divergences using TempEst. Evolutionary hypotheses using strict and relaxed-clock models were tested using BEAST for a variety of demographic models and assuming a general time reversible substitution model. Model selection was performed by estimating Bayes Factors using path sampling and stepping-stone sampling. The selected model was then used for applying discrete trait models comparing different scenarios of transmission between locations (i.e. bidirectional symmetric/asymmetric or unidirectional). Our preliminary phylodynamic inference results indicate that the origin of this subtype was in Europe and dates back to 1990 (HPD 95%: 1984–2001). We report an exponential growth rate of 0.362 per year, which corresponds to a doubling time of 1.43 years. Our results suggest that this subtype was introduced to the US in the year 2000 (HPD 95%: 1994–2006). Phylodynamic analysis suggests that the recent increase in isolation of Salmonella 4,[5],12:i:- from different sources in the USA may be due to the exponential expansion of an emerging clone which originated in Europe and then expanded to the USA. The emergence and expansion of this serovar is of great public health importance due to the high prevalence of multidrug resistance traits found in USA isolates from this group and especially due to the presence of plasmid-mediated resistance genes for quinolones and extended spectrum cephalosporins, key antimicrobials used for the treatment of invasive Salmonella infections.

Making Steppingstones out of Stumbling Blocks: A Bayesian Model Evidence Estimator with Application to Groundwater Transport Model Selection

Bayesian model evidence (BME) is a measure of the average fit of a model to observation data given all the parameter values that the model can assume. By accounting for the trade-off between goodness-of-fit and model complexity, BME is used for model selection and model averaging purposes. For strict Bayesian computation, the theoretically unbiased Monte Carlo based numerical estimators are preferred over semi-analytical solutions. This study examines five BME numerical estimators and asks how accurate estimation of the BME is important for penalizing model complexity. The limiting cases for numerical BME estimators are the prior sampling arithmetic mean estimator (AM) and the posterior sampling harmonic mean (HM) estimator, which are straightforward to implement, yet they result in underestimation and overestimation, respectively. We also consider the path sampling methods of thermodynamic integration (TI) and steppingstone sampling (SS) that sample multiple intermediate distributions that link the prior and the posterior. Although TI and SS are theoretically unbiased estimators, they could have a bias in practice arising from numerical implementation. For example, sampling errors of some intermediate distributions can introduce bias. We propose a variant of SS, namely the multiple one-steppingstone sampling (MOSS) that is less sensitive to sampling errors. We evaluate these five estimators using a groundwater transport model selection problem. SS and MOSS give the least biased BME estimation at an efficient computational cost. If the estimated BME has a bias that covariates with the true BME, this would not be a problem because we are interested in BME ratios and not their absolute values. On the contrary, the results show that BME estimation bias can be a function of model complexity. Thus, biased BME estimation results in inaccurate penalization of more complex models, which changes the model ranking. This was less observed with SS and MOSS as with the three other methods.

Stepping-stone sampling algorithm for calculating the evidence of gravitational wave models

Bayesian statistical inference has become increasingly important for the analysis of observations from the Advanced LIGO and Advanced Virgo gravitational wave detectors. To this end, iterative simulation techniques, in particular nested sampling and parallel tempering, have been implemented in the software library LALInference to sample from the posterior distribution of waveform parameters of compact binary coalescence events. Nested sampling was mainly developed to calculate the marginal likelihood of a model but can produce posterior samples as a byproduct. Thermodynamic integration is employed to calculate the evidence using samples generated by parallel tempering but has been found to be computationally demanding. Here we propose the stepping-stone sampling algorithm, originally proposed by Xie et al. (2011) in phylogenetics and a special case of path sampling, as an alternative to thermodynamic integration. The stepping-stone sampling algorithm is also based on samples from the power posteriors of parallel tempering but has superior performance as fewer temperature steps and thus computational resources are needed to achieve the same accuracy. We demonstrate its performance and computational costs in comparison to thermodynamic integration and nested sampling in a simulation study and a case study of computing the marginal likelihood of a binary black hole signal model applied to simulated data from the Advanced LIGO and Advanced Virgo gravitational wave detectors. To deal with the inadequate methods currently employed to estimate the standard errors of evidence estimates based on power posterior techniques, we propose a novel block bootstrap approach and show its potential in our simulation study and LIGO application.

Thermodynamic integration and steppingstone sampling methods for estimating Bayes factors: A tutorial

One of the more principled methods of performing model selection is via Bayes factors. However, calculating Bayes factors requires marginal likelihoods, which are integrals over the entire parameter space, making estimation of Bayes factors for models with more than a few parameters a significant computational challenge. Here, we provide a tutorial review of two Monte Carlo techniques rarely used in psychology that efficiently compute marginal likelihoods: thermodynamic integration (Friel & Pettitt, 2008; Lartillot & Philippe, 2006) and steppingstone sampling (Xie, Lewis, Fan, Kuo, & Chen, 2011). The methods are general and can be easily implemented in existing MCMC code; we provide both the details for implementation and associated R code for the interested reader. While Bayesian toolkits implementing standard statistical analyses (e.g., JASP Team, 2017; Morey & Rouder, 2015) often compute Bayes factors for the researcher, those using Bayesian approaches to evaluate cognitive models are usually left to compute Bayes factors for themselves. Here, we provide examples of the methods by computing marginal likelihoods for a moderately complex model of choice response time, the Linear Ballistic Accumulator model (Brown & Heathcote, 2008), and compare them to findings of Evans and Brown (2017), who used a brute force technique. We then present a derivation of TI and SS within a hierarchical framework, provide results of a model recovery case study using hierarchical models, and show an application to empirical data. A companion R package is available at the Open Science Framework: https://osf.io/jpnb4.

Bayesian model evidence as a practical alternative to deviance information criterion.

While model evidence is considered by Bayesian statisticians as a gold standard for model selection (the ratio in model evidence between two models giving the Bayes factor), its calculation is often viewed as too computationally demanding for many applications. By contrast, the widely used deviance information criterion (DIC), a different measure that balances model accuracy against complexity, is commonly considered a much faster alternative. However, recent advances in computational tools for efficient multi-temperature Markov chain Monte Carlo algorithms, such as steppingstone sampling (SS) and thermodynamic integration schemes, enable efficient calculation of the Bayesian model evidence. This paper compares both the capability (i.e. ability to select the true model) and speed (i.e. CPU time to achieve a given accuracy) of DIC with model evidence calculated using SS. Three important model classes are considered: linear regression models, mixed models and compartmental models widely used in epidemiology. While DIC was found to correctly identify the true model when applied to linear regression models, it led to incorrect model choice in the other two cases. On the other hand, model evidence led to correct model choice in all cases considered. Importantly, and perhaps surprisingly, DIC and model evidence were found to run at similar computational speeds, a result reinforced by analytically derived expressions.

Thermodynamic Bayesian model comparison

Thermodynamics have been shown to have direct applications in Bayesian model evaluation. Within a tempered transitions scheme, the Boltzmann---Gibbs distribution pertaining to different Hamiltonians is implemented to create a path which links the distributions of interest at the endpoints. As illustrated here, an optimal temperature exists along the path which directly provides the free energy, which in this context corresponds to the marginal likelihood and/or Bayes factor. Estimators which have been developed under this framework are organised here using a unifying approach, in parallel with their stepping-stone sampling counterparts. New estimators are presented and the use of compound paths is introduced. As a byproduct, it is shown how the thermodynamic integral allows for the estimation of probability distribution divergences and measures of statistical entropy. A geometric approach is employed here to illustrate the importance of the choice of the path in terms of the corresponding estimator's error (path-related variance), which provides a more intuitive approach in tuning the error sources.

Hidden diversity of marine borderline lichens and a new order of fungi: Collemopsidiales (Dothideomyceta)

The fungal genus Collemopsidium comprises species that develop so-called borderline lichen symbioses with algae or cyanobacteria. Together with morphologically similar pyrenocarpous fungi it has been assigned to the family Xanthopyreniaceae. The adscription of this family to higher taxonomic ranks remain uncertain. Using sequence data of five nuclear genomic regions (nuLSU, nuSSU, tef1-α, rpb1 and rpb2) and one mitochondrial locus (mtSSU) we found that the studied representatives of this family are placed in the Dothideomyceta, yet relationships with the classes Dothideomycetes and Arthoniomycetes remain uncertain. We describe the new order Collemopsidiales to accommodate the genus Collemopsidium (paraphyletic as currently understood) and the lichenicolous genus Zwackhiomyces. Using five fungal fossils as calibrations points, we infer an age of c. 230 Mya for the crown of Collemopsidiales. Based on two molecular markers, we also provide insight into the global diversity of marine species of the genus Collemopsidium. According to the species delimitation algorithm GMYC, c. 26 putative species exist, far more than the six species recognized hitherto. We have confirmed this result by comparing the two alternative species models by means of Bayes factors, using path sampling and stepping-stone sampling algorithms to estimate the marginal likelihood of each model. Finally, our observations suggest rock-boring ability evolved in parallel in the different lineages within this group of fungi.

Estimating Divergence Times and Substitution Rates in Rhizobia.

Accurate estimation of divergence times of soil bacteria that form nitrogen-fixing associations with most leguminous plants is challenging because of a limited fossil record and complexities associated with molecular clocks and phylogenetic diversity of root nodule bacteria, collectively called rhizobia. To overcome the lack of fossil record in bacteria, divergence times of host legumes were used to calibrate molecular clocks and perform phylogenetic analyses in rhizobia. The 16S rRNA gene and intergenic spacer region remain among the favored molecular markers to reconstruct the timescale of rhizobia. We evaluate the performance of the random local clock model and the classical uncorrelated lognormal relaxed clock model, in combination with four tree models (coalescent constant size, birth–death, birth–death incomplete sampling, and Yule processes) on rhizobial divergence time estimates. Bayes factor tests based on the marginal likelihoods estimated from the stepping-stone sampling analyses strongly favored the random local clock model in combination with Yule process. Our results on the divergence time estimation from 16S rRNA gene and intergenic spacer region sequences are compatible with age estimates based on the conserved core genes but significantly older than those obtained from symbiotic genes, such as nodIJ genes. This difference may be due to the accelerated evolutionary rates of symbiotic genes compared to those of other genomic regions not directly implicated in nodulation processes.

Genealogical Working Distributions for Bayesian Model Testing with Phylogenetic Uncertainty.

Marginal likelihood estimates to compare models using Bayes factors frequently accompany Bayesian phylogenetic inference. Approaches to estimate marginal likelihoods have garnered increased attention over the past decade. In particular, the introduction of path sampling (PS) and stepping-stone sampling (SS) into Bayesian phylogenetics has tremendously improved the accuracy of model selection. These sampling techniques are now used to evaluate complex evolutionary and population genetic models on empirical data sets, but considerable computational demands hamper their widespread adoption. Further, when very diffuse, but proper priors are specified for model parameters, numerical issues complicate the exploration of the priors, a necessary step in marginal likelihood estimation using PS or SS. To avoid such instabilities, generalized SS (GSS) has recently been proposed, introducing the concept of "working distributions" to facilitate--or shorten--the integration process that underlies marginal likelihood estimation. However, the need to fix the tree topology currently limits GSS in a coalescent-based framework. Here, we extend GSS by relaxing the fixed underlying tree topology assumption. To this purpose, we introduce a "working" distribution on the space of genealogies, which enables estimating marginal likelihoods while accommodating phylogenetic uncertainty. We propose two different "working" distributions that help GSS to outperform PS and SS in terms of accuracy when comparing demographic and evolutionary models applied to synthetic data and real-world examples. Further, we show that the use of very diffuse priors can lead to a considerable overestimation in marginal likelihood when using PS and SS, while still retrieving the correct marginal likelihood using both GSS approaches. The methods used in this article are available in BEAST, a powerful user-friendly software package to perform Bayesian evolutionary analyses.

Origin and diversification of living cycads: a cautionary tale on the impact of the branching process prior in Bayesian molecular dating

BackgroundBayesian relaxed-clock dating has significantly influenced our understanding of the timeline of biotic evolution. This approach requires the use of priors on the branching process, yet little is known about their impact on divergence time estimates. We investigated the effect of branching priors using the iconic cycads. We conducted phylogenetic estimations for 237 cycad species using three genes and two calibration strategies incorporating up to six fossil constraints to (i) test the impact of two different branching process priors on age estimates, (ii) assess which branching prior better fits the data, (iii) investigate branching prior impacts on diversification analyses, and (iv) provide insights into the diversification history of cycads.ResultsUsing Bayes factors, we compared divergence time estimates and the inferred dynamics of diversification when using Yule versus birth-death priors. Bayes factors were calculated with marginal likelihood estimated with stepping-stone sampling. We found striking differences in age estimates and diversification dynamics depending on prior choice. Dating with the Yule prior suggested that extant cycad genera diversified in the Paleogene and with two diversification rate shifts. In contrast, dating with the birth-death prior yielded Neogene diversifications, and four rate shifts, one for each of the four richest genera. Nonetheless, dating with the two priors provided similar age estimates for the divergence of cycads from Ginkgo (Carboniferous) and their crown age (Permian). Of these, Bayes factors clearly supported the birth-death prior.ConclusionsThese results suggest the choice of the branching process prior can have a drastic influence on our understanding of evolutionary radiations. Therefore, all dating analyses must involve a model selection process using Bayes factors to select between a Yule or birth-death prior, in particular on ancient clades with a potential pattern of high extinction. We also provide new insights into the history of cycad diversification because we found (i) periods of extinction along the long branches of the genera consistent with fossil data, and (ii) high diversification rates within the Miocene genus radiations.Electronic supplementary materialThe online version of this article (doi:10.1186/s12862-015-0347-8) contains supplementary material, which is available to authorized users.

Exploring new dating approaches for parasites: The worldwide Apodanthaceae (Cucurbitales) as an example

Gene trees of holoparasitic plants usually show distinctly longer branch lengths than seen in photosynthetic closest relatives. Such substitution rate jumps have made it difficult to infer the absolute divergence times of parasites. An additional problem is that parasite clades often lack a fossil record. Using nuclear and mitochondrial DNA sequences of Apodanthaceae, a worldwide family of endoparasites living inside Fabaceae and Salicaceae, we compared several dating approaches: (i) an uncorrelated lognormal (UCLN) model calibrated with outgroup fossils, (ii) ages of host lineages as a maximal age in an UCLN model, (iii) user-assigned local clocks, and (iv) outgroup-fossil-calibrated random local clocks (RLC) with varying prior probabilities on the number of permitted rate changes (RLCu and RLCp models), a variable that has never been explored. The resulting dated phylogenies include all 10 species of the family, three in Australia, one in Iran, one in Africa, and the remainder in the Americas. All clock models infer a drastic rate jump between nonparasitic outgroups and Apodanthaceae, but since they distribute the rate heterogeneity differently, they result in much-different age estimates. Bayes factors using path and stepping-stone sampling indicated that the RLCp model fit poorly, while for matR, topologically unconstrained RLCu and UCLN models did not differ significantly and for 18S, the UCLN model was preferred. Under the equally well fitting models, the Apodanthaceae appear to be a relatively old clade, with a stem age falling between 65 and 81my, the divergence of Apodanthes from Pilostyles between 36 and 57my ago, and the crown age of the Australian clade 8–18my ago. In our study system, host-age calibrations did not yield well-constrained results, but they may work better in other parasite clades. For small data sets where statistical convergence can be reached even with complex models, random local clocks should be explored as an alternative to the exclusive reliance on UCLN clocks.

Bayesian tests of topology hypotheses with an example from diving beetles.

We review Bayesian approaches to model testing in general and to the assessment of topological hypotheses in particular. We show that the standard way of setting up Bayes factor tests of the monophyly of a group, or the placement of a sample sequence in a known reference tree, can be misleading. The reason for this is related to the well-known dependency of Bayes factors on model-specific priors. Specifically, when testing tree hypotheses it is important that each hypothesis is associated with an appropriate tree space in the prior. This can be achieved by using appropriately constrained searches or by filtering trees in the posterior sample, but in a more elaborate way than typically implemented. If it is difficult to find the appropriate tree sets to be contrasted, then the posterior model odds may be more informative than the Bayes factor. We illustrate the recommended techniques using an empirical test case addressing the issue of whether two genera of diving beetles (Coleoptera: Dytiscidae), Suphrodytes and Hydroporus, should be synonymized. Our refined Bayes factor tests, in contrast to standard analyses, show that there is strong support for Suphrodytes nesting inside Hydroporus, and the genera are therefore synonymized. [Bayes factor; Coleoptera; Dytiscidae; marginal likelihood; model testing; posterior odds; reversible-jump MCMC; stepping-stone sampling.]

Bayesian evolutionary model testing in the phylogenomics era: matching model complexity with computational efficiency

The advent of new sequencing technologies has led to increasing amounts of data being available to perform phylogenetic analyses, with genomic data giving rise to the field of phylogenomics. High-performance computing is becoming an indispensable research tool to fit complex evolutionary models, which take into account specific genomic properties, to large datasets. Here, we perform an extensive Bayesian phylogenetic model selection study, comparing codon and nucleotide substitution models, including codon position partitioning for nucleotide data as well gene-specific substitution models for both data types. For the best fitting partitioned models, we also compare independent partitioning with standard diffuse prior specification to conditional partitioning via hierarchical prior specification. To compare the different models, we use state-of-the-art marginal likelihood estimation techniques, including path sampling and stepping-stone sampling. We show that a full codon model best describes the features of a whole mitochondrial genome dataset, consisting of 12 protein-coding genes, but only when each gene is allowed to evolve under a separate codon model. However, when using hierarchical prior specification for the partition-specific parameters instead of independent diffuse priors, codon position partitioned nucleotide models can still outperform standard codon models. We demonstrate the feasibility of fitting such a combination of complex models using the BEAGLE library for BEAST in combination with recent graphics cards. We argue that development and use of such models needs to be accompanied by state-of-the-art marginal likelihood estimators because the more traditional and computationally less demanding estimators do not offer adequate accuracy.

Species history and divergence times of viviparous and oviparous Chinese toad-headed sand lizards (Phrynocephalus) on the Qinghai-Tibetan Plateau

The Qinghai-Tibetan Plateau (QTP) is an important biogeographical area and has recently become a focus for biodiversity studies. Phyrnocephalus lizards form a widespread Eurasian group with oviparous and viviparous reproductive modes, but two previous mtDNA studies of species from around the QTP have provided different phylogenetic hypotheses. We analysed three loci (mtDNA, RAG-1, AME) from all recognised Chinese Phrynocephalus species to reconstruct the speciation history of the group and to estimate species divergence times. The effects of mtDNA partitioning strategy on phylogenetic inference were examined. Bayes factor comparisons of marginal likelihoods (mLs) estimated using stepping-stone sampling revealed that partitioning strategy had a major impact on mL. Nevertheless, it had a negligible effect on the inferred tree topology. The impact of hard-bound uniform or equivalent soft-bound gamma speciation time calibration priors as well as the use of a fixed topology (as opposed to integration over all possible species histories) on divergence time estimation were also assessed, and found to have little impact on posterior estimates. All three gene trees and the species tree supported the hypothesis that the Chinese species form oviparous and viviparous sister clades. This was in agreement with an early mtDNA study but differed from a subsequent reanalysis of the mtDNA data. Inclusion of mtDNA from more widely distributed Phrynocephalus (from previous studies) indicates that the oviparous P. interscapularis from Central Asia lies outside the clade of Chinese viviparous and oviparous species, but that other Asian oviparous species lie within the Chinese oviparous clade. The median of the posterior on the divergence time of Chinese oviparous and viviparous species was 9.7Ma ago (95% interval: 7.2–13.0Ma ago), which coincides with major uplifting of the QTP and indicates that viviparity evolved when this clade became restricted to regions of high elevation. We also found that cladogenesis within the viviparous clade began around 5Ma ago whereas those in the oviparous clade began around 8.6Ma ago. We establish more robust estimates of divergence times and relationships within this important group and so provide improved insights into the origins of Phrynocephalus diversity across the QTP.

Make the most of your samples: Bayes factor estimators for high-dimensional models of sequence evolution

BackgroundAccurate model comparison requires extensive computation times, especially for parameter-rich models of sequence evolution. In the Bayesian framework, model selection is typically performed through the evaluation of a Bayes factor, the ratio of two marginal likelihoods (one for each model). Recently introduced techniques to estimate (log) marginal likelihoods, such as path sampling and stepping-stone sampling, offer increased accuracy over the traditional harmonic mean estimator at an increased computational cost. Most often, each model’s marginal likelihood will be estimated individually, which leads the resulting Bayes factor to suffer from errors associated with each of these independent estimation processes.ResultsWe here assess the original ‘model-switch’ path sampling approach for direct Bayes factor estimation in phylogenetics, as well as an extension that uses more samples, to construct a direct path between two competing models, thereby eliminating the need to calculate each model’s marginal likelihood independently. Further, we provide a competing Bayes factor estimator using an adaptation of the recently introduced stepping-stone sampling algorithm and set out to determine appropriate settings for accurately calculating such Bayes factors, with context-dependent evolutionary models as an example. While we show that modest efforts are required to roughly identify the increase in model fit, only drastically increased computation times ensure the accuracy needed to detect more subtle details of the evolutionary process.ConclusionsWe show that our adaptation of stepping-stone sampling for direct Bayes factor calculation outperforms the original path sampling approach as well as an extension that exploits more samples. Our proposed approach for Bayes factor estimation also has preferable statistical properties over the use of individual marginal likelihood estimates for both models under comparison. Assuming a sigmoid function to determine the path between two competing models, we provide evidence that a single well-chosen sigmoid shape value requires less computational efforts in order to approximate the true value of the (log) Bayes factor compared to the original approach. We show that the (log) Bayes factors calculated using path sampling and stepping-stone sampling differ drastically from those estimated using either of the harmonic mean estimators, supporting earlier claims that the latter systematically overestimate the performance of high-dimensional models, which we show can lead to erroneous conclusions. Based on our results, we argue that highly accurate estimation of differences in model fit for high-dimensional models requires much more computational effort than suggested in recent studies on marginal likelihood estimation.

Delayed colonisation of Acacia by thrips and the timing of host-conservatism and behavioural specialisation.

BackgroundRepeated colonisation of novel host-plants is believed to be an essential component of the evolutionary success of phytophagous insects. The relative timing between the origin of an insect lineage and the plant clade they eat or reproduce on is important for understanding how host-range expansion can lead to resource specialisation and speciation. Path and stepping-stone sampling are used in a Bayesian approach to test divergence timing between the origin of Acacia and colonisation by thrips. The evolution of host-plant conservatism and ecological specialisation is discussed.ResultsResults indicated very strong support for a model describing the origin of the common ancestor of Acacia thrips subsequent to that of Acacia. A current estimate puts the origin of Acacia at approximately 6 million years before the common ancestor of Acacia thrips, and 15 million years before the origin of a gall-inducing clade. The evolution of host conservatism and resource specialisation resulted in a phylogenetically under-dispersed pattern of host-use by several thrips lineages.ConclusionsThrips colonised a diversity of Acacia species over a protracted period as Australia experienced aridification. Host conservatism evolved on phenotypically and environmentally suitable host lineages. Ecological specialisation resulted from habitat selection and selection on thrips behavior that promoted primary and secondary host associations. These findings suggest that delayed and repeated colonisation is characterised by cycles of oligo- or poly-phagy. This results in a cumulation of lineages that each evolve host conservatism on different and potentially transient host-related traits, and facilitates both ecological and resource specialisation.

Improving the Accuracy of Demographic and Molecular Clock Model Comparison While Accommodating Phylogenetic Uncertainty

Recent developments in marginal likelihood estimation for model selection in the field of Bayesian phylogenetics and molecular evolution have emphasized the poor performance of the harmonic mean estimator (HME). Although these studies have shown the merits of new approaches applied to standard normally distributed examples and small real-world data sets, not much is currently known concerning the performance and computational issues of these methods when fitting complex evolutionary and population genetic models to empirical real-world data sets. Further, these approaches have not yet seen widespread application in the field due to the lack of implementations of these computationally demanding techniques in commonly used phylogenetic packages. We here investigate the performance of some of these new marginal likelihood estimators, specifically, path sampling (PS) and stepping-stone (SS) sampling for comparing models of demographic change and relaxed molecular clocks, using synthetic data and real-world examples for which unexpected inferences were made using the HME. Given the drastically increased computational demands of PS and SS sampling, we also investigate a posterior simulation-based analogue of Akaike's information criterion (AIC) through Markov chain Monte Carlo (MCMC), a model comparison approach that shares with the HME the appealing feature of having a low computational overhead over the original MCMC analysis. We confirm that the HME systematically overestimates the marginal likelihood and fails to yield reliable model classification and show that the AICM performs better and may be a useful initial evaluation of model choice but that it is also, to a lesser degree, unreliable. We show that PS and SS sampling substantially outperform these estimators and adjust the conclusions made concerning previous analyses for the three real-world data sets that we reanalyzed. The methods used in this article are now available in BEAST, a powerful user-friendly software package to perform Bayesian evolutionary analyses.