Modeling Assay-Specific Observation Error

Why Assay Error Belongs in the Calibration Model

A QSP model predicts quantities such as a biomarker concentration or a response level, while an assay reports an observed value. Those values need not agree exactly, even when the biological model is appropriate: measurement variability can be small in absolute terms, grow with signal magnitude, or act multiplicatively. The observation-error model states how the observed assay output varies around its model prediction. It therefore determines which discrepancies count as large during calibration.

The choice is consequential. A difference of 2 units is substantial when the prediction is 3, but modest when it is 100. An error model that treats those differences as equally surprising can distort the fit. In a multiple-assay calibration, the model can allow different error scales for different assays; this tutorial focuses on the form of error for continuous measurements, not on dependence among repeated observations.

Definition: Let \(y_i\) be a continuous observation from an assay and \(m_i=f_i(\theta)\) the corresponding QSP model prediction at parameter vector \(\theta\). The observation-error model specifies the conditional distribution of \(y_i\) given \(m_i\). The raw residual is \(r_i=y_i-m_i\). For assay \(a\), its error-scale parameter is denoted \(\sigma_a\).

The basic assumptions used below are that predictions are defined on the same scale and in the same units as the observations, the error model is appropriate for that assay, and observations are conditionally independent given the predictions and error parameters. Independence is a working assumption for these derivations, not a claim that all assay data are inherently independent. If an assay has a meaningful zero or can produce negative readings, that also affects which models are plausible.

Three Models and the Residuals They Imply

With additive normal error, the observation differs from the prediction by a normally distributed amount with constant standard deviation. The raw residual is the relevant error-scale quantity: observations with the same absolute discrepancy have the same variance, regardless of the prediction.

y_i\mid m_i \sim N(m_i,\sigma_a^2), \qquad r_i=y_i-m_i, \qquad z_i=\frac{r_i}{\sigma_a}.

Here \(z_i\) is the standardized residual. This model is often reasonable when assay variability is roughly constant in the reported units across the relevant measurement range. It permits negative observations, which may be acceptable for a background-corrected signal but is usually not appropriate for a strictly positive concentration.

With proportional normal error, the standard deviation grows in proportion to the prediction. A convenient form for positive predictions is:

y_i\mid m_i \sim N(m_i,\sigma_a^2 m_i^2), \qquad z_i=\frac{y_i-m_i}{\sigma_a m_i}.

In this parameterization, \(\sigma_a\) is a relative error scale: the standard deviation is \(\sigma_a m_i\), and the coefficient of variation is \(\sigma_a\). The relative residual \((y_i-m_i)/m_i\) expresses the discrepancy as a fraction of the prediction; dividing it by \(\sigma_a\) gives the standardized residual. This model can also generate negative observations, so it is an approximation rather than a natural positive-valued distribution when the coefficient of variation is not small.

With lognormal error, positive observations vary multiplicatively around a positive prediction. In the parameterization below, \(m_i\) is the conditional median, not the arithmetic mean:

\log(y_i)\mid m_i \sim N(\log(m_i),\sigma_a^2), \qquad u_i=\log\left(\frac{y_i}{m_i}\right), \qquad z_i=\frac{u_i}{\sigma_a}.

The log residual \(u_i\) is positive when the observation exceeds the prediction and negative when it falls below. Equal multiplicative discrepancies have equal absolute log residuals: for example, ratios of 1.2 and \(1/1.2\) give log residuals of equal magnitude and opposite sign. The lognormal model requires positive observations and positive predictions. Its conditional arithmetic mean is \(m_i\exp(\sigma_a^2/2)\), so if the intended prediction represents the arithmetic mean, the log-scale location must instead be \(\log(m_i)-\sigma_a^2/2\).

Likelihoods: Why the Error Model Changes Calibration

For independent normal errors with observation-specific variance \(v_i\), the negative log-likelihood, apart from no omitted data-dependent terms, is:

-\ell(\theta,\{v_i\}) = \frac{1}{2}\sum_i \left[ \log(2\pi v_i) + \frac{(y_i-m_i)^2}{v_i} \right].

This expression includes both a squared, variance-scaled residual and a term for the variance itself. Under additive error, \(v_i=\sigma_a^2\); under proportional normal error, \(v_i=\sigma_a^2m_i^2\). Thus, proportional error does not merely divide every raw residual by the same constant: its variance changes with the prediction. If the prediction depends on \(\theta\), the variance term also depends on \(\theta\).

