Heteroscedastic Biomarker Residuals

Why allow residual variances to differ?

A longitudinal biomarker may be more variable in one treatment group than another, or its short-term variability may change over follow-up. A model with one residual variance for every observation cannot represent either pattern. If that assumption is implausible, the fitted model may misrepresent the covariance of the data and the precision of estimated effects.

The important distinction is between residual variation and subject-to-subject variation. A subject-specific random intercept, for example, describes persistent differences among subjects. Residual variation describes departures from a subject’s modeled trajectory at particular observations. The earlier material on marginal covariance decomposition gives the corresponding covariance distinction; here the focus is on allowing the residual part to have unequal variances.

Definition: In a Gaussian longitudinal mixed model, let $Y_{ij}$ be the biomarker for subject $i$ at observation $j$, with fixed-effect covariate vector $x_{ij}$, random-effect covariate vector $z_{ij}$, fixed-effect vector $\beta$, and subject random effects $b_i$. Write $$ Y_{ij}=x_{ij}^{T}\beta+z_{ij}^{T}b_i+\epsilon_{ij}. $$ Assume $b_i\sim N(0,D)$, and that the residual vector $\epsilon_i$, conditional on $b_i$, has mean zero and covariance matrix $R_i$. Subjects are independent, and the random effects are independent of the residuals. A homoscedastic model uses the same residual variance $\sigma^2$ for every observation. A heteroscedastic model allows $\operatorname{Var}(\epsilon_{ij})$ to depend on a specified group or observation characteristic.

For a simple model with conditionally independent residuals, $R_i$ is diagonal. If residual correlation is also modeled, a useful representation is $$ R_i=A_i C_i A_i, $$ where $A_i$ is diagonal with entries $\sigma_{ij}=\sqrt{\operatorname{Var}(\epsilon_{ij})}$, and $C_i$ is a residual correlation matrix with diagonal entries 1. This separates the choice of residual variances from the choice of residual correlations. Whatever correlation structure is used, $R_i$ must be a valid covariance matrix: its variances must be positive and it must be positive definite.

Specifying the variance pattern

A variance model should reflect a pattern that is scientifically plausible and can be estimated from the available observations. Two common choices are group-specific variances and time-specific variances. The variances are positive parameters; in computation, they are often represented through their logarithms to ensure positivity.

Definition: With $G$ treatment groups, a group-specific residual variance model sets $\operatorname{Var}(\epsilon_{ij})=\sigma^2_{g(i)}$, where $g(i)$ identifies subject $i$'s group. With $K$ follow-up categories, a time-specific model sets $\operatorname{Var}(\epsilon_{ij})=\sigma^2_{k(j)}$, where $k(j)$ identifies the observation’s time category. A constant-variance model is the special case in which all these variances are equal.

For instance, a two-group model can estimate $\sigma^2_1$ and $\sigma^2_2$, while a model with early, middle, and late visits can estimate $\sigma^2_1,\sigma^2_2,\sigma^2_3$. If follow-up time is continuous, one may instead specify a variance function, such as $$ \operatorname{Var}(\epsilon_{ij})=\sigma^2\exp(2\gamma t_{ij}), $$ where $t_{ij}$ is time and $\gamma$ controls the change in residual log standard deviation per unit of time. This assumes a smooth exponential pattern rather than a separate free variance at every visit. It is not automatically preferable: its form should be defensible and supported by the observation pattern.

Worked Example: Comparing Residual Variability by Treatment

Suppose a hypothetical biomarker model has fitted the mean trajectory and subject effects for two treatment groups. For illustration, the resulting conditional residuals, in biomarker units, are $-2,-1,1,2$ in Group A and $-6,-3,3,6$ in Group B. Treat these residuals as a small descriptive illustration, not as a full variance-component analysis.

The residual means in both groups are zero. The sums of squared residuals are $$ \sum e_{Aj}^{2}=(-2)^2+(-1)^2+1^2+2^2=10, \qquad \sum e_{Bj}^{2}=(-6)^2+(-3)^2+3^2+6^2=90. $$

Using the within-group maximum-likelihood variance calculation for this simple illustration, divide each sum by its group’s four observations: $$ \widehat{\sigma}^{2}_{A}=10/4=2.5, \qquad \widehat{\sigma}^{2}_{B}=90/4=22.5. $$ The estimated variance ratio is $22.5/2.5=9$, and the estimated residual standard deviations are $\sqrt{2.5}=1.581$ and $\sqrt{22.5}=4.743$ biomarker units. Thus the estimated residual standard deviation in Group B is three times that in Group A, since $\sqrt{9}=3$. These calculations illustrate interpretation; in an actual mixed model, variance estimates account for fitted fixed effects, random effects, and the chosen residual correlation structure.

Worked Example: Residual Variance Across Follow-Up

Consider hypothetical conditional residuals at two visits, after accounting for the fitted trajectory and subject effects. At the early visit they are $-1,1,-2,2$; at the late visit they are $-2,2,-4,4$. The residual sums of squares are $$ SS_{\text{early}}=1+1+4+4=10, \qquad SS_{\text{late}}=4+4+16+16=40. $$

Dividing by four observations at each visit gives illustrative variance estimates of $10/4=2.5$ and $40/4=10$ biomarker-units squared. The estimated late-to-early variance ratio is $10/2.5=4$; the corresponding standard-deviation ratio is $\sqrt{4}=2$. The late-visit residual spread is therefore estimated to be twice the early-visit spread. This pattern motivates considering a time-dependent variance, but these small calculations alone do not establish that such a model fits better: repeated observations, uncertainty in the fitted effects, and residual correlation must be handled in the full model.

What changes in estimation?

Under the model, the marginal covariance matrix for subject $i$'s observations is $$ V_i=Z_i D Z_i^T+R_i, $$ where $Z_i$ stacks the random-effect covariates. Changing the residual variance structure changes $R_i$, and therefore changes the covariance used to estimate fixed effects and variance parameters. It does not, by itself, change the fixed-effect mean formula. Because observations with different modeled variances have different precision, the model’s likelihood does not treat their residual deviations as equally informative.

The role of variance in the likelihood can be seen in a simpler conditional setting. Suppose independent Gaussian residuals $e_j$ have specified variances $v_j>0$. Their log likelihood, apart from no omitted data-dependent terms, is $$ \ell=-\frac{1}{2}\sum_j\left\{\log(2\pi v_j)+\frac{e_j^2}{v_j}\right\}. $$ The squared residual is scaled by its variance, and the log-variance term also contributes. Thus one cannot obtain a valid heteroscedastic likelihood merely by dividing residuals by their standard deviations and then ignoring the variance parameters.

Justification: For independent $e_j\sim N(0,v_j)$, the density of one residual is $(2\pi v_j)^{-1/2}\exp[-e_j^2/(2v_j)]$. Taking its logarithm gives $-\frac12\{\log(2\pi v_j)+e_j^2/v_j\}$. Independence makes the joint log likelihood the sum of these terms. The factor $1/v_j$ shows how residual scale affects the quadratic contribution, while $\log(v_j)$ accounts for the changing density scale.

In mixed models, covariance parameters are usually estimated jointly or iteratively with fixed effects under the selected likelihood or restricted likelihood. When comparing residual variance structures by restricted likelihood, the fixed-effects design must be the same. A likelihood comparison is evidence about relative fit under the candidate models, not proof that one variance pattern is true. Model complexity, stability of estimates, and whether the proposed grouping or time pattern makes sense should also matter.

Worked Example: Comparing Constant and Group-Specific Variances

Suppose a hypothetical two-treatment mixed-model analysis uses the same fixed-effects design and random-effects structure in both fits. The constant-residual-variance model has restricted deviance $-2\ell_R=214.6$; the model with separate residual variances by treatment has $-2\ell_R=207.8$. The latter adds one variance parameter.

The likelihood-ratio statistic is $$ \Lambda=214.6-207.8=6.8. $$ Under the usual regularity conditions for comparing these nested covariance models, the null hypothesis of equal positive variances gives an approximate $\chi^2$ reference distribution with $df=1$. The resulting upper-tail p-value is approximately $0.0091$. At significance level $\alpha=0.05$, this is evidence against the constant-variance restriction in favor of the group-specific model, conditional on the rest of the model being correctly specified.

The comparison does not say that the probability the null hypothesis is true is $0.0091$. Nor does it quantify the practical importance of the difference: the estimated variance ratio and its uncertainty are needed for that. The approximation also depends on regularity and a valid nested comparison; if the fitting method or model setup does not support that reference distribution, the reported test should not be used mechanically.

When is constant variance implausible?

A constant residual variance is questionable when residual spread shows a reproducible pattern after accounting for the fitted mean and subject-level structure. Examples include systematically wider residual scatter in one treatment group, increasing or decreasing spread over follow-up, or a clear relationship between residual magnitude and fitted biomarker level. Such patterns should be judged against the observation design: a visit with few observations or an unusual mixture of subjects can produce an unstable apparent spread.

A practical assessment uses the scientific context and focused summaries or plots of residuals by group, visit, or fitted value. The next tutorial develops residual diagnostics in detail; here the key point is to look for structured differences in spread rather than treating every large residual as evidence of heteroscedasticity. Model-based standardized residuals divide residuals by their modeled standard deviations, but they are not automatically independent or identically distributed when the model includes random effects or residual correlation.

Common pitfalls and limitations

Confusing residual variation with random-effect variation. A group with more variable subject trajectories may need a different random-effects covariance, not necessarily a larger within-subject residual variance. The two sources enter different parts of $V_i$.
Interpreting a variance ratio as a standard-deviation ratio. If one residual variance is four times another, its standard deviation is twice as large, not four times as large.
Fitting too many separate variances. Estimating one variance per visit can be unstable when some visits have few observations or substantial missingness. A parsimonious variance function or fewer scientifically meaningful categories may be more reliable.
Changing the variance structure while ignoring correlation. Unequal variances and within-subject correlation are distinct features, but both shape $R_i$. A misspecified correlation pattern can complicate interpretation of estimated residual spreads.
Treating model selection as certainty. A small p-value or improved likelihood does not establish that the chosen structure is the true data-generating mechanism. Estimates can be imprecise, and conclusions may depend on the candidate models.

Heteroscedastic models also have limits. They cannot repair a misspecified mean trajectory, unmodeled dependence, or informative missingness. Group- or time-specific variance estimates describe conditional residual variation under the chosen model, not every source of observed biomarker variability. A useful analysis therefore selects a modest, interpretable variance structure and checks whether its estimates are supported by the data.

Check Your Understanding

Answer each question, then compare with the solutions below.

In a model with residual variances $4$ and $25$ biomarker-units squared for two groups, what is the variance ratio and what is the ratio of residual standard deviations?
Why does a treatment-specific residual variance not mean that the treatment necessarily changes the average biomarker trajectory?
For independent Gaussian residuals with variance $v_j$, identify the two contributions of $v_j$ to the log likelihood.
What covariance component changes when residual variances vary by visit: $Z_iDZ_i^T$, $R_i$, or both?
Give one reason that separate residual variances at every visit may be a poor choice.

Solutions. (1) The variance ratio is $25/4=6.25$; the standard-deviation ratio is $\sqrt{6.25}=2.5$. (2) The residual variance specifies dispersion around the modeled mean; the mean trajectory is specified by the fixed effects. A variance difference alone is not a mean difference. (3) The log likelihood includes $-\frac12\log(v_j)$ and $-\frac12 e_j^2/v_j$, in addition to the constant $-\frac12\log(2\pi)$. (4) $R_i$ changes; the random-effects contribution $Z_iDZ_i^T$ is not directly changed by that residual variance specification. (5) With few observations at some visits, the individual variance estimates may be unstable or poorly supported.