FOCE Conditional Likelihood

Why Linearize at an Individual Estimate?

In a nonlinear mixed-effects model, the effect of a subject’s random effects on the observations can vary substantially across the random-effects distribution. First-order estimation by linearization approximates each subject’s prediction near zero random effects. That approximation can be poor when the observations for a subject support a substantially nonzero individual effect. First-order conditional estimation (FOCE) instead linearizes near that subject’s empirical Bayes estimate: the conditional mode of the random effects given the subject’s data and the current population parameters.

The distinction is about the expansion point, not about replacing integration with a simple plug-in estimate. FOCE uses the conditional mode to build a local linear model, then integrates the random effects in that linearized model. The result is an approximation to the marginal likelihood, not merely the likelihood conditional on the estimated random effects.

Model and assumptions: For subject \(i\), let \(y_i\) be the vector of observations and let \(f_i(\theta,\eta_i)\) be the model prediction, where \(\theta\) denotes population parameters and \(\eta_i\sim N(0,\Omega)\) is a vector of random effects. Assume \(y_i\mid\eta_i\sim N(f_i(\theta,\eta_i),R_i)\), with residual error independent of \(\eta_i\) and \(R_i\) not depending on \(\eta_i\). Assume \(R_i\) and \(\Omega\) are positive definite, and \(f_i\) is differentiable near the conditional mode. These assumptions describe the non-interaction residual-error setting.

Given \(\theta\), the conditional mode \(\hat{\eta}_i\) maximizes the joint density of \(y_i\) and \(\eta_i\), or equivalently minimizes the negative log density up to constants:

Q_i(\eta_i)= \{y_i-f_i(\theta,\eta_i)\}^{\mathsf T}R_i^{-1} \{y_i-f_i(\theta,\eta_i)\} +\eta_i^{\mathsf T}\Omega^{-1}\eta_i .

At an interior local mode, the gradient is zero; if the Hessian exists, it is negative semidefinite for the conditional log density (equivalently, positive semidefinite for \(Q_i\)). It need not be unique. In practice, the mode is found numerically for the current population parameters; a converged optimizer alone does not establish that it found a suitable mode.

Deriving the FOCE Approximation

Write \(\hat f_i=f_i(\theta,\hat{\eta}_i)\), and let \(G_i\) be the Jacobian matrix of \(f_i\) with respect to \(\eta_i\), evaluated at \(\hat{\eta}_i\). A first-order Taylor expansion around that point is

f_i(\theta,\eta_i)\ \approx\ \hat f_i+G_i(\eta_i-\hat{\eta}_i) =a_i+G_i\eta_i, \qquad a_i=\hat f_i-G_i\hat{\eta}_i .

The intercept \(a_i\) matters: this conditional linearization is not generally the same as using \(f_i(\theta,0)\) as the intercept. Holding the expansion point, \(a_i\), and \(G_i\) fixed, the linearized observation model is Gaussian. Integrating its random effects gives

\int \phi\!\left(y_i;\,a_i+G_i\eta_i,R_i\right) \phi\!\left(\eta_i;\,0,\Omega\right)\,d\eta_i = \phi\!\left(y_i;\,a_i,V_i\right), \qquad V_i=R_i+G_i\Omega G_i^{\mathsf T}.

Here \(\phi(\cdot;\mu,\Sigma)\) denotes a normal density with mean \(\mu\) and covariance matrix \(\Sigma\). The identity follows because, under the fixed linearized model, \(G_i\eta_i\) is normal with covariance \(G_i\Omega G_i^{\mathsf T}\), independent of the residual error. The sum of these independent Gaussian terms has covariance \(V_i\).

There is an important qualification. In FOCE, \(\hat{\eta}_i\), \(a_i\), and \(G_i\) are selected using the observed \(y_i\). Thus the normal density above, evaluated at the observed data with these quantities, is a pointwise approximation to the marginal likelihood. It does not, in general, define a normalized sampling distribution for \(y_i\) when the expansion point is selected from \(y_i\). The fixed-expansion-point Gaussian integral justifies the approximation at that point; it does not remove its data dependence.

