Efficient Simulation and Likelihood Methods for Non-Neutral Multi-Allele Models

2.1. Neutral K-allele model

In a population with K allelic types, evolving for a long time under the neutral model with parent-independent mutation, the stationary distribution of allele frequencies ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x_1, x_2, \ldots x_K$$ \end{document} ), where x_i denotes the frequency of the ith allelic type, is given by the Dirichlet distribution (Wright, 1949): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f_{N} ({\bf x} ; \theta \nu) = \frac {1} {b (\theta \nu)} x_1^{\theta \nu_1 - 1} x_2^{\theta \nu_2 - 1} \ldots x_K^{\theta \nu_K - 1} I \left\{\sum_{i = 1}^K x_i = 1 \right\}. \tag {1} \end{align*} \end{document}

Here, θ = 4Nμ is the population-scaled mutation parameter, where N is the effective population size, μ is the mutation rate per gene per generation, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\nu_i \left(\sum\nolimits_{i = 1}^K \nu_i = 1 \right)$$ \end{document} is the fraction of the mutation rate corresponding to allele i, and I{A} is the indicator function taking a value of 1 on set A and 0 otherwise. The normalizing constant b(θ ν ) is a well-known ratio of gamma functions and is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} b (\theta \nu) = \frac {\prod_{i = 1}^K \Gamma (\theta \nu_i)} {\Gamma (\theta \sum_{i = 1}^K \nu_i)}. \end{align*} \end{document}

The sampling distribution L_N(n) of a sample of size n under the multinomial sampling from distribution given in equation (1) is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} L_N ({\bf n}) = \frac {n! (\theta \nu_1)_{(n_1)} (\theta \nu_2)_{(n_2)} \cdots (\theta \nu_K)_{(n_K)}} {n_1! n_2! \cdots n_K! (\theta)_{(n)}}, \end{align*} \end{document}

where, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\bf n} (n_1, n_2, \ldots, n_K)$$ \end{document} , n_i is the number of copies of ith allele in the sample, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$n = \sum\nolimits_{i = 1}^K n_i$$ \end{document} is the sample size, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$z_{(n)} = z (z + 1) \cdots (z + n - 1)$$ \end{document} is the rising factorial. The multinomial sampling probabilities in this case are also given by b(θ ν  + n)/b(θ ν ). The Dirichlet distribution is conjugate under multinomial sampling (Kotz, 2000), meaning that if the prior probabilities on the population frequencies is Dirichlet distributed, then the posterior probabilities of the population frequencies given a sample is also Dirichlet distributed. The conditional distribution of the population frequencies given a sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\bf n} = (n_1, n_2, \cdots, n_K)$$ \end{document} from the distribution in equation (1) is Dirichlet with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f_{N \mid {\bf n}} ({\bf x} ; \theta \nu) = \frac {1} {b (\theta \nu + {\bf n})} x_1^{\theta \nu_1 + n_1 - 1} x_2^{\theta \nu_2 + n_2 - 1} \ldots x_K^{\theta \nu_K + n_K - 1} I \left\{\sum_{i = 1}^K x_i = 1 \right\} . \tag {2} \end{align*} \end{document}

2.2. K-allele model with selection

The stationary distribution of allele frequencies in the presence of selection is related to the stationary distribution under neutrality given in equation (1) in the following way. For a population with allele frequencies x, we define the mean fitness \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar{\sigma} ({\bf x})$$ \end{document} by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \bar{\sigma} ({\bf x}) = \sum_{i, j = 1}^K \sigma (A_i, A_j) x_i x_j, \end{align*} \end{document}

