Residual Error Models Across Software Syntaxes

Why Residual Error Syntax Needs Care

A structural model supplies a prediction for an observation; a residual error model describes how an observed value may differ from that prediction. The same verbal label—such as “combined error”—does not guarantee that two software implementations use the same equation. One implementation may add two independent error terms, while another may specify the standard deviation of a single normal error directly. Those constructions agree only when their scales are translated correctly.

This matters when comparing or transferring models among NONMEM, Monolix, nlmixr2, and mrgsolve. The goal is not to make the syntax look alike. It is to verify that each version gives the observation the same conditional mean, variance, and units. The notation below uses \(f_{ij}\) for the model-predicted observation for individual \(i\) at observation \(j\), and \(Y_{ij}\) for the corresponding observed value. For a concentration model, both are measured in concentration units.

Observation Equations and Assumptions

A common additive model is \(Y_{ij}=f_{ij}+\epsilon_{a,ij}\), where \(\epsilon_{a,ij}\) has conditional mean zero and standard deviation \(\sigma_a\). If the errors are normal, then \(Y_{ij}\mid f_{ij}\sim N(f_{ij},\sigma_a^2)\). The parameter \(\sigma_a\) has the same units as \(Y_{ij}\), such as mg/L for a concentration.

A common proportional model is \(Y_{ij}=f_{ij}(1+\epsilon_{p,ij})\), where \(\epsilon_{p,ij}\) has conditional mean zero and standard deviation \(\sigma_p\). Under a normal error assumption, \(Y_{ij}\mid f_{ij}\sim N(f_{ij},(f_{ij}\sigma_p)^2)\). The parameter \(\sigma_p\) is dimensionless: it is a relative error scale, not a concentration scale.

One common combined model adds independent additive and proportional normal errors:

Assuming the two errors are independent, the conditional mean is \(f_{ij}\), and the conditional variance is the sum of their variances:

Thus the combined model can also be expressed as a single normal observation with standard deviation \(\sqrt{\sigma_a^2+(f_{ij}\sigma_p)^2}\). This equivalence depends on the errors being independent and normal. A different combined convention, such as using \(\sigma_a+|f_{ij}|\sigma_p\) as the standard deviation, is a different model and should not be substituted without checking the intended equation.

Worked Example: Calculating a Combined Error Scale

Translating the Equations Across Four Tools

The following snippets illustrate the same normal-error assumptions and the independent-error combined variance above. They are schematic model fragments: the prediction \(f\) must be replaced with the variable or expression used in the particular structural model, and parameter declarations must match the surrounding model file. Always compare the implemented equation, not just parameter names.

NONMEM. In a typical control stream, residual normal random variables are represented by EPS terms and their variances are specified in the SIGMA matrix. An additive model can be written in the error code as Y = F + EPS(1); for an additive variance \(\sigma_a^2\), the corresponding diagonal SIGMA value is \(\sigma_a^2\). A proportional model can be written as Y = F + F*EPS(1), with SIGMA value \(\sigma_p^2\). For the independent combined model, use separate EPS terms, for example Y = F + EPS(1) + F*EPS(2), with diagonal SIGMA values \(\sigma_a^2\) and \(\sigma_p^2\). Here EPS(1) has observation units and EPS(2) is dimensionless. In NONMEM, the SIGMA entries are variances, not standard deviations.

Monolix. In an observation model, the prediction is linked to the observed variable and the error model is specified explicitly. For example, the schematic form y = {type = continuous, prediction = f, errorModel = constant(a)} represents a constant additive standard deviation \(a\); errorModel = proportional(b) represents a proportional scale \(b\). For the combined variance used here, the Monolix form is errorModel = combined2(a,b), where \(a\) is the additive scale and \(b\) the proportional scale. Verify that the selected model option denotes \(SD=\sqrt{a^2+(bf)^2}\) in the software version and model specification being used.

