Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Wald and Likelihood-Ratio Tests for Pseudobulk Models

Differential Expression With Pseudobulk and Mixed Models · Tutorial 485 of 700

Wald and Likelihood-Ratio Tests for Pseudobulk Models

Compare Wald and likelihood-ratio tests for pseudobulk condition effects, work through their statistics and degrees of freedom, and identify when their conclusions may be unreliable.

Advanced 10 min read

What You'll Learn

  • State scalar and joint null hypotheses for condition effects in a pseudobulk model.
  • Calculate Wald statistics and likelihood-ratio statistics from hypothetical model outputs.
  • Determine degrees of freedom from the number of tested model parameters.
  • Explain the asymptotic assumptions behind each test and its p-value.
  • Diagnose warning signs such as poor convergence, influential samples, and unstable estimates.
  • Choose between a targeted contrast test and a nested-model comparison.

Why Two Tests for a Condition Effect?

For each gene, a pseudobulk analysis models counts aggregated by biological sample and cell type. The model estimates condition effects while accounting for the design and library-size offsets. Once a condition effect has been estimated, the inferential question is whether the data provide evidence against a specified null value for that effect. A Wald test and a likelihood-ratio test answer that question in different ways: the Wald test uses an estimate and its uncertainty, while the likelihood-ratio test compares the fit of nested models.

Both procedures rely on large-sample approximations, and neither removes the need to check whether the fitted model is appropriate. In particular, a small p-value does not establish that condition caused a change in expression. It indicates incompatibility between the observed data and a specified null model, under the assumptions used to calculate the test.

Definition: For gene \(g\) and pseudobulk sample \(j\), a negative-binomial generalized linear model may be written as \(Y_{gj}\sim NB(\mu_{gj},\phi_g)\), with \(\operatorname{Var}(Y_{gj})=\mu_{gj}+\phi_g\mu_{gj}^{2}\), and $$ \log(\mu_{gj})=o_{gj}+\mathbf{x}_j^{T}\boldsymbol{\beta}_g. $$ Here \(Y_{gj}\) is the aggregated count, \(\mu_{gj}\) its mean, \(\phi_g\) the dispersion, \(o_{gj}\) the known library-size offset on the log scale, \(\mathbf{x}_j\) the design vector, and \(\boldsymbol{\beta}_g\) the gene-specific coefficient vector. The tests below assess one or more restrictions on \(\boldsymbol{\beta}_g\).

The model assumes that, conditional on its mean and dispersion, each count follows the stated negative-binomial distribution. It also assumes that the biological samples included in the model are independent, unless their dependence is explicitly represented in the model. The design must identify the coefficients being tested: exact confounding or a rank-deficient design prevents some effects from being estimated. Dispersion may be estimated or moderated across genes, as in the earlier material on empirical Bayes shrinkage of gene dispersions. The test must use a likelihood and uncertainty calculation consistent with the dispersion treatment used by the fitting procedure.

Wald Tests: Estimate Divided by Uncertainty

A Wald test is suited to a specific coefficient or contrast. Suppose \(c^{T}\boldsymbol{\beta}_g\) is the log-scale condition contrast of interest, where \(c\) is a specified contrast vector. The null hypothesis is that this contrast equals a value \(\delta_0\), usually zero. If the fitted contrast estimate is \(\hat{\delta}_g=c^{T}\hat{\boldsymbol{\beta}}_g\), and its estimated standard error is \(SE(\hat{\delta}_g)\), the scalar Wald statistic is

$$ Z=\frac{\hat{\delta}_g-\delta_0}{SE(\hat{\delta}_g)}. $$

Under the null hypothesis and regularity conditions, this statistic is approximately standard normal. For a two-sided test, the p-value is \(2P(Z_{\mathrm{std}}\geq |z_{\mathrm{obs}}|)\), where \(Z_{\mathrm{std}}\sim N(0,1)\). Equivalently, \(Z^2\) is approximately chi-square with one degree of freedom. A Wald test is not restricted to one coefficient: for \(q\) restrictions collected in \(C\boldsymbol{\beta}_g=\mathbf{d}\), its quadratic form is approximately chi-square with \(q\) degrees of freedom.

Worked Example: A One-Degree-of-Freedom Wald Test

Suppose a hypothetical study has 12 independent donor-level pseudobulk samples from one cell type, with six donors in each of two conditions. A fitted negative-binomial model includes an intercept, a condition term, and the appropriate offsets. Assume the design is full rank, the model converges, and the diagnostics do not show a major mean-variance or sample-fit problem. The estimated treatment-versus-control log-fold change for a gene is \(0.50\), with standard error \(0.20\). Test \(H_0:\delta_g=0\) against a two-sided alternative.

