Improving marginal likelihood estimation for Bayesian phylogenetic model selection

Abstract

The marginal likelihood is commonly used for comparing different evolutionary models in Bayesian phylogenetics and is the central quantity used in computing Bayes Factors for comparing model fit. A popular method for estimating marginal likelihoods, the harmonic mean (HM) method, can be easily computed from the output of a Markov chain Monte Carlo analysis but often greatly overestimates the marginal likelihood. The thermodynamic integration (TI) method is much more accurate than the HM method but requires more computation. In this paper, we introduce a new method, steppingstone sampling (SS), which uses importance sampling to estimate each ratio in a series (the "stepping stones") bridging the posterior and prior distributions. We compare the performance of the SS approach to the TI and HM methods in simulation and using real data. We conclude that the greatly increased accuracy of the SS and TI methods argues for their use instead of the HM method, despite the extra computation needed.

FIGURE 1.

Three possible priors (dotted lines) for an analysis in which the likelihood is indicated by a solid line. a) Flat prior. b) Informative prior that strongly overlaps the likelihood. c) Informative prior with almost no overlap with likelihood.1

FIGURE 2.

Bias of the TI method for different numbers and distributions of β values in the standard normal example. K is the number of β intervals, and β values are evenly spaced quantiles from a Beta(α,1.0) distribution. Bias is defined as .

FIGURE 3.

RMSE of the TI and SS methods for different numbers and distributions of β values in the standard normal example. K is the number of β intervals, and β values are evenly spaced quantiles from a Beta(α,1.0) distribution.

FIGURE 4.

The topology assumed for the green plant rbcL example. The tree topology and branch lengths were estimated by maximum likelihood using the GTR+G model.

FIGURE 5.

Log marginal likelihood for three estimation methods as a function of the number of β intervals, K, for the green plant Ribulose Bisphosphate Carboxylase/Oxygenase large subunit (rbcL) example. β values are evenly spaced quantiles from a Beta(0.3,1.0) distribution. Error bars represent ±1 standard error based on 30 independent MCMC analyses.