Likelihood Uncertainty in Haplotype-Frequency Estimates

Why Uncertainty Requires the Likelihood Shape

An algorithm can return a single estimated frequency even when the data do not pin that frequency down precisely. This matters for unphased two-locus genotypes: a double heterozygote is compatible with two different haplotype pairings, and the observed genotype does not reveal which pairing occurred. The likelihood must account for both possibilities. Its shape near and away from its maximum indicates how strongly the observed genotypes distinguish among candidate haplotype frequencies.

The focus here is the observed-data likelihood and its profile, not just the point at which an estimation algorithm stops. A sharply peaked profile suggests that nearby frequency values fit the data less well; a broad or irregular profile signals weaker resolution. That assessment depends on the genetic sampling model as well as the observed genotype counts.

Model and Observed-Data Likelihood

Consider two biallelic loci, with alleles \(A/a\) at the first locus and \(B/b\) at the second. Denote the haplotype frequencies by \(f_{AB}\), \(f_{Ab}\), \(f_{aB}\), and \(f_{ab}\). They are nonnegative and sum to one. A diploid individual carries an unordered pair of haplotypes; the pair determines the two-locus genotype, but an unphased record may not identify the pair.

Model assumptions: Each sampled individual receives two haplotypes drawn independently from a population with fixed haplotype frequencies, as under random mating with Hardy–Weinberg proportions. Sampled individuals are independent. Genotypes are recorded without error, and the sampled population is represented by one set of haplotype frequencies. The likelihood below is conditional on these assumptions; departures such as population mixture or genotyping error can change its shape and interpretation.

For two distinct haplotypes \(h_i\) and \(h_j\), the probability of the unordered diplotype is \(2f_i f_j\), because either haplotype can be inherited from either parent. For a pair of identical haplotypes \(h_i/h_i\), the probability is \(f_i^2\). An observed unphased genotype \(g\) can be compatible with one or more such pairs. Its probability is the sum of the probabilities of all compatible pairs.

Definition: Let \(C(g)\) be the set of unordered haplotype pairs compatible with genotype \(g\). For a pair \(h_i/h_j\), define its probability weight as \(w_{ij}=f_i^2\) if \(i=j\), and \(w_{ij}=2f_i f_j\) if \(i\ne j\). Then the observed-genotype probability is \(P(g\mid\mathbf f)=\sum_{(i,j)\in C(g)}w_{ij}\), where \(\mathbf f\) is the vector of four haplotype frequencies.

If \(n_g\) individuals have observed genotype \(g\), the genotype-count likelihood is multinomial. Up to the multinomial coefficient, its log likelihood is a sum of count-weighted log probabilities. The coefficient does not depend on the frequencies, so it does not change the maximizing estimates or likelihood ratios.

L(\mathbf f)=\frac{n!}{\prod_g n_g!}\prod_g P(g\mid\mathbf f)^{n_g}, \qquad \ell(\mathbf f)=\sum_g n_g\log P(g\mid\mathbf f)+\text{constant}, \qquad \sum_h f_h=1,\quad f_h\ge 0.

Worked Example: Building a Likelihood Contribution

Suppose a hypothetical sample of six individuals contains three double heterozygotes \(Aa,Bb\), two genotypes \(AA,Bb\), and one genotype \(aa,bb\). Evaluate the likelihood at the candidate frequencies \(f_{AB}=0.40\), \(f_{Ab}=0.10\), \(f_{aB}=0.20\), and \(f_{ab}=0.30\). These values are nonnegative and sum to 1.00, as required.

For \(Aa,Bb\), the compatible pairs are \(AB/ab\) and \(Ab/aB\). Thus its probability is \(2(0.40)(0.30)+2(0.10)(0.20)=0.24+0.04=0.28\). For \(AA,Bb\), the pair must be \(AB/Ab\), giving probability \(2(0.40)(0.10)=0.08\). For \(aa,bb\), both haplotypes must be \(ab\), giving probability \((0.30)^2=0.09\).

The multinomial coefficient is \(6!/(3!2!1!)=60\). Therefore,

L(\mathbf f)=60(0.28)^3(0.08)^2(0.09) =60(0.021952)(0.0064)(0.09) \approx 0.0007587.

This is the likelihood of the stated genotype-count pattern at these candidate frequencies. It is not, by itself, a probability that the frequencies are correct: to assess relative support, compare likelihoods across candidate frequencies or examine a profile.

How Phase Ambiguity Changes Curvature

At a double heterozygote, the observed probability includes both phase configurations:

P(Aa,Bb\mid\mathbf f)=2f_{AB}f_{ab}+2f_{Ab}f_{aB}.

The two terms are added because either configuration produces the same unphased genotype. Treating one phase as known would replace this sum with a single term and would generally use the wrong likelihood. The sum can also make different frequency combinations fit similarly well: an increase in one compatible-pair probability can partly compensate for a decrease in the other.

The curvature of the log likelihood describes how quickly its slope changes as frequencies move. In a one-parameter setting, a large negative second derivative near the maximum indicates a relatively sharp peak; a small-magnitude second derivative indicates a flatter peak. With multiple frequency parameters, curvature is described by the Hessian matrix of second derivatives. Weak curvature in some direction means the data distinguish poorly among frequency combinations in that direction. A phase-ambiguous genotype may contribute information, but its contribution is through a sum of alternatives rather than a directly observed phase.

Worked Example: Comparing the Two Phases

For the frequencies in the previous example, compute how the model divides the probability of a double heterozygote between its two compatible configurations. The coupling configuration \(AB/ab\) has probability \(2(0.40)(0.30)=0.24\). The repulsion configuration \(Ab/aB\) has probability \(2(0.10)(0.20)=0.04\). Their sum, \(0.28\), is the probability of the observed unphased genotype.

Conditional on observing this genotype and these specified frequencies, the model assigns the coupling configuration the fraction \(0.24/0.28\approx0.8571\) of the probability and the repulsion configuration \(0.04/0.28\approx0.1429\). These are model-based conditional phase probabilities at fixed frequencies, not evidence that the phase was directly observed.

For five independent double heterozygotes at these same frequencies, the probability of the observed unphased pattern is \((0.28)^5\approx0.001721\). If one incorrectly treated all five as known coupling phase, the corresponding probability would instead be \((0.24)^5\approx0.0007963\). The latter calculation discards compatible configurations and is not the observed-data likelihood.

Profile Likelihood for a Frequency

A likelihood profile focuses on one frequency parameter while allowing the other frequencies to adjust. Write the parameter of interest as \(\psi\), and let \(\boldsymbol\eta\) denote nuisance parameters, such as the remaining independent haplotype frequencies. The profile log likelihood at \(\psi\) is the largest log likelihood possible at that value of \(\psi\), after maximizing over the nuisance parameters subject to the frequency constraints.

Definition: The profile log likelihood is \(\ell_p(\psi)=\max_{\boldsymbol\eta}\ell(\psi,\boldsymbol\eta)\), where the maximization is over values satisfying the model constraints. A likelihood-ratio profile compares each value with the maximum over all parameters: \(2\{\ell_p(\hat\psi)-\ell_p(\psi)\}\). Large values mean that the candidate value fits worse relative to the best-fitting value; small values mean that it remains comparatively compatible with the observed data.

The profile is useful because it accounts for the fact that nuisance frequencies can change when \(\psi\) changes. Simply holding all other frequencies fixed can exaggerate or understate uncertainty. A broad profile indicates that a range of values has similar fit. Asymmetry can show that uncertainty differs on the two sides of the estimate; a boundary maximum or multiple peaks calls for particular care. A numerical optimizer’s single answer does not reveal these features.

Under regular conditions, a likelihood-ratio statistic for one parameter is often compared with a chi-square distribution with one degree of freedom. This is an approximation, not a universal rule: boundary estimates, weak identification, or model violations can make the usual reference unreliable. The profile itself remains a useful description of relative fit, but converting it into a formal interval or error rate requires appropriate inferential conditions.

Worked Example: Reading Curvature in a One-Parameter Model

Consider a hypothetical dataset of 90 individuals. Suppose the two allele frequencies are fixed at 0.50 at each locus, and let \(D\) vary the haplotype frequencies as \(f_{AB}=f_{ab}=0.25+D\) and \(f_{Ab}=f_{aB}=0.25-D\), with \(-0.25\le D\le0.25\). This restriction gives a transparent one-parameter illustration; it is not a claim that marginal frequencies are always known in practice.

The observed counts are 30 double heterozygotes; 10 in each of the four classes that are heterozygous at one locus and homozygous at the other; and 5 in each of the four classes homozygous at both loci. Under the stated frequency restriction, the probabilities are \(0.25+4D^2\) for a double heterozygote, \(0.125-2D^2\) for each single-heterozygote class, and either \((0.25+D)^2\) or \((0.25-D)^2\) for the four double-homozygote classes. Thus, omitting the count-only multinomial coefficient, the log likelihood is

\ell(D)=30\log(0.25+4D^2)+40\log(0.125-2D^2)+10\log\{(0.25+D)^2\}+10\log\{(0.25-D)^2\}.

