Introduction
A likelihood describes how compatible observed genetic data are with different values of a model parameter. The maximum-likelihood estimation problem is to find the parameter value, among those allowed by the model, that makes the observed data most likely. Sometimes the answer follows from a closed-form formula. In other genetic models, including models with uncertain genotypes or several interacting parameters, it must be found numerically.
The computational task is only part of the work. An estimate must satisfy the model's parameter constraints, the search must have converged to a relevant solution, and that solution must be checked against other possible maxima, including boundary values. A software optimizer's output is a candidate until these checks support it.
From a Genetic Likelihood to an Estimate
Start with the likelihood for the observed data, using the genetic model and assumptions appropriate to those data. If individuals or families contribute independent observations under the model, their likelihood contributions multiply; the log-likelihood then adds the log contributions. This can make optimization easier, particularly when the product contains many small probabilities.
Next, specify the allowed values of every parameter. An allele frequency lies in \([0,1]\), and a penetrance probability also lies in \([0,1]\). A variance must be positive. The parameter space matters: an unconstrained calculation can produce a value that is mathematically stationary but biologically or statistically invalid. If the optimum may lie at a boundary, that boundary must be included in the search.
For a single parameter with a tractable likelihood, differentiation or a known closed-form result may identify a candidate estimate. For more complicated models, a numerical optimizer evaluates the likelihood at a sequence of parameter values and searches for a larger value. Either approach still requires checking that the result is a maximum over the allowed space, not merely a point where the optimizer stopped.
Write the likelihood for the observed genotypes, phenotypes, read counts, or family data, and identify the parameter or parameters to estimate.
State the parameter space explicitly, including closed boundaries such as probabilities equal to zero or one when the model permits them.
Use a closed-form expression when available; otherwise optimize the likelihood or log-likelihood numerically while respecting the constraints.
Check convergence and parameter validity, then compare relevant boundaries and other candidate solutions to determine whether the result is a maximum over the full parameter space.
Worked Example: Allele-Frequency Estimation Under Hardy–Weinberg Equilibrium
Worked Example: Allele-Frequency Estimation Under Hardy–Weinberg Equilibrium
Suppose a hypothetical sample has 18 \(AA\), 24 \(Aa\), and 8 \(aa\) individuals at a biallelic locus. Assume Hardy–Weinberg equilibrium, with \(p\) the frequency of allele \(A\) and \(1-p\) the frequency of allele \(a\). Estimate \(p\).
For genotype counts under this model, the closed-form maximum-likelihood estimate is the observed fraction of allele copies that are \(A\). The sample contains \(2(18)+24=60\) copies of \(A\) among \(2(18+24+8)=100\) total allele copies. Thus $$ \hat{p}=\frac{2n_{AA}+n_{Aa}}{2n} =\frac{2(18)+24}{2(50)} =\frac{60}{100} =0.600. $$
This estimate is in the allowed range \([0,1]\). In context, the fitted model estimates the \(A\)-allele frequency in the sampled population to be 0.600, with the \(a\)-allele frequency estimated as \(1-0.600=0.400\). This is a point estimate under the stated equilibrium model; it does not by itself establish that the model assumptions hold.
Closed-Form Estimation for Genotype-Specific Penetrance
A simple closed-form result arises when the observations are independent binary outcomes with a common probability. For \(r\) observed cases among \(m\) individuals in a genotype group, the likelihood for that group's penetrance \(\theta\) is proportional to \(\theta^r(1-\theta)^{m-r}\), for \(0\leq\theta\leq1\). If both cases and non-cases are observed, maximizing this likelihood gives \(\hat{\theta}=r/m\). The same calculation can be applied separately to groups when their penetrances are distinct parameters.
Worked Example: Estimating Two Genotype-Specific Penetrances
In a hypothetical study, 7 of 20 individuals with genotype \(G\) have a binary phenotype, while 3 of 15 individuals without \(G\) have it. Assume independent Bernoulli outcomes within and between the groups, with penetrances \(\theta_G\) and \(\theta_0\), respectively.
The likelihood, omitting factors that do not depend on the parameters, is $$ L(\theta_G,\theta_0) \propto \theta_G^7(1-\theta_G)^{13} \theta_0^3(1-\theta_0)^{12}, \qquad 0\leq\theta_G,\theta_0\leq1. $$
For the \(G\) group, the log-likelihood contribution is \(7\log(\theta_G)+13\log(1-\theta_G)\). Its derivative is \(7/\theta_G-13/(1-\theta_G)\); setting this equal to zero gives \(7(1-\theta_G)=13\theta_G\), so \(\hat{\theta}_G=7/20=0.350\). For the group without \(G\), the same calculation gives \(\hat{\theta}_0=3/15=0.200\).
Both estimates lie within their probability constraints. Each group's log-likelihood has second derivative \(-r/\theta^2-(m-r)/(1-\theta)^2\), which is negative throughout the interior when \(0<r<m\). Thus each interior stationary point is a maximum for its group's likelihood. The fitted likelihood has log value $$ 7\log(0.350)+13\log(0.650)+3\log(0.200)+12\log(0.800) \approx -20.455. $$ The estimates describe phenotype probabilities within these genotype groups under the assumed model. They do not show that genotype causes the phenotype.
Numerical Optimization and Maximum Checks
A numerical method is useful when the likelihood has no convenient closed-form maximizer. It repeatedly updates a parameter value in an attempt to increase the likelihood, stopping according to a numerical criterion such as a small change in the parameter or objective. The method may operate on the likelihood or, more commonly for products of probabilities, the log-likelihood. Since the logarithm is increasing, these objectives have the same maximizer wherever the likelihood is positive.
A stopping rule does not prove that the optimizer found the global maximum. It may have found a local maximum, stopped near a flat region, or failed because of poor scaling or an unsuitable starting value. For a bounded one-dimensional parameter, evaluate the objective at any interior candidate and at both endpoints. For several parameters, use starts at multiple plausible values, inspect the reported convergence, and compare solutions. A transformation can help enforce constraints—for example, mapping an unconstrained real value to a probability between zero and one—but the boundary cases still need attention if they are permitted and relevant.
Worked Example: Checking a Bounded Numerical Maximum
Suppose a hypothetical genetic model with uncertain genotype information yields the following likelihood as a function of an allele-frequency parameter \(p\), up to a fixed positive multiplier: $$ L(p)=0.1+0.64p-1.11p^2+0.40p^3+0.07p^4, \qquad 0\leq p\leq1. $$ A bounded numerical search returns a stationary point near \(p=0.3676\). Check whether it is a maximum over the full allowed interval.
The derivative is \(L'(p)=0.64-2.22p+1.20p^2+0.28p^3\). A numerical root search gives \(p\approx0.3676\). To check the shape over the interval, note that \(L''(p)=-2.22+2.40p+0.84p^2\), which is strictly increasing on \([0,1]\), since its derivative is \(2.40+1.68p>0\). It changes sign once, near \(p=0.736\). Therefore \(L'\) first decreases and then increases. Since \(L'(0)=0.64>0\), while \(L'(1)=-0.10<0\), and the minimum of \(L'\) is negative, \(L'\) has exactly one root in \([0,1]\). It is positive before that root and negative after it, so the likelihood increases and then decreases: the root is the unique interior maximum.
Substitution provides an independent numerical check. At \(p=0.3676\), \(p^2=0.13512976\), \(p^3\approx0.04967370\), and \(p^4\approx0.01826005\). Hence $$ \begin{aligned} L(0.3676) &=0.1+0.64(0.3676)-1.11(0.13512976)\\ &\quad+0.40(0.04967370)+0.07(0.01826005)\\ &\approx0.2064. \end{aligned} $$
The endpoints give \(L(0)=0.1000\) and \(L(1)=0.1000\), both smaller than the interior value. The parameter is within \([0,1]\), and the derivative changes from positive to negative at the candidate. These checks support \(\hat{p}\approx0.3676\) as the global maximum on the stated interval. The likelihood value is not itself a probability that the estimate is correct; it is the objective evaluated at that parameter value.
Why the Checks Matter
Genetic likelihoods often combine several sources of information: genotype counts, uncertain genotype calls, phenotype outcomes, or transmissions within families. Summing over an unobserved genotype can make the likelihood nonlinear, and several parameters may be constrained simultaneously. In such cases, a numerical solution can depend on starting values or lie close to a boundary. Careful checking helps distinguish a valid model estimate from an artifact of the search.
A useful final record of an optimization includes the fitted parameter values, their constraints, the maximized log-likelihood, the method and starting values used, and the checks made for convergence and competing maxima. If estimates repeatedly differ across starting points, or a parameter lands at a boundary, that is a reason to investigate the likelihood and model rather than report a single unexplained software result. An estimate at a boundary can be the genuine maximum; it should not be discarded merely because it is unusual.
Common Mistakes
- Treating any optimizer output as the MLE. A reported stopping point may be a local solution or a boundary-incomplete search. Check the parameter space, convergence, and competing candidates.
- Ignoring parameter constraints. A probability estimate below zero or above one is invalid for the stated model. Optimize over the allowed parameter space rather than accepting an unconstrained result.
- Checking only an interior stationary point. A maximum may occur at \(p=0\), \(p=1\), or another parameter boundary. Compare boundary likelihoods whenever those values are allowed.
- Confusing a high likelihood with certainty. The likelihood compares parameter values for the observed data; it is not the probability that a parameter value is true.
- Interpreting a fitted genetic association causally. A likelihood estimate describes the fitted model and data. Causal claims require additional design and assumptions.
Check Your Understanding
Try these questions using the definitions and checks from this tutorial.
- In a likelihood problem, why does maximizing the log-likelihood give the same estimate as maximizing the likelihood?
- A sample has 10 \(AA\), 20 \(Aa\), and 10 \(aa\) individuals under a Hardy–Weinberg model. What is the closed-form estimate of the \(A\)-allele frequency?
- For a genotype group with 4 affected individuals among 10, what is the Bernoulli maximum-likelihood estimate of its penetrance? State the estimate as a probability.
- A bounded optimizer reports an interior candidate with a larger likelihood than one starting point. What additional checks are needed before calling it the global maximum?
- Why should an estimate at a permitted boundary not automatically be treated as an optimization failure?