How does limma-trend work? An empirical partially Bayes perspective

Abstract: In high-throughput biology, it is common to fit thousands of linear regressions -- one per gene, protein, or other unit -- with very few samples per unit. Limma-trend, one of the most widely used methods in this setting, improves power by shrinking variance estimates parametrically toward a fitted curve (the trend) relating variance to a unit-level summary (e.g., average intensity, peptide count), before computing p-values and applying the Benjamini-Hochberg procedure to control the false discovery rate (FDR). We study limma-trend through the lens of empirical partially Bayes inference, a paradigm in which a prior is posited and estimated for the nuisance parameters while parameters of interest remain fixed. From this perspective, limma-trend computes approximate partially Bayes p-values that condition on the residual sample variance and the unit-level summary. The same framework explains why MAnorm2, a popular variant for ChIP-seq, can sometimes fail to control FDR. We then derive a nonparametric generalization of limma-trend that estimates the residual variance prior using nonparametric maximum likelihood. Under dense signals, this procedure asymptotically controls the FDR -- even when the trend is misspecified or inconsistently estimated. To allow the full shape of the conditional variance distribution to depend on the unit-level summary, we develop a second procedure that learns it directly.

• The paper establishes limma-trend as empirical partially Bayes inference, showing that trend-informed variance shrinkage produces conditionally calibrated p-values and identifying an orthogonality condition needed for valid FDR control.
• The paper introduces Reg-NPMLE and JT-NPMLE, which retain asymptotic FDR control under trend misspecification while achieving minimax or near-parametric convergence rates and competitive statistical power.
• The paper demonstrates practical benefits across RNA-seq, ChIP-seq, and proteomics, while showing that violating the design condition can inflate FDR to 12.6% for MAnorm2 and over 50% for MAP in simulations.

Overview and motivation

