Estimating Allele Frequency Under an Inbreeding Model

Why Estimate Allele Frequency and Inbreeding Together?

Genotype counts contain information about both how common an allele is and how often the two alleles occur together in a homozygote. Under random mating, genotype probabilities are determined by allele frequency alone. An inbreeding model adds a second parameter that allows the heterozygote probability to be lower—or, in a broader statistical model, higher—than that random-mating value. Estimating the two parameters together separates these aspects of the genotype distribution.

The observations here are diploid genotypes at one biallelic locus. Let $A$ be the allele of interest, with population frequency $p$, and let $a$ have frequency $q=1-p$. Write the genotype counts as $n_{AA}$, $n_{Aa}$, and $n_{aa}$, with total sample size $n=n_{AA}+n_{Aa}+n_{aa}$. The model assumes individuals are sampled independently from a population with a common $p$ and a common inbreeding coefficient $F$, and that genotype calls are correct. Population structure, related sampling, or genotype error can violate these assumptions and affect the estimates.

Definition: The biallelic inbreeding model assigns genotype probabilities $$ P(AA)=p^2+Fpq,\qquad P(Aa)=2pq(1-F),\qquad P(aa)=q^2+Fpq, $$ where $q=1-p$. The parameter $F$ measures departure of heterozygosity from its random-mating value: the model's heterozygote probability is $1-F$ times $2pq$. In the usual inbreeding interpretation, $F\geq 0$; the broader probability model can also allow negative $F$, provided all three probabilities are nonnegative.

For $0<p<1$, the probabilities must be nonnegative. In particular, $F\leq 1$ is required for the heterozygote probability to be nonnegative, while the homozygote probabilities impose lower bounds on $F$. Specifically, $F\geq -p/q$ and $F\geq -q/p$. Thus the feasible range is $\max(-p/q,-q/p)\leq F\leq1$. At $p=0$ or $p=1$, only the corresponding homozygote can occur and $F$ cannot be identified from genotype probabilities.

Likelihood and the Estimates

Under independent sampling, the three genotype counts have a multinomial distribution with the probabilities above. Ignoring the multinomial coefficient, which does not depend on $p$ or $F$, the likelihood and log-likelihood are:

L(p,F)\propto (p^2+Fpq)^{n_{AA}} \{2pq(1-F)\}^{n_{Aa}} (q^2+Fpq)^{n_{aa}}, \qquad \ell(p,F)=\sum_{g\in\{AA,Aa,aa\}} n_g\log P(g).

The observed genotype proportions provide a direct route to the maximum-likelihood estimates when the fitted proportions correspond to feasible parameter values. Define $r_{AA}=n_{AA}/n$, $r_{Aa}=n_{Aa}/n$, and $r_{aa}=n_{aa}/n$. The allele frequency represented by these proportions is $r_{AA}+r_{Aa}/2$: each $AA$ contributes two copies of $A$, and each $Aa$ contributes one.

Result: If all three observed genotype proportions are positive and the resulting parameter values are feasible, the maximum-likelihood estimates are $$ \hat{p}=\frac{2n_{AA}+n_{Aa}}{2n}, \qquad \hat{q}=1-\hat{p}, \qquad \hat{F}=1-\frac{n_{Aa}/n}{2\hat{p}\hat{q}}. $$ The estimate $\hat{p}$ is the observed allele-copy proportion, and $\hat{F}$ compares observed heterozygosity with the heterozygosity $2\hat{p}\hat{q}$ implied by random mating.

To justify the result, first consider the three probabilities as free category probabilities $(\pi_{AA},\pi_{Aa},\pi_{aa})$ summing to one. Their multinomial likelihood is maximized at the observed proportions $(r_{AA},r_{Aa},r_{aa})$. In the interior, the inbreeding model can reproduce any such three proportions: set $p=r_{AA}+r_{Aa}/2$, then solve the heterozygote equation for $F$. The remaining homozygote probability follows because the probabilities sum to one. Therefore, when this solution is feasible, it reaches the multinomial maximum.

p=r_{AA}+\frac{r_{Aa}}{2}, \qquad F=1-\frac{r_{Aa}}{2pq}.

This result also explains why estimating $F$ does not change the allele-copy estimate in the unconstrained model: the allele frequency is the genotype distribution's marginal frequency of $A$. If $F$ is restricted to be nonnegative and the unconstrained estimate is negative, the unconstrained solution is not allowed. The maximum then lies at the boundary $F=0$, and $p$ must be estimated under that restriction.