FOCE likelihood contribution: Under the stated assumptions, the subject’s approximate log-likelihood contribution is the log of the Gaussian density evaluated at the observed \(y_i\), using the conditionally linearized \(a_i\) and \(V_i\). The approximation includes integration over the random effects; using only \(p(y_i\mid\hat{\eta}_i)\) would be a conditional plug-in likelihood and would omit that integration.

A useful equivalent expression shows how the conditional residual and random-effects penalty enter. Let \(\hat r_i=y_i-\hat f_i\), and define \(S_i=\Omega^{-1}+G_i^{\mathsf T}R_i^{-1}G_i\). The mode condition for the linearized model is

G_i^{\mathsf T}R_i^{-1}\hat r_i=\Omega^{-1}\hat{\eta}_i .

This condition also follows from the gradient of the original joint log density at an interior conditional mode. For the linearized model, the quadratic form in the integrand is a quadratic function of \(\eta_i\). Its gradient vanishes at \(\hat{\eta}_i\), so completing the square gives the following identity:

\begin{aligned} &\{y_i-a_i-G_i\eta_i\}^{\mathsf T}R_i^{-1} \{y_i-a_i-G_i\eta_i\} +\eta_i^{\mathsf T}\Omega^{-1}\eta_i\\ &\qquad = \hat r_i^{\mathsf T}R_i^{-1}\hat r_i +\hat{\eta}_i^{\mathsf T}\Omega^{-1}\hat{\eta}_i +(\eta_i-\hat{\eta}_i)^{\mathsf T}S_i(\eta_i-\hat{\eta}_i). \end{aligned}

Integrating the final centered quadratic contributes the determinant term. The matrix determinant lemma gives \(|V_i|=|R_i||\Omega||S_i|\). Consequently, the negative twice approximate log likelihood, apart from no omitted constants, is

-2\ell_{i,\mathrm{FOCE}}= m_i\log(2\pi)+\log|R_i|+\log|\Omega|+\log|S_i| +\hat r_i^{\mathsf T}R_i^{-1}\hat r_i +\hat{\eta}_i^{\mathsf T}\Omega^{-1}\hat{\eta}_i ,

where \(m_i\) is the number of observations for subject \(i\). This form makes clear why evaluating only the residual sum of squares at the conditional mode is not the FOCE likelihood: the random-effects penalty and the determinant terms are also part of the integrated approximation.

Worked Example: Conditional Linearization in a Nonlinear Model

Consider one hypothetical subject with two observations, residual covariance \(R_i=0.04I\), and one random effect \(\eta_i\sim N(0,0.25)\). Let \(f_1(\eta)=1+\eta+0.5\eta^2\) and \(f_2(\eta)=2+0.5\eta+0.25\eta^2\). Suppose the observations are \(y_1=1.24\) and \(y_2=2.123333\). Check whether \(\hat{\eta}=0.2\) is a conditional mode, then calculate the FOCE covariance and quadratic contribution.

At \(\eta=0.2\), the predictions are \(f_1=1.22\) and \(f_2=2.11\), so the residuals are \(0.02\) and approximately \(0.013333\). The derivatives are \(f'_1=1.2\) and \(f'_2=0.6\). The derivative of the conditional log density with respect to \(\eta\) is

\frac{(y_1-f_1)f'_1+(y_2-f_2)f'_2}{0.04} -\frac{\eta}{0.25} = \frac{0.02(1.2)+0.013333(0.6)}{0.04}-0.8 \approx 0.

With the stated rounded observation, the score is approximately zero; taking \(y_2=2.123333333\ldots\) makes it exactly zero at \(\eta=0.2\). The second derivative of the log density there is

\frac{-(1.2)^2+0.02(1)}{0.04} +\frac{-(0.6)^2+0.013333333\ldots(0.5)}{0.04} -\frac{1}{0.25} \approx -48.33<0 .

Thus this is a local conditional mode (for the exactly specified observation, with the rounded values giving the same result to the displayed precision). The Jacobian is \(G=(1.2,0.6)^{\mathsf T}\), and \(a=(1.22,2.11)^{\mathsf T}-0.2G=(0.98,1.99)^{\mathsf T}\). The covariance is

V_i=0.04I+0.25 \begin{pmatrix}1.2\\0.6\end{pmatrix} \begin{pmatrix}1.2&0.6\end{pmatrix} = \begin{pmatrix}0.40&0.18\\0.18&0.13\end{pmatrix}.