where σ(A_i, A_j) denotes the parameter for the strength of selection of a genotype with allelic types A_i and A_j. We let Σ be the symmetric matrix whose (i, j)th element is σ(A_i, A_j) = σ_ij. Then, the stationary distribution of allele frequencies in the K-allele model with selection (Wright, 1949) can be expressed as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f_{S} ({\bf x} ; \Sigma, \theta \nu) & = \frac {e^{{\bf x}^{\prime} \Sigma {\bf x}}} {c (\Sigma, \theta \nu)} x_1^{\theta \nu_1 - 1} x_2^{\theta \nu_2 - 1} \ldots x_K^{\theta \nu_K - 1} I \left\{ \sum_{i = 1}^K x_i = 1 \right\} \\ & = \frac {e^{\sum_{i = 1}^K \sum_{j = 1}^K \sigma_{ij} x_i x_j}} {c (\Sigma, \theta \nu)} x_1^{\theta \nu_1 - 1} x_2^{\theta \nu_2 - 1} \ldots x_K^{\theta \nu_K - 1} I \left\{\sum_{i = 1}^K x_i = 1 \right\}, & (3) \end{align*} \end{document}

where c(Σ, θ ν ) is the normalizing constant and x′ denotes the transpose of the vector x. The normalizing constant is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} c (\Sigma, \theta \nu) = \int e^{\sum_{i = 1}^K \sum_{j = 1}^K \sigma_{ij} x_i x_j} x_1^{\theta \nu_1 - 1} x_2^{\theta \nu_2 - 1} \ldots x_{K - 1}^{\theta \nu_{K - 1} - 1} x_{K}^{\theta \nu_K - 1} dx_K \cdots dx_2 dx_1, \end{align*} \end{document}

subject to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sum\nolimits_{i = 1}^K x_i = 1$$ \end{document} , where we used the integral sign for multiple integration over the allele-frequency space. An argument in Joyce (1994) uses equation (2) to show that the relationship between the sampling distribution of a sample of size n under selection, L_S(n), and under neutrality, L_N(n), is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} L_S ({\bf n}) & = L_N ({\bf n}) \frac {E_N (e^{{\bf x}^{\prime} \Sigma {\bf x}} \mid {\bf n})} {E_N (e^{{\bf x}^{\prime} \Sigma {\bf x}})} =\, L_N ({\bf n}) \frac {c (\Sigma, \theta \nu + {\bf n})} {c (\Sigma, \theta \nu)} \frac {b (\theta \nu)} {b (\theta \nu + {\bf n})} \\ & = \frac {c (\Sigma, \theta \nu + {\bf n})} {c (\Sigma, \theta \nu)}, & (4) \end{align*} \end{document}

where E_N(·) is the expectation with respect to the distribution under neutrality given by equation (1) and the last equality follows by the fact that L_N(n) = b(θ ν  + n)/b(θ ν ).

Under non-neutral models we describe, the likelihoods are intractable because the normalizing constants c(Σ, θ ν ) and c(Σ, θ ν  + n) in equation (4) are difficult to compute. Monte Carlo integration, based on generating many independent copies of x _i under the neutral model and approximating the expectation \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$E_N (e^{{\bf x}^{\prime} \Sigma {\bf x}})$$ \end{document} in equation (4) by the average \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac {1} {M} \sum_{i = 1}^M e^{{\bf x}_i^{\prime} \Sigma {\bf x}_i} \tag {5} \end{align*} \end{document}

is inefficient when selection is strong, because the allele-frequency patterns under the non-neutral model are substantially different than the allele-frequency patterns under the neutral model. Importance sampling improves the estimation of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$E_N (e^{{\bf x}^{\prime} \Sigma {\bf x}})$$ \end{document} to some extent, although obtaining many such estimates remains computationally infeasible.

3.1. Computation of the normalizing constant when Σ is diagonal

