Residuals Marginalized Over Random Effects

Why Marginalize Over Random Effects?

A generalized linear mixed model (GLMM) describes outcomes conditional on cluster-specific random effects. A residual calculated after substituting an estimated random effect asks whether an observation is surprising for that particular cluster, given its estimated latent effect. A different question is whether the model predicts outcomes well for the population of clusters represented by the random-effects distribution. To ask the second question, calculate outcome probabilities after averaging over that distribution.

This distinction matters when clusters are sparse. An estimated random effect for a cluster may be imprecise, especially when it has few observations. Conditional residuals can therefore depend strongly on noisy cluster-specific estimates. Marginal residuals avoid conditioning on those estimates, but they do not reveal every within-cluster discrepancy: integrating over random effects targets population-level fit.

Definition: For observation $i$ in cluster $j$, let $Y_{ij}$ have conditional distribution $P(Y_{ij}=y\mid b_j)$, where $b_j$ is the cluster random effect. The marginal outcome distribution averages over the assumed random-effects distribution: $$ P_M(Y_{ij}=y)=\int P(Y_{ij}=y\mid b)\,dF(b). $$ Here $F$ is the model-specified distribution of the random effects. The subscript $M$ denotes marginalization over random effects; this probability is not conditional on an estimated $b_j$.

For a common GLMM specification, the conditional mean is related to the linear predictor by $g(\mu_{ij}(b))=\mathbf{x}_{ij}^{T}\boldsymbol{\beta}+\mathbf{z}_{ij}^{T}b$, and the random effects have a specified distribution, often $b_j\sim N(0,D)$. The marginal mean is

\mu_{ij}^{M}=E_b\{E(Y_{ij}\mid b)\} =\int \mu_{ij}(b)\,dF(b).

In general, this is not the conditional mean evaluated at $b=0$, nor is it obtained simply by applying the inverse link to the fixed-effects linear predictor. Nonlinear links make the integration consequential. Numerical quadrature or simulation is often needed.

Marginal Residual Definitions

Once the marginal distribution is available, several residual summaries can be formed. A marginal response residual subtracts the marginal mean from the observed outcome. A marginal Pearson residual additionally scales by the marginal standard deviation. For a discrete outcome, a marginal randomized quantile residual uses the entire marginal cumulative distribution rather than only its mean and variance.

Definition: Let $\mu_{ij}^{M}=E_M(Y_{ij})$ and $V_{ij}^{M}=\operatorname{Var}_M(Y_{ij})$. The marginal response and Pearson residuals are $$ e_{ij}^{M}=y_{ij}-\mu_{ij}^{M}, \qquad r_{ij}^{M}=\frac{y_{ij}-\mu_{ij}^{M}}{\sqrt{V_{ij}^{M}}}. $$ For a discrete outcome, define $F_{ij}^{M}(y)=P_M(Y_{ij}\le y)$. Draw $U_{ij}$ uniformly on the interval $[F_{ij}^{M}(y_{ij}^{-}),F_{ij}^{M}(y_{ij})]$, where $F_{ij}^{M}(y_{ij}^{-})=P_M(Y_{ij}<y_{ij})$; for an integer-valued outcome, this lower endpoint is $F_{ij}^{M}(y_{ij}-1)$. The marginal randomized quantile residual is $q_{ij}^{M}=\Phi^{-1}(U_{ij})$, where $\Phi$ is the standard normal cumulative distribution function.

For counts with support beginning at zero, the lower endpoint is $F_{ij}^{M}(-1)=0$. For binary outcomes, the interval for an observed zero is $[0,P_M(Y=0)]$, while the interval for an observed one is $[P_M(Y=0),1]$. Randomization fills the gaps that would otherwise occur between discrete cumulative probabilities.

Under a correctly specified marginal distribution with known parameters, randomized probability integral transforms are uniform on $(0,1)$; applying $\Phi^{-1}$ therefore gives standard normal marginal residuals. This is a statement about each observation's marginal distribution. Outcomes in the same cluster can remain dependent because they share random effects, so the residuals need not be independent or jointly standard normal.

Worked Marginal Residual Calculations

