Fast Bayesian Choice of Phylogenetic Models: Prospecting Data Augmentation–Based Thermodynamic Integration

With the increasing number of substitution models being proposed over recent years, there is a need for accurate and fast methods for performing model selection in phylogenetics. However, the most popular computational approach used in Bayesian phylogenetic contexts—the Bayes factor approximated with the harmonic mean estimator (HME) (Newton and Raftery 1994)—has proved unreliable (Lartillot and Philippe 2006, Fan et al. 2011, Xie et al. 2011). Here, we discuss recent advances in reliable computational methods based on thermodynamic integration principles for Bayesian model selection and emphasize the potential of data augmentation–based methods for producing fast and accurate results. The HME to the marginal likelihood was proposed as a method requiring nothing more than the output of the MCMC system for sampling from the posterior in Equation (3) (Newton and Raftery 1994). The popularity of the HME, as exemplified in the recent phylogenetics literature (e.g., Bleidorn et al. 2007, Brown and Lemmon 2007, Alekseyenko et al. 2008, Praz et al. 2008, Hampl et al. 2009), is undoubtedly due to its methodological simplicity and to its availability in a number of user-friendly programs, such as MrBayes (Ronquist and Huelsenbeck 2003) and BEAST (Drummond and Rambaut 2007). This continued practice is unfortunate, however, given the problems with the HME. Lartillot and Philippe (2006), for instance, convincingly showed that the HME and its so-called stabilized version are unreliable. They proposed the use of another MCMC-based method, called thermodynamic integration (TI), and demonstrated its reliability in analytical contexts. Moreover, Lartillot and Philippe (2006) showed that Bayes factors computed with the HME and TI can produce drastically different results in amino acid substitution model comparisons, to the point of leading to opposite conclusions in model ranking. Unfortunately, the reliability of TI comes with the cost of greater computing time. Lartillot and Philippe (2006) suggested that TI could require a computational time about one order of magnitude greater than that needed for a plain MCMC run under the most complex of the models being evaluated. It seems likely that TI has not been widely adopted partly due to this increased computing cost. Recent efforts have been made by Xie et al. (2011) and Fan et al. (2011) to speed up TI-based calculations with an approach they call stepping-stone (SS). Their basic idea takes advantage of both TI-inspired schemes and importance sampling, in a manner showing great computational promise. In particular, their approach enables reliable quantitative Bayesian model ranking wherever substitution models are tractable in the full pruning-based (also called peeling-based or dynamic programming–based) Bayesian MCMC framework. The pursuit of general estimation approaches (such as TI and SS) for performing Bayesian model ranking is illustrative of a tendency among MCMC developers and practitioners: that of striving for computational systems that are independent of the particular modeling context. However, with regard to the means of computing Bayes factors, so long as accurate methods are used (thus excluding the HME), generality is not in itself necessary. For now, given the technical challenge of the calculation, it might be warranted to use less general approaches in some modeling contexts, as we illustrate below. Following Lartillot and Philippe (2006), Fan et al. (2011), and Xie et al. (2011), we demonstrate the issue at stake when using the HME approach to model ranking with a set of codon models described in detail in Appendix 1. The codon substitution models of interest to us here build upon the approach proposed by Muse and Gaut (1994). We write MG to refer to a codon substitution formulation based only on a set of nucleotide-level mutational parameters (see Appendix 1). As presented by Yang and Nielsen (2008), the models of interest are based on the mutation–selection principle and parameterize a substitution process in terms of mutational propensities combined with a fixation factor that distinguishes between the types of (nonsynonymous) events. Based on the model nomenclature explained in Appendix 1, MG-MutSelYN is defined as a global codon substitution model, where all model parameters are common to all codon sites in the alignment. At the other end of the spectrum is the model referred to as MG-MutSelHB, where each codon site has its own fixation factors, whereas all sites share the same mutational parameters (inspired by Halpern and Bruno 1998). Between these two extremes, the models denoted as MG-MutSelC20, MG-MutSelC40, and MG-MutSelC60 are empirical mixtures of 20, 40, and 60 components, which all share the same mutational parameters, but have distinct parameters governing fixation factors (as presented in Rodrigue et al. 2010). These mixture models provide a simple, crude, yet convenient compromise between the global homogeneous MG-MutSelYN model and the site-specific MG-MutSelHB model. The importance of using reliable methods for computing Bayes factors in this mutation–selection modeling framework is immediately apparent from a simple simulation experiment. Using the alignment of 17 vertebrate globin genes taken from Yang et al. (2000) (and subsequently studied repeatedly by many authors), along with the tree topology used therein, we ran a plain MCMC sampling under the MG model. For 12 random draws from the posterior, we simulated new data sets conditional on those parameters (under the same MG model). For the resulting replicate data sets, we computed the natural log Bayes factors contrasting the MG-MutSelYN model to the MG model using both the HME approach and a TI model switch approach (see Rodrigue et al. 2008a). We repeated the calculations 10 times for each replicate, and the means and standard deviations of the natural log Bayes factors are shown in Table 1. Under these simulation conditions, based on the simple MG model, we expect that an appropriate model selection approach should favor the genuine MG model, thereby disfavoring the over-parameterized MG-MutSelYN model. Our results reveal that the Bayes factors computed with HME would lead one to believe that the MG and MG-MutSelYN are nearly equivalent (with natural log Bayes factors relatively close to 0). In contrast, with the TI calculations, the log Bayes factors are systematically negative, in strong favor of the simpler correct model. Although the standard deviations obtained are greater for the HME than for the TI calculations, the problem with the HME does not appear to be one of a trivial Monte Carlo error. Mean and standard deviation (from 10 distinct calculations) of the natural logarithm of the Bayes factor in favor of the MG-MutSelYN model over the MG model on simulated data Notes: A negative value corresponds to a Bayes factor in favor of the simpler MG model, which was used to simulate the data. Results are reported for the HME and TI. Mean and standard deviation (from 10 distinct calculations) of the natural logarithm of the Bayes factor in favor of the MG-MutSelYN model over the MG model on simulated data Notes: A negative value corresponds to a Bayes factor in favor of the simpler MG model, which was used to simulate the data. Results are reported for the HME and TI. The estimation of Bayes factors with TI in the mutation–selection codon substitution context is tractable under a model based on a single substitution matrix such as MG-MutSelYN (Rodrigue et al. 2008b), albeit still requiring weeks of computation for small single-gene data sets on typical desktop computers. The same TI-based estimation of Bayes factor under the site heterogeneous MG-MutSelC20, MG-MutSelC40, or MG-MutSelC60 models pushes the limits of tractability, requiring months of computing time on a modern desktop computer (Rodrigue et al. 2010). Although the computational demand could be reduced by using SS instead of TI, it would likely still not allow for the calculation of Bayes factors involving the MG-MutSelHB model, which specifies as many codon substitution matrices as there are codon sites in the alignments. Indeed, the MG-MutSelHB model remains essentially intractable within a full pruning-based MCMC sampling system, and the applicability of such rich mutation–selection models has only recently become feasible thanks to another methodological development receiving increasing attention in Bayesian phylogenetic contexts: data augmentation (e.g., Lartillot 2006, Mateiu and Rannala 2006, de Koning et al. 2010). A difficulty for Bayesian model ranking with approaches employing the data augmentation scheme described above is that the best-performing TI method (the model switch approaches, Lartillot and Philippe 2006, Rodrigue et al. 2008a) and the best-performing SS method (the generalized stepping-stone, Fan et al. 2011) rely on sampling from distributions involving two distinctly defined substitution processes, via a “morphing” between the posteriors under the two models of interest. However, generating detailed substitution realizations along phylogenetic trees conditional on a morphing of posteriors implicating two substitution processes is not immediately obvious. In the terminology of Xie et al. (2011) and Fan et al. (2011), detailed data augmentation of the sort described above is not directly compatible with sampling from “power posteriors.” There may be several different approaches to exploiting data augmentation within TI or SS calculations, depending on the specific type of sampling system. Previous works using data augmentation within a TI sampler have done so because no other alternatives were available (e.g., Rodrigue et al. 2006, Baele et al. 2008). We emphasize here that such a methodology is relevant even when the full pruning based sampler is technically possible. As an example, in Appendix 2, we describe one way of employing detailed data augmentation within a TI-MCMC framework for mutation–selection models. The basic idea is to perform a morphing of the substitution process itself. By using this method to analyze the same globin data set as above, we computed the natural log Bayes factors contrasting the mutation–selection models described in Appendix 1 against the MG model (used as a reference model). The approach yields precise results (Table 2), in far less computing time than the full pruning-based model switch TI. For instance, the full pruning-based TI performed in Rodrigue et al. (2010) to compute natural log Bayes factors of the MG-MutSelC20, MG-MutSelC40, and MG-MutSelC60 models against MG required nearly 4 months of computing time per run (at least two runs were performed to assess precision) and returned values of 236, 235, and 269 all within about one log unit. Here, on the other hand, with data augmentation and TI, the results for the same model comparisons (see Table 2) required less than 12 hours per run. This considerable speedup allowed us to perform 10 runs (5 bidirectional runs), which returned results all within half a log unit. Note that these 10 runs were performed to assess the precision of the calculations. However, in normal practice with the methods, performing only a few (say, one bidirectional pair of calculations) would likely be sufficient. Mean and standard deviation (from 10 distinct calculations) of the natural logarithm of the Bayes factor for the models considered, with MG used as a reference Note: Results are reported for the HME and data augmentation–based thermodynamic integration (DA-TI). Mean and standard deviation (from 10 distinct calculations) of the natural logarithm of the Bayes factor for the models considered, with MG used as a reference Note: Results are reported for the HME and data augmentation–based thermodynamic integration (DA-TI). The Bayes factors computed with the HME and TI methods are once again quite different. Most strikingly, calculations done with the HME would lead one to elect MG-MutSelHB as the most appropriate model, whereas TI calculations give favor to the MG-MutSelC60 model. These results corroborate previous studies (Lartillot and Philippe 2006, Fan et al. 2011, Xie et al. 2011), indicating that Bayes factors computed using the HME tend to mislead one into selecting the most parameter-rich model. Note that this tendency is not a bias of the Bayes factor itself, but rather points to a serious issue in the way it is approximated with the HME. Data augmentation–based posterior sampling is still relatively new to phylogenetic modeling contexts, and new methods are being proposed regularly (e.g., Pedersen and Jensen 2001, Nielsen 2002, Robinson et al. 2003, Rodrigue et al. 2008b, de Koning et al. 2010). Reliable methods for Bayesian model ranking, that is, excluding HME, have received less attention. Methods such as reversible jump sampling (Green 1995) or the Dickie-Savage ratio (Verdinelli and Wasserman 1995) have rarely been used in phylogenetics beyond the seminal papers introducing them, respectively, by Huelsenbeck et al. (2004) and Suchard et al. (2001). Each method poses its own set of technical challenges and limitations, and the basic thermodynamic method (Lartillot and Philippe 2006) is no exception, where the computational costs have likely discouraged many. Recent work by Xie et al. (2011) and Fan et al. (2011), however, could help revive this area. The method we present in Appendix 2 also shows promise and is illustrative of a context-specific design of TI for the mutation–selection codon models. Our method, however, has a number of limitations. For instance, additional computer code needs to be written to accommodate the morphing between the particular substitution models of interest. Also, the data augmentation approach employed here has the weakness of working only on a fixed topology. This is mainly because no MH kernels for performing joint moves on topology and data augmentation have yet been proposed. One potential alternative would be to revert to full pruning-based moves on the tree topology, with calculations dispatched to multiple processors or to the processing units present on most graphic cards (Suchard and Rambaut 2009)—provided that the hardware improves to allow for double precision calculations. Many possible sampling approaches can readily be envisaged, and future studies could be based on exploiting the advantages of various existing ideas, such as approximating Bayes factors via combinations of TI-based approaches (Lartillot and Philippe 2006; Fan et al. 2011; Xie et al. 2011), Savage-Dickie ratios (Suchard et al. 2001), and Laplace approximations (e.g., Rodrigue et al. 2007). As we have explored here, using data-augmentation ideas within such methodologies could be computationally fruitful, and approaches other than the full Gibbs updating employed here (e.g., de Koning et al. 2010) should also be investigated in such contexts. Altogether, through advances along these lines, there is hope to eventually achieve a more efficient and reliable computational framework for comparing hypotheses under a given model, as well as model ranking in Bayesian phylogenomic analysis. Supplementary material, including data files and online-only appendices, can be found at http://www.sysbio.oxfordjournals.org/. This work was supported by the Natural Sciences and Engineering Research Council of Canada (N.R., S.A.-B.), the Canada Foundation for Innovation (S.A.-B.), and the Quebec Centre for Biodiversity Science (N.R.). 𝒜: a and b are synonymous and differ only at codon position c; ℬ: a and b are nonsynonymous and differ only at codon position c; ϱ = (ϱlm)1 ≤ l,m ≤ 4 is a set of symmetrical nucleotide relative exchangeability parameters, with the arbitrary constraint ∑1 ≤ l < m ≤ 4ϱlm = 1; , represents a set of global nucleotide equilibrium propensities; β is a parameter to be used for the TI: when β = 0, the model amounts to a simple model without selection, which we refer to as MG (because it takes inspiration from Muse and Gaut 1994); when β = 1, we have the mutation–selection model of interest, denoted MG-MutSelYN (because it includes the mutation–selection approach used by Yang and Nielsen 2008); S is the scaled selection coefficient (scaled by twice the effective population size) associated with an amino acid replacing mutation. As described in Rodrigue et al. (2010), S = ln(φf(b)/φf(a)), where f(a) returns an index (ranging from 1 to 20) of the amino acid encoded by codon a, and φ = (φn)1 ≤ n ≤ 20, with , is a set amino acid preference parameters. Note that when S or β (or their product) approaches 0, the fixation factor approaches 1, but for numerical stability, we use the third-order Taylor approximation of the exponential term in the denominator in such cases. This leads to a fixation factor approximated as . We also use empirical mixtures models, designated as MG-MutSelC20, MG-MutSelC40, and MG-MutSelC60, constructed from 20, 40 and 60 different φ vectors, fixed to the values obtained by Le et al. (2008). The free parameters introduced by these models only consist of weights associated with the components of the mixture. However, we use a sampling system based on an allocation variable, specifying a configuration in which each site is affiliated to one of the 20, 40, or 60 different components (e.g., see Le et al. 2008). Performing MH updates to the allocation variable implicitly integrates away the weight parameters. Finally, we also explore a model in which each site has its own φ vector. Because it is inspired from Halpern and Bruno (1998), we refer to it as MG-MutSelHB. We used the same priors described in Rodrigue et al. (2010). For the MG-MutSelHB model, not included in the latter study, we used independent flat Dirichlet priors on each vector (we did not explore a flexible hyperprior structure controlling the overall mean and variance of site-specific variables). zij is the number of substitution events for site i between nodes jup and j; tijk is the timing of the kth event along branch j at site i, with tij0 = 0; sijk − 1 and sijk are the codon states before and after event k along branch j for site i, with sij0 = sijup and sijzij = sij; is the rate away from codon state a; and λj is the length of branch j. Again for reasons of numerical stability in cases where βS approaches 0 (in particular at one end of the TI where β is very small), we use third-order Taylor approximations. We used sampling methods described elsewhere (e.g., Rodrigue et al. 2008b, Rodrigue et al. 2009, Rodrigue et al. 2010) and used K = 10,000 (δβ = 0.0001). Moreover, we performed bidirectional integrations, which consist in running two independent calculations, one with β going from 0 to 1 and the other with β going from 1 to 0 (our 10-fold calculations using the proposed data augmentation–based TI actually consist of five bidirectional calculations). We wish to thank Nicolas Lartillot for discussions and comments on the manuscript, as well as Mark Holder, David Posada, Vladimir Minin, and an anonymous reviewer.

Fetched live from OpenAlex and de-inverted. Abstracts are not stored in this database: the inverted indexes are 8.6 GB of the frame’s 9.3 GB of text, and the host has 13 GB free.