We let Σ be a diagonal matrix with elements \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sigma = (\sigma_1, \sigma_2, \cdots, \sigma_K), \sigma_i < 0, \forall i$$ \end{document} , in the diagonal and 0 elsewhere. We implement the numerical methods with σ_i <  0 because this type of selection matrix characterizes balancing selection, a selection scheme of biological importance as we discuss in our examples in Section 3.2. The normalizing constant for the stationary distribution of allele frequencies under the non-neutral model under this parameterization is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} c (\sigma, \theta \nu) = & \int_0^1 x_1^{\theta \nu_1 - 1} \int_0^{1 - x_1} x_2^{\theta \nu_2 - 1} \cdots \int_0^{1 - \sum_{i = 1}^{K - 3}x_i} x_{K - 2}^{\theta \nu_{K - 2} - 1} \\ & \int_0^{1 - \sum_{i = 1}^{K - 2}x_i} x_{K - 1}^{\theta \nu_{K - 1} - 1} \int_0^{1 - \sum_{i = 1}^{K - 1}x_i} x_K^{\theta \nu_K - 1} \\ & \times e^{- \sum_{i = 1}^{K} \sigma_i x_i^2} dx_K \cdots dx_2 dx_1. \end{align*} \end{document}

The special structure of the integrand allow computing c( σ , θ ν ) as a sequence of one-dimensional integrals. In order to demonstrate this, we first separate the inner exponential term into factors which are distributed to the respective outer integrals, and we rewrite c( σ , θ ν ) as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} c (\sigma, \theta \nu) = & \int_0^1 x_1^{\alpha_1}e^{- \sigma_1 x_1^2} \int_0^{1 - x_1} x_2^{\alpha_2}e^{- \sigma_2 x_2^2} \cdots \int_0^{1 - \sum_{i = 1}^{K - 3} x_i} x_{K - 2}^{\alpha_{K - 2}}e^{- \sigma_{K - 2} x_{K - 2}^2} \\ & \int_0^{1 - \sum_{i = 1}^{K - 2} x_i} x_{K - 1}^{\alpha_{K - 1}}e^{- \sigma_{K - 1} x_{K - 1}^2}g_K \left(1 - \sum_{i = 1}^{K - 1} x_i \right) dx_{K - 1} \cdots dx_2 dx_1, & (6) \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$g_K (y) \equiv y^{\alpha K}e^{- \sigma Ky^2}$$ \end{document} and α_i ≡ θν_i − 1. We note that the constant c( σ , θ ν ) is expressed as a (K − 1)-dimensional integral in equation (6). This is possible since the frequency of the Kth allele given the frequencies \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x_1, x_2, \ldots, x_{K - 1}$$ \end{document} , can be determined using the constraint \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sum\nolimits_{i = 1}^K x_i = 1$$ \end{document} on the allele frequencies, which allows us to write \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x_K = 1 - \sum\nolimits_{i = 1}^{K - 1} x_i$$ \end{document} . Now we consider the integral with respect to x_K−1 on the right side of equation (6) which we rewrite as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \int_0^{1 - \sum_{i = 1}^{K - 2} x_i} x_{K - 1}^{\alpha_{K - 1}} e^{- \sigma_{K - 1} x_{K - 1}^2} g_K \left(1 - \sum_{i = 1}^{K - 2} x_i - x_{K - 1} \right) dx_{K - 1}, \end{align*} \end{document}

and we note that this last expression depends on the variables \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$x_1, x_2, \cdots, x_{K - 2}$$ \end{document} only through the term \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$1 - \sum\nolimits_{i = 1}^{K - 2} x_i$$ \end{document} . We define \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$y_{K - 1} = 1 - \sum\nolimits_{i = 1}^{K - 2}x_i$$ \end{document} and denote the last integral by the function g_K−1(y_K−1). We have \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} g_{K - 1} \left(y_{K - 1} \right) = \int_0^{y_{K - 1}} x_{K - 1}^{\alpha_{K - 1}} e^{- \sigma_{K - 1} x_{K - 1}^2} g_K \left(y_{K - 1} - x_{K - 1} \right) dx_{K - 1}. \tag{7} \end{align*} \end{document}

