Worked FOCE Example With Proportional Error

Why Work Through a Proportional-Error FOCE Calculation?

FOCE requires two linked calculations: first find the conditional mode of the random effect for the subject, then linearize the individual prediction around that mode to construct an approximate marginal contribution. With proportional residual error, the residual variance depends on the prediction, so the conditional objective itself changes with the random effect. The interaction and non-interaction formulations handle that dependence differently.

The example here uses one observation and one random effect so each step can be checked directly. The same logic extends to several observations, with vector predictions and covariance matrices, but the scalar version makes it especially clear why the mode, linearized mean, and variance must all be calculated using a consistent formulation.

Definition: Let \(Y\mid\eta\sim N(f(\eta),R(\eta))\), where \(f(\eta)>0\) is the prediction and \(\eta\sim N(0,\omega^2)\). The conditional mode \(\hat\eta\) maximizes the conditional posterior density of \(\eta\) given the observation, or equivalently minimizes its negative log density. In FOCE, the prediction is then linearized around \(\hat\eta\) to approximate the subject's marginal likelihood contribution.

Model and FOCE Calculation

Suppose residual error is proportional, with coefficient of variation \(\sigma\). In this tutorial's model, the conditional observation variance is \(R(\eta)=\sigma^2 f(\eta)^2\). The normal log density includes both the standardized squared residual and the logarithm of this variance. That logarithm matters when finding a conditional mode because \(R\) changes with \(\eta\).

For the interaction formulation, the conditional mode is found using this eta-dependent variance. The FOCE approximation retains the residual variance evaluated at that mode. For the non-interaction formulation, residual variance is held at its typical-subject value \(R(0)\) while finding the mode and forming the FOCE contribution. These are distinct approximations; they should not be mixed within a single calculation.

q_I(\eta)=\frac{1}{2}\left\{\log R(\eta)+\frac{[y-f(\eta)]^2}{R(\eta)}\right\}+\frac{\eta^2}{2\omega^2}

Terms constant in \(\eta\) have been omitted from this negative log posterior. The interaction mode \(\hat\eta_I\) minimizes \(q_I\). In the non-interaction calculation, put \(R_0=R(0)\) and use

q_N(\eta)=\frac{[y-f(\eta)]^2}{2R_0}+\frac{\eta^2}{2\omega^2}.

To obtain the FOCE contribution, linearize the prediction at the chosen mode: \(f(\eta)\approx f(\hat\eta)+f'(\hat\eta)(\eta-\hat\eta)\). Integrating this linearized normal model over \(\eta\sim N(0,\omega^2)\) gives the approximate mean \(m\) and variance \(V\) below. Here \(R_*\) is \(R(\hat\eta_I)\) for interaction, and \(R_0\) for non-interaction.

m=f(\hat\eta)-f'(\hat\eta)\hat\eta,\qquad V=[f'(\hat\eta)]^2\omega^2+R_*