For lognormal error with median prediction \(m_i\), the density is defined for \(y_i>0\). Its negative log-likelihood is:

-\ell(\theta,\sigma_a) = \sum_i \left[ \log(y_i)+\log(\sigma_a)+\frac{1}{2}\log(2\pi) + \frac{\{\log(y_i/m_i)\}^2}{2\sigma_a^2} \right].

The \(\log(y_i)\) term arises because the normal density is specified for \(\log(y_i)\) while the likelihood is for \(y_i\). When comparing candidate model predictions for the same positive observations, this term is constant in \(\theta\), but it remains part of the density. For any of these models, minimizing a weighted sum of squared residuals is equivalent to maximizing the likelihood only when the omitted terms are constant with respect to the parameters being estimated.

Justification for the additive-error scale estimate: Suppose the predictions are fixed and \(n\) independent residuals follow \(N(0,\sigma_a^2)\). The negative log-likelihood is \(n\log(\sigma_a)+\sum_i r_i^2/(2\sigma_a^2)\), up to a constant. Differentiating with respect to \(\sigma_a\) and setting the derivative to zero gives \(n/\sigma_a-\sum_i r_i^2/\sigma_a^3=0\). Therefore the maximum-likelihood estimate is \(\hat{\sigma}_a=\sqrt{\sum_i r_i^2/n}\). This uses divisor \(n\), not \(n-1\), because it is a likelihood estimate with the predictions treated as fixed.

Worked Residual Calculations

Worked Example: Constant Absolute Error

Suppose a hypothetical assay reports concentrations in ng/mL. For three samples, a fitted QSP model predicts 10, 20, and 40 ng/mL, while the observed values are 12, 17, and 44 ng/mL. Assume independent additive normal errors with a common standard deviation for this assay.

The raw residuals, observation minus prediction, are \(2\), \(-3\), and \(4\) ng/mL. Their squared sum is \(2^2+(-3)^2+4^2=29\) \((\mathrm{ng/mL})^2\). Treating these predictions as fixed, the maximum-likelihood error scale is:

\hat{\sigma}_a=\sqrt{\frac{29}{3}}=3.109\ \mathrm{ng/mL}.

Using this fitted scale, the standardized residuals are \(2/3.109=0.643\), \(-3/3.109=-0.965\), and \(4/3.109=1.286\). Under the stated model, none is more than 1.3 fitted standard deviations from its prediction. This small example cannot establish that the error is truly constant across the measurement range; it only shows the calculations under that assumption.

Worked Example: Error That Scales With the Prediction

A second hypothetical assay predicts 5, 20, and 50 response units, with observed values 6, 18, and 55 response units. Suppose the assay is modeled with proportional normal error, so its standard deviation at each prediction is \(\sigma_a m_i\).

The relative residuals are \((6-5)/5=0.20\), \((18-20)/20=-0.10\), and \((55-50)/50=0.10\). With predictions treated as fixed, the maximum-likelihood estimate of the relative error scale is the root mean square of these relative residuals:

\hat{\sigma}_a = \sqrt{\frac{0.20^2+(-0.10)^2+0.10^2}{3}} = \sqrt{0.0200} = 0.1414.

The fitted standard deviations are \(0.1414(5)=0.707\), \(0.1414(20)=2.828\), and \(0.1414(50)=7.071\) response units. The standardized residuals are therefore \(1.414\), \(-0.707\), and \(0.707\). The same absolute discrepancy would be judged differently at different predictions: a 1-unit discrepancy is large at a prediction of 5 but small at 50.

Worked Example: Multiplicative Error on the Log Scale

A hypothetical positive-valued assay has model predictions of 10, 12, and 20 mg/L and observations of 8, 12, and 18 mg/L. Assume lognormal error with each prediction representing the conditional median. The observations and predictions are all positive, as required.

The log residuals are \(\log(8/10)=-0.2231\), \(\log(12/12)=0\), and \(\log(18/20)=-0.1054\). The maximum-likelihood estimate of the log-scale standard deviation, with these predictions treated as fixed, is:

\hat{\sigma}_a = \sqrt{\frac{(-0.2231)^2+0^2+(-0.1054)^2}{3}} = \sqrt{0.02030} = 0.1425.

The standardized log residuals are approximately \(-1.566\), \(0\), and \(-0.740\). The first observation is 20% below its prediction, and the third is 10% below; their log residuals are not simply the raw differences in mg/L. Under this model, a prediction of 10 mg/L is a median. The corresponding arithmetic mean is \(10\exp(0.1425^2/2)=10.10\) mg/L, not exactly 10 mg/L.