Substituting the right hand side of equation 7 into equation (6), we get \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} c (\sigma, \theta \nu) = & \int_0^1 x_1^{\alpha_1}e^{- \sigma_1 x_1^2} \int_0^{1 - x_1} x_2^{\alpha_2}e^{- \sigma_2 x_2^2} \cdots \\ & \int_0^{1 - \sum_{i = 1}^{K - 3} x_i} x_{K - 2}^{\alpha_{K - 2}}e^{- \sigma_{K - 2} x_{K - 2}^2}g_{K - 1} \left(1 - \sum_{i = 1}^{K - 2} x_i \right) dx_{K - 2} \cdots dx_2 dx_1. \end{align*} \end{document}

Continuing in this manner, we define the successive integrals by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} g_i (y_{i}) = \int_0^{y_{i}} t^{\alpha_i} e^{- \sigma_i t^2}g_{i + 1} (y_{i} - t) dt, \tag{8} \end{align*} \end{document}

for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$y_i = 1 - \sum\nolimits_{j = 1}^{i - 1} x_j y_1 = 1$$ \end{document} , and i = K − 1, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$K - 2, \ldots, 1$$ \end{document} . The required normalizing constant c( σ , θ ν ) is given by g₁(1).

For numerical computations, we introduce the same mesh k/m, for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$k = 0, 1, \ldots, m$$ \end{document} , for each of the y_i variables. Then the g_i(y_i) values on this mesh are approximated using a trapezoidal integration rule, and g₁(1) is approximated by computing g_i(y_i) mesh approximations for i = K − 1, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$K - 2, \ldots, 2$$ \end{document} . The resulting g₁(1) approximation can be accurately and efficiently computed in O(m²) time (Genz and Joyce, 2003), except in some cases, where the α_i values produce end-point singularities for the g_i(y_i) integrals. In those cases, we use standard numerical methods for “subtracting out the singularities” (David and Rabinowitz, 1984), to eliminate the singular integrand behaviors and produce rapidly converging approximations to g₁(1).

3.2. Generating samples when Σ is a diagonal matrix

We define the cumulative distribution functions F_i(·; z) with parameter z, for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$i = 1, 2, \cdots, K$$ \end{document} as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} F_i (y;z) = \begin{cases} 0 &\quad y \leq 0 \\ \\ \frac {\int_0^y t^{\alpha_i} \exp (- \sigma_i t^2) g_{i + 1} (z - t) dt} {g_i (z)} &\quad 0 \leq y \leq z \\ \\ 1 & \quad y > z \end{cases} \end{align*} \end{document}

where g_i(y) is defined by equation (8). We note that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$P (X_i \leq y \mid X_{i - 1}, \cdots, X_1) = F_i (y;1 - X_1 - \cdots - X_{i - 1})$$ \end{document} . Since F_i(·;z) is a strictly increasing function for 0 < y < z, it has a well-defined inverse that we denote by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F_i^{- 1} (\cdot ;z)$$ \end{document} . Given these definitions, we now generate a sequence of random variables \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\bf X} = (X_1, X_2, \cdots, X_K)$$ \end{document} with joint probability density function given by equation (3) when the matrix of selection parameters is diagonal. The method, which is an application of inversion sampling (i.e., “the CDF method”) is as follows. We let U₁ be a uniformly distributed random number on [0, 1]. That is, U₁ ∼UNIF[0, 1]. We define \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X_1 = F_1^{- 1} (U_1, 1)$$ \end{document} . Now we define iteratively

1. Generate U_i ∼ UNIF \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$[0, 1 - X_1 - X_2 - \cdots - X_{i - 1}]$$ \end{document} ,

If the normalizing constant is calculated using the method described in Section 3.1, then approximations for the g_i(y) functions have also been computed, and we can generate independent random vectors of allele frequencies using the above algorithm. Each vector X is a realization of the (random) population allele frequencies under selection. To simulate a sample of size n with allele counts given by n and distribution given by equation (3), it is sufficient to draw a multinomial sample conditional on the given population frequencies X.

4.1. Bayesian inference under non-neutral models

