Model-Based Clustering with Ancestry Proportions

Why Model Ancestry as a Mixture?

Principal components summarize major directions of genetic variation, but they do not directly specify a probability model for how an individual's genotypes relate to ancestral allele-frequency profiles. A model-based mixture gives a more explicit description: at each locus, an individual's allele copies are associated probabilistically with \(K\) components, and each component has its own allele frequencies. The individual's ancestry proportions determine how strongly each component contributes.

The model is useful because it links two different quantities: population-level allele-frequency profiles and individual-level mixture proportions. These are estimated together from genotype data by maximizing a likelihood, often using the expectation-maximization algorithm. The resulting components are statistical parts of a model; their interpretation depends on the data and assumptions, not on an automatic identification of discrete historical populations.

Definition: In an admixture model with \(K\) components, individual \(i\) has ancestry proportions \(q_{i1},\ldots,q_{iK}\), where \(q_{ik}\geq 0\) and \(\sum_{k=1}^{K}q_{ik}=1\). At locus \(\ell\), component \(k\) has allele-frequency profile \(f_{k\ell a}\), the frequency of allele \(a\). Each allele copy is assigned probabilistically to a component according to the individual's proportions, and its allele is then drawn according to that component's frequency profile.

Model and Genotype Likelihood

Begin with diploid individuals and biallelic loci. Call one allele at locus \(\ell\) the designated allele, and let \(f_{k\ell}\) be its frequency in component \(k\). The other allele has frequency \(1-f_{k\ell}\). Let \(x_{i\ell}\in\{0,1,2\}\) be the number of designated-allele copies observed for individual \(i\). For any one allele copy from that individual, the probability of observing the designated allele is the mixture

p_{i\ell}=\sum_{k=1}^{K}q_{ik}f_{k\ell}.

Under the model, the two allele copies at a locus are conditionally independent given the parameters, and copies at different loci are treated as independent. Thus the genotype count has a binomial probability with two trials and success probability \(p_{i\ell}\). The likelihood across \(n\) individuals and \(L\) loci is

L(\mathbf q,\mathbf f) =\prod_{i=1}^{n}\prod_{\ell=1}^{L} \binom{2}{x_{i\ell}}p_{i\ell}^{x_{i\ell}} (1-p_{i\ell})^{2-x_{i\ell}}.

Here \(\mathbf q\) denotes all ancestry proportions and \(\mathbf f\) all component-specific frequencies. The likelihood is the probability of the observed genotypes as a function of these unknown parameters. A maximum-likelihood estimate is a set of parameter values that maximizes it, subject to the constraints that each ancestry vector sums to one and each allele frequency lies between zero and one. The binomial coefficient does not depend on the parameters, so it can be omitted when maximizing the likelihood.

This formulation assumes that the loci included in the likelihood can be treated as independent. It also assumes that the specified \(K\) profiles and individual mixture proportions adequately describe the allele-frequency patterns in the data. These are modeling assumptions, not facts established by computing a likelihood. For a locus with more than two alleles, replace \(f_{k\ell}\) with a frequency for each allele, with the frequencies summing to one within each component and locus.

Worked Example: A Genotype Probability Under a Mixture

Suppose a hypothetical individual has \(q_{i1}=0.60\) and \(q_{i2}=0.40\). At one biallelic locus, the designated-allele frequencies are \(f_{1\ell}=0.80\) and \(f_{2\ell}=0.30\). Find the probability of the designated allele in a copy, and the probability of observing a heterozygous genotype.

The copy-level probability is the ancestry-weighted frequency:

p_{i\ell}=(0.60)(0.80)+(0.40)(0.30)=0.48+0.12=0.60.

A heterozygote has \(x_{i\ell}=1\). Its probability is therefore

P(X_{i\ell}=1)=\binom{2}{1}(0.60)^1(1-0.60)^1 =2(0.60)(0.40)=0.48.

The value \(0.60\) is a probability for one allele copy, whereas \(0.48\) is the probability of the diploid heterozygous genotype. They are not interchangeable. The calculation uses the model's conditional-independence assumption for the two copies.

Latent Ancestry Assignments

The observed-data likelihood sums over the possible component assignments of each allele copy. The EM algorithm makes these hidden assignments explicit as latent variables. For an observed allele copy \(c\) from individual \(i\) at locus \(\ell\), let \(a_{i\ell c}\) denote its allele. If it is the designated allele, the component-specific frequency is \(f_{k\ell}\); if it is the other allele, the frequency is \(1-f_{k\ell}\). The conditional probability that this copy came from component \(k\), given the observed allele and current parameter values, is called its responsibility:

r_{i\ell ck} = \frac{q_{ik}f_{k\ell a_{i\ell c}}} {\sum_{h=1}^{K}q_{ih}f_{h\ell a_{i\ell c}}}.

