Tutorials › Population Pharmacokinetics › Worked Residual Model Selection Across Concentrations

Random-Effect Structure and Error Models · Tutorial 164 of 599

Worked Residual Model Selection Across Concentrations

Learn to compare how competing residual-error models predict concentration variability across the measurement range and interpret the resulting model-selection evidence.

Intermediate 11 min read

What You'll Learn

  • Calculate predicted standard deviations under additive, proportional, and combined error models.
  • Estimate a single residual-scale parameter by maximum likelihood for additive and proportional models.
  • Compare fitted error models using the normal likelihood and AIC.
  • Use replicate residuals at different concentration levels to assess variance patterns.
  • Recognize when sparse concentration coverage or a poor mean model undermines error-model selection.

Why Concentration Range Matters

A residual-error model describes the remaining difference between an observed concentration and the concentration predicted by the structural model, after accounting for the modeled individual effects. Its choice matters because it determines how much variability the analysis expects at each predicted concentration. A constant standard deviation may be plausible near the lower end of the range but underestimate variability at high concentrations; a proportional standard deviation may do the reverse if measurement noise includes a substantial concentration-independent component.

The question is not simply whether residuals are “large” or “small.” It is whether their spread changes with the predicted concentration in a way that a candidate variance function can represent. The earlier material on replicate concentrations established why observations at different concentration levels help reveal residual variance. Here, the emphasis is on calculating and comparing the implications of competing models across that range.

Definition: Let \(y_i\) be an observed concentration, \(f_i\) the corresponding model-predicted concentration, and \(e_i=y_i-f_i\) the residual. Under a normal residual model, \(y_i\mid f_i\) has mean \(f_i\) and variance \(v_i\). An additive model has \(v_i=\sigma_a^2\); a proportional model has \(v_i=(\sigma_p f_i)^2\); and a combined model has \(v_i=\sigma_a^2+(\sigma_p f_i)^2\). The corresponding standard deviation is \(\sqrt{v_i}\). Here, \(\sigma_a\) has concentration units, while \(\sigma_p\) is dimensionless.

These are variance functions, not claims that every observed residual must increase smoothly with concentration. Under each model, residuals vary randomly around zero; the function specifies their expected variance at a given prediction. The combined model treats the additive and proportional contributions as independent variance components, so their variances, not their standard deviations, add.

Assumptions and Predicted Standard Deviations

For the normal likelihood calculations below, assume that, conditional on the model predictions and any modeled individual effects, residuals have mean zero, follow normal distributions with the specified variances, and are independent across observations. When comparing candidate error models, hold the structural mean model and the data fixed; otherwise a better likelihood could reflect a change in the mean model rather than the residual variance function. A candidate model must also assign positive variance to each observation.

Predicted standard deviations make the candidate models directly interpretable. For the proportional model, standard deviation is the same fraction of the prediction at every level. For the additive model, it is constant in concentration units. The combined model has a nonzero additive contribution at low concentrations and increasing proportional contribution as predictions rise.

Worked Example: Calculating Spread Across the Range

Suppose the additive standard deviation is \(0.40\) mg/L and the proportional standard deviation fraction is \(0.15\). Calculate predicted standard deviations at predictions of 1, 5, and 20 mg/L under all three models. These are hypothetical parameter values used to illustrate the variance functions, not estimates from a study.

At \(f=1\) mg/L, the additive model gives \(0.40\) mg/L and the proportional model gives \(0.15(1)=0.15\) mg/L. The combined model gives \(\sqrt{0.40^2+(0.15\times1)^2}=\sqrt{0.1825}=0.4272\) mg/L.

At \(f=5\) mg/L, the predictions are \(0.40\) mg/L for the additive model and \(0.15(5)=0.75\) mg/L for the proportional model. The combined standard deviation is \(\sqrt{0.40^2+0.75^2}=\sqrt{0.7225}=0.8500\) mg/L.

At \(f=20\) mg/L, the additive prediction remains \(0.40\) mg/L, while the proportional prediction is \(0.15(20)=3.00\) mg/L. The combined prediction is \(\sqrt{0.40^2+3.00^2}=\sqrt{9.16}=3.027\) mg/L.

The additive component is most consequential at low concentrations, while the proportional component dominates at high concentrations. These calculations describe the models’ implied spreads; they do not, by themselves, show which model fits data best.

Fitting and Comparing Candidate Models

For independent normal residuals with variances \(v_i\), the log likelihood, apart from any parameters in the structural model, is \(\ell=-\frac{1}{2}\sum_i[\log(2\pi)+\log(v_i)+e_i^2/v_i]\). The first variance-dependent term penalizes models that assign large variances, while the standardized squared residual term penalizes models that assign variances too small to explain the residuals. Both matter: comparing only residual sums of squares ignores how the candidate models scale those residuals.