By Bayes' Theorem, under the non-neutral model with diagonal matrix of selection parameters, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Sigma = {\rm diag} (\sigma, \sigma, \ldots, \sigma)$$ \end{document} , the joint posterior distribution of selection and mutation parameters can be written as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \pi (\sigma, \theta \nu \mid {\bf x}) = \frac {f ({\bf x} \mid \sigma, \theta \nu) \pi (\sigma, \theta \nu)} {\int \int f ({\bf x} \mid \sigma, \theta \nu) \pi (\sigma, \theta \nu) d \sigma d \theta \nu}, \tag {9} \end{align*} \end{document}

where f (x|σ, θν) is the likelihood and π(σ, θν) is the joint prior. For simplicity, we assume that mutation is non-preferential among the K alleles and we fix \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\nu = (1 / K, 1 / K, \ldots 1 / K)$$ \end{document} although use of other schemes does not add difficulty to the analysis. Further, we assume the prior independence of selection and mutation parameters and write π(σ, θν) = π(σ)π(θν). When sampling the joint posterior distribution in equation (9), the normalizing constant in the denominator need not be computed, since most stochastic samplers require the posterior distribution to be known up to a constant free of parameters. However, this is not the normalizing constant computed by our methods in Section 3.1. In equation (9), the parameter-dependent normalizing constant computed by our methods is hidden in the likelihood f (x|σ, θ ν ). We let, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\bf x}_o = ({x_1}_o, {x_2}_o, \ldots, {x_K}_o)$$ \end{document} and n _o to denote the observed population frequencies and sample frequencies respectively. If the sample frequencies are taken as proxy to population frequencies and the population distribution is used, the likelihood f (x _o |σ, θ ν ) can be explicitly written as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f ({\bf x}_o \mid \sigma, \theta \nu) = \frac {e^{- \sigma \sum_{i = 1}^K {{x_i}_o}^2}} {c (\sigma, \theta \nu)} {{x_1}_o}^{\theta / K - 1} {{x_2}_o}^{\theta / K - 1} \ldots {{x_K}_o}^{\theta / K - 1}, \tag {10} \end{align*} \end{document}

and if the sampling distribution is used, it can be written as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f ({\bf n}_o \mid \sigma, \theta \nu) = L_N ({\bf n}_o) \frac {c (\sigma, \theta \nu + {\bf n}_o)} {c (\sigma, \theta \nu)} \frac {b (\theta \nu)} {b (\theta \nu + {\bf n}_o)} = \frac {c (\sigma, \theta \nu + {\bf n}_o)} {c (\sigma, \theta \nu)} . \tag {11} \end{align*} \end{document}

We now define the statistics \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} F = \sum_{i = 1}^{K}x_{i}^2, \quad G = - \sum_{i = 1}^{K} \log x_i \tag{12} \end{align*} \end{document}

and express the density of the non-neutral model with diagonal matrix of selection parameters as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f ({\bf x} \mid \sigma, \theta \nu) & = \frac {e^{- \sigma \sum_{i = 1}^K {x_i}^2}} {c (\sigma, \theta \nu)} {x_1}^{\theta / K - 1} {x_2}^{\theta / K - 1} \ldots {x_K}^{\theta / K - 1} \\ & = \frac {e^{- \sigma F}} {c (\sigma, \theta \nu)} e^{- \log \prod_{i = 1}^K x_i^{\theta / K - 1}} \\ & = \frac {e^{- \sigma F + (\theta / K - 1) G}} {c (\sigma, \theta \nu)} . \end{align*} \end{document}

The right-hand side of the last equality depends on the allele frequencies x only through the functions of F and G, and the Fisher-Neyman factorization theorem implies that F and G are jointly sufficient for the parameters σ and θ. To sample the joint posterior distribution of the selection and mutation parameters we use an approximate Bayesian computation algorithm with the likelihood based on equation (10):

Algorithm–ABC (with rejection):

1. Simulate a value θ* from the prior distribution π(θ) and a value σ* from the prior distribution π(σ)

