Serially Correlated Residuals in Dense Sampling

Why Model Residual Correlation?

In dense pharmacokinetic sampling, successive concentration measurements are often close in time. Even after the structural model and subject-specific parameters determine the predicted concentration profile, neighboring observations may tend to lie on the same side of that profile. This pattern can arise from short-term processes not represented in the mean model. If the remaining errors are correlated, treating each observation as independent gives the wrong joint probability for the subject’s concentration record.

Serial correlation is a statement about the residuals after conditioning on the predicted concentrations and subject effects. It does not replace the structural model, and it does not by itself describe differences between subjects. The model below is a discrete-time AR(1) covariance: residual correlation decreases geometrically with the number of sampling intervals between observations.

Definition: For subject \(i\), let \(y_{ij}\) be the observed concentration at sample \(j\), and let \(f_{ij}(\eta_i,\theta)\) be its model-predicted concentration, given subject random effects \(\eta_i\) and population parameters \(\theta\). The residual is \(e_{ij}=y_{ij}-f_{ij}(\eta_i,\theta)\). For \(n_i\) equally spaced samples, an AR(1) residual model assumes that, conditional on \(\eta_i\), the residual vector has mean zero and covariance \(R_i=D_i C_i D_i\), where \(D_i=\operatorname{diag}(s_{i1},\ldots,s_{in_i})\) contains residual standard deviations and \((C_i)_{jk}=\rho^{|j-k|}\). Here \(-1<\rho<1\), and \(\rho\) is the correlation between residuals at adjacent samples. Residual vectors from different subjects are assumed independent conditional on their random effects.

The standard deviation \(s_{ij}\) can be constant or can come from an already specified residual variance function. For constant residual standard deviation \(s\), \(D_i=sI\), so \(\operatorname{Cov}(e_{ij},e_{ik}\mid\eta_i)=s^2\rho^{|j-k|}\). Thus \(\rho=0\) gives independent residuals, while positive \(\rho\) means nearby residuals tend to have the same sign. Negative values are mathematically allowed, though the pattern they imply may be less plausible for concentration errors that vary smoothly over time.

The model requires a correctly specified conditional mean and residual variance, a multivariate Gaussian residual distribution, and a common correlation parameter for successive equally spaced samples under the stated AR(1) structure. The observations must be ordered by time. The lag \(|j-k|\) counts sample intervals, not elapsed time in an arbitrary unit; this model does not automatically accommodate unequal gaps by treating them as equal. These assumptions should be considered alongside the usual NLME assumptions for the subject effects.

Constructing the Covariance Matrix

Worked Example: Building a Three-Sample Covariance

Suppose a hypothetical subject has three equally spaced concentration samples. The predicted concentrations are \(10\), \(6\), and \(3\) mg/L, while the observations are \(11\), \(5\), and \(4\) mg/L. The residuals are therefore \(1\), \(-1\), and \(1\) mg/L. Assume a constant residual standard deviation \(s=2\) mg/L and adjacent residual correlation \(\rho=0.5\).

The correlation at lag 1 is \(0.5\), and at lag 2 it is \(\rho^2=0.25\). Multiplying the correlation matrix by \(s^2=4\;(\mathrm{mg/L})^2\) gives

R_i=4 \begin{pmatrix} 1 & 0.5 & 0.25\\ 0.5 & 1 & 0.5\\ 0.25 & 0.5 & 1 \end{pmatrix} = \begin{pmatrix} 4 & 2 & 1\\ 2 & 4 & 2\\ 1 & 2 & 4 \end{pmatrix} \;(\mathrm{mg/L})^2.

For example, the model covariance between the first and third residuals is \(4(0.25)=1\;(\mathrm{mg/L})^2\). The observed residuals are one possible realization; their particular signs do not define the covariance. The covariance describes the distribution of residual vectors the model assigns before observing them.

How the Covariance Changes the Likelihood

For a Gaussian residual vector \(e_i\) with mean zero and positive-definite covariance \(R_i\), the conditional log likelihood contribution from subject \(i\), given the subject effects, is

\ell_i = -\frac{1}{2} \left[ n_i\log(2\pi) +\log|R_i| +e_i^{\mathsf T}R_i^{-1}e_i \right].

