Worked Empirical Bayes Analysis of Genomic Tests

From Z-Statistics to a Fitted Mixture

A genomic analysis can produce thousands of test statistics, but a local false discovery rate for one test depends on more than its z-statistic alone. It also depends on the estimated distribution of statistics under the null and on how much probability the fitted model assigns to the null group. A worked calculation makes those dependencies visible: a modest change in a component density or mixture weight can change the posterior null probability.

This tutorial works through a two-group normal mixture. It treats the z-statistics as observations from a population of tests and uses an empirical Bayes model to estimate the groups from their overall distribution. The fitted posterior quantities are conditional on that model; they do not prove that any particular gene is null or non-null.

Definition: Let \(Z\) denote a test statistic drawn from a mixture with null proportion \(\pi_0\). Its density is \(f(z)=\pi_0 f_0(z)+(1-\pi_0)f_1(z)\), where \(f_0\) is the null density and \(f_1\) is the alternative density. The local false discovery rate at \(z\) is \(\operatorname{lfdr}(z)=P(\text{null}\mid Z=z)=\frac{\pi_0 f_0(z)}{\pi_0 f_0(z)+(1-\pi_0)f_1(z)}\).

For the examples, use a two-normal model:

Z \sim \pi_0 N(\mu_0,\sigma_0^2)+(1-\pi_0)N(\mu_1,\sigma_1^2).

Here \(N(\mu,\sigma^2)\) denotes a normal distribution with mean \(\mu\) and variance \(\sigma^2\). In a theoretical-null analysis, \(f_0\) is fixed at the standard normal density, with \(\mu_0=0\) and \(\sigma_0=1\). In an empirical-null analysis, \(\mu_0\) and \(\sigma_0\) are estimated, allowing for a systematic shift or excess spread in the central test-statistic distribution. The alternative distribution may also be estimated; it need not be a standard normal.

A basic two-group fit assumes the test statistics are exchangeable draws from the same mixture, that the component densities are adequately represented by the chosen model, and that the mixture is identifiable enough for estimation. In practice, dependence among genomic tests can affect uncertainty and calibration, and different parameter combinations can fit similarly. A fitted mixture is therefore a modeling approximation, not a substitute for checking the data’s distribution.

Estimating the Mixture With EM

The EM algorithm alternates between estimating each observation’s probability of belonging to each component and updating the component parameters using those probabilities as weights. Let \(z_1,\ldots,z_m\) be the test statistics, and let \(r_i\) be the current posterior probability that observation \(i\) is null. Given current parameter values,

r_i=\frac{\pi_0 f_0(z_i)} {\pi_0 f_0(z_i)+(1-\pi_0)f_1(z_i)}.

The M-step updates the null proportion to the average null responsibility and updates each normal component using its weighted mean and variance. For the null component, for example,

\pi_0^{\mathrm{new}}=\frac{1}{m}\sum_{i=1}^m r_i,\qquad \mu_0^{\mathrm{new}}=\frac{\sum_{i=1}^m r_i z_i}{\sum_{i=1}^m r_i},\qquad (\sigma_0^2)^{\mathrm{new}}=\frac{\sum_{i=1}^m r_i(z_i-\mu_0^{\mathrm{new}})^2}{\sum_{i=1}^m r_i}.

The alternative updates use weights \(1-r_i\) in the same formulas. Repeat the E- and M-steps until the log-likelihood or parameter estimates change by less than a chosen tolerance. The fit can depend on initialization; multiple starts and substantive constraints or checks may be needed to avoid a poor or label-switched solution. The null component is conventionally identified by its role as the component centered near the calibrated null distribution, not merely by its label in software.

Worked Example: Reading a Compact EM Fit

Consider a hypothetical genomic analysis with \(m=1{,}000\) test statistics. The displayed values below are a small subset of that panel; the fit uses all 1,000 statistics. After EM convergence under a two-normal model, suppose the fitted parameters are \(\hat{\pi}_0=0.5428\), \(\hat{\mu}_0=0\), \(\hat{\sigma}_0=0.88\), \(\hat{\mu}_1=2.00\), and \(\hat{\sigma}_1=1.00\). These are illustrative fitted values, not estimates from a real study.

The implied null proportion is 0.5428, so the fitted alternative proportion is \(1-0.5428=0.4572\). The empirical null is \(N(0,0.88^2)\): its mean remains zero, while its standard deviation is below one. This fitted spread is a feature of this hypothetical panel, not a general expectation for genomic z-statistics.

For a single EM update, suppose an observation \(z_i=0.8\) has current null and alternative densities \(f_0(0.8)=0.2999\) and \(f_1(0.8)=0.1942\). Its current null responsibility is \(\frac{0.5428(0.2999)}{0.5428(0.2999)+0.4572(0.1942)}\). The null-weighted contribution is \(0.5428(0.2999)=0.16278572\); the alternative-weighted contribution is \(0.4572(0.1942)=0.08878824\). Thus the responsibility is approximately \(0.16278572/0.25157396=0.6471\). In an M-step this observation contributes 0.6471 to the total null weight and 0.3529 to the alternative weight. Its contribution alone does not determine either component’s updated mean or variance; those use the weighted contributions of all 1,000 observations.

This illustrates why one test cannot identify the mixture parameters. Estimation pools information across the panel, whereas the posterior null probability for a particular test combines the fitted weights with the two densities at that test’s statistic.

Calculating and Interpreting Local FDR