where F_o and G_o are the statistics computed from the observed data, || · || is the Euclidean distance, and δ is a suitable tolerance parameter. The posterior sample obtained by Algorithm–ABC is from the correct posterior distribution as δ approaches zero. We also note that, since the summary statistics are sufficient for the parameters under this non-neutral model, there is no approximation associated with reducing the full-data likelihood to a likelihood based on the summary statistics. An example for the performance of inference on the joint distribution of the parameters is given in Figure 2, using the likelihood based on the population frequencies (i.e., equation (10)). We simulated the “true” population frequencies used as test data, from a model with 10 replications, from a locus with K = 10 alleles, under a selection parameter σ = 10 and a mutation parameter θ = 2. We generated 10⁵ data sets under the model with 10⁵ parameter values (σ, θ), drawn from their prior distribution (Step 2 in Algorithm-ABC above), and 1% of the parameters which generated the summary statistics (F*, G*) closest to the observed summary statistics (F_o, G_o) in the Euclidean sense, are taken as a sample from the posterior distribution.

OPEN IN VIEWER

FIG. 2.

Posterior distributions (top) and boxplots (bottom) of selection and mutation parameters. The “true” population frequencies are simulated under a K = 10 allele model with 10 replicates (e.g., 10 loci) and a diagonal selection matrix where all elements are equal (i.e., symmetric balancing selection model). The “true” selection parameter is σ = 10 and the mutation parameter is θ = 2. Top left and right histograms are posterior samples of selection and mutation parameters respectively, both obtained by approximate Bayesian computation as explained in Algorithm-ABC in the text. Bottom left and right series of box plots are from posterior samples for 30 independent runs of “true” population frequencies, obtained by the same procedure as above.

4.2. The distribution of the r-largest allele frequencies

In models considered thus far, we assumed that all the allelic types present in the population are observed in the sample, treating K as fixed and known. In many practical problems, this assumption might not be reasonable. For example, in highly-polymorphic immune system genes in humans such as Human Leukocyte Antigen (HLA), some of the low-frequency alleles in the population will likely be missing in the sample. The stationary distribution of allele frequencies for a diagonal matrix of selection parameters for various cases of missing data and when K is unkown is studied in detail by Watterson (1977). Here, we consider the distribution of the r-largest allele frequencies in a population with K total alleles (r < K), to investigate the effect of missing low-frequency alleles in the sample when performing inference on the strength of balancing selection.

The distribution of the r-largest allele frequencies \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\bf x} = (x_1, x_2, \ldots, x_r), x_1 > x_2 > \cdots > x_r > \cdots > x_K$$ \end{document} , in a population of K total types can be written as (Watterson, 1977), \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} & f_{S}^{(K, r)} ({\bf x} ; \sigma, \theta \nu) = \frac {K!} {(K - r) !} \times e^{- \sigma \sum_{i = 1}^r x^2_i} x_1^{\theta \nu_1 - 1} x_2^{\theta \nu_2 - 1} \ldots x_r^{\theta \nu_r - 1} \times \frac {c^{(K, r)} (\sigma, \theta \nu)} {c (\sigma, \theta \nu)}, & (13) \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} c^{(K, r)} (\sigma, \theta \nu) = & \int_0^{1 - \sum_{i = 1}^{r}x_i}x_{r + 1}^{\theta \nu_{r + 1} - 1} \int_0^{1 - \sum_{i = 1}^{r + 1}x_i} x_{r + 2}^{\theta \nu_{r + 2} - 1} \cdots & (14) \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \int_0^{1 - \sum_{i = 1}^{K - 3}x_i} x_{K - 2}^{\theta \nu_{K - 2} + n_{K - 2} - 1} \int_0^{1 - \sum_{i = 1}^{K - 2}x_i}x_{K - 1}^{\theta \nu_{K - 1} + n_{K - 1} - 1} & (15) \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \int_0^{1 - \sum_{i = 1}^{K - 1}x_i}x_K^{\theta \nu_K + n_K - 1} \times e^{- \sigma \sum_{i = r + 1}^{K} x_i^2} dx_K \cdots dx_2 dx_{r + 1}. & (16) \end{align*} \end{document}