Set \(t=D^2\). The expression becomes \(30\log(0.25+4t)+40\log(0.125-2t)+20\log(0.0625-t)\). Its derivative is

\frac{d\ell}{dt}=\frac{120}{0.25+4t}-\frac{80}{0.125-2t}-\frac{20}{0.0625-t}.

At \(t=0\), this derivative is \(120/0.25-80/0.125-20/0.0625=480-640-320=-480\). Moreover, the second derivative is \(-480/(0.25+4t)^2-160/(0.125-2t)^2-20/(0.0625-t)^2\), which is negative wherever these probabilities are positive. The log likelihood therefore decreases as \(t\) increases from zero, so its maximum occurs at \(D=0\). Because the likelihood depends on \(D^2\), it is symmetric: the candidate values \(D\) and \(-D\) have equal likelihood in this example.

To see how quickly fit declines, compare \(D=0.05\) with the maximum. At \(D=0.05\), the three distinct probabilities in the expression are \(0.26\), \(0.12\), and \(0.06\), compared with \(0.25\), \(0.125\), and \(0.0625\) at zero. Therefore,

\ell(0.05)-\ell(0) =30\log(0.26/0.25)+60\log(0.12/0.125) \approx 1.1766-2.4493=-1.2727.

The likelihood-ratio comparison is \(2\{\ell(0)-\ell(0.05)\}\approx2.545\). A value of \(D=0.05\) thus fits less well than zero, but it is not separated from the maximum by a large likelihood drop on this scale. For \(D=0.10\), the probabilities are \(0.29\), \(0.105\), and \(0.0525\), giving

\ell(0.10)-\ell(0) =30\log(0.29/0.25)+60\log(0.105/0.125) \approx 30(0.1484)+60(-0.1744) \approx -6.01.

This larger displacement has a substantially worse fit. The example shows how a likelihood profile distinguishes a modest change from a change that the data oppose more strongly. Its symmetry also means the profile alone gives equal support to positive and negative \(D\) values here; it does not select one sign.

Using Profiles Without Overreading Them

In an unrestricted analysis, the profile for a haplotype frequency is found by fixing that frequency at a sequence of feasible values and maximizing the observed-data likelihood over the other frequencies at each value. The resulting curve should be checked over the full feasible range, not only near the estimate. In particular, a profile can reveal a broad shoulder, a boundary, or more than one region of good fit that a point estimate conceals.

There is a practical distinction between optimization uncertainty and statistical uncertainty. Optimization uncertainty concerns whether a numerical routine found the intended maximum. Statistical uncertainty concerns how much the observed genotypes distinguish that maximum from alternatives. Checking convergence from multiple starting values can help diagnose the first; examining the profile addresses the second. Neither check repairs an incorrect sampling model.

Specify the genotype model.
List the compatible haplotype pairs for each observed genotype and state the assumptions that justify their probabilities.

Construct the observed-data likelihood.
Sum over compatible pairs within each genotype probability, then combine those probabilities using the observed genotype counts.

Profile the target frequency.
At each feasible target value, maximize over the remaining frequencies while respecting nonnegativity and the sum-to-one constraint.

Interpret the shape.
Compare likelihood drops, check for asymmetry or multiple peaks, and distinguish relative-fit evidence from a formal confidence statement.

A common mistake is to say, “The algorithm estimated this frequency, so its uncertainty must be small.” That conclusion does not follow. The algorithm reports an optimizer, not the width or shape of the likelihood around it. The correction is to inspect the profile: a wide range of nearly comparable likelihoods indicates weak precision, even if the maximizing value is reported to several decimal places. Conversely, a flat or symmetric profile can arise from phase ambiguity or model structure, not simply from a software failure.

Key takeaway: For unphased two-locus genotypes, add the probabilities of all compatible haplotype pairs before forming the likelihood. Then use a profile—not just a point estimate—to see how strongly the observed genotypes support alternative haplotype frequencies.

Check Your Understanding

Answer each question by referring to the observed-data likelihood and its assumptions.

Why does the probability of a double heterozygote contain two terms, and why are those terms added?
Under random pairing, what is the probability of the unordered diplotype \(AB/ab\)? How does it differ from the probability of \(AB/AB\)?
In a profile likelihood, what happens to nuisance frequencies when the target frequency is fixed at a candidate value?
What does a broad likelihood profile suggest, and why is that different from saying the optimizer failed?
In the worked one-parameter model, why do \(D=0.05\) and \(D=-0.05\) have the same likelihood?