A useful special case is estimating one common scale parameter. If a model has \(v_i=\theta^2 g_i^2\), where the \(g_i\) are known positive values, then its log likelihood depends on \(\theta\) through \(n\log(\theta^2)+\sum_i e_i^2/(\theta^2g_i^2)\). When \(\sum_i e_i^2/g_i^2>0\), differentiating with respect to \(\theta^2\), setting the derivative to zero, and solving gives the maximum-likelihood estimate below. If the sum is zero, all residuals are zero and the likelihood is unbounded as \(\theta^2\to0^+\), so no finite maximum-likelihood estimate exists over positive variances.

$$ \widehat{\theta^2} = \frac{1}{n}\sum_{i=1}^{n}\frac{e_i^2}{g_i^2}. $$

For additive error, \(g_i=1\) and \(\theta=\sigma_a\); for proportional error, \(g_i=f_i\) and \(\theta=\sigma_p\). This result assumes the predictions are fixed for the calculation and the normal variance model is appropriate. In a population analysis, structural and random-effect parameters may also be estimated, so the full likelihood must account for them.

Worked Example: Likelihood Comparison at Low and High Predictions

Consider four hypothetical observations with fixed predictions \(2,2,20,20\) mg/L and residuals \(-0.20,+0.20,-4.0,+4.0\) mg/L. Compare additive and proportional normal residual models. Both models have one residual-scale parameter; any parameters common to the structural model are omitted because they do not change the AIC difference.

For the additive model, the maximum-likelihood variance is the mean squared residual: \((2(0.20^2)+2(4.0^2))/4=(0.0800+32.00)/4=8.020\) \((\mathrm{mg/L})^2\). Thus \(\hat{\sigma}_a=\sqrt{8.020}=2.832\) mg/L. The sum of squared residuals divided by the fitted variance is \(32.08/8.02=4.000\), and the sum of log variances is \(4\log(8.02)\). Therefore, \(-2\ell=4\log(2\pi)+4\log(8.02)+4=19.679\), approximately.

For the proportional model, the scale estimate is \(\hat{\sigma}_p^2=[2(0.20/2)^2+2(4.0/20)^2]/4=[2(0.0100)+2(0.0400)]/4=0.1000/4=0.02500\). Thus \(\hat{\sigma}_p=\sqrt{0.02500}=0.1581\). The fitted variances are \(0.02500(2^2)=0.1000\) \((\mathrm{mg/L})^2\) at 2 mg/L and \(0.02500(20^2)=10.00\) \((\mathrm{mg/L})^2\) at 20 mg/L. The sum of log variances is \(2\log(0.10)+2\log(10)=0\). The standardized squared-residual sum is \(2(0.20^2/0.10)+2(4.0^2/10)=0.8000+3.200=4.000\). Hence \(-2\ell=4\log(2\pi)+4=11.352\).

AIC is \(-2\ell+2k\), where \(k\) is the number of estimated parameters. The two models have the same parameter count here, so their AIC difference based on the unrounded likelihoods is approximately \(8.328\), favoring the proportional model (subtracting the displayed rounded values gives \(19.679-11.352=8.327\)). The residual-to-prediction magnitudes are \(0.20/2=0.10\) at the low predictions and \(4.0/20=0.20\) at the high predictions; they are not equal. The likelihood favors the proportional model among these two candidates, but this small data set does not establish constant relative residual magnitudes or prove the proportional model is adequate.

The calculation illustrates why a model can be favored without matching every simple visual summary exactly. Model selection balances the fit of all observations with the variances the model assigns, rather than requiring identical residual ratios across a few concentration levels. It also illustrates why arithmetic should be checked: the additive sum of squares is \(32.08\), not \(32.16\), and the proportional scale estimate uses a sum of \(0.10\), not \(0.16\).

Using Replicates to Assess a Combined Variance Function

Replicate observations at a common prediction can provide information about the residual spread at that level, provided the replicates share the relevant modeled conditions. If residual variance changes with prediction, the combined model can be examined by relating variance to squared prediction. It has the linear form \(v_i=A+P f_i^2\), where \(A=\sigma_a^2\) and \(P=\sigma_p^2\). A meaningful combined variance function requires \(A\geq0\) and \(P\geq0\).

Worked Example: Estimating Combined Components From Replicate Spreads

Suppose replicate residuals at a predicted concentration of 2 mg/L are \(-0.6062,+0.6062,-0.6062,+0.6062\) mg/L. At a prediction of 10 mg/L, suppose they are \(-2.598,+2.598,-2.598,+2.598\) mg/L. The values are hypothetical and chosen to make the calculations transparent. Calculate the sample variance at each level, then use the two variance equations to estimate the combined-model components by moments.

At 2 mg/L, the residual mean is zero and the sum of squared deviations is approximately \(4(0.6062^2)=1.470\). The sample variance, using denominator \(4-1=3\), is \(1.470/3=0.4900\) \((\mathrm{mg/L})^2\). At 10 mg/L, the sum of squared deviations is approximately \(4(2.598^2)=27.00\), giving sample variance \(27.00/3=9.000\) \((\mathrm{mg/L})^2\).