Choosing and Checking an Error Model

The models encode different assumptions about how measurement variability behaves. The table summarizes their residual scales and practical restrictions.

Error model	Residual used for calibration	Assay variability assumption	Main restriction
Additive normal	\(y_i-m_i\), scaled by \(\sigma_a\)	Constant absolute standard deviation	Can predict negative observations
Proportional normal	\((y_i-m_i)/m_i\), scaled by \(\sigma_a\)	Standard deviation proportional to prediction	Requires positive predictions; can predict negative observations
Lognormal	\(\log(y_i/m_i)\), scaled by \(\sigma_a\)	Variation is multiplicative on the original scale	Requires positive observations and predictions

A sensible choice reflects assay behavior and the range of measurements, not merely which model produces the smallest calibration objective. Useful checks include plotting raw residuals against predictions, plotting standardized residuals against predictions, and examining residuals separately by assay. For additive error, the raw residual spread should be roughly stable over the range. For proportional error, relative residuals should have roughly stable spread. For lognormal error, log residuals should not show a systematic trend or a marked change in spread with fitted value.

A funnel-shaped raw-residual plot, with spread increasing as predictions rise, is evidence against a constant-variance additive model. Conversely, if relative residuals become extremely variable near zero predictions, a purely proportional model may be unsuitable. Large systematic runs of positive or negative residuals may indicate a biased prediction pattern or another model discrepancy, not just the wrong variance scale. A few small illustrative observations cannot reliably distinguish these patterns; diagnostics require adequate data across the assay’s working range.

The models are not interchangeable. Lognormal error keeps simulated observations positive and treats reciprocal multiplicative deviations symmetrically on the log scale. Proportional normal error also scales variability with the prediction but retains a symmetric normal distribution on the original scale. Additive error is often natural near a measurement background where absolute noise is approximately constant. A combined additive-plus-proportional variance can be useful when both background noise and signal-dependent noise are scientifically justified, but it introduces additional assumptions and parameters and is not needed simply because one model fails to fit.

Common Pitfalls and Limitations

Comparing raw residuals across different magnitudes. A 5-unit residual is not automatically equally surprising at predictions of 10 and 100. Use residuals standardized according to the chosen error model.
Using lognormal error for nonpositive data. The logarithm is undefined for zero or negative observations. Do not add an arbitrary constant without a measurement-based justification; that changes the model.
Calling a lognormal prediction the mean without checking the parameterization. In the formulation here, \(m_i\) is the median and the mean is \(m_i\exp(\sigma_a^2/2)\).
Dropping variance terms from the likelihood. When the variance depends on the prediction or an estimated error scale, a residual sum of squares alone generally is not the full likelihood objective.
Interpreting a fitted error scale as biological variability. The error model describes conditional observation variability as specified; unexplained biological variation or model misspecification may also contribute to residuals.

These error models are conditional descriptions, not guarantees that every discrepancy is measurement error. They may be inadequate when an assay has a background offset, a detection boundary, batch effects, or dependence among observations. They also do not correct a systematically biased QSP prediction. A useful calibration therefore checks both the residual scale and residual structure, and reports the measurement scale and the error-model parameterization unambiguously.

Key takeaway: Additive error standardizes absolute differences, proportional error standardizes relative differences, and lognormal error standardizes log ratios. Choose the model whose assumptions match the assay over the observed range, then interpret calibration residuals on that model’s scale.

Check Your Understanding

Answer each question, then compare with the solutions below.

An assay prediction is 25 units and its observation is 30 units. What are the raw and relative residuals?
For proportional normal error with \(\sigma_a=0.12\) and prediction \(m_i=40\) units, what is the observation standard deviation?
For lognormal error with median prediction 50 units and \(\sigma_a=0.20\), what is the conditional arithmetic mean?
Why is minimizing only the sum of squared raw residuals not generally appropriate for proportional normal error?

Solutions. (1) The raw residual is \(30-25=5\) units; the relative residual is \(5/25=0.20\), or 20%. (2) The standard deviation is \(\sigma_a m_i=0.12(40)=4.8\) units. (3) The mean is \(50\exp(0.20^2/2)=50\exp(0.02)=51.01\) units. (4) The variance changes with the prediction, so residuals must be scaled by their variances, and the likelihood also includes variance terms that can depend on the model predictions.