An alternative marginal likelihood estimator for phylogenetic models

Bayesian phylogenetic methods are generating noticeable enthusiasm in the field of molecular systematics. Many phylogenetic models are often at stake and different approaches are used to compare them within a Bayesian framework. The Bayes factor, defined as the ratio of the marginal likelihoods of two competing models, plays a key role in Bayesian model selection. We focus on an alternative estimator of the marginal likelihood whose computation is still a challenging problem. Several computational solutions have been proposed none of which can be considered outperforming the others simultaneously in terms of simplicity of implementation, computational burden and precision of the estimates. Practitioners and researchers, often led by available software, have privileged so far the simplicity of the harmonic mean estimator (HM) and the arithmetic mean estimator (AM). However it is known that the resulting estimates of the Bayesian evidence in favor of one model are biased and often inaccurate up to having an infinite variance so that the reliability of the corresponding conclusions is doubtful. Our new implementation of the generalized harmonic mean (GHM) idea recycles MCMC simulations from the posterior, shares the computational simplicity of the original HM estimator, but, unlike it, overcomes the infinite variance issue. The alternative estimator is applied to simulated phylogenetic data and produces fully satisfactory results outperforming those simple estimators currently provided by most of the publicly available software.

💡 Research Summary

Bayesian phylogenetics has become a cornerstone of modern molecular systematics, largely because it provides a coherent framework for estimating evolutionary trees while simultaneously quantifying uncertainty. Central to Bayesian model comparison is the Bayes factor, which is defined as the ratio of marginal likelihoods (also called Bayesian evidence) for two competing models. Unfortunately, marginal likelihoods are high‑dimensional integrals that cannot be evaluated analytically for realistic phylogenetic models, and their estimation remains a computational bottleneck.

Historically, the community has relied on two very simple estimators: the harmonic mean (HM) estimator and the arithmetic mean (AM) estimator. The HM estimator uses only posterior samples from a standard MCMC run and computes the inverse of the average of the reciprocal likelihoods. Its appeal lies in its trivial implementation, but it suffers from an infinite‑variance problem: when the likelihood surface has heavy tails, a few posterior draws with extremely low likelihood dominate the sum, leading to wildly unstable estimates. The AM estimator, which averages the likelihood over prior draws, is equally problematic because it is highly sensitive to the choice of prior and often requires a prohibitively large number of prior samples to achieve reasonable precision.

More sophisticated approaches such as thermodynamic integration, stepping‑stone sampling, and path sampling have been proposed to overcome these deficiencies. While they generally produce more accurate evidence estimates, they demand additional MCMC chains at intermediate power posteriors, careful temperature schedules, and substantial programming effort. Consequently, many practitioners continue to use HM or AM simply because they are already built into popular software packages such as MrBayes and BEAST.

In this context, the authors introduce a Generalized Harmonic Mean (GHM) estimator that retains the computational elegance of the original HM while eliminating its infinite‑variance pathology. The key idea is to replace the implicit proposal distribution of the HM (which is the posterior itself) with a deliberately broader “importance” distribution (h(\theta)). The marginal likelihood is then estimated as

\