The distribution of the r-largest allele frequencies has a normalizing constant that is a ratio of two integrals. The first integral, which appears in the denominator of equation (13) is the normalizing constant of the full model with K-alleles. The second integral, which appears in the numerator of equation (13) is a (K − r)-dimensional integral and is a sum over the frequencies of unobserved alleles. The methods of Section 3.1 can be modified to compute the integral in the numerator and hence the normalizing constant in equation (13). Application of our numerical methods, once for the integral in the numerator and once for the integral in the denominator, allows us to evaluate likelihoods in r-largest allele frequencies based on equation (13).

To demonstrate the use of the distribution of r-largest allele frequencies, we consider the following idealized example, where sample frequencies are reasonable proxy for population frequencies. We simulate a vector of population frequencies under the model with K = 10, σ = 100, and θ = 5 to represent the “true” state of the population at a locus. We now consider three analyses. First, we assume that r = 7 largest-frequency alleles are observed in a sample with the true frequencies as in the population and analyze these frequencies using a K = 7 allele model. This is a case where three alleles are missing in the data set, since the true K is 10, and we are unaware that we are missing the three smallest-frequency alleles, using a model with K = 7 to analyze the data. Second, we assume that r = 7 largest-frequency alleles are observed in a sample with the true frequencies as in the population and analyze these data using equation (13) where K = 10. This is a case where again three alleles are missing in the data set, since the true K is 10, but this time we know that we are missing the three smallest-frequency alleles and use a model with K = 10 to analyze the data. Third, we assume all K = 10 alleles are observed in a sample with the true frequencies as in the population and analyze these data with a model with K = 10 alleles. The posterior distribution obtained from the third analysis represents the ideal situation, in which there are no missing alleles in the sample and the sample frequencies are proxy for the population frequencies. Therefore, the posterior distribution obtained from the third procedure is the gold standard for the posterior distributions obtained from the first and second analyses. We obtain posterior samples for all three scenarios using an MCMC approach with 100,000 iterations and 10,000 burn in steps.

The posterior distribution of the selection and mutation parameters obtained using the first analysis results in biased estimates for selection and mutation parameters (Fig. 3). Low-frequency alleles carry valuable information on mutation and consequently the posterior distribution of the mutation parameter is substantially affected. Using the distribution of the r-largest allele frequencies results in improved estimation, as the model compensates the missing alleles. However, a good estimate of K is required in this case, close to its true value. These results suggest that in a K-allele balancing selection model, it might be important to capture the low-frequency alleles in the sample to obtain good estimates of the strength of selection, especially when K is unknown and no good estimates for K exist.

OPEN IN VIEWER

FIG. 3.

Kernel density estimates of posterior distributions of selection (top) and mutation (bottom) parameters obtained analyzing allele-frequencies under three non-neutral models. The models differ in their assumptions on the number of allelic types in the population, K. The “true” population has K = 10 alleles but only 7 largest-frequency alleles are observed in the sample. Green plots are obtained by analyzing these 7 frequencies under the assumption that r = 7 largest-frequency alleles out of K = 10 total alleles are observed in the data. Red plots are obtained by analyzing the 7 frequencies under the assumption that these are all the existing allelic types in the population. Hence, the second model corresponds to using a K = 7 allele model, whereas the frequencies are actually generated under a model with K = 10. Using the second model, the mutation parameter is overestimated (bottom, red plot) because of missing three small-frequency alleles. Consequently, estimates for the selection parameter are not accurate. Blue plots are obtained by assuming ideal conditions where all 10 allelic types are observed and the sample frequencies are proxy for the population frequencies. Hence, the blue plots are obtained under the correct model when there is no missing data and inference in this case is the gold standard for the other two models.