Serial Correlation In Within-Subject Residuals

Why Within-Subject Residual Correlation Matters

A pharmacokinetic model predicts a concentration for each observation, conditional on the subject’s covariates and random effects. Residuals are the differences between observed concentrations and those predictions. A common starting assumption is that residuals from different observations are conditionally independent. That assumption can be inadequate when nearby observations from the same subject tend to lie on the same side of the model prediction: an assay batch, a short-lived physiological fluctuation, or a local structural mismatch can make successive deviations related.

Serial correlation represents dependence among a subject’s residuals ordered over time, after conditioning on the model and the subject’s random effects. It is not a substitute for specifying the pharmacokinetic mean model, and it does not automatically explain why the residuals are related. Its role is to describe a pattern in their conditional covariance. This distinction matters because unmodelled serial dependence can make a model treat repeated evidence as more independent than it really is.

The qualifier conditional is essential. Even if residuals are conditionally independent, two concentrations from one subject can be associated marginally because they share random effects. Conversely, a serial residual covariance describes dependence that remains after conditioning on those effects. The earlier material on random effects and residual error provides the distinction between these sources of variation; here the focus is how to express the residual dependence itself.

Covariance Models for Repeated Concentrations

Suppose the conditional residual variance at observation \(j\) for subject \(i\) is \(\sigma_{ij}^{2}\). A general covariance specification is

where \(R_i\) is a correlation matrix: it is symmetric, has ones on the diagonal, and is positive semidefinite. In matrix notation, with \(D_i=\operatorname{diag}(\sigma_{i1},\ldots,\sigma_{in_i})\), the within-subject covariance matrix is \(\Sigma_i=D_iR_iD_i\). This form allows residual scale to vary with predicted concentration while expressing serial dependence separately through \(R_i\). For constant residual standard deviation \(\sigma\), it reduces to \(\Sigma_i=\sigma^2R_i\).

The AR(1) index is the observation order, not automatically elapsed time. It is most directly interpretable when observations are equally spaced. If gaps vary, using \(\rho^{|j-k|}\) makes correlation depend on the number of sampling steps rather than their actual duration. The exponential model instead uses elapsed time. For \(\lambda>0\), its correlation decreases as the time separation increases; at \(\lambda=0\), it is identically one for every separation, the limiting degenerate case rather than a decreasing-correlation model.

Both models assume a particular covariance pattern. The AR(1) form permits negative adjacent correlation when \(\rho<0\), with the sign alternating at successive lags. The exponential form has nonnegative correlation and a monotone decay when \(\lambda>0\). Neither form asserts that concentrations themselves follow an autoregression: the correlation applies to residuals after subtracting the pharmacokinetic predictions.

Why the AR(1) Covariance Has This Form

One way to justify the equally spaced AR(1) covariance is by a stationary Gaussian recursion. Let \(u_{ij}\) be independent standard normal innovations, and let the initial value \(z_{i1}\) be standard normal and independent of all innovations. For \(-1<\rho<1\), define

The independence condition on the initial value is part of the construction. Since \(z_{i1}\) has mean zero and variance one, each subsequent value has mean zero and variance

For \(k<j\), the innovation \(u_{ij}\) is independent of \(z_{ik}\), so the recursion gives \(\operatorname{Cov}(z_{ij},z_{ik})=\rho\operatorname{Cov}(z_{i,j-1},z_{ik})\). Applying this relation repeatedly yields \(\operatorname{Cov}(z_{ij},z_{ik})=\rho^{j-k}\). Because both variances are one, this is also the correlation. If the initial value were merely specified to have mean zero and variance one but not to be independent of the innovations, the covariance calculation would not follow; if it were not normal, the sequence would not generally be Gaussian. Those conditions should not be left implicit when using this recursive construction.

For irregular gaps, the exponential covariance has a related conditional transition representation. If adjacent times are separated by \(\Delta_{ij}=t_{ij}-t_{i,j-1}>0\), set \(a_{ij}=\exp(-\lambda\Delta_{ij})\), and use \(z_{ij}=a_{ij}z_{i,j-1}+\sqrt{1-a_{ij}^{2}}u_{ij}\), with a standard normal initial value independent of the independent standard normal innovations. This preserves unit variance. Multiplying the transition coefficients across intervening gaps gives correlation \(\exp[-\lambda(t_{ij}-t_{ik})]\) for \(j>k\).

Implications Compared With Independent Residuals

The uncorrelated-residual case is \(R_i=I\), the identity matrix. Under a conditionally multivariate normal residual model, this also implies conditional independence; otherwise, independence requires a separate assumption. Under a serial model, a run of positive residuals is less surprising than it would be under independence when \(\rho>0\) or when the exponential correlation is appreciable. The covariance also changes the joint contribution of a subject’s observations to the likelihood. Under a multivariate normal residual model, the conditional log-likelihood contribution, apart from constants, is

where \(\mathbf{e}_i\) is the vector of that subject’s residuals. Independence uses a diagonal \(\Sigma_i\); serial correlation uses a structured non-diagonal matrix. It is therefore not generally correct to fit each observation as if it were independent and then adjust only a standard error afterward: the covariance structure changes how the joint residual pattern is evaluated.

Assumptions, Diagnostics, and Limitations

A Gaussian serial residual model commonly assumes conditional mean zero, the specified covariance matrix, and multivariate normality of each subject’s residual vector conditional on its random effects and covariates. It also usually assumes residual vectors are independent across subjects. The covariance matrix must be valid and positive semidefinite; the AR(1) condition \(-1<\rho<1\) ensures the standard stationary model is nondegenerate, and the exponential form with \(\lambda\geq0\) is a valid covariance model. If residual scale depends on predicted concentration, the chosen \(\sigma_{ij}\) must also be appropriate; correlation modeling does not repair a misspecified variance function or mean model.

Diagnostics should examine residuals in time order within subjects, not only pooled residuals. Useful warning signs include long runs of residuals with the same sign, a residual autocorrelation pattern that persists after fitting the proposed covariance, and a fitted correlation estimate at or near a boundary. Check that the assumed sampling-time structure is sensible: a lag-one model can be misleading for highly irregular gaps, while a single exponential decay may fail if correlation does not decrease smoothly with time. Also inspect whether apparent residual runs are better explained by a systematic concentration-time misspecification or an omitted subject effect. A covariance model can absorb patterns statistically without identifying their biological source.

Serial correlation is most distinguishable when subjects contribute several appropriately timed observations, and when sampling schedules provide information across multiple time separations. With few observations per subject, the data may offer little information to separate correlation from residual scale or other model components. Replicate concentrations at the same time can help characterize residual variation, as discussed in the earlier material on replicate concentrations, but same-time replicates alone do not establish how correlation decays across time.

The principal alternative is conditionally independent residual error, which is simpler and may be adequate when residual dependence is negligible. A structured serial covariance is more economical than estimating an unrestricted covariance for every pair of sampling times, but its parsimony comes from assumptions such as constant lag decay or exponential time decay. Neither alternative should be selected merely because it produces a preferred result. Compare whether the covariance assumptions match the design, inspect conditional residual patterns, and evaluate model fit using methods appropriate to the fitted likelihood and model complexity. A serial residual model also does not automatically account for outliers or changing residual variance; those are distinct features requiring their own assessment.