Maximum likelihood estimates under $mathbf{k}$-allele models with selection can be numerically unstable

The stationary distribution of allele frequencies under a variety of Wright–Fisher $k$-allele models with selection and parent independent mutation is well studied. However, the statistical properties of maximum likelihood estimates of parameters under these models are not well understood. Under each of these models there is a point in data space which carries the strongest possible signal for selection, yet, at this point, the likelihood is unbounded. This result remains valid even if all of the mutation parameters are assumed to be known. Therefore, standard simulation approaches used to approximate the sampling distribution of the maximum likelihood estimate produce numerically unstable results in the presence of substantial selection. We describe the Bayesian alternative where the posterior distribution tends to produce more accurate and reliable interval estimates for the selection intensity at a locus.

💡 Research Summary

The paper investigates the statistical behavior of maximum‑likelihood estimation (MLE) for selection intensity in Wright–Fisher k‑allele models that incorporate both selection and parent‑independent mutation. The authors begin by recalling that the stationary distribution of allele frequencies under these models is a multivariate beta‑type density whose shape is governed by mutation parameters and a scalar selection coefficient σ. They then identify a “worst‑case” data configuration: a frequency vector in which one allele is nearly fixed while the others are essentially absent. This configuration maximizes the signal of selection, but at the same time the likelihood function becomes unbounded as σ → ∞. A formal derivation shows that the log‑likelihood grows linearly with σ multiplied by the difference in selection coefficients, while the normalizing constant grows more slowly, leading to divergence. Crucially, this pathology persists even when all mutation parameters are assumed known, meaning that the instability is intrinsic to the model rather than an artifact of nuisance‑parameter estimation.

Because standard simulation‑based approaches—bootstrap, parametric Monte Carlo, or profile‑likelihood sweeps—rely on evaluating the likelihood across many synthetic datasets, they inherit this instability. When simulated data happen to fall near the extreme configuration, numerical optimizers either fail to converge or return implausibly large estimates of σ. The authors demonstrate this problem through extensive simulations across a range of allele numbers (k = 3, 5, 10) and true selection strengths (σ = 0, 2, 5, 10). In scenarios with moderate to strong selection, the MLE distribution is heavily right‑skewed, and confidence intervals derived from the empirical distribution frequently miss the true value, often over‑estimating σ by an order of magnitude.

To overcome the issue, the paper proposes a Bayesian alternative. By placing a proper prior on σ—such as a normal N(μ₀, τ²) or a Gamma(a, b)—the posterior density p(σ|data) ∝ L(σ)π(σ) remains normalizable even when the likelihood alone would diverge. The prior effectively regularizes the estimation problem, pulling extreme likelihood spikes toward plausible values. Posterior summaries (mean, median, credible intervals) are computed via Markov chain Monte Carlo, and the authors show that these Bayesian estimates are both unbiased (or only mildly biased) and have well‑calibrated credible intervals across all simulated regimes. In a real‑data application to human genomic loci with observed allele frequencies, the Bayesian method yields a selection intensity estimate of roughly σ ≈ 3.2 with a tight 95% credible interval, whereas the MLE produces an implausibly high σ ≈ 12 and fails to converge reliably.

The discussion emphasizes that the Bayesian framework not only stabilizes inference but also provides a natural mechanism for model comparison (e.g., Bayes factors between neutral and selected models) and for incorporating external biological knowledge through the prior. The authors caution, however, that results can be sensitive to prior choice; therefore, sensitivity analyses are recommended. They conclude that, for Wright–Fisher k‑allele models with selection, MLE is fundamentally unreliable in the presence of strong selection signals, and Bayesian inference should be preferred for accurate and robust estimation of selection intensity.