The determinant term measures the volume of the covariance distribution, while the quadratic term measures how far the observed residual vector lies from zero in the covariance-scaled geometry. Under independent residuals, \(R_i\) is diagonal, so the quadratic term is a sum of individually scaled squared residuals. With serial correlation, off-diagonal covariance terms enter both the determinant and the quadratic form. Consequently, the joint likelihood is not the product of independent single-observation likelihoods.

For constant \(s\) and equally spaced AR(1) residuals, the determinant and quadratic form have convenient expressions:

|R_i|=s^{2n_i}(1-\rho^2)^{n_i-1}, \qquad e_i^{\mathsf T}R_i^{-1}e_i = \frac{e_{i1}^2}{s^2} + \sum_{j=2}^{n_i} \frac{(e_{ij}-\rho e_{i,j-1})^2}{s^2(1-\rho^2)}.

These expressions rely on the stationary AR(1) Gaussian model and its equally spaced observations. The quadratic form can be read as a sum of squared innovations: the first residual is scaled by \(s\), and each later residual is compared with the value predicted from the preceding residual, \(\rho e_{i,j-1}\).

Justification: Under the stationary AR(1) model, \(e_{i1}\sim N(0,s^2)\), and for \(j\geq2\), \(e_{ij}\mid e_{i,j-1}\sim N(\rho e_{i,j-1},s^2(1-\rho^2))\). The joint density factors into the initial density and these conditional densities. Its log density has quadratic terms \(e_{i1}^2/s^2\) and \((e_{ij}-\rho e_{i,j-1})^2/[s^2(1-\rho^2)]\), giving the stated quadratic form. The same factorization gives the determinant: the initial variance contributes \(s^2\), and each of the \(n_i-1\) conditional variances contributes \(s^2(1-\rho^2)\). Their product is \(s^{2n_i}(1-\rho^2)^{n_i-1}\).

Worked Example: Comparing Conditional Log Likelihoods

Consider three hypothetical equally spaced residuals \(e=(1,1,1)^{\mathsf T}\) mg/L. Compare the independent model with the AR(1) model at fixed values \(s=1\) mg/L and \(\rho=0.5\). This calculation compares likelihoods at the specified parameter values; it does not establish which model fits a data set better after parameter estimation.

Under independence, \(R=I\), so \(\log|R|=0\) and \(e^{\mathsf T}R^{-1}e=1^2+1^2+1^2=3\). Therefore,

\ell_{\mathrm{ind}} = -\frac{1}{2}\left[3\log(2\pi)+0+3\right] =-4.257.

For the AR(1) model, the determinant is \((1-\rho^2)^2=0.75^2=0.5625\), so \(\log|R|=\log(0.5625)=-0.5754\). The quadratic form is

e^{\mathsf T}R^{-1}e = 1^2+ \frac{(1-0.5\times1)^2+(1-0.5\times1)^2}{1-0.5^2} = 1+\frac{0.25+0.25}{0.75} =1.667.

Thus

\ell_{\mathrm{AR(1)}} = -\frac{1}{2}\left[3\log(2\pi)-0.5754+1.667\right] =-3.302.

At these fixed parameter values, the AR(1) log likelihood is higher by \(-3.302-(-4.257)=0.954\). A run of same-sign residuals is more compatible with positive serial correlation than with independence. In an actual model comparison, parameters would be estimated and the extra correlation parameter would need to be accounted for; this fixed-value calculation alone is not a model-selection conclusion.

Innovations and Practical Checks

The conditional formulation also suggests a way to inspect the assumed correlation. With constant \(s\), define standardized innovations by

u_{i1}=\frac{e_{i1}}{s}, \qquad u_{ij}=\frac{e_{ij}-\rho e_{i,j-1}}{s\sqrt{1-\rho^2}},\quad j\geq2.

If the conditional Gaussian AR(1) model is correct and its parameters are known, these innovations have a standard normal distribution and are mutually independent. In practice, parameters are estimated, so the innovations are not perfectly independent in a finite sample. Still, persistent patterns in them can indicate that the chosen correlation has not adequately captured the residual dependence.

Worked Example: Transforming a Residual Sequence

Suppose a hypothetical subject has residuals \(e=(1.0,1.2,0.9,0.8)\) mg/L from four equally spaced samples. Let \(s=0.5\) mg/L and \(\rho=0.6\). Since \(\sqrt{1-\rho^2}=\sqrt{0.64}=0.8\), first standardize the residuals: \(e/s=(2.0,2.4,1.8,1.6)\).

