The EM Algorithm for Two-Group Mixtures

Why Use EM for a Two-Group Mixture?

In a two-group model for genomic test statistics, each observed statistic is treated as coming from either a null distribution or an alternative distribution. The group membership is not observed. This makes direct likelihood maximization harder than it would be if the memberships were known, because each observation's density is a weighted sum of two component densities. The expectation-maximization (EM) algorithm handles this by alternating between estimating the unobserved memberships probabilistically and updating the model parameters using those probabilities as weights.

EM is an optimization procedure for the chosen mixture model; it does not decide whether the null and alternative distributions are scientifically appropriate. Its updates depend on the model and on which parameters are fixed or estimated. Here the central derivation is for a two-component mixture, followed by a concrete normal-mixture example.

Model and Assumptions

The likelihood used below treats the \(z_i\)'s as independent draws from this common mixture, so the observed-data log likelihood is the sum of the log marginal densities. This independence assumption may be an approximation for correlated genomic tests. In that setting, the algorithm can still maximize the product of marginal mixture densities as a working likelihood, but the resulting fit does not account for dependence.

The component densities must be valid densities, the mixing proportions must be nonnegative and sum to one, and any parameter constraints in the chosen component models must be respected. The null and alternative are assigned distinct meanings by the model; if their distributions are both unconstrained, exchanging their labels can describe the same mixture. Anchoring the null distribution to a specified or separately estimated null model gives the labels their intended interpretation.

The E-Step and M-Step

Introduce an indicator \(I(C_i=c)\), equal to one when observation \(i\) belongs to component \(c\) and zero otherwise. If the labels were observed, the complete-data log likelihood would be

Because the labels are missing, EM replaces each indicator with its conditional expectation under the current parameters. At iteration \(t\), define the responsibility \(\tau_i^{(t)}\) as the posterior probability that observation \(i\) is null, and \(\gamma_i^{(t)}\) as its posterior probability of being alternative:

These are the E-step quantities. They are probabilities, not hard classifications. The M-step maximizes the expected complete-data log likelihood, conditional on the responsibilities just computed:

If both mixing proportions are estimated without additional constraints, maximizing their part of \(Q\), subject to \(\pi_0+\pi_1=1\), gives

The component parameters are updated by maximizing their responsibility-weighted log densities. In particular, \(\theta_1^{(t+1)}\) maximizes \(\sum_i\gamma_i^{(t)}\log f_1(z_i;\theta_1)\). The corresponding null-parameter update uses weights \(\tau_i^{(t)}\), if null parameters are being estimated. The exact formulas depend on the component model; parameters fixed by the analysis are not updated.

Why the Observed Likelihood Does Not Decrease

The EM updates have a useful likelihood property when each M-step maximizes \(Q\) (or increases it). For each observation, let \(q_i(c)\) be any probabilities over the two components. Jensen's inequality applied to the log function gives a lower bound for the log marginal density:

Set \(q_i(0)=\tau_i^{(t)}\) and \(q_i(1)=\gamma_i^{(t)}\), using the posterior probabilities at the current parameter values. Summing the bound over observations gives \(Q(\theta\mid\theta^{(t)})+\sum_i H(q_i)\), where \(H(q_i)=-\sum_c q_i(c)\log q_i(c)\). The entropy terms do not depend on the candidate \(\theta\). At \(\theta=\theta^{(t)}\), the bound is exact because the \(q_i(c)\)'s are the current posterior probabilities. Therefore, if the M-step increases \(Q\), the observed log likelihood cannot decrease. This is a monotonicity result, not a guarantee that EM finds the global maximum.

A Full Numerical Iteration

Consider four hypothetical test statistics \(z=(-0.2,0.4,1.2,2.5)\). For this example, assume independent observations, a fixed null density \(N(0,1)\), and an alternative density \(N(\mu,1)\). Both variances are fixed at one. Start with \(\pi_1^{(0)}=0.20\), \(\pi_0^{(0)}=0.80\), and \(\mu^{(0)}=1\). These assumptions specify the model for this calculation; in particular, they do not establish that a normal alternative is suitable for other data.

Under these densities, the ratio of the alternative density to the null density is

With \(\mu^{(0)}=1\), the alternative responsibility can be calculated from this ratio \(r_i\) as \(\gamma_i^{(0)}=0.20r_i/(0.80+0.20r_i)\). For example, at \(z=-0.2\), \(r=\exp(-0.7)\approx0.4966\), so the alternative responsibility is \(0.20(0.4966)/[0.80+0.20(0.4966)]\approx0.1104\).

The M-step averages the alternative responsibilities and uses them as weights for the alternative mean:

The updated alternative proportion is about \(0.3196\), and the alternative mean is about \(1.623\) statistic units. These are one-iteration updates, not necessarily final estimates. A subsequent E-step must recompute all four responsibilities using the updated \(\pi_1\) and \(\mu\).

Implementation, Diagnostics, and Alternatives

A practical EM fit repeats the E- and M-steps until a declared stopping rule is met. Common checks include the change in observed log likelihood, changes in the parameter estimates, and whether all responsibilities and mixture proportions remain within their allowed ranges. Likelihood should not decrease beyond small numerical tolerance when exact EM updates are used. Log-density calculations or stable log-sum-exp arithmetic can help avoid underflow when component densities are very small.

Large changes from different starting values, very small effective component weights, or estimates repeatedly reaching parameter constraints can indicate a difficult or weakly determined fit. These are prompts to inspect the model and initialization, not automatic proof that the implementation is wrong. For a two-component model with both components free, label switching can also make nominally different solutions equivalent after swapping labels.

Direct numerical maximization of the observed likelihood is an alternative: it optimizes the sum of log mixture densities without introducing explicit latent labels. It can be useful for checking an EM solution, but mixture likelihoods may be difficult to optimize and require constraints to keep parameters valid. EM often has simple updates and a nondecreasing likelihood, but can converge slowly near a maximum and can settle at a local rather than global maximum. Hard assignment, which assigns each observation wholly to its most likely group, is not the EM E-step: it discards uncertainty and generally does not maximize the same expected complete-data likelihood.

The method's conclusions remain conditional on the chosen mixture densities and independence working likelihood. A converged fit does not establish that the null is calibrated, that the alternative distribution captures the data, or that dependence among genomic tests is negligible. Inspect the fitted model and compare plausible initializations; EM's monotonicity is a property of optimization, not a diagnostic certificate of scientific adequacy.

Common Mistake