Tutorials › High-Dimensional Genomic Data Analysis › Maximum Likelihood for Two-Group Models

Empirical Bayes and Local FDR · Tutorial 180 of 795

Maximum Likelihood for Two-Group Models

Learn to write and optimize a two-group mixture likelihood, interpret its estimates, and recognize when dependence or model misspecification makes those estimates unreliable.

Advanced 12 min read

What You'll Learn

  • Write the likelihood and log likelihood for independent test statistics under a two-group mixture.
  • Estimate the null proportion when the component densities are specified.
  • Estimate an alternative-component parameter by maximizing the mixture likelihood.
  • Distinguish full-likelihood fitting from central matching and null-only fitting.
  • Diagnose numerical and model-fitting problems that can distort estimates.

Why Write a Likelihood for Test Statistics?

A two-group model treats each test statistic as arising from one of two latent groups: a null group and a nonnull group. The group membership is unknown for an individual test, so the observed statistic has a mixture density. Maximum likelihood estimates the mixture parameters by choosing values that make the observed collection of statistics plausible under that density. This makes explicit which assumptions connect the data to estimates of the null proportion and component distributions.

The introductory material on two-group mixture models establishes the basic mixture representation and its connection to posterior null probabilities. Here the focus is the likelihood: how it is constructed, what is optimized, and what the resulting estimates do and do not tell us. We use \(Z_1,\ldots,Z_m\) for the observed test statistics, where \(m\) is the number of tests.

Definition: In a parametric two-group model, each test statistic has density \(f(z;\theta)=\pi_0 f_0(z;\eta_0)+(1-\pi_0)f_1(z;\eta_1)\). Here \(0\leq\pi_0\leq1\) is the null proportion, \(f_0\) is the null-group density with parameter \(\eta_0\), and \(f_1\) is the alternative-group density with parameter \(\eta_1\). The parameters may be partly or fully specified in advance. The group label is latent: it is not observed for each test.

For the ordinary likelihood below, assume the test statistics are independent and each has the same marginal mixture density \(f\). Also assume the proposed densities are valid probability densities, the parameter values lie in their allowed parameter spaces, and \(f(Z_i;\theta)>0\) for every observed value at the parameter values being considered. Under these assumptions, the joint density of the observed statistics is the product of their marginal densities.

$$ L(\theta;z_1,\ldots,z_m)=\prod_{i=1}^{m} \left\{\pi_0 f_0(z_i;\eta_0)+(1-\pi_0)f_1(z_i;\eta_1)\right\}, \qquad \ell(\theta)=\sum_{i=1}^{m}\log f(z_i;\theta). $$

The maximum-likelihood estimate (MLE) is any allowed parameter value \(\hat{\theta}\) that maximizes \(L\), or equivalently \(\ell\), since the logarithm is strictly increasing. The log form is used in calculation because it replaces a potentially tiny product by a sum. Each observation contributes the logarithm of its mixture density, not a log density selected after assigning the observation to a group. The latent group is not known, so the contribution is a sum of the two weighted densities inside the logarithm.

Worked Likelihood With Specified Components

Worked Example: Evaluating a Two-Group Likelihood

Suppose three hypothetical test statistics are \(-1,0,2\). For illustration, take \(f_0\) to be the standard normal density and \(f_1\) to be the normal density with mean 0 and variance 4. Both densities are fully specified, and consider evaluating the likelihood at \(\pi_0=0.5\). Assume the three statistics are independent under this model.

The densities are \(f_0(z)=\phi(z;0,1)\) and \(f_1(z)=\phi(z;0,4)\). Their values, rounded to six decimal places, are \(f_0(-1)=0.241971\), \(f_1(-1)=0.176033\); \(f_0(0)=0.398942\), \(f_1(0)=0.199471\); and \(f_0(2)=0.053991\), \(f_1(2)=0.120985\). Each density integrates to one, the mixing proportion is between zero and one, and independence is being stipulated for this hypothetical calculation.

At \(\pi_0=0.5\), the mixture density at each observation is the average of the component densities:

$$ f(-1)=0.5(0.241971)+0.5(0.176033)=0.209002, \quad f(0)=0.5(0.398942)+0.5(0.199471)=0.299207, $$
$$ f(2)=0.5(0.053991)+0.5(0.120985)=0.087488. $$