High-throughput biological assays routinely require fitting one linear regression per unit—gene, protein, or genomic interval—with degrees of freedom $K-p$ often in the single digits. The resulting per-unit residual variance estimates are too noisy for reliable $t$-tests, and limma-trend addresses this by shrinking variances toward a fitted mean-variance trend before computing p-values. Despite its ubiquity across RNA-seq, proteomics, ChIP-seq, methylation, cytometry, and lipidomics pipelines, the procedure has lacked a formal statistical account. This paper supplies one by interpreting limma-trend as an instance of empirical partially Bayes inference (Cox's framework): a prior is posited and estimated for the nuisance parameter $\sigma_i^2$ while the parameter of interest $\theta_i$ remains fixed. The paper then develops two nonparametric generalizations with asymptotic false discovery rate (FDR) guarantees, and identifies a previously unstated design condition whose violation explains FDR inflation in MAnorm2.

The partially Bayes interpretation of limma-trend

For each unit $i$, ordinary least squares yields a contrast estimate $Z_i$, residual variance $S_i^2$, and side-information $M_i$. Limma-trend posits that $\sigma_i^2 \mid M_i \sim \xi_0^2(M_i) G$, decomposing $\log(\sigma_i^2)$ into a trend component plus independent residual heterogeneity. The oracle partially Bayes p-value integrates the Gaussian tail probability over the posterior of $t$0 given $t$1:

$t$2

These p-values are conditionally uniform under the null given both $t$3 and $t$4—a property the plain $t$5-test lacks, since null units with by-chance small $t$6 produce inflated statistics. Parametric limma-trend is recovered as the special case where $t$7 is a scaled inverse-$t$8 distribution; the authors prove this equivalence explicitly.

A key structural issue arises because the default summary $t$9 (average intensity) is computed from the same outcomes used for inference. The paper resolves this via an orthogonality condition on the design: $\sigma_i^2$0 and $\sigma_i^2$1. This assumption—which the authors state has not been spelled out before—is satisfied by standard two-sample and treatment-effect designs but violated by MAnorm2's use of the unweighted group-mean average $\sigma_i^2$2 under unbalanced group sizes. Simulations confirm the consequence: MAnorm2 inflates FDR to 12.6% at nominal 5% in an unbalanced-design setting, while MAP—a method treating $\sigma_i^2$3 as a deterministic function of $\sigma_i^2$4—inflates FDR to over 50% when residual heterogeneity is present.

Nonparametric generalization: Reg-NPMLE

The first proposal replaces the parametric prior with the nonparametric maximum likelihood estimator (NPMLE) applied to trend-adjusted variances $\sigma_i^2$5, solved as a convex program via MOSEK. The main convergence theorem shows estimated p-values approach the oracle p-values at rate $\sigma_i^2$6, which matches the minimax rate of trend estimation up to logarithmic factors—for parametric or Hölder-smooth trends this is $\sigma_i^2$7 or $\sigma_i^2$8 respectively. Under a criticality (dense signal) condition, Benjamini–Hochberg applied to these p-values controls FDR asymptotically.

The most consequential result concerns misspecification: even if the fitted trend $\sigma_i^2$9 converges to some function $\theta_i$0 unrelated to the true data-generating process, FDR control persists. The mechanism is that the nonparametric prior on residual heterogeneity absorbs trend error—the misspecified p-values remain uniform conditional on the trend-adjusted variance $\theta_i$1. This robustness is strictly stronger than what deterministic plug-in methods like MAP can offer, and it implies practitioners need not worry about parametric spline trend fits being inconsistent. The guarantee does have a boundary: conditional calibration given $\theta_i$2 itself (e.g., within each peptide-count stratum) is lost under misspecification, which matters in applications where stratum-level error control is desired.

Joint estimation: JT-NPMLE

The second proposal dispenses with explicit trend fitting entirely. Treating $\theta_i$3 jointly as nuisance parameters drawn from an unknown bivariate prior $\theta_i$4, it estimates $\theta_i$5 by NPMLE from the marginal distribution of $\theta_i$6. The analysis rests on an Eddington–Tweedie representation expressing the oracle p-values as functionals of the marginal densities $\theta_i$7 and $\theta_i$8. Remarkably, despite the nonparametric bivariate estimation, the estimated p-values converge to the oracle at a near-parametric rate $\theta_i$9, and BH again achieves asymptotic FDR control. A discrete variant—fitting a separate NPMLE within each level of a discrete $i$0—serves as a methodological successor to DEqMS for proteomics.

Numerical evidence

Across four simulation settings ($i$1 units, 90% null), the proposed procedures are the only feasible methods besides the oracle to control FDR throughout, while matching or exceeding competitor power. In Setting 4, where the trend perfectly determines the variance (MAP's ideal scenario), Reg-NPMLE matches MAP's power (59.3% vs. 59.3%) while controlling FDR—showing the nonparametric residual prior costs essentially nothing even when unnecessary. Real-data analyses reinforce the pattern:

• RNA-seq (CD4⁺ T cells, melanoma): untrended methods make zero discoveries; trended methods make roughly 76–133.
• ChIP-seq (H3K4me3, lymphoblastoid cells): Untrended-Inv$i$2 yields 197 discoveries versus 4,620+ for trended methods, reflecting clear parametric misspecification.
• Proteomics (A431 miRNA mimics): discrete JT-NPMLE attains the most discoveries across all six contrasts. Here the well-specified model visibly fails—low-peptide-count proteins show greater residual variability—and untrended methods concentrate discoveries disproportionately among low-$i$3 proteins, exactly the failure mode DEqMS was designed to prevent.

Limitations and open questions

Several caveats are stated plainly in the paper. The theory ignores discretization error and data-driven support selection in the NPMLE sieve approximations. The criticality condition requires dense signals; behavior under sparse signal regimes is not established. The Gaussian model is applied to log-transformed count data out of convention, and the authors explicitly defer studying the impact of this approximation to future work. Under trend misspecification, calibration conditional on $i$4 is not guaranteed, so stratum-level type-I error claims are unavailable. Finally, the frequentist (compound decision) analysis shows FDR control holds at level $i$5 rather than $i$6, precluding null-proportion-adaptive procedures such as Storey's method, and p-values satisfy only asymptotic compound validity rather than conditional uniformity.

Conclusion

This paper provides the first rigorous statistical foundation for limma-trend, recasting it as approximate partially Bayes inference with a trend-dependent nuisance prior. The framework yields practical dividends: a checkable orthogonality condition that diagnoses MAnorm2's FDR inflation, two NPMLE-based procedures with asymptotic FDR guarantees robust to trend misspecification, and near-parametric convergence rates enabled by Eddington–Tweedie identities. The empirical results indicate the nonparametric priors deliver power comparable to aggressive plug-in methods without sacrificing error control, positioning them as principled replacements for DEqMS-style heuristics in proteomics and for parametric limma-trend more broadly.