nlmixr2. The residual formula is attached to the observation variable. A schematic model statement is DV ~ add(add.sd) + prop(prop.sd) for an additive-plus-proportional model; the component parameters are standard-deviation scales. The corresponding single-component forms are DV ~ add(add.sd) and DV ~ prop(prop.sd). Check the fitted model’s parameterization and prediction convention when translating an existing model, particularly if the source uses a different combined-error formula.

mrgsolve. mrgsolve can implement a residual equation for simulation, but it is not a population parameter estimation engine. For reproducible simulations, one option is to provide standardized residual draws as input columns and apply the equation in the model. A schematic fragment is:

Here EADD and EPROP are externally supplied standard-normal draws, SIGA is in observation units, and SIGP is dimensionless. With independent standard-normal inputs, this produces the same conditional variance as the combined equation. Supplying the draws explicitly makes the random inputs inspectable and reusable; setting them to zero would produce predictions without residual noise, not a residual-error simulation.

Worked Example: Matching Four Implementations

Derivation and Model Checks

The combined variance follows directly from variance rules. Write \(Y=f+\epsilon_a+f\epsilon_p\), treating \(f\) as fixed when conditioning on the model prediction. Since both error terms have mean zero, \(E(Y\mid f)=f\). For independent components, the covariance is zero, so

The square root gives the standard deviation. This derivation also shows why an additive and proportional term should not automatically be combined by adding their standard deviations: independent variances add, whereas standard deviations do not. If the errors are dependent, the covariance term need not be zero and the target equation must state how that dependence is handled.

Useful implementation checks include evaluating the implied conditional standard deviation at several predictions, confirming parameter units, and checking that each software parameter is a variance or standard deviation as expected. At \(f=0\), the proportional contribution vanishes and the combined model has standard deviation \(\sigma_a\). For predictions near zero, this means proportional error alone can imply a very small observation scale; a combined model may be preferable when a nonzero absolute scale is scientifically justified. If a prediction can be negative, the variance still uses \(f^2\), but a proportional normal model can generate observations of either sign. Whether that is acceptable depends on the measured quantity and the model’s intended observation distribution.

Diagnostics, Alternatives, and Limitations

After fitting or simulating, inspect residual patterns against predictions and time, and compare residual spread across the prediction range. A roughly constant spread is consistent with additive error; spread increasing approximately in proportion to the magnitude of predictions is consistent with proportional error. A combined model allows both patterns. These are diagnostic clues, not proof that a particular model is correct. Residual trends can also reflect a misspecified structural model or data problems.

Common signs of a translation failure include a different spread at large predictions, unexplained differences in the objective or likelihood, and residuals that become systematically too large or too small across the prediction range. Check for a parameter entered as a standard deviation where the software expects a variance, or vice versa. Also confirm whether a software option defines a combined scale as a sum or as a square root of summed squares. A simulation comparison is informative only when the same predictions, parameter values, and random standardized draws are used.

These normal-error forms can produce negative observations, particularly when predictions are small relative to the error scale. That may be unsuitable for a strictly positive measurement. Changing to a transformation or another observation distribution changes the model and requires its own assumptions; it is not merely a syntax adjustment. Likewise, mrgsolve’s explicit simulation equation does not estimate residual parameters, so estimation and simulation responsibilities should not be confused when reproducing a workflow.

Common Mistakes

• Assigning concentration units to the proportional scale. The parameter \(\sigma_p\) is dimensionless; only the product \(f\sigma_p\) has observation units.
• Entering a standard deviation as a NONMEM SIGMA value. A diagonal SIGMA entry is a variance, so a standard deviation of \(0.50\) corresponds to a variance of \(0.25\).
• Adding independent standard deviations. Independent variances add; the combined standard deviation is the square root of their sum.
• Assuming all “combined” options mean the same thing. Verify the exact equation and implied variance in the selected software option.
• Treating simulation code as an estimation model. An mrgsolve observation equation can generate simulated noisy observations, but does not by itself fit the residual model.