Worked Example: Estimating Both Parameters From Genotype Counts

Suppose a hypothetical sample has $n_{AA}=48$, $n_{Aa}=40$, and $n_{aa}=12$, for $n=100$. The counts are all positive, so an interior fit can be considered.

Allele frequency. There are $2(48)+40=136$ copies of $A$ among $2(100)=200$ allele copies:

\hat{p}=\frac{2(48)+40}{2(100)}=\frac{136}{200}=0.6800, \qquad \hat{q}=1-0.6800=0.3200.

Inbreeding coefficient. The observed heterozygote proportion is $40/100=0.4000$, while the fitted random-mating heterozygote probability is $2(0.6800)(0.3200)=0.4352$. Thus

\hat{F}=1-\frac{0.4000}{0.4352} =1-0.9191 \approx 0.08088.

As a check, the fitted probabilities are $0.6800^2+(0.08088)(0.6800)(0.3200)\approx0.4800$, $2(0.6800)(0.3200)(1-0.08088)\approx0.4000$, and $0.3200^2+(0.08088)(0.6800)(0.3200)\approx0.1200$. They reproduce the observed proportions, up to rounding. In this hypothetical sample, the fitted heterozygosity is about 8.1% below its random-mating value at the estimated allele frequency.

Negative Estimates and Parameter Constraints

An estimate below zero is not automatically an arithmetic error. It indicates more heterozygotes than expected under random mating at the estimated allele frequency. Such a value is feasible in the broader genotype-probability model if the three probabilities remain nonnegative. It does not have the usual interpretation of positive inbreeding, however. If the scientific model specifically requires $F\geq0$, the estimate must be obtained subject to that constraint rather than reported as an allowed negative value.

Worked Example: Excess Heterozygosity and a Nonnegative Constraint

Consider hypothetical counts $n_{AA}=20$, $n_{Aa}=60$, and $n_{aa}=20$, with $n=100$. The allele-copy estimate is

\hat{p}=\frac{2(20)+60}{200}=0.5000, \qquad \hat{q}=0.5000.

The observed heterozygote proportion is $60/100=0.6000$, compared with $2\hat{p}\hat{q}=2(0.5000)(0.5000)=0.5000$ under random mating. Therefore,

\hat{F}=1-\frac{0.6000}{0.5000}=-0.2000.

This is feasible in the unrestricted probability model: the fitted probabilities are $0.25-0.05=0.20$ for each homozygote and $0.50(1.20)=0.60$ for the heterozygote. They equal the observed proportions. If the model instead imposes $F\geq0$, the unconstrained fit is excluded. For any fixed $F=0$, the likelihood is maximized at the allele-copy proportion $p=0.5000$; consequently the constrained maximum is $(\hat{p},\hat{F})=(0.5000,0)$. The constraint changes the fitted genotype probabilities to $0.25,0.50,0.25$, rather than reproducing the observed proportions.

Quantifying Uncertainty With Likelihood Information

Point estimates do not show how much the estimates would vary across repeated samples. For an interior fit, a standard large-sample approximation uses the Fisher information matrix. Let $\boldsymbol{\pi}(p,F)$ be the vector of the three genotype probabilities. For one sampled individual, its information matrix is the sum of the outer products of probability derivatives, each weighted by the reciprocal of its category probability. For $n$ independent individuals, information is multiplied by $n$. The inverse information matrix approximates the covariance matrix of the maximum-likelihood estimates.

I(p,F)=n\sum_{g\in\{AA,Aa,aa\}} \frac{1}{P(g)} \begin{pmatrix} \frac{\partial P(g)}{\partial p}\\[2pt] \frac{\partial P(g)}{\partial F} \end{pmatrix} \begin{pmatrix} \frac{\partial P(g)}{\partial p} & \frac{\partial P(g)}{\partial F} \end{pmatrix}, \qquad \widehat{\operatorname{Cov}} \begin{pmatrix}\hat{p}\\\hat{F}\end{pmatrix} \approx I(\hat{p},\hat{F})^{-1}.