Set these sample variances equal to the combined variance function at the two predictions: \(0.4900=A+4P\) and \(9.000=A+100P\). Subtracting the first equation from the second gives \(8.510=96P\), so \(P=8.510/96=0.08865\). Then \(A=0.4900-4(0.08865)=0.1354\). Both estimates are nonnegative, giving \(\hat{\sigma}_a=\sqrt{0.1354}=0.3680\) mg/L and \(\hat{\sigma}_p=\sqrt{0.08865}=0.2977\).

The implied standard deviations are \(\sqrt{0.1354+0.08865(4)}=0.7000\) mg/L at 2 mg/L and \(\sqrt{0.1354+0.08865(100)}=3.000\) mg/L at 10 mg/L. They match the two sample standard deviations by construction. With only two concentration levels, however, this moment calculation cannot test whether the combined function describes intermediate levels well. It is an illustration of the variance components, not compelling evidence that the combined model should be selected.

Diagnostics, Alternatives, and Limitations

After fitting a candidate model, inspect residuals on the scale implied by that model. For a normal variance function, standardized residuals are \(e_i/\sqrt{\hat v_i}\). If the variance function is adequate, their spread should not show a systematic trend with predicted concentration. Also inspect residuals against prediction and time, and examine replicate spreads where available. A trend with prediction can indicate a poor variance function, but can also result from a misspecified structural mean, unmodeled correlation, or influential observations. Residual plots are diagnostic evidence, not a substitute for a sound model.

The principal alternatives have different implications and trade-offs:

ModelStandard deviation at prediction \(f\)Useful whenPotential failure
Additive\(\sigma_a\)Absolute measurement variability is roughly stable.Can imply implausibly large relative error at low concentrations or understate high-level spread.
Proportional\(\sigma_p f\)Spread grows approximately in proportion to the prediction.Spread tends to zero as the prediction tends to zero, which may be unrealistic where an absolute error floor exists.
Combined\(\sqrt{\sigma_a^2+(\sigma_p f)^2}\)Both a low-level absolute component and a concentration-dependent component are plausible.Extra flexibility can be weakly identified without enough data across the range.

Use likelihood-based criteria such as AIC only for models fitted to the same observations with comparable likelihood definitions. The earlier material on comparing residual error models with replicate data motivates such comparisons; neither an AIC difference nor a visually appealing plot guarantees that the chosen model is scientifically adequate. In particular, a complicated variance function may fit a small sample better simply because it has more flexibility. Conversely, a small difference in AIC is not strong evidence for a practically important improvement.

Selection can fail when predictions cover only a narrow range, when there are few replicates, or when one region of the range has little information. The additive and proportional components can then trade off: a slightly larger constant component may be offset by a smaller proportional component. A poor mean model can also create structured residuals that look like changing variance. Finally, the normal independent-residual assumption may not hold. Dependence or heavy tails require different modeling choices; choosing a different variance function alone does not repair those problems.

Common Mistakes

  • Adding standard deviations in the combined model. The model adds variances: \(\sigma_a^2+(\sigma_p f)^2\). Take the square root only after adding them.
  • Choosing proportional error because high-level residuals are larger. Larger absolute residuals alone are not decisive. Check their scale relative to predictions and compare fitted models using the same data and mean model.
  • Interpreting a favored model as proven correct. AIC compares candidate models in the set examined; diagnostics and scientific plausibility remain necessary.
  • Using a negative variance component. A negative estimate of \(A\) or \(P\) cannot define a valid combined variance function. It signals that unconstrained moment estimates do not support that component; constrained fitting and broader diagnostics are needed.
Key takeaway: Residual-model selection across concentrations is a comparison of predicted variance functions. Calculate the implied standard deviations, fit candidates on the same data, inspect standardized residual patterns, and treat model-selection evidence as conditional on the mean model, assumptions, and concentration coverage.

Check Your Understanding

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

  1. At a prediction of 8 mg/L, what standard deviation does a proportional model with \(\sigma_p=0.12\) predict?
  2. A combined model has \(\sigma_a=0.30\) mg/L and \(\sigma_p=0.10\). Calculate its standard deviation at a prediction of 4 mg/L.
  3. Why is the combined variance \(\sigma_a^2+(\sigma_p f)^2\), rather than the sum of the two standard deviations?
  4. Two candidate models have the same number of fitted parameters. Their \(-2\ell\) values are 24.0 for additive and 20.5 for proportional error. Which has the lower AIC, and by how much?
  5. What does a systematic increase in standardized-residual spread with prediction suggest, and what other modeling problem could produce a similar pattern?

Solutions. 1. The predicted standard deviation is \(0.12(8)=0.960\) mg/L. 2. The variance is \(0.30^2+(0.10\times4)^2=0.0900+0.1600=0.2500\) \((\mathrm{mg/L})^2\), so the standard deviation is \(0.5000\) mg/L. 3. The combined model represents independent variance contributions, which add; the square root of their sum gives the total standard deviation. 4. With equal parameter counts, the AIC difference is the difference in \(-2\ell\): \(24.0-20.5=3.5\). The proportional model has AIC lower by 3.5. 5. The candidate variance function may be inadequate. A misspecified structural mean can also create a pattern that resembles changing residual spread, so the mean model should be checked as well.

Clinical Biostats