Marginal Residuals For Longitudinal Submodels

Why Examine Marginal Residuals?

A longitudinal model can fit each subject's measurements reasonably well after accounting for that subject's estimated random effects and still describe the population mean trajectory poorly. Conditional residuals compare observations with subject-specific fitted values. Marginal residuals instead compare observations with the model's mean after averaging over the random-effects distribution. That difference makes marginal residuals useful for asking whether the model's population-level trajectory is adequate.

The distinction is important in joint models. The longitudinal submodel describes repeated marker values, while shared random effects may also connect the marker process to the event process. Here the diagnostic target is the population-level mean in the longitudinal submodel, obtained by integrating over the random effects. It is not a subject-specific fit, and it is not automatically the mean among subjects who remain event-free at a particular time.

Definition and Gaussian Model

For subject $i$ at measurement occasion $j$, consider a Gaussian longitudinal submodel

Y_{ij}=x_{ij}^{T}\beta+z_{ij}^{T}b_i+\epsilon_{ij}, \qquad b_i\sim N(0,D), \qquad \epsilon_{ij}\sim N(0,\sigma^2).

Here $x_{ij}$ and $z_{ij}$ are known covariate vectors, $\beta$ is the fixed-effects parameter, and $b_i$ is a vector of subject-specific random effects. Assume the random effects are independent across subjects, the measurement errors are independent across occasions and subjects, and $b_i$ and the errors are independent. More generally, the within-subject measurement errors can have covariance matrix $R_i$; the independent-error model has $R_i=\sigma^2 I$.

Definition: The marginal fitted mean at occasion $j$ is the conditional mean averaged over the random-effects distribution, $$ \mu_{ij}^{\mathrm{marg}}=E(Y_{ij}\mid x_{ij},z_{ij}) =\int E(Y_{ij}\mid b_i,x_{ij},z_{ij})\,f(b_i)\,db_i =x_{ij}^{T}\beta. $$ The marginal residual is $r_{ij}^{\mathrm{marg}}=Y_{ij}-\mu_{ij}^{\mathrm{marg}}$. In an analysis, replace unknown parameters by their fitted estimates.

The integral uses the population distribution of $b_i$, whose mean is zero. It does not substitute an estimated subject effect such as an empirical Bayes estimate. In a nonlinear model, the marginal mean on the response scale is generally the integral of the inverse link over the random-effects distribution; it need not equal the inverse link evaluated at the mean linear predictor. The calculations below focus on the Gaussian marker model, for which the integral has the simple form shown above.

What Marginal Residuals Contain

A marginal residual includes both the subject's deviation from the population trajectory and measurement error. In the Gaussian model it can be written as

r_{ij}^{\mathrm{marg}}=z_{ij}^{T}b_i+\epsilon_{ij}.

Thus, a nonzero marginal residual is not by itself evidence that the model is wrong: subjects are expected to vary around the population mean. Integrating over random effects also induces dependence among residuals from the same subject. For two occasions $j$ and $k$,

\operatorname{Var}(r_{ij}^{\mathrm{marg}})=z_{ij}^{T}Dz_{ij}+\sigma^2, \qquad \operatorname{Cov}(r_{ij}^{\mathrm{marg}},r_{ik}^{\mathrm{marg}}) =z_{ij}^{T}Dz_{ik}\quad (j\ne k),

under independent measurement errors. In vector form, if $Z_i$ is the matrix of random-effects design rows, the marginal covariance of the residual vector is $Z_iDZ_i^T+R_i$. Consequently, marginal residuals should not be treated as independent observations merely because they have been standardized one at a time.

Worked Example: Computing a Marginal Residual

Suppose the fitted population trajectory is $\widehat{\mu}(t)=10+0.5t$, where time $t$ is measured in months and the marker is measured in units. For one subject, a measurement at month 2 is 12.4 units. The fitted random-intercept and random-slope variances are 4 and 0.09, respectively, with zero covariance; the measurement-error variance is 1.

The population-level fitted mean is $10+0.5(2)=11$ units. Therefore the marginal residual is