Here \(f_{k\ell a}\) means the frequency of the observed allele \(a\) in component \(k\). The numerator is the joint probability of choosing component \(k\) and then observing that allele; the denominator sums this quantity over all components. For each allele copy, the responsibilities are nonnegative and sum to one across \(k\). They are conditional allocation probabilities under the current fitted parameter values, not observed ancestry labels.

Worked Example: Computing Responsibilities at Two Loci

Consider a hypothetical individual with \(q_{i1}=q_{i2}=0.50\), heterozygous at each of two biallelic loci. At locus 1, the designated-allele frequencies are \(0.20\) and \(0.70\) in components 1 and 2. At locus 2, they are \(0.30\) and \(0.80\). Compute the responsibilities for each allele type.

At locus 1, for a designated-allele copy, the component weights are \((0.50)(0.20)=0.10\) and \((0.50)(0.70)=0.35\). Their sum is \(0.45\), so the responsibilities are

r_{\text{designated},1}=\frac{0.10}{0.45}=0.2222, \qquad r_{\text{designated},2}=\frac{0.35}{0.45}=0.7778.

For the other allele at locus 1, the frequencies are \(0.80\) and \(0.30\). The weights are \(0.40\) and \(0.15\), summing to \(0.55\). Thus the responsibilities are \(0.40/0.55=0.7273\) for component 1 and \(0.15/0.55=0.2727\) for component 2.

At locus 2, a designated copy has weights \((0.50)(0.30)=0.15\) and \((0.50)(0.80)=0.40\), summing to \(0.55\). Its responsibilities are \(0.15/0.55=0.2727\) for component 1 and \(0.40/0.55=0.7273\) for component 2. For the other allele, the frequencies are \(0.70\) and \(0.20\), giving weights \(0.35\) and \(0.10\), with sum \(0.45\). The responsibilities are \(0.35/0.45=0.7778\) and \(0.10/0.45=0.2222\), respectively.

As a check, the two copy probabilities at locus 2 are \(0.55\) for the designated allele and \(0.45\) for the other allele. Their product, multiplied by two possible copy orders, gives heterozygote probability \(2(0.55)(0.45)=0.495\). The weights used in each responsibility calculation sum to the corresponding copy probability.

Likelihood-Based Estimation With EM

Directly maximizing the observed-data likelihood can be difficult because each copy's component of origin is unknown. EM alternates between estimating the conditional distribution of those assignments and maximizing the expected complete-data log-likelihood. In the E-step, calculate responsibilities using the current parameters. In the M-step, treat these responsibilities as fractional expected counts and update the ancestry proportions and allele frequencies.

Initialize parameters.
Choose valid starting proportions \(q_{ik}\) and frequencies \(f_{k\ell}\), respecting their sum and range constraints.

Compute responsibilities.
For every observed allele copy, calculate \(r_{i\ell ck}\) from the current proportions and frequencies.

Update proportions and frequencies.
Use expected ancestry counts for each individual to update \(q_{ik}\), and expected allele and total-copy counts within each component and locus to update \(f_{k\ell}\).

Iterate.
Recalculate responsibilities and update parameters until the likelihood or parameter estimates change by less than a chosen tolerance.

For diploid data, write \(r_{i\ell ck}\) for the responsibility of component \(k\) for copy \(c\). The proportion update is the expected number of copies assigned to component \(k\) in individual \(i\), divided by the \(2L\) copies observed for that individual. The frequency update is the expected number of designated copies assigned to component \(k\) at locus \(\ell\), divided by the expected total number of copies assigned there:

q_{ik}^{\text{new}} =\frac{\sum_{\ell=1}^{L}\sum_{c=1}^{2}r_{i\ell ck}}{2L}, \qquad f_{k\ell}^{\text{new}} =\frac{\sum_{i=1}^{n}\sum_{c=1}^{2} r_{i\ell ck}\,\mathbb{1}(a_{i\ell c}=\text{designated})} {\sum_{i=1}^{n}\sum_{c=1}^{2}r_{i\ell ck}}.

The indicator \(\mathbb{1}(\cdot)\) equals one when its condition holds and zero otherwise. The proportion updates sum to one because, for every copy, the responsibilities over components sum to one. The frequency update is a weighted allele fraction: each copy contributes to component \(k\) in proportion to its responsibility.

Worked Example: One EM Update for Two Individuals

Consider a hypothetical diploid sample with two individuals and one biallelic locus. Individual A has genotype \(x=2\), and individual B has \(x=0\). Start with \(q_{A1}=q_{A2}=q_{B1}=q_{B2}=0.50\), and frequencies \(f_{1}=0.80\), \(f_{2}=0.20\). Perform one EM update.

For each designated copy in A, the weights are \((0.50)(0.80)=0.40\) and \((0.50)(0.20)=0.10\). The responsibilities are \(0.40/0.50=0.80\) for component 1 and \(0.10/0.50=0.20\) for component 2. A has two such copies, so its expected ancestry counts are \(1.60\) and \(0.40\), giving \(q_{A1}^{\text{new}}=1.60/2=0.80\) and \(q_{A2}^{\text{new}}=0.40/2=0.20\).