For a single observation, the resulting log contribution is \(\ell_{\mathrm{FOCE}}=-\tfrac12\{\log(2\pi V)+(y-m)^2/V\}\). The mode is not substituted directly as the mean of the marginal distribution: the linearized mean includes the adjustment \(-f'(\hat\eta)\hat\eta\).

Worked Calculation: Interaction Mode

Worked Example: Finding the Interaction Conditional Mode

Consider the hypothetical observation \(y=120\) concentration units and prediction \(f(\eta)=100e^\eta\) concentration units. Let \(\eta\sim N(0,0.09)\), so \(\omega^2=0.09\), and let the proportional error coefficient be \(\sigma=0.20\). Thus \(R(\eta)=0.04f(\eta)^2\), in squared concentration units. All parameters and data here are illustrative.

For this model, define \(r=y/f(\eta)\). Since \(f'(\eta)=f(\eta)\), differentiating the interaction objective gives

q_I'(\eta)=-\frac{(r-1)r}{0.04}+1+\frac{\eta}{0.09}.

The first term comes from differentiating the standardized residual, the \(+1\) comes from \(\tfrac12\log R(\eta)\), and the last term is the normal random-effect penalty. To locate the zero, at \(\eta=0.105\) the prediction is \(100e^{0.105}\approx111.071\), so \(r\approx1.080389\). Substitution gives

q_I'(0.105)\approx-\frac{(1.080389-1)(1.080389)}{0.04}+1+\frac{0.105}{0.09} \approx -0.00463.

The derivative is slightly negative. The second derivative is \(q_I''(\eta)=r(2r-1)/0.04+1/0.09\), which is approximately \(42.5\) at this point. A Newton step for solving \(q_I'(\eta)=0\) is therefore \(\eta_{\mathrm{new}}=0.105-q_I'(0.105)/q_I''(0.105)\approx0.105-(-0.004618)/42.5=0.105109\). Numerically solving the score equation gives \(\hat\eta_I\approx0.10511\). The positive Newton adjustment is consistent with a negative score and positive curvature.

At this mode, \(f(\hat\eta_I)\approx111.08\) concentration units and \(R(\hat\eta_I)=0.04f(\hat\eta_I)^2\approx493.6\) squared concentration units. The random-effect penalty is \(\hat\eta_I^2/(2\omega^2)\approx0.0614\). These are mode-dependent quantities, so they should be recomputed if the mode is updated.

Worked Calculation: Non-Interaction Mode

Worked Example: Holding the Residual Variance Fixed

Use the same hypothetical observation, prediction, and random-effect distribution, but use the non-interaction formulation. At the typical value \(\eta=0\), \(f(0)=100\) concentration units, so the fixed residual variance is \(R_0=0.04(100)^2=400\) squared concentration units.

Differentiating \(q_N\) gives the score

q_N'(\eta)=-\frac{[y-f(\eta)]f(\eta)}{400}+\frac{\eta}{0.09}.

At \(\eta=0.13705\), \(f(\eta)=100e^{0.13705}\approx114.68855\) concentration units. The two score terms are approximately \(-[(120-114.68855)(114.68855)]/400\approx-1.5229\) and \(0.13705/0.09\approx1.5228\). Their sum is about \(-0.00013\), so the score is near zero and \(\hat\eta_N\approx0.13705\). The two modes differ because the interaction objective accounts for changing residual variance, whereas this objective uses \(R_0\) throughout.

For a curvature check, \(q_N''(\eta)=f(\eta)[2f(\eta)-y]/400+1/0.09\), which is positive at this mode. In fact it is positive for every positive \(f\) in this example: the minimum of \(f(2f-120)/400\) over \(f>0\) is \(-4.5\), while \(1/0.09\approx11.11\), so the total curvature is bounded below by \(11.11-4.5\approx6.61>0\). Thus \(q_N\) is strictly convex here and its stationary point is its unique minimum.

Worked Calculation: FOCE Contributions

Worked Example: Comparing the Two FOCE Contributions

Continue with the hypothetical data above. For the exponential prediction, \(f'(\eta)=f(\eta)\), so the linearized mean is \(m=f(\hat\eta)(1-\hat\eta)\). For interaction, use \(\hat\eta_I\approx0.10511\), \(f(\hat\eta_I)\approx111.08\), and \(R_*=R(\hat\eta_I)\approx493.6\).

The interaction mean is \(m_I\approx111.08(1-0.10511)\approx99.41\) concentration units. Its variance is \(V_I\approx(111.08)^2(0.09)+493.6\approx1604.1\) squared concentration units. The squared standardized difference is \((120-99.41)^2/1604.1\approx0.2644\). Therefore,

\ell_{\mathrm{FOCE},I}\approx-\frac12\{\log(2\pi\cdot1604.1)+0.2644\}\approx-4.741.

For non-interaction, use \(\hat\eta_N\approx0.13705\), \(f(\hat\eta_N)\approx114.69\), and \(R_0=400\). Then \(m_N\approx114.69(1-0.13705)\approx98.97\) concentration units, and \(V_N\approx(114.69)^2(0.09)+400\approx1583.8\) squared concentration units. The squared standardized difference is \((120-98.97)^2/1583.8\approx0.2792\), so

\ell_{\mathrm{FOCE},N}\approx-\frac12\{\log(2\pi\cdot1583.8)+0.2792\}\approx-4.742.

The two approximate log contributions are close in this example, but their closeness is not guaranteed. The residual-variance convention changes the conditional mode and the variance used in the marginal approximation; differences may be larger with more influential random effects or more strongly varying predictions.

Why the Formula Has This Form

The linearized prediction can be written \(f(\eta)\approx a+b\eta\), where \(b=f'(\hat\eta)\) and \(a=f(\hat\eta)-b\hat\eta\). Because \(\eta\sim N(0,\omega^2)\), the linearized prediction has mean \(a\) and variance \(b^2\omega^2\). Adding independent residual error with variance \(R_*\) gives mean \(a\) and total variance \(b^2\omega^2+R_*\). This yields the FOCE mean and variance used above.

This argument also explains a frequent error: using \(f(\hat\eta)\) as the marginal mean while adding the random-effect variance term. Under the linearized model, the random effect still has population mean zero, not conditional mean \(\hat\eta\); the intercept of the linearized function is \(f(\hat\eta)-f'(\hat\eta)\hat\eta\). The mode determines where the slope is evaluated, but does not replace the zero-mean distribution in the integration.

A second important detail is the normalizing term \(\log R(\eta)\) in the interaction conditional objective. A normal density with variance \(R(\eta)\) is proportional to \(R(\eta)^{-1/2}\exp\{-[y-f(\eta)]^2/[2R(\eta)]\}\). Omitting the variance factor changes the conditional score and therefore the mode. In the non-interaction objective, \(R_0\) is constant with respect to \(\eta\), so its log term can be omitted when locating the mode.

Assumptions, Diagnostics, and Limitations

These calculations assume a positive prediction, a normal random effect with mean zero and variance \(\omega^2\), and a conditionally normal observation with the stated proportional variance. They also treat the population parameters as fixed during the subject-level mode calculation. For several independent residuals, the conditional objective is the sum of their log-variance and standardized-residual terms, together with the random-effect penalty; the FOCE covariance then uses the prediction gradient and the residual covariance matrix.

Several checks help detect mistakes or numerical trouble:

Check the score. Substitute the reported mode into \(q'(\hat\eta)\). It should be near zero at the precision claimed. A score with the wrong sign can reveal a Newton-step sign error.
Check curvature and minimization. A positive second derivative at a stationary point supports a local minimum. If optimization gives a boundary value, multiple candidate roots, or unstable modes under small changes, investigate the objective rather than assuming the result is valid.
Check units and scale. \(f\), \(m\), and \(y\) have concentration units; \(R\) and \(V\) have squared concentration units; and the standardized squared difference is dimensionless.
Check the variance convention. Interaction uses \(R(\hat\eta_I)\); non-interaction uses \(R_0\). Using one convention for the mode and the other for the FOCE contribution produces an internally inconsistent calculation.
Check the prediction and error model. Proportional-only error is unsuitable when predictions can be zero or negative, and may be inadequate when variability does not grow in proportion to prediction. A different residual model may then be needed.

FOCE is a local linearization, not an exact evaluation of the integral over random effects. Strong curvature near the mode, sparse observations, or a conditional distribution that is far from locally well represented by the linearized model can weaken the approximation. First-order (FO) methods instead linearize around the population mean of the random effect, typically zero, and do not adapt the expansion point to the subject's data. FOCE's conditional mode can improve that local fit, but it does not remove approximation error. Direct numerical integration is an alternative when a reliable marginal likelihood is needed and the integration is computationally feasible.

Common pitfall: “The residual variance is proportional, so I can use the same conditional mode and merely change the final variance.” That is incorrect for the interaction formulation: \(R(\eta)\) appears in both the standardized residual and the log-variance term, so it changes the score equation and the mode. Solve the mode under the selected formulation, then use that same formulation consistently in the FOCE contribution.

Key takeaway: In proportional-error FOCE, distinguish the mode-finding objective from the final integrated approximation, but keep their variance convention consistent. Verify the score at the mode, then calculate the linearized mean and total variance before evaluating the normal log contribution.

Check Your Understanding

Use the definitions and calculations above. The solutions are included so each numerical step can be checked.

In the interaction objective, why does differentiating the log-variance term contribute \(+1\) for \(f(\eta)=100e^\eta\) and \(R(\eta)=0.04f(\eta)^2\)?
For the non-interaction calculation, what fixed residual variance is used if \(f(0)=80\) concentration units and \(\sigma=0.15\)?
For a linearized prediction with \(f(\hat\eta)=50\), \(f'(\hat\eta)=20\), and \(\hat\eta=0.3\), calculate the linearized mean.
At a proposed mode, the score is \(-0.02\) and the curvature is \(4\). What is one Newton update for the mode, and in which direction does it move?

Solutions. (1) Since \(\log R=\log(0.04)+2\log f\) and \(f'/f=1\), the derivative of \(\tfrac12\log R\) is \(\tfrac12(2f'/f)=1\). (2) \(R_0=(0.15\cdot80)^2=12^2=144\) squared concentration units. (3) \(m=f(\hat\eta)-f'(\hat\eta)\hat\eta=50-(20)(0.3)=44\), in the prediction's units. (4) The Newton update is \(\eta_{\mathrm{new}}=\eta-(-0.02)/4=\eta+0.005\); it moves upward because the score is negative and curvature positive.