Once parameters are fitted, local FDR is a direct application of Bayes’ rule. The numerator is the joint density contribution from the null group; the denominator is the total mixture density. Densities are not probabilities of observing an exact continuous value, but their ratio gives the conditional component probability at that value.

Worked Example: Local FDR at a Z-Statistic of 0.8

Use the fitted values from the previous example. The null density is approximately 0.2999 at \(z=0.8\), and the alternative density is approximately 0.1942. The weighted null contribution is \(0.5428(0.2999)=0.16278572\), and the weighted alternative contribution is \(0.4572(0.1942)=0.08878824\).

The total mixture density is \(0.16278572+0.08878824=0.25157396\). Therefore,

\widehat{\operatorname{lfdr}}(0.8) =\frac{0.5428(0.2999)} {0.5428(0.2999)+0.4572(0.1942)} =\frac{0.16278572}{0.25157396} \approx 0.6471.

Under this fitted model, a test with \(z=0.8\) has an estimated posterior null probability of about 0.6471, or 64.71%. This is not a statement that there is a 64.71% chance the null hypothesis is true in a frequentist sense; it is a model-based conditional probability, given the observed statistic and the fitted mixture.

Worked Example: Comparing a More Extreme Statistic

For \(z=2.0\), keep the same fitted parameters. The null density is \(\frac{1}{0.88\sqrt{2\pi}}\exp[-\frac12(2/0.88)^2]\), approximately 0.0342609. The alternative density under \(N(2,1)\) is \(\frac{1}{\sqrt{2\pi}}\exp[-\frac12(2-2)^2]\), approximately 0.3989423.

The weighted null contribution is \(0.5428(0.0342609)\approx0.0185968\); the weighted alternative contribution is \(0.4572(0.3989423)\approx0.1823964\). Thus \(\widehat{\operatorname{lfdr}}(2.0)\approx0.0185968/(0.0185968+0.1823964)=0.09252\).

The estimated null probability is about 9.25% at \(z=2.0\), much lower than at \(z=0.8\) in this fitted model. The comparison follows from the relative component densities and mixture weights; it is not a universal rule that every statistic of magnitude 2 has this local FDR.

Worked Example: Averaging Local FDRs for a Small Set

Suppose three hypothetical tests have fitted local FDR values 0.6471, 0.0925, and 0.1200. If these three tests are selected as a set, the average fitted local FDR is \((0.6471+0.0925+0.1200)/3=0.8596/3=0.2865\).

Under the fitted model, this average estimates the expected fraction of null tests in this particular set, conditional on the observed statistics. It is not the probability that all three are null, nor does it guarantee that exactly 28.65% are null. The interpretation depends on the model’s calibration and on the chosen set.

Why the Empirical Null and Diagnostics Matter

If the null distribution is fixed incorrectly, the numerator of local FDR can be too large or too small across a broad range of statistics. An empirical null can accommodate central shifts or spread changes that a standard normal null cannot. But adding estimated null parameters also adds modeling choices and can make the fit less stable. A two-group model is useful when its components provide a reasonable description of the test-statistic distribution, not simply because it returns a number for every test.

Useful checks include comparing the fitted mixture density with a histogram or density estimate, inspecting quantile-quantile behavior in the central region, checking whether fitted component parameters are scientifically plausible, and repeating the fit from different initial values. A fitted density that misses the center or tails, a null component displaced implausibly far from zero, or estimates that change markedly with initialization are warnings. Dependence among tests also means that a large number of statistics need not provide as much independent information as the same number of independent draws.

The earlier material on local FDR from z-scores develops the posterior-density calculation; this worked analysis emphasizes how estimated mixture parameters enter it. A standard-null calculation fixes \(f_0\) at \(N(0,1)\), which is simpler and avoids estimating the null mean and variance. An empirical-null mixture estimates those quantities from the observed panel and may better reflect systematic departures from the theoretical null. Neither choice is automatically superior: conclusions can be sensitive to the null model, alternative density, and fit constraints.

A common mistake is to report the null density alone as the local FDR, or to multiply the component densities by the wrong mixture weights. The correction is to form both weighted contributions, add them for the mixture density, and divide the weighted null contribution by that total. In particular, with the rounded inputs above, \(0.5428(0.2999)=0.16278572\), which rounds to 0.1628, and the resulting local FDR rounds to 0.6471. Intermediate rounding should be delayed so that the final ratio is not shifted by avoidable arithmetic error.

Key takeaway: A local FDR is the fitted null contribution divided by the total fitted mixture density at the observed statistic. Its interpretation is conditional on the estimated model, so report the model and check whether its empirical null and mixture fit are credible.

Check Your Understanding

Answer each question using the definitions and calculations above. Solutions are included for self-checking.

In a two-group mixture, what are the numerator and denominator of the local FDR?
For a fitted null proportion of 0.60, what is the fitted alternative proportion?
Why does an empirical null differ from a theoretical standard normal null?
Using the worked \(z=0.8\) calculation, what is the estimated local FDR, and what does it mean?
Name two diagnostics that could flag a poor mixture fit.

Solutions. (1) The numerator is \(\pi_0f_0(z)\), the null contribution; the denominator is \(\pi_0f_0(z)+(1-\pi_0)f_1(z)\), the total mixture density. (2) \(1-0.60=0.40\). (3) A theoretical null fixes mean 0 and standard deviation 1, whereas an empirical null estimates its parameters from the observed test statistics. (4) It is approximately 0.6471, a model-based posterior null probability of about 64.71% conditional on \(z=0.8\). (5) Examples include poor agreement between fitted density and observed distribution, implausible null parameters, and substantial changes across EM initializations.