r^{\mathrm{marg}}=12.4-11=1.4\text{ units}.

At month 2, the random-effects design vector is $(1,2)^T$. The marginal residual variance is

(1,2) \begin{pmatrix}4&0\\0&0.09\end{pmatrix} \begin{pmatrix}1\\2\end{pmatrix}+1 =4+4(0.09)+1=5.36\text{ units}^2.

The marginal standard deviation is $\sqrt{5.36}=2.315$ units, so the residual standardized by its marginal standard deviation is $1.4/2.315=0.605$. This is a modest positive residual relative to the model's marginal spread. It does not establish a problem at this occasion, and it does not account for the dependence with this subject's other residuals.

Derivation: Marginal Mean and Dependence

The marginal formulas follow by taking expectations and covariances over both sources of variation. Since $E(b_i)=0$ and $E(\epsilon_{ij})=0$,

E(Y_{ij}\mid x_{ij},z_{ij}) =x_{ij}^{T}\beta+z_{ij}^{T}E(b_i)+E(\epsilon_{ij}) =x_{ij}^{T}\beta.

Subtracting this expectation leaves $z_{ij}^{T}b_i+\epsilon_{ij}$. Independence of $b_i$ and $\epsilon_{ij}$ gives the variance as the sum of their variances:

\operatorname{Var}(z_{ij}^{T}b_i+\epsilon_{ij}) =z_{ij}^{T}Dz_{ij}+\sigma^2.

For distinct occasions, the measurement errors have covariance zero, while both residuals include the same $b_i$. Therefore

\operatorname{Cov}(z_{ij}^{T}b_i+\epsilon_{ij},z_{ik}^{T}b_i+\epsilon_{ik}) =z_{ij}^{T}Dz_{ik}.

This derivation explains why marginal residuals are suitable for assessing the population mean but are not independent replicates of model error. Their dependence is part of the model, not automatically a failure of it.

Worked Example: A Pattern Across Measurement Times

Suppose a fitted linear population trajectory for a marker is $\widehat{\mu}(t)=5+t$, with time in years and marker values in units. Across a group of subjects, the observed average marker values at years 0, 1, and 2 are 5.1, 6.2, and 8.1 units. For this illustration, treat these averages as summaries of the observed data, not as exact population means.

Subtract the fitted population mean at each time:

\begin{array}{c|c|c|c} t\ (\text{years}) & \widehat{\mu}(t)\ (\text{units}) & \text{Observed average} & \text{Average marginal residual}\\ 0 & 5 & 5.1 & 5.1-5=0.1\\ 1 & 6 & 6.2 & 6.2-6=0.2\\ 2 & 7 & 8.1 & 8.1-7=1.1 \end{array}

The positive residual averages increase markedly by year 2. This pattern flags a possible mismatch between the fitted linear population trajectory and the observed mean pattern, such as insufficient curvature. It is not proof of misspecification: sampling variability, parameter estimation, changing composition of the observed subjects, or other features of the design could contribute. A useful diagnostic question is whether the pattern is systematic across subjects and supported at times with adequate data, rather than driven by a few observations.

Using the Residuals and Recognizing Failure Modes

A practical diagnostic begins by calculating $r_{ij}^{\mathrm{marg}}=y_{ij}-x_{ij}^{T}\widehat{\beta}$ for each observed measurement. Plot residuals against time and important covariates, and inspect summaries by time or covariate group. A persistent curve or directional trend can suggest that the population mean structure is too simple. A cloud centered near zero without a systematic pattern is more consistent with an adequate mean structure, although it does not certify the whole joint model.

A marginal standardization divides by the model-based marginal standard deviation, $\sqrt{z_{ij}^{T}\widehat{D}z_{ij}+\widehat{\sigma}^2}$, in the independent-error case. This places residuals on a comparable scale when their marginal variances differ. It does not remove within-subject correlation. Nor does an estimated standardized residual have exactly a standard normal distribution in finite samples: the parameters are estimated, and observations may influence their own fitted values.