Worked Example: A Binary Outcome Averaged Over a Random Intercept

Consider a hypothetical binary-outcome GLMM with a random intercept taking values $-1,0,1$ with probabilities $0.25,0.50,0.25$. The fixed linear predictor is zero, and the conditional success probability is the logistic function of the random intercept. For one observation, suppose $y=1$.

The conditional probabilities at the three support points are $\operatorname{logit}^{-1}(-1)=0.2689$, $\operatorname{logit}^{-1}(0)=0.5000$, and $\operatorname{logit}^{-1}(1)=0.7311$. Thus the marginal success probability is

$$ p^M=0.25(0.2689)+0.50(0.5000)+0.25(0.7311)=0.5000. $$

A Bernoulli variable with success probability $0.5000$ has marginal variance $p^M(1-p^M)=0.2500$. Therefore $e^M=1-0.5000=0.5000$, and $r^M=0.5000/\sqrt{0.2500}=1.000$. This residual says that the observed success is one marginal standard deviation above the marginal mean.

If instead one substituted the random intercept $b=1$, the conditional residual would be $1-0.7311=0.2689$, with conditional Pearson scaling $0.2689/\sqrt{0.7311(1-0.7311)}\approx0.607$. These numbers answer different questions: the first averages over the assumed cluster population; the second conditions on a particular random-effect value.

Worked Example: A Randomized Quantile Residual for a Count

Suppose a hypothetical Poisson random-intercept model is evaluated by a three-point numerical integration rule with weights $1/6, 2/3, 1/6$. At the three random-effect support points, the conditional probabilities of zero counts are $0.050000, 0.100000,$ and $0.260472$. For an observation with $y=0$, calculate its marginal zero probability and illustrate one randomized quantile residual.

The integrated probability of zero is the weighted average:

P_M(Y=0)=\frac{1}{6}(0.050000)+\frac{2}{3}(0.100000)+\frac{1}{6}(0.260472) =\frac{0.050000+0.400000+0.260472}{6} =0.118412.

Since the observation is zero, its marginal randomized interval is $[0,0.118412]$. Suppose the randomization produces $U=0.01635$, which lies in that interval. Then

q^M=\Phi^{-1}(0.01635)\approx -2.1357.

This realized residual is in the lower tail of the standard normal scale. The value depends on the random draw within the valid interval; a different draw for the same observation could produce a different residual. The integrated probability $0.118412$, however, is fixed by the model and the specified integration rule.

Worked Example: Marginal Variance for a Clustered Binomial Count

Suppose a hypothetical observation is a count of successes among four trials. Its cluster-specific success probability $p$ takes values $0.2,0.5,0.8$ with probabilities $0.25,0.50,0.25$, and $Y\mid p\sim\operatorname{Binomial}(4,p)$. If the observed count is $y=4$, calculate the marginal mean, variance, and Pearson residual.

The marginal mean is $E(Y)=4E(p)$. Here $E(p)=0.25(0.2)+0.50(0.5)+0.25(0.8)=0.5000$, so $\mu^M=2.000$. The variance decomposition gives

V^M=E\{\operatorname{Var}(Y\mid p)\}+\operatorname{Var}\{E(Y\mid p)\}.

First, $E\{\operatorname{Var}(Y\mid p)\}=4E\{p(1-p)\}=4[0.25(0.16)+0.50(0.25)+0.25(0.16)]=0.8200$. Next, $E(p^2)=0.25(0.04)+0.50(0.25)+0.25(0.64)=0.2950$, so $\operatorname{Var}(p)=0.2950-0.5000^2=0.04500$. Since $E(Y\mid p)=4p$, its variance is $16(0.04500)=0.7200$. Thus $V^M=0.8200+0.7200=1.540$, and

r^M=\frac{4-2.000}{\sqrt{1.540}}\approx1.612.

Using only the average success probability in a binomial variance would give $4(0.5)(0.5)=1.000$, omitting the $0.7200$ contribution from variation in $p$ across clusters. Marginalizing the distribution, not merely substituting a marginal mean into a conditional variance formula, is essential.

Why the Variance Decomposition Matters