B has two copies of the other allele. Its component weights are \((0.50)(1-0.80)=0.10\) and \((0.50)(1-0.20)=0.40\). Thus each copy has responsibilities \(0.20\) and \(0.80\). The expected counts are \(0.40\) and \(1.60\), so \(q_{B1}^{\text{new}}=0.20\) and \(q_{B2}^{\text{new}}=0.80\).

For the frequency update, the expected number of designated copies assigned to component 1 is \(1.60\), all from A; for component 2 it is \(0.40\). The expected total copy counts at the locus are \(1.60+0.40=2.00\) for component 1 and \(0.40+1.60=2.00\) for component 2. Therefore

f_{1}^{\text{new}}=\frac{1.60}{2.00}=0.80, \qquad f_{2}^{\text{new}}=\frac{0.40}{2.00}=0.20.

The frequencies remain unchanged in this iteration. The observed-data likelihood, omitting no parameter-dependent terms, increases from \(0.50^2(0.50)^2=0.0625\) initially to \((0.68)^2(0.32)^2=0.0473\) after the update? This comparison must use genotype probabilities: initially each homozygote has probability \(0.50^2=0.25\), so the joint likelihood is \(0.25^2=0.0625\). After the update, A's copy probability is \(0.80(0.80)+0.20(0.20)=0.68\), and B's designated-copy probability is \(0.20(0.80)+0.80(0.20)=0.32\); its genotype probability is \((1-0.32)^2=0.68^2\). Thus the joint likelihood is \(0.68^2\times0.68^2=0.2138\), greater than \(0.0625\).

Why the EM Updates Make Sense

Each EM iteration does not simply assign every allele copy to the most likely component. Instead, it allocates each copy fractionally across components according to its responsibility. This retains uncertainty about component of origin and gives the M-step expected counts appropriate for maximum-likelihood updates.

The nondecreasing-likelihood property follows from the EM construction. Let \(Z\) denote all latent copy assignments and let \(r(Z)\) be their conditional distribution given the observed genotypes under the current parameters \(\theta\). For proposed parameters \(\theta'\), Jensen's inequality gives

\log P(X\mid\theta') =\log\sum_Z r(Z)\frac{P(X,Z\mid\theta')}{r(Z)} \geq \sum_Z r(Z)\log\frac{P(X,Z\mid\theta')}{r(Z)}.

The right-hand side is the expected complete-data log-likelihood under the current assignments plus their entropy. At \(\theta'=\theta\), this bound equals \(\log P(X\mid\theta)\), because \(r(Z)\) is the posterior assignment distribution under \(\theta\). The M-step chooses \(\theta'\) to maximize the expected complete-data log-likelihood, so the bound cannot decrease. Consequently, the observed-data log-likelihood cannot decrease in an exact EM iteration. This does not guarantee that the algorithm finds the global maximum: results can depend on starting values, and convergence to a local maximum or a stationary point is possible.

Interpretation and Common Pitfalls

An estimated \(q_{ik}\) is the fitted fraction of individual \(i\)'s allele-copy ancestry attributed to component \(k\) under the specified model. It is a model-based summary, not a direct observation of genealogical ancestry or a statement that an individual belongs to a single discrete group. Likewise, an estimated \(f_{k\ell}\) is a component-specific allele frequency within this fitted representation.

Confusing proportions with hard assignments. A value such as \(q_{i1}=0.60\) does not mean that exactly 60% of the individual's observed loci have been assigned to component 1. It is a fitted mixture proportion.
Using the wrong allele frequency in a responsibility. A copy of the other allele uses \(1-f_{k\ell}\), not \(f_{k\ell}\). Normalize the component weights using the frequencies of the allele actually observed.
Forgetting the denominator in the frequency update. The denominator is the expected total number of copies assigned to that component at that locus, not the number of individuals in the sample.
Reading the likelihood as proof that the model is true. Maximizing likelihood finds the best-fitting parameter values within the specified model. It does not establish that the independence assumptions or component interpretation are correct.

Key takeaway: The admixture likelihood combines an individual's ancestry proportions with component-specific allele frequencies to give genotype probabilities. EM estimates the parameters by repeatedly allocating expected ancestry counts and updating the proportions and frequencies; its output must be interpreted in the context of the model assumptions.

Check Your Understanding

Try these questions using the model and update equations above.

With two components, \(q_{i1}=0.25\), \(q_{i2}=0.75\), and designated-allele frequencies \(0.60\) and \(0.20\), what is the designated-allele probability for one copy?
Why does the genotype probability for a heterozygote include a factor of two?
For a copy of the other allele, which frequency enters the responsibility numerator for component \(k\)?
In the M-step, how is the updated individual ancestry proportion calculated from responsibilities across loci?
Does a nondecreasing likelihood guarantee that EM has found the global maximum? Explain briefly.