The Wald statistic is $$ Z=\frac{0.50-0}{0.20}=2.50. $$

The two-sided standard-normal p-value is approximately \(2P(Z_{\mathrm{std}}\geq 2.50)=0.0124\). Thus, at \(\alpha=0.05\), this test provides evidence against a zero log-fold change for this gene. The estimated fold change is \(\exp(0.50)=1.649\), or about 1.65 times the modeled expression in the treatment condition relative to control, for the specified contrast and model. The p-value is not the probability that the null hypothesis is true.

For a \(q\)-dimensional restriction, let \(\hat{\boldsymbol{\delta}}\) be the estimated vector of restrictions and let \(\widehat{V}\) be its estimated covariance matrix. The multi-parameter Wald statistic is

$$ W=(\hat{\boldsymbol{\delta}}-\mathbf{d})^{T}\widehat{V}^{-1}(\hat{\boldsymbol{\delta}}-\mathbf{d})\ \overset{H_0}{\approx}\ \chi^2_q. $$

The degrees of freedom equal the number of independent restrictions, provided the covariance matrix is nonsingular and the restrictions are estimable. For example, a condition with three levels represented by two coefficients can be tested jointly by testing that both coefficients are zero; this is a two-degree-of-freedom test. A single contrast between two specified conditions is instead one degree of freedom.

Likelihood-Ratio Tests: Compare Nested Models

A likelihood-ratio test compares a full model, which includes the effects under investigation, with a reduced model that imposes the null hypothesis. The models must be nested: every parameter configuration allowed by the reduced model must also be allowed by the full model. Both models must use the same observations and compatible offsets and likelihood specification.

Definition: If \(\ell_F\) and \(\ell_R\) are the maximized log-likelihoods under the full and reduced models, respectively, the likelihood-ratio statistic is $$ D=2(\ell_F-\ell_R). $$ Under the null hypothesis, regularity conditions, and a correctly specified model, \(D\) is approximately chi-square with \(q\) degrees of freedom, where \(q\) is the difference in the number of freely estimated parameters in the full and reduced models.

The reduced model cannot fit better at its maximum likelihood than the full model when both are fitted to the same data under nested restrictions, so \(\ell_F\geq\ell_R\) and \(D\geq0\), apart from numerical tolerance. When dispersion is a nuisance parameter, the likelihood comparison must handle it consistently and account for its estimation according to the fitting procedure. A likelihood-ratio calculation from incompatible likelihoods—for example, fits using different observations—is not a valid nested-model test.

Worked Example: Testing a Three-Level Condition Jointly

Consider a hypothetical experiment with 18 independent donor-level pseudobulk samples and a condition variable with three levels. For one gene, the full model includes two condition coefficients, in addition to the intercept and offsets; the reduced model omits both condition coefficients. Suppose both models converge, use the same samples and offset specification, and the design is full rank. The maximized log-likelihoods are \(\ell_F=-120.4\) and \(\ell_R=-124.0\).

The null hypothesis is that both condition coefficients are zero. The full model adds \(q=2\) parameters, so $$ D=2\{(-120.4)-(-124.0)\}=2(3.6)=7.20. $$

The reference distribution is \(\chi^2_2\), and the upper-tail probability is \(P(\chi^2_2\geq7.20)=\exp(-7.20/2)=0.0273\). At \(\alpha=0.05\), the test provides evidence that at least one of the two condition coefficients differs from zero. It does not identify which condition level differs or give the direction of a particular contrast; those require examining suitable estimates or contrasts.

Why the Reference Distributions Arise

The Wald approximation follows from the asymptotic normality of the maximum-likelihood estimator. Under regularity conditions and the null hypothesis, the standardized estimation error for a scalar contrast is approximately \(N(0,1)\). Squaring a standard normal variable gives a chi-square variable with one degree of freedom. For \(q\) independent standardized restrictions, the sum of their squared components gives the chi-square reference with \(q\) degrees of freedom.

The likelihood-ratio approximation can also be understood from the local shape of the log-likelihood. Near its maximum, a regular log-likelihood is approximately quadratic. Under the null, the difference between the unconstrained maximum and the maximum subject to \(q\) restrictions is therefore approximately the squared length of \(q\) standard-normal components. Twice that log-likelihood difference has an approximate \(\chi^2_q\) distribution. This argument requires identifiable parameters, a sufficiently informative sample, an interior null value, and a model that is regular near the fitted values. It can be unreliable when these conditions fail.

Worked Example: Comparing One-Parameter Tests