The calculation in the binomial example follows the law of total variance. For any outcome with finite variance, conditional on a random effect $b$,

\operatorname{Var}_M(Y) =E_b\{\operatorname{Var}(Y\mid b)\} +\operatorname{Var}_b\{E(Y\mid b)\}.

To justify it, write $Y-E(Y)=[Y-E(Y\mid b)]+[E(Y\mid b)-E(Y)]$, square, and take expectations. The cross-product has expectation zero because $E[Y-E(Y\mid b)\mid b]=0$. The first squared term averages the conditional variance; the second is the variance of the conditional mean. Consequently, marginal Pearson residuals include both ordinary outcome variation and heterogeneity induced by random effects.

For a vector of outcomes in the same cluster, the same random effect also induces marginal association. Individual marginal residuals use each observation's integrated distribution; they do not remove that dependence. A check that groups observations by cluster or examines within-cluster residual patterns addresses different features of fit.

Assumptions, Diagnostics, and Limitations

Marginal residuals are meaningful only relative to a specified model. The conditional outcome family and link must be appropriate, the random-effects distribution and covariance structure must be adequate for the intended population, and the fixed-effects predictors and offsets must be correctly formed. The integration must also be accurate enough: sparse outcomes and nonlinear links can make a coarse quadrature approximation consequential. Estimated fixed effects and covariance parameters are normally substituted, so the displayed residuals also have parameter-estimation uncertainty that their simple formulas do not show.

When marginal residuals show systematic departures from their reference behavior, possible explanations include a misspecified mean structure, an inadequate random-effects distribution, omitted heterogeneity, or an incorrect conditional outcome distribution. A pattern is a diagnostic signal, not a unique identification of the cause. Residual magnitude should be considered alongside the outcome scale, sample size, and the uncertainty introduced by estimating model parameters.

Marginal residuals have two useful limitations to keep in view. First, integrating over random effects can conceal a cluster-specific lack of fit: positive and negative conditional discrepancies may average out at the population level. Second, marginal residuals within a cluster are dependent, so procedures that treat them as independent observations can give misleading uncertainty. For formal calibration, simulation from the fitted model can reproduce the clustering and parameter-fitting process; that is distinct from treating the residuals as independent standard normal data.

Common pitfall: Using an estimated cluster random effect to calculate a residual and then describing it as population-marginal is incorrect. That calculation is conditional on the estimated effect. For a marginal residual, integrate the conditional outcome probabilities over the random-effects distribution $F$, then calculate the residual from the resulting marginal distribution. Conversely, a marginal residual should not be claimed to diagnose every within-cluster pattern.

The main alternatives are complementary rather than interchangeable. Conditional residuals use estimated random effects and focus on fit within observed clusters, but can be unstable when cluster information is limited. Marginal residuals use the random-effects distribution and target population-level predictions, but can mask local cluster misfit. Simulation-based residual checks can reflect discrete outcomes and clustered dependence when simulation and refitting are set up appropriately, at the cost of additional computation and simulation variability. Selecting among them depends on whether the scientific diagnostic question concerns a new or average cluster, or the fit within particular observed clusters.

Check Your Understanding

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

For a binary observation with marginal success probability $0.30$ and observed value $y=0$, state the randomized quantile interval.
Why is the marginal variance of a clustered binomial count generally not obtained by inserting the marginal success probability into the conditional binomial variance?
A residual calculated using an estimated cluster random intercept is described as a population-level marginal residual. What is wrong, and what should be done instead?
Do correctly specified marginal randomized quantile residuals have to be independent within a cluster? Explain.

Solutions. (1) The zero outcome has interval $[0,P_M(Y=0)]=[0,0.70]$, because $P_M(Y=0)=1-0.30$. (2) Random effects make the conditional mean vary across clusters, adding $\operatorname{Var}_b\{E(Y\mid b)\}$ to the average conditional variance. (3) Substituting an estimated random intercept conditions on it. Instead, integrate the conditional outcome probabilities over the random-effects distribution before computing the residual. (4) No. Outcomes in one cluster share random effects, so their marginal residuals can remain dependent even when each residual has the correct marginal reference distribution.