The derivatives needed in this calculation are $\partial P(AA)/\partial p=2p+F(1-2p)$, $\partial P(Aa)/\partial p=2(1-F)(1-2p)$, and $\partial P(aa)/\partial p=-2q+F(1-2p)$. The derivatives with respect to $F$ are $pq$, $-2pq$, and $pq$, respectively. Estimated standard errors are the square roots of the diagonal entries of the estimated covariance matrix. The off-diagonal entry records estimated covariance between $\hat{p}$ and $\hat{F}$.

Worked Example: Information-Based Standard Errors

Return to a hypothetical sample of $n=250$ with counts $133,84,33$ for $AA,Aa,aa$, respectively. The fitted proportions are $0.532,0.336,0.132$. They give $\hat{p}=0.532+0.336/2=0.7000$ and $\hat{F}=1-0.336/\{2(0.7000)(0.3000)\}=0.2000$. These interior estimates reproduce the three observed proportions.

At these fitted values, the probability derivatives with respect to $p$, in genotype order $AA,Aa,aa$, are $1.32,-0.64,-0.68$. The derivatives with respect to $F$ are $0.21,-0.42,0.21$. For example, the information entry for $p$ per individual is

\frac{1.32^2}{0.532} +\frac{(-0.64)^2}{0.336} +\frac{(-0.68)^2}{0.132} =3.275188+1.219048+3.503030 =7.997266.

Using the same formula for the other entries gives the information per individual and the total information:

\frac{I(\hat{p},\hat{F})}{250} \approx \begin{pmatrix} 7.9973 & 0.2392\\ 0.2392 & 0.9420 \end{pmatrix}, \qquad I(\hat{p},\hat{F})^{-1} \approx \begin{pmatrix} 0.0005040 & -0.0001280\\ -0.0001280 & 0.004279 \end{pmatrix}.

Thus $SE(\hat{p})=\sqrt{0.0005040}\approx0.02245$, and $SE(\hat{F})=\sqrt{0.004279}\approx0.06541$. For orientation, a normal-approximation 95% interval for $p$ is $0.7000\pm1.96(0.02245)$, or approximately $0.6560$ to $0.7440$. The corresponding approximation for $F$ is $0.2000\pm1.96(0.06541)$, or approximately $0.0718$ to $0.3282$. These are approximate intervals based on large-sample likelihood theory; they do not assign a 95% probability to a fixed parameter lying in the calculated interval.

Interpretation and Practical Cautions

The estimated $p$ describes the frequency of allele $A$ under the assumed sampling model. The estimated $F$ describes the relative heterozygote deficit or excess at that frequency. A positive estimate can reflect inbreeding, but a deficit can also arise from population mixture, nonrandom sampling, or other departures from the model. The genotype counts alone do not establish a causal explanation.

Information-based standard errors and normal intervals work best when the sample is sufficiently large and the estimates are well inside their feasible parameter region. Near a boundary—such as $F=0$ under a nonnegative constraint—or with sparse genotype counts, the sampling distribution may be asymmetric and the normal approximation may be poor. In those cases, likelihood-based intervals or a carefully designed parametric bootstrap can better reflect the constraint and finite-sample uncertainty. If any observed genotype category has count zero, the simple interior formulas may land at a boundary; the constrained likelihood should then be maximized directly.

Common pitfall: Treating every negative $\hat{F}$ as an impossible result, or silently replacing it with zero, changes the model and potentially the maximum-likelihood fit. First decide whether the model permits negative $F$. If it does, check that all genotype probabilities are valid and interpret the estimate as heterozygote excess. If it requires $F\geq0$, maximize the likelihood under that constraint and report that the fit is on the boundary when appropriate.

The core calculation is therefore compact: estimate $p$ from allele copies, compare observed heterozygosity with $2\hat{p}\hat{q}$ to estimate $F$, and use the joint likelihood information—not a one-parameter calculation—to quantify uncertainty in both estimates.

Check Your Understanding

Use the model and procedures in this tutorial to answer the following questions.

For counts $n_{AA}=30$, $n_{Aa}=40$, and $n_{aa}=30$, calculate $\hat{p}$ and $\hat{F}$.
Why does the allele-frequency estimate use one-half of the heterozygote count?
What does a negative estimate of $F$ say about observed heterozygosity relative to $2\hat{p}\hat{q}$?
Why can the observed genotype proportions be used to obtain both maximum-likelihood estimates when the fitted parameters are interior and feasible?
What does the off-diagonal entry of the estimated covariance matrix describe?
Name one reason a normal-approximation interval may be unreliable when the fitted $F$ is near a boundary.