Using the exact \(y_2=2.123333333\ldots\), \(y_i-a_i=(0.26,0.133333333\ldots)^{\mathsf T}\). The determinant is \(0.40(0.13)-0.18^2=0.0196\), and the Gaussian quadratic is

(y_i-a_i)^{\mathsf T}V_i^{-1}(y_i-a_i) = \frac{0.13(0.26)^2-2(0.18)(0.26)(0.133333333\ldots)+0.40(0.133333333\ldots)^2}{0.0196} \approx 0.1744.

This agrees with the conditional residual-plus-penalty form: \(0.02^2/0.04+0.013333333\ldots^2/0.04+0.2^2/0.25\approx0.1744\). The calculation is internally consistent, and the determinant term remains necessary for the likelihood.

How FOCE Differs from Linearization at Zero

First-order estimation by linearization expands \(f_i\) at \(\eta_i=0\), using the prediction \(f_i(\theta,0)\) and the Jacobian there. FOCE uses the conditional-mode prediction and Jacobian, and adjusts the intercept to \(a_i=\hat f_i-G_i\hat{\eta}_i\). For a nonlinear model, both the slope and the resulting intercept can differ. For a model linear in \(\eta_i\), the slope is constant and the adjusted intercept equals \(f_i(\theta,0)\); the two linearized likelihoods then coincide.

Worked Example: When FO and FOCE Coincide

For a hypothetical subject, suppose \(f(\eta)=(1+\eta,\ 2+0.5\eta)^{\mathsf T}\), \(R=0.09I\), \(\Omega=0.16\), and \(y=(1.3,\ 2.2)^{\mathsf T}\). The model is linear in \(\eta\), so its Jacobian is \(G=(1,0.5)^{\mathsf T}\) at both zero and the conditional mode.

The conditional-mode equation gives

\hat{\eta} = \frac{G^{\mathsf T}R^{-1}(y-f(0))} {\Omega^{-1}+G^{\mathsf T}R^{-1}G} = \frac{0.4/0.09}{1/0.16+1.25/0.09} = \frac{32}{145} \approx0.2207.

The conditional residuals are \(0.3-\hat{\eta}\approx0.07931\) and \(0.2-0.5\hat{\eta}\approx0.08966\). The conditional intercept is \(f(\hat{\eta})-G\hat{\eta}=(1,2)^{\mathsf T}=f(0)\), and the covariance is

V=0.09I+0.16GG^{\mathsf T} = \begin{pmatrix}0.25&0.08\\0.08&0.13\end{pmatrix}.

The determinant is \(0.25(0.13)-0.08^2=0.0261\). With \(y-f(0)=(0.3,0.2)^{\mathsf T}\), the quadratic is

\frac{0.13(0.3)^2-2(0.08)(0.3)(0.2)+0.25(0.2)^2}{0.0261} = \frac{0.0121}{0.0261} \approx0.4636.

Because the model is linear in the random effect, expansion at zero is exact; shifting the expansion point changes neither the intercept nor the Jacobian. FO and FOCE therefore give the same approximate marginal likelihood in this example.

Diagnostics, Failure Modes, and Limitations

FOCE is most useful when the conditional mode is meaningfully informed by the subject’s data and the prediction is locally close to linear over the random-effects region that contributes to the integral. The approximation can be unreliable when that local description is poor or when the conditional mode is difficult to determine.

Check the mode, not just optimizer status. Inspect the score near the reported mode and the curvature of the conditional log density. A substantial nonzero score means the point is not an interior mode. A flat-curvature direction is a diagnostic warning requiring further assessment; it does not by itself establish weak identification or rule out a local maximum. If the Hessian exists, a direction of positive curvature in the conditional log density rules out a local maximum.
Look for sensitivity to starting values. Different solutions can indicate multiple conditional modes. If plausible starts lead to different modes or materially different FOCE contributions, a single mode may not represent the integral well.
Assess local linearity. Large curvature of \(f_i\) near \(\hat{\eta}_i\), or predictions that change strongly across plausible random-effect values, warns that a first-order expansion may be inadequate.
Watch sparse subjects. With little information per subject, the data may not locate \(\hat{\eta}_i\) reliably. The mode can then be strongly influenced by the random-effects distribution or be unstable to small data changes.
Verify the error assumptions. The derivation above assumes Gaussian residuals with covariance independent of \(\eta_i\). It should not be transferred unchanged to a model in which residual variability depends on random effects.

