Why Replicate Concentrations Help
Residual-error models describe how observed concentrations vary around the model-predicted concentrations. Choosing an unsuitable model can make residuals appear systematically too large in one concentration range and too small in another. Replicate concentrations are useful because repeated measurements at the same prediction level can reveal how within-level variability changes with concentration. The introductory material on using replicate concentrations to identify residual variance established this role; here, the focus is comparing candidate variance functions and their likelihood-based evidence.
Let \(C_{rj}\) be replicate \(j\) at prediction level \(r\), and let \(F_r\) be the model-predicted concentration for that level. For this comparison, assume the \(F_r\) values are treated as fixed, the residuals have conditional mean zero, and replicate residuals are independent and normally distributed. These assumptions make the likelihood calculations below explicit; in an application, dependence or a misspecified mean model can undermine the comparison.
In formulas, the candidate models are:
Here \(N(\mu,\sigma^2)\) uses variance as its second argument. The parameters \(\sigma_a^2\), \(\sigma_p^2\), and \(\sigma^2\) are nonnegative; the power parameter \(p\) controls how the standard deviation changes with prediction. In the proportional model, \(\sigma_p\) is dimensionless when \(C\) and \(F\) use the same concentration units. In the power model, the units of \(\sigma\) depend on \(p\), so that \(\sigma F^p\) has concentration units. The power model requires positive prediction levels for this interpretation and for the log calculations used to estimate \(p\).
An additive model has constant standard deviation across predictions. A proportional model has standard deviation \(\sigma_pF_r\), so doubling the prediction quadruples the variance. A combined model allows a nonzero variance floor as well as a concentration-dependent component. A power model permits an intermediate or steeper variance trend, depending on \(p\). These descriptions concern residual variation around \(F_r\); they do not describe variation in predictions between people.
Estimating Variance Patterns From Replicates
At a prediction level with \(n_r\) replicates, define residuals \(e_{rj}=C_{rj}-F_r\). If the prediction is the correct conditional mean and is treated as known, the empirical maximum-likelihood variance at that level is the residual sum of squares divided by the number of replicates:
This is not the usual sample variance with a denominator of \(n_r-1\). It is the maximum-likelihood estimate when the mean \(F_r\) is specified. If the mean is instead estimated from the same replicates, the variance calculation must account for that estimation. Residual plots and replicate summaries are informative only when the groups being compared genuinely share a prediction level and relevant conditions, such as assay and occasion.
Worked Example: Reading the Variance Pattern
Suppose four hypothetical replicate concentrations are observed at each of two fixed predictions. At \(F_1=5\) mg/L, the observations are 3, 7, 3, and 7 mg/L. At \(F_2=10\) mg/L, they are 7, 13, 7, and 13 mg/L. The predictions are assumed to be the correct conditional means, and the replicates are assumed independent within each level.
At 5 mg/L, the residuals are \(-2,2,-2,2\) mg/L. Thus,
At 10 mg/L, the residuals are \(-3,3,-3,3\) mg/L, giving
The estimated variance rises by a factor of \(9/4=2.25\), while the prediction doubles. A strictly additive model would predict no variance increase; a strictly proportional model would predict a variance increase by a factor of \(2^2=4\). The pattern is between those simple expectations. Two levels alone do not establish which model will generalize, but they suggest evaluating models that allow a nonconstant variance.
Likelihood Comparison and a Short Derivation
Under the independent normal assumptions, the log likelihood for a candidate variance function \(V_r\) is
The first term penalizes large variances through the log variance; the second penalizes residuals that are large relative to the proposed variance. Therefore, a model does not improve simply by inflating its variance everywhere. It must balance the two terms across prediction levels.
For a single shared variance \(V\), the log likelihood is \(-\frac12[n\log(2\pi V)+\mathrm{SSE}/V]\), where \(n\) is the total number of observations and \(\mathrm{SSE}\) is the total residual sum of squares. Differentiating with respect to \(V\), setting the derivative to zero, and solving gives the maximum-likelihood estimate:
When \(\mathrm{SSE}>0\), the second derivative is negative at \(\widehat V=\mathrm{SSE}/n\), so this stationary point is a maximum. If \(\mathrm{SSE}=0\), then \(\ell=-\frac n2\log(2\pi V)\) increases without bound as \(V\to0^+\), so there is no finite maximum over positive variances. For a variance function with multiple parameters, the same likelihood principle applies, but the parameters are generally estimated jointly by numerical optimization.
Worked Example: Fitting and Comparing Four Models
Use the hypothetical data above, with four replicates at each prediction, total \(n=8\), and group residual sums of squares 16 and 36 \((\text{mg/L})^2\). Assume the predictions are fixed and all eight residuals are independent normal observations with mean zero under each candidate model.
Additive. The shared variance estimate is \( (16+36)/8=6.5\;(\text{mg/L})^2\). The fitted variance at both levels is 6.5.
Proportional. Write the variance as \(qF_r^2\), where \(q=\sigma_p^2\). Maximizing the likelihood gives \(\hat q=\sum e_{rj}^2/\sum F_r^2=52/[4(5^2)+4(10^2)]=52/500=0.104\). The fitted variances are \(0.104(25)=2.6\) and \(0.104(100)=10.4\;(\text{mg/L})^2\).
Combined. The two empirical variances satisfy \(v_r=A+BF_r^2\) for \(A=\sigma_a^2\) and \(B=\sigma_p^2\). Solving \(A+25B=4\) and \(A+100B=9\) gives \(B=5/75=0.06667\) and \(A=4-25(0.06667)=2.333\). Both estimates are nonnegative, so this combined variance is valid and fits both empirical variances exactly.
Power. A power model has variance \(\sigma^2F_r^{2p}\). With two distinct positive prediction levels and strictly positive empirical variances \(v_1,v_2\), it can match both by solving \(v_r=\sigma^2F_r^{2p}\). Taking the ratio gives
Then \(\hat\sigma=\sqrt{4}/5^{0.5850}\approx0.780\), with units chosen so that \(\hat\sigma F^{\hat p}\) is in mg/L. This model also matches both empirical variances. The positivity condition matters: if either empirical variance is zero and the other is positive, no finite \(\sigma,p\) can match both at positive prediction levels.
To compare likelihoods, use \(-2\ell=8\log(2\pi)+\sum_r[n_r\log(V_r)+\mathrm{SSE}_r/V_r]\). The resulting values are approximately 37.68 for additive, 37.51 for proportional, and 37.04 for both combined and power. Thus the two two-parameter models attain the highest likelihood here. However, AIC adds a penalty for estimated parameters: \(\mathrm{AIC}=-2\ell+2k\), where \(k\) is the number of estimated parameters. Counting only residual-model parameters, the AIC values are 39.68 for additive, 39.51 for proportional, and 41.04 for each two-parameter model. The proportional model has the lowest AIC in this small example, but the differences are modest and are not proof that it is the true data-generating model.
Diagnostics, Assumptions, and Failure Modes
A candidate model should be checked against the pattern of residuals, not selected by a criterion alone. For each prediction level, compare the observed replicate spread with the model-implied standard deviation \(\sqrt{V_r}\). Residuals divided by \(\sqrt{V_r}\) should have roughly similar spread across prediction levels if the variance function is adequate. Also inspect residuals against predictions and time: a trend in their mean may indicate a structural-model problem, while systematic changes in spread may indicate a residual-variance problem.
For a numerical diagnostic, suppose another hypothetical data set has four replicates at each of \(F=2,5,10\) mg/L, with residuals alternating between \(-1\) and \(1\), \(-2\) and \(2\), and \(-3\) and \(3\) mg/L, respectively. Its within-level variances are 1, 4, and 9 \((\text{mg/L})^2\). Consider the candidate combined variance \(V(F)=1+0.08F^2\). It predicts variances 1.32, 3, and 9. The observed residual sums of squares divided by the corresponding model variances and replicate counts are \(4/(4\cdot1.32)=0.758\), \(16/(4\cdot3)=1.333\), and \(36/(4\cdot9)=1.000\). These values show no clear monotonic increase or decrease in scaled dispersion, although four replicates per level are too few to establish adequacy conclusively.
The likelihood comparison assumes the mean predictions are appropriate and that residuals are independent and normally distributed with the specified variances. Replicates can help distinguish residual variation from subject-level variability when the design and model support that separation, but they do not automatically remove confounding. If replicate measurements are correlated, treating them as independent overstates the information. Likewise, an omitted time trend or an incorrect structural prediction can be misread as a variance pattern. The earlier material on serial correlation and separating interindividual variability from residual error is relevant when these features are plausible.
Information criteria compare relative fit while penalizing complexity; they do not validate the assumptions or provide a universal test of which error model is correct. A likelihood-ratio comparison also requires care: candidate models must be nested for the usual likelihood-ratio framework, and boundary cases may alter its reference distribution. In practice, consider the variance pattern, residual diagnostics, numerical stability, scientific plausibility, and information criteria together.
Comparing the Models in Practice
| Model | Variance pattern | Useful indication | Important limitation |
|---|---|---|---|
| Additive | Constant across predictions | Spread is roughly constant in concentration units | May understate high-concentration variability |
| Proportional | Increases as \(F^2\) | Relative spread is roughly constant | Variance approaches zero as prediction approaches zero |
| Combined | Constant component plus \(F^2\) component | Both low-level and concentration-dependent spread are plausible | Extra parameter may be weakly identified without enough prediction range |
| Power | Proportional to \(F^{2p}\); increases for \(p>0\), is constant for \(p=0\), and decreases for \(p<0\) | A flexible monotone power trend is plausible | May be difficult to distinguish from combined variance with sparse levels |
With only two distinct positive prediction levels, ordered \(F_2>F_1\), and strictly positive empirical variances, a two-parameter power variance can interpolate those two variances. A combined variance can do so only if the fitted components are nonnegative: \(B=(v_2-v_1)/(F_2^2-F_1^2)\ge0\) and \(A=(F_2^2v_1-F_1^2v_2)/(F_2^2-F_1^2)\ge0\). This does not mean they make equivalent predictions between or beyond the observed levels. More prediction levels, adequate replication, and a broad prediction range help reveal whether their different variance curves matter. If an empirical variance is zero, it should not be treated as precise evidence for a zero population variance; with few replicates, zero can occur by chance.
Common Mistake
A common error is to choose the model with the largest likelihood and conclude that it is automatically best. For nested models, the maximized likelihood cannot decrease when the larger model's parameter space contains the smaller model's; candidate variance models that are not nested can have either a higher or lower maximized likelihood and should be compared using appropriate criteria and diagnostics. Another error is to call a larger variance at high concentrations “proportional error” without checking its rate of increase: proportional variance scales with \(F^2\), whereas a general power model scales with \(F^{2p}\). Finally, do not use the ordinary sample variance with denominator \(n_r-1\) as though it were the maximum-likelihood variance when the prediction is a fixed mean; that estimate is based on a different calculation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Under a proportional error model, what happens to the variance if the predicted concentration doubles?
- At a fixed prediction, four replicate residuals are \(1,-1,1,-1\) mg/L. What is the maximum-likelihood variance estimate?
- Two distinct positive prediction levels have strictly positive empirical variances. Can a two-parameter power model match both? What changes if one empirical variance is zero and the other is positive?
- Why might the model with the largest likelihood not have the lowest AIC?
- Give one reason why a changing residual spread across time might not be solved by changing only the residual variance function.
Solutions. (1) The variance is multiplied by \(2^2=4\). (2) The residual sum of squares is \(1+1+1+1=4\;(\text{mg/L})^2\), so the maximum-likelihood estimate is \(4/4=1\;(\text{mg/L})^2\). (3) Yes: solve the variance ratio for \(p\), then solve for \(\sigma\). If one variance is zero and the other positive, no finite power-model parameters can match both at positive prediction levels. (4) AIC penalizes additional estimated parameters using \(2k\), so a higher likelihood can be outweighed by the complexity penalty. (5) A time trend in mean residuals, serial correlation, or a misspecified structural prediction can create a pattern that a variance function alone cannot correct.