Common pitfall: “The conditional residuals are small, so the population trajectory must fit.” This reasoning is wrong because conditional residuals subtract estimated subject effects, which can absorb subject-specific departures from the population curve. To assess the population-level mean, compare observations with the mean integrated over the random-effects distribution; then interpret the resulting residual patterns while accounting for their dependence.

In a joint model, event-related dropout deserves particular attention. The marginal mean defined here averages over the specified random-effects distribution for the longitudinal submodel. If subjects with certain random effects experience events earlier, the subjects still observed at later times may have a different random-effects composition. Their observed average can therefore depart from the unconditional marginal mean even when the fitted longitudinal submodel is appropriate. Distinguish this survivor-composition issue from a misspecified population trajectory, and be explicit about the target population and time range of the diagnostic.

How This Diagnostic Differs From Alternatives

Diagnostic quantity	Reference fit	Useful for	Main caution
Marginal residual	Mean integrated over random effects	Assessing population-level mean patterns	Contains between-subject variation and is correlated within subject
Conditional residual	Subject-specific fit using an estimated random effect	Assessing deviations after accounting for subject effects	Can conceal population mean misspecification
Subject-level residual summary	A summary across a subject's repeated residuals	Comparing residual behavior across subjects	Can hide occasion-specific patterns when aggregated

These approaches answer different questions rather than competing to identify one universally best residual. Marginal residuals are especially informative about the fixed population mean, while conditional residuals focus on fit after subject-specific effects. A marginal residual plot alone cannot establish correct random-effects distributions, residual covariance, event submodel, or joint association structure. It is one diagnostic component, to be considered alongside model assumptions and other checks.

Worked Example: Why Standardization Does Not Remove Dependence

Suppose a random-intercept model has random-intercept variance 1 unit squared and independent measurement-error variance 0.25 unit squared. The random-intercept design value is 1 at both occasions. For two measurements from one subject, the marginal variance at each occasion is $1+0.25=1.25$ unit squared, and their covariance is 1 unit squared.

The model-implied correlation is therefore

\operatorname{Corr}(r_{i1}^{\mathrm{marg}},r_{i2}^{\mathrm{marg}}) =\frac{1}{\sqrt{1.25(1.25)}}=\frac{1}{1.25}=0.800.

If the observed marginal residuals are 1 and $-1$ units, their individually standardized values are $1/\sqrt{1.25}=0.894$ and $-1/\sqrt{1.25}=-0.894$. The standardized values remain correlated under the model; dividing by each marginal standard deviation has not made them independent. The opposing signs in this particular pair may occur, but the pair must not be analyzed as two independent pieces of evidence.

Check Your Understanding

Answer each question, then compare with the solutions below.

In the Gaussian model, why is the marginal fitted mean $x_{ij}^{T}\beta$ rather than $x_{ij}^{T}\beta+z_{ij}^{T}\widehat b_i$?
A fitted marginal mean is 7.5 units and the observed measurement is 8.0 units. Calculate the marginal residual and interpret its sign.
With a random-intercept variance of 2 and independent measurement-error variance of 0.5, calculate the variance and covariance of marginal residuals at two distinct occasions.
Why can a small conditional residual fail to demonstrate that the population mean trajectory is adequate?
What extra caution applies when late measurements are observed only for subjects who have not yet experienced the event?

Solutions. (1) The marginal mean averages over $b_i$, and $E(b_i)=0$, so the random-effects contribution vanishes. The subject-specific term uses an estimated individual effect instead of integrating over its population distribution. (2) $8.0-7.5=0.5$ units; the observation is above the fitted population mean by 0.5 units. (3) Each marginal variance is $2+0.5=2.5$ units squared, and the covariance is 2 units squared because the same random intercept contributes to both occasions. (4) The estimated subject effect can absorb departures from the population trajectory, so conditional fit does not by itself assess the integrated mean. (5) Event-related dropout can change the random-effects composition of those still observed, making their average differ from the unconditional marginal mean even without longitudinal mean misspecification.