For a hypothetical gene and a single treatment-versus-control coefficient, suppose the fitted estimate is \(0.60\) with standard error \(0.30\). A Wald test of \(H_0:\delta_g=0\) gives $$ Z=\frac{0.60}{0.30}=2.00, $$ with a two-sided p-value of approximately \(0.0455\). Now suppose the full-versus-reduced model likelihood comparison gives \(D=3.10\) for the same one-parameter null. Its p-value is \(P(\chi^2_1\geq3.10)\), approximately \(0.0782\).

These hypothetical results place the Wald p-value just below \(0.05\) and the likelihood-ratio p-value above it. They are not required to agree exactly: the Wald test uses a local estimate-and-standard-error approximation, while the likelihood-ratio test compares fitted likelihoods. A discrepancy near a decision threshold is a reason to check the model fits and the stability of the estimate, not to select whichever result is more favorable. The reported evidence should identify which test was used.

Choosing a Test and Checking Its Reliability

QuestionWald testLikelihood-ratio test
What is compared?Estimate of a specified coefficient or contrast with its standard errorMaximized likelihoods for nested full and reduced models
Natural useA targeted, often one-dimensional contrastA joint test of one or more terms added to the model
Degrees of freedomOne for a scalar contrast; \(q\) for \(q\) independent restrictionsNumber of parameters added in the full model, \(q\)
Key warningCan be unstable when estimates or standard errors are poorly behavedRequires nested models fitted to the same data with a compatible likelihood

For a single restriction, the Wald and likelihood-ratio tests often give similar results in well-behaved, adequately informative samples. They are not interchangeable formulas, and finite-sample results can differ. Use a Wald test when the scientific question is a particular contrast and its estimate and standard error are of direct interest. Use a likelihood-ratio test when the question concerns whether a set of terms improves fit as a group. A joint test is especially useful when no one coefficient alone represents the question.

Before interpreting either test, examine whether the fitting procedure converged, whether the design identifies the tested parameters, and whether there are influential samples or extreme coefficient estimates. Residual diagnostics for donor-aware models are relevant here: systematic residual patterns or a few samples dominating the fit weaken confidence in the model-based reference distribution. Also check that offsets and samples match between the full and reduced fits for an LRT. A singular covariance estimate can make a multi-parameter Wald statistic undefined or unstable.

Both tests can perform poorly with few independent biological samples, sparse counts, or parameters near a boundary. In such cases the asymptotic chi-square or normal approximation may be inaccurate even if the software returns a p-value. A failed or unstable fit is not evidence for the null; it is a warning that the proposed inference is not reliable. Conversely, a converged fit does not by itself verify the distributional assumptions or the independence of samples.

Finally, these are gene-level tests. A p-value for one gene does not account for having tested many genes; the multiplicity problem is a separate part of differential-expression inference. Report the tested contrast or restrictions, the test used, its degrees of freedom, the effect estimate where relevant, and the model context. A p-value alone cannot convey the magnitude or direction of a condition effect.

Key takeaway: A Wald test asks whether an estimated contrast is far from its null value relative to its standard error. A likelihood-ratio test asks whether the full model fits sufficiently better than its nested null model. Their degrees of freedom follow the number of independent restrictions, and their conclusions depend on model fit and asymptotic assumptions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

  1. A scalar condition contrast is estimated as \(-0.40\) with standard error \(0.20\). What is the Wald statistic, and what is the approximate two-sided p-value?
  2. A full model adds three estimable condition-related coefficients to a reduced model. What is the likelihood-ratio reference distribution, and what does rejection establish?
  3. Why is it invalid to compare full and reduced models fitted to different pseudobulk samples using the usual likelihood-ratio test?
  4. A joint Wald test involves two independent restrictions. What are its degrees of freedom, and what does rejection say about the individual coefficients?
  5. A Wald p-value is 0.04 and an LRT p-value is 0.08 for the same one-parameter null. What should be checked and reported?

Solutions. (1) \(Z=(-0.40)/0.20=-2.00\); the two-sided p-value is approximately \(0.0455\). (2) Use \(\chi^2_3\). Rejection is evidence that at least one of the three restrictions is false; it does not identify which coefficient differs. (3) The usual LRT relies on a nested comparison of likelihoods for the same observations and compatible model specification. Different samples do not satisfy that comparison. (4) The reference distribution is \(\chi^2_2\). Rejection indicates that the two restrictions are not jointly supported, but does not show that each coefficient is individually nonzero. (5) Check convergence, design identifiability, influential samples, model diagnostics, and whether the same data and likelihood assumptions were used. Report the chosen test and its result transparently; do not choose a test solely because its p-value crosses a threshold.

Clinical Biostats