Thus \(L=(0.209002)(0.299207)(0.087488)\), approximately \(0.005471\), and \(\ell=\log(0.209002)+\log(0.299207)+\log(0.087488)\approx-5.208\). This is the log likelihood at the specified parameter value, not evidence by itself that \(0.5\) is the MLE. To obtain an estimate, the likelihood must be compared or optimized across allowed parameter values.

Estimating the Null Proportion

When both component densities are specified, only \(\pi_0\) remains to be estimated. The log likelihood is a sum of logarithms of linear functions of \(\pi_0\). Its derivative indicates whether increasing the null proportion locally increases or decreases the likelihood:

$$ \ell'(\pi_0)= \sum_{i=1}^{m} \frac{f_0(z_i)-f_1(z_i)} {\pi_0 f_0(z_i)+(1-\pi_0)f_1(z_i)}. $$

An interior MLE satisfies \(\ell'(\hat{\pi}_0)=0\). If the derivative does not cross zero on \([0,1]\), the maximum can occur at a boundary. For fixed component densities, each term in the log likelihood has second derivative \(-[f_0(z_i)-f_1(z_i)]^2/f(z_i;\pi_0)^2\), which is nonpositive. Therefore the log likelihood is concave in \(\pi_0\); a root of the derivative is a global maximum (though it need not be unique in degenerate cases where the densities do not distinguish the groups).

Worked Example: Finding the Null-Proportion MLE

Use the three statistics and component densities from the previous example. For compactness, let \(r_i=f_1(z_i)/f_0(z_i)\). The ratios at \(-1,0,2\) are approximately \(0.7275,0.5000,2.241\), respectively. Dividing numerator and denominator in each score term by \(f_0(z_i)\) gives

$$ \ell'(\pi_0)= \sum_{i=1}^{3}\frac{1-r_i}{\pi_0+(1-\pi_0)r_i}. $$

At \(\pi_0=0.5\), the score is approximately \((0.2725/0.86375)+(0.5000/0.7500)+(-1.241/1.620)\), or about \(0.216\), so the likelihood is still increasing there. Solving the score equation numerically on \([0,1]\) gives \(\hat{\pi}_0\approx0.685\). At this value, the three score contributions are approximately \(0.298\), \(0.593\), and \(-0.892\), which sum to approximately zero, allowing for rounding. Concavity makes this interior solution a global maximum for this specified model.

The estimate says that the fitted model assigns a mixing weight of about 0.684 to its null component. It does not mean that exactly 68.4% of these three particular tests are null, nor does it identify which tests are null. That distinction matters because the model estimates a population mixture proportion from the distribution of the observed statistics.

Estimating a Component Parameter

The same likelihood can estimate parameters within a component, provided its density is specified up to those parameters. For example, an alternative distribution might have an unknown mean \(\eta_1=\mu_1\). The likelihood then has to account for the uncertainty about group membership at every observed statistic. It is not generally valid to select a subset of observations as “alternative” first and calculate a mean as if those labels were known.

When a closed-form maximizer is unavailable, numerical maximization of the log likelihood is standard. The parameter space must be respected; for example, a standard deviation must be positive and \(\pi_0\) must remain in \([0,1]\). An optimizer’s output should be checked by evaluating the likelihood near the reported solution, checking constraints, and, where practical, using multiple starting values.

Worked Example: Estimating an Alternative Mean

Consider hypothetical statistics \(1,2,3\), with the null density fixed as \(N(0,1)\). Let the alternative density be \(N(\mu_1,1)\), and fix \(\pi_0=0.5\) to isolate the estimation of the alternative mean. Assume independence. The log likelihood as a function of \(\mu_1\) is

$$ \ell(\mu_1)= \sum_{z\in\{1,2,3\}}\log\left\{ 0.5\phi(z;0,1)+0.5\phi(z;\mu_1,1) \right\}. $$

For a normal density with variance 1, its derivative with respect to its mean is \((z-\mu_1)\phi(z;\mu_1,1)\). Consequently, differentiating each mixture term gives the score

$$ \ell'(\mu_1)=\sum_{i=1}^{3} w_i(\mu_1)(z_i-\mu_1), \qquad w_i(\mu_1)= \frac{\phi(z_i;\mu_1,1)} {\phi(z_i;0,1)+\phi(z_i;\mu_1,1)}. $$

For instance, at \(\mu_1=2.24\), the alternative-to-null density ratios for \(z=1,2,3\) are approximately \(0.764,7.18,67.4\), using the ratio \(\exp(\mu_1z-\mu_1^2/2)\). Thus the corresponding weights are about \(0.433,0.878,0.985\). Substitution gives a score of approximately \(0.433(1-2.24)+0.878(2-2.24)+0.985(3-2.24)\approx0.001\). Numerical maximization gives \(\hat{\mu}_1\) about \(2.24\). The log likelihood there is about \(-4.94\), compared with about \(-9.76\) at \(\mu_1=0\), where the two component densities coincide and all three statistics have the standard normal density.

This fit is conditional on the specified null, equal component variances, fixed mixing proportion, and independence assumption. It is not a general estimate of a biological effect size; it is the fitted mean of the alternative component on the test-statistic scale under this particular model.

Joint Fitting, Diagnostics, and Alternatives

In a general parametric fit, the unknown parameter vector may include \(\pi_0\), parameters of the null density, and parameters of the alternative density. The objective remains the same sum of log mixture densities, but the parameters must be optimized jointly unless some are fixed. A numerical optimizer can return a plausible-looking value even when the model is poorly supported, so optimization and model adequacy are separate questions.

Useful checks include verifying that every evaluated density is finite and positive, confirming that estimates satisfy parameter constraints, and plotting the fitted mixture against the observed test-statistic distribution. A solution that changes substantially with starting values, lies at an unexpected boundary, or has a nearly flat likelihood deserves careful investigation. A flat likelihood indicates that quite different parameter values explain the data similarly; reporting many decimal places in that situation would overstate precision. The later material on mixture identifiability treats such issues in greater depth.

The independence assumption also needs attention. If test statistics are dependent, the product likelihood is not the true joint likelihood, even if each statistic has the stated marginal mixture density. The fitted values may still be used as a working or composite-likelihood estimate in some settings, but usual likelihood-based uncertainty calculations cannot automatically be justified by treating the tests as independent. Dependence among genomic tests is a substantive feature to diagnose, not something repaired by taking logarithms.

Full likelihood fitting differs from central matching, which estimates null parameters by matching the distribution in a chosen central region. Central matching can be useful when the center is thought to be dominated by null statistics, but it depends on the region and matching criterion. Full likelihood uses all observations and requires a credible specification of both component densities. A null-only fit instead estimates a null distribution without fitting the entire two-group mixture; it does not, by itself, estimate the same joint set of mixture parameters. These are different targets and assumptions, not interchangeable computational recipes.

The parametric approach can be efficient when its component distributions are well chosen, but it can be sensitive to misspecification. For example, a heavy-tailed null distribution may be mistaken for a larger alternative contribution if the fitted null is too narrow. Conversely, a flexible alternative may absorb features that the model intended to attribute to the null. The likelihood maximizer is the best fit within the specified model family; maximizing a misspecified likelihood does not make its assumptions true.

Common Mistakes

  • Taking logs before adding components. The contribution is \(\log\{\pi_0f_0(z)+(1-\pi_0)f_1(z)\}\), not a weighted average of \(\log f_0(z)\) and \(\log f_1(z)\).
  • Assigning each statistic to a group first. Group membership is latent. The observed-data likelihood sums the weighted component densities rather than treating an uncertain label as known.
  • Calling a likelihood value a probability that a parameter is correct. A likelihood compares parameter values for observed data; it is not a posterior probability of the parameter or of a null hypothesis.
  • Ignoring dependence or model fit. A numerical maximum does not validate the independence assumption or the component distributions.
Key takeaway: Construct each observation’s weighted mixture density first, multiply those contributions under independence, and maximize the resulting log likelihood over the allowed parameter space. Interpret the estimates conditionally on the component models and sampling assumptions that produced them.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. For independent observations, why does the likelihood contribution for one statistic involve a sum of two weighted densities rather than one selected component density?
  2. If \(f_0(z_i)=0.20\), \(f_1(z_i)=0.10\), and \(\pi_0=0.60\), calculate that observation’s mixture density.
  3. With fixed component densities, what does a zero interior score \(\ell'(\hat{\pi}_0)=0\) establish when the log likelihood is concave in \(\pi_0\)?
  4. Name one reason an optimizer’s reported maximum might not justify trusting the fitted null proportion.

Solutions. (1) The group label is latent, so the density of the observed statistic is the probability-weighted sum over the two possible groups. (2) \(0.60(0.20)+0.40(0.10)=0.12+0.04=0.160\). (3) It identifies a global maximum over the allowed interval, subject to possible non-uniqueness in degenerate cases. (4) Examples include dependence that invalidates the product likelihood, a poor component-distribution fit, a boundary estimate, or sensitivity to starting values.

Clinical Biostats