Worked Example: Sparse Data and Multiple Conditional Modes

Consider a hypothetical one-observation model \(y=\eta^2+\epsilon\), with \(y=1\), \(\eta\sim N(0,1)\), and \(\epsilon\sim N(0,0.04)\). The conditional log-density score is \(2\eta(1-\eta^2)/0.04-\eta\). Its roots include \(\eta=0\) and \(\eta=\pm\sqrt{0.98}\approx\pm0.9899\).

At zero, the second derivative of the log density is \(2(1)/0.04-1=49>0\), so zero is not a local maximum. At either nonzero root, the residual is \(1-0.98=0.02\), and the second derivative is

\frac{-(2\eta)^2+0.02(2)}{0.04}-1 = \frac{-4(0.98)+0.04}{0.04}-1 =-98<0.

There are two symmetric conditional modes. At the positive mode, \(G=2\sqrt{0.98}\approx1.9799\), \(\hat f=0.98\), and \(a=\hat f-G\hat{\eta}=0.98-1.96=-0.98\). The FOCE variance is \(V=0.04+G^2=3.96\), and \(y-a=1.98\). The Gaussian quadratic is \(1.98^2/3.96=0.99\), exactly matching the residual-plus-penalty value \(0.02^2/0.04+0.98=0.99\). The negative mode gives the same approximation by symmetry.

By contrast, linearizing at zero gives \(f(0)=0\), \(G(0)=0\), and variance \(0.04\), producing quadratic \(1^2/0.04=25\). This large discrepancy illustrates how the expansion point can matter in a nonlinear model. The tied modes also illustrate a diagnostic limitation: the data do not distinguish the sign of \(\eta\), so the individual mode is not uniquely determined even though the two local FOCE calculations agree.

FOCE’s principal alternative in this comparison is linearization at zero random effects, which is simpler and avoids subject-specific mode searches but can be inaccurate when individual effects are far from zero. Exact integration of the nonlinear subject likelihood avoids the local Taylor approximation, but is generally more demanding. FOCE occupies a middle ground: it uses subject-specific information while retaining an analytically integrated local Gaussian approximation. It can still be biased when the local expansion is poor, and numerical stability does not by itself establish approximation accuracy.

Common mistake: “FOCE evaluates the likelihood at the empirical Bayes estimate, so I can use \(p(y_i\mid\hat{\eta}_i)\) as the subject contribution.” This omits the random-effects integration and its determinant contribution. FOCE instead linearizes around \(\hat{\eta}_i\) and evaluates the resulting integrated Gaussian approximation. Also, because the expansion point depends on \(y_i\), that pointwise approximation should not be described as a normalized sampling distribution for \(y_i\).

Check Your Understanding

Answer each question, then compare with the solutions below.

What makes an empirical Bayes estimate the expansion point in FOCE, and what condition must hold at an interior conditional mode?
Why is the conditional FOCE intercept \(a_i\) generally not \(f_i(\theta,0)\)?
For a linear model in \(\eta_i\), explain why FO and FOCE can coincide.
Why does evaluating only the residual likelihood at \(\hat{\eta}_i\) fail to give the FOCE likelihood?
What does it mean to call the FOCE Gaussian expression a pointwise approximation rather than a normalized sampling distribution?

Solutions. (1) It is the mode of the joint density of the subject’s observations and random effects, conditional on the current population parameters; the conditional log-density score must be zero and its curvature must indicate a local maximum. (2) The intercept is \(\hat f_i-G_i\hat{\eta}_i\), which adjusts the local line so it passes through the prediction at the conditional mode; nonlinearity means this need not equal the prediction at zero. (3) A linear model has a constant Jacobian, and the adjusted intercept equals its value at zero, so the two linearizations are the same. (4) The integrated likelihood also includes the random-effects penalty and determinant terms, as well as the conditional residual contribution. (5) The expansion point and linearized coefficients were chosen using the observed data, so the Gaussian density at those data is an approximation to the marginal likelihood value, not generally a sampling density normalized over all possible observations.

What You'll Learn