The first innovation is \(u_{i1}=2.0\). The others are

u_{i2}=\frac{2.4-0.6(2.0)}{0.8}=1.50,\qquad u_{i3}=\frac{1.8-0.6(2.4)}{0.8}=0.450,\qquad u_{i4}=\frac{1.6-0.6(1.8)}{0.8}=0.650.

The original residuals are all positive and fairly close together. The transformed values show what remains after accounting for the model’s expected dependence on the previous residual. They are not proof that the model is correct: across subjects and profiles, innovations should be checked for remaining time patterns, changing spread, and unusually influential observations.

Choosing and Interpreting the Structure

An independent residual model is the main alternative. It is simpler and is appropriate when conditional residuals are plausibly uncorrelated at the sampling resolution. The AR(1) model is more appropriate when adjacent residuals show dependence that decays with lag. A fully unstructured residual covariance could represent a different correlation at every pair of sample positions, but it uses many more covariance parameters and may be unstable with limited data. AR(1) imposes a parsimonious geometric pattern instead.

These choices concern residual dependence, not all dependence within an individual’s observations. Random effects can induce covariance in predicted concentrations across samples; the residual covariance describes remaining variation around those predictions. A positive residual correlation should therefore not automatically be interpreted as evidence for a particular random effect, nor should adding random effects be treated as a guaranteed substitute for an appropriate residual covariance. The components have different roles in the model.

Useful checks include plotting conditional residuals in time order, examining whether neighboring residuals tend to share a sign, and inspecting innovations after fitting the model. Problems are suggested when residual patterns persist after the AR(1) transformation, when the estimated \(\rho\) is close to its allowable boundary, or when estimates change substantially with small changes to the sampling subset or variance assumptions. A large estimated correlation can also be compensating for an inadequate mean trajectory or another omitted source of structure; correlation should not be used to conceal such misspecification.

Dense sampling does not itself guarantee that \(\rho\) is estimable. There must be enough informative within-subject residual sequences, and the fitted model must distinguish residual correlation from subject-level variation and other sources of dependence. Moreover, an AR(1) correlation is tied to a fixed sampling interval: its lag-one parameter has the same meaning for each adjacent pair in this model. Do not interpret it as a general correlation rule for arbitrary time gaps.

Common Pitfall

Common mistake: “The residuals are correlated, so I can keep the independent likelihood and simply report a positive correlation afterward.” This is incorrect because the joint density depends on both \(\log|R_i|\) and \(e_i^{\mathsf T}R_i^{-1}e_i\). The correction is to use the specified covariance inside the likelihood itself. A related mistake is to call \(\rho\) the correlation between all observations: it is the adjacent-sample correlation, and the correlation at lag \(k\) is \(\rho^k\) under this model.

Key takeaway: An AR(1) residual model replaces independent within-subject errors with a joint covariance whose off-diagonal entries follow \(\rho^{|j-k|}\). The determinant and covariance-scaled residual quadratic form both change, so the likelihood must be calculated jointly across a subject’s samples.

Check Your Understanding

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

For four equally spaced samples with \(s=3\) mg/L and \(\rho=0.4\), what is the covariance between the first and fourth residuals?
In the AR(1) model, what does \(\rho=0\) imply, and what does \(\rho\) represent when it is nonzero?
Why does a serially correlated residual model require changing the likelihood rather than only changing how residuals are described afterward?
A subject’s samples are separated by unequal time gaps. Why is it not automatically appropriate to apply the lag-index AR(1) formula as though each gap were identical?

Solutions. (1) The lag is 3, so the covariance is \(s^2\rho^3=9(0.4^3)=0.576\;(\mathrm{mg/L})^2\). (2) At \(\rho=0\), all distinct residuals are uncorrelated and, under the joint Gaussian assumption, independent. Otherwise, \(\rho\) is the correlation between adjacent equally spaced residuals; lag-\(k\) correlation is \(\rho^k\). (3) The joint density contains the covariance determinant and inverse, so both its normalization and residual quadratic form depend on correlation. (4) The formula defines correlation by the number of equal sampling intervals, not by elapsed time; unequal gaps do not share the same interval interpretation.