Why FOCE Requires an Integrated Subject Contribution
The marginal likelihood for a subject averages that subject’s likelihood over the population distribution of random effects. In a nonlinear model, this average is generally not available as a closed-form integral. FOCE constructs an approximation by expanding the individual prediction around the subject’s conditional random-effect mode, then integrating the resulting local Gaussian model. The integration matters: evaluating the conditional likelihood only at the mode does not give the FOCE subject contribution.
This derivation assumes normally distributed random effects and normally distributed residuals whose covariance does not depend on the random effects. The residual covariance may differ between subjects, but it is held fixed as the random effect varies. The approximation is local: the nonlinear prediction is replaced by its first-order Taylor expansion at the conditional mode.
Write \(J_i\) for the matrix of first derivatives of \(f_i\) with respect to \(\eta_i\), evaluated at \(\hat\eta_i\). Its dimensions are the number of observations for subject \(i\) by the number of random effects. The mode is assumed to be an interior minimizer. Because \(R_i\) is fixed with respect to \(\eta_i\), differentiating \(q_i\) gives the score condition below.
This condition is a useful check on mode optimization, but it does not itself determine the FOCE likelihood. The likelihood comes from integrating the locally linearized conditional model against the population distribution of \(\eta_i\).
Deriving the FOCE Subject Contribution
At the conditional mode, the first-order Taylor approximation to the prediction is \(f_i(\eta_i)\approx f_i(\hat\eta_i)+J_i(\eta_i-\hat\eta_i)\). Conditional on \(\eta_i\), the approximate observation model is therefore Gaussian. Since \(\eta_i\) is also Gaussian, integrating this local model gives a Gaussian distribution for \(y_i\). The steps below show its mean and covariance.
Replace the nonlinear prediction by \(f_i(\hat\eta_i)+J_i(\eta_i-\hat\eta_i)\) in the neighborhood of the conditional mode.
Rewrite the linear prediction as \(f_i(\hat\eta_i)-J_i\hat\eta_i+J_i\eta_i\). Its population mean over \(\eta_i\) is \(f_i(\hat\eta_i)-J_i\hat\eta_i\).
The random component \(J_i\eta_i\) has covariance \(J_i\Omega J_i'\). Adding independent residual covariance \(R_i\) gives total covariance \(R_i+J_i\Omega J_i'\).
Use the resulting multivariate normal density, including its determinant and quadratic terms, as the FOCE approximation to the marginal subject likelihood.
Here \(n_i\) is the number of observations for subject \(i\), and \(\phi_{n_i}(y;\mu,S)\) denotes the \(n_i\)-dimensional normal density at \(y\), with mean \(\mu\) and covariance \(S\). Equivalently, defining \(z_i=r_i(\hat\eta_i)+J_i\hat\eta_i\), the log contribution is
The vector \(z_i\) is the difference between the observations and the mean of the integrated local model. The determinant term is essential: it accounts for the volume of the Gaussian distribution as its covariance changes. Omitting it, or replacing the expression by the conditional log likelihood at \(\hat\eta_i\), does not produce this normalized FOCE contribution.
This construction uses only a first-order expansion of the prediction and then integrates that approximation. It is not a claim that the original nonlinear marginal density is exactly Gaussian. In particular, the approximation depends on the chosen expansion point, which is why locating the conditional mode and checking the local behavior are important.
Worked Calculations
Worked Example: A Scalar Nonlinear Prediction and Its FOCE Contribution
Consider one hypothetical subject with \(y=1.9375\), prediction \(f(\eta)=1+1.5\eta+\eta^2\), residual variance \(R=4\), and random-effect variance \(\Omega=1\). Take the proposed conditional mode to be \(\hat\eta=0.25\). Then \(f(\hat\eta)=1.4375\), \(J=f'(\hat\eta)=2\), and \(r(\hat\eta)=0.5\).
First check that this is the conditional mode, rather than assuming it from the value alone. The score condition is \(Jr/R=2(0.5)/4=0.25=\Omega^{-1}\hat\eta=1(0.25)\). To verify a minimum, the second derivative of \(q\) is \(q''(\eta)=1.09375+2.25\eta+1.5\eta^2\). This quadratic is positive for all \(\eta\): its leading coefficient is positive and its discriminant is \(2.25^2-4(1.5)(1.09375)=-1.5<0\). Thus \(q\) is strictly convex, and the score-zero point is its unique global minimum.
The FOCE mean and variance are \(\mu^{\mathrm{FOCE}}=f(\hat\eta)-J\hat\eta=1.4375-2(0.25)=0.9375\) and \(S=R+J^2\Omega=4+2^2(1)=8\). The observation’s difference from the integrated mean is \(z=y-\mu^{\mathrm{FOCE}}=1.9375-0.9375=1\), also equal to \(r(\hat\eta)+J\hat\eta=0.5+0.5=1\).
Therefore the approximate log contribution is \(-\tfrac12[\log(2\pi)+\log(8)+1^2/8]=-\tfrac12(4.042319)=-2.021\). The corresponding approximate density is about \(0.1324\). These are contributions for the stated hypothetical model and data, not estimates from a real study.
Worked Example: A Two-Observation Subject Contribution
For a hypothetical subject, suppose \(R=\begin{pmatrix}1&0\\0&4\end{pmatrix}\), \(\Omega=1\), \(J=\begin{pmatrix}1\\2\end{pmatrix}\), \(\hat\eta=0.20\), and \(r(\hat\eta)=\begin{pmatrix}0.08\\0.24\end{pmatrix}\). These local quantities satisfy the mode score condition: \(J'R^{-1}r=0.08+2(0.24/4)=0.20=\Omega^{-1}\hat\eta\).
The covariance of the integrated local model is \(S=R+J\Omega J'=\begin{pmatrix}2&2\\2&8\end{pmatrix}\), with determinant \(|S|=16-4=12\). Its inverse is \(S^{-1}=\tfrac1{12}\begin{pmatrix}8&-2\\-2&2\end{pmatrix}\). The residual from the FOCE mean is \(z=r(\hat\eta)+J\hat\eta=\begin{pmatrix}0.28\\0.64\end{pmatrix}\).
Substitution into the quadratic term gives \(z'S^{-1}z=[8(0.28^2)-4(0.28)(0.64)+2(0.64^2)]/12 =(0.6272-0.7168+0.8192)/12=0.0608\). Thus \(\log\tilde L^{\mathrm{FOCE}}=-\tfrac12[2\log(2\pi)+\log(12)+0.0608]\approx-3.111\). Both the determinant and the quadratic term use the full two-observation covariance, including its off-diagonal entries.
Worked Example: Comparing Linearization at the Mode and at Zero
Use the scalar data and model from the first example to compare the FOCE expansion point with the population mean \(\eta=0\). At zero, \(f(0)=1\) and \(f'(0)=1.5\). Linearizing there gives an approximate marginal mean of \(1\) and variance \(4+1.5^2=6.25\). The observed difference from this mean is \(1.9375-1=0.9375\).
The corresponding first-order log contribution is \(-\tfrac12[\log(2\pi)+\log(6.25)+0.9375^2/6.25]\). Since \(0.9375^2/6.25=0.140625\), this is \(-\tfrac12(1.837877+1.832581+0.140625)=-1.906\). The FOCE value from the mode expansion was \(-2.021\). The two approximations differ because they linearize the nonlinear prediction at different points. This one calculation does not establish which approximation is more accurate for other data or models.
Assumptions, Diagnostics, and Alternatives
The formula requires a valid Gaussian random-effect distribution, positive-definite \(\Omega\), positive-definite residual covariance \(R_i\), and a differentiable prediction at the conditional mode. The residual covariance must not change with \(\eta_i\) for the stated derivation. If it does, the conditional density’s covariance also varies across random effects, so the derivation above no longer applies as written. The mode should be an interior minimum of the penalized conditional objective, not merely a point where an optimizer stopped.
In practice, check the conditional-mode gradient, the objective value, and convergence from reasonable starting values. The score residual \(J_i'R_i^{-1}r_i(\hat\eta_i)-\Omega^{-1}\hat\eta_i\) should be small relative to the scales of its components. Check that \(S_i=R_i+J_i\Omega J_i'\) is numerically positive definite and that its determinant and linear solves are computed stably. A poorly conditioned covariance or markedly different modes from different starts can signal weak information, numerical trouble, or an unsuitable local approximation.
A further diagnostic is to compare FOCE with a more direct numerical evaluation of the subject integral for selected subjects, especially when predictions are strongly curved near the conditional modes. Such comparisons can reveal approximation sensitivity; they do not turn one method into a universal benchmark in every computational setting. For modest random-effect dimension, numerical integration may be a useful check but can become expensive as dimension increases.
| Approach | Expansion or calculation | Practical distinction |
|---|---|---|
| FOCE | Linearizes predictions at each subject’s conditional mode, then integrates the local Gaussian model. | Uses subject-specific derivatives and accounts for the mode’s displacement from zero. |
| First-order approximation | Linearizes predictions at the population mean random effect, zero. | Simpler, but can be sensitive when the conditional modes are far from zero or predictions are curved. |
| Numerical integration | Numerically evaluates the nonlinear subject integral without replacing the prediction by this local linear expansion. | Can serve as a check, but computational cost can grow rapidly with random-effect dimension. |
A common mistake is to calculate the conditional likelihood at \(\hat\eta_i\) and call it the FOCE contribution. That omits integration over the population random-effect distribution and misses the FOCE covariance \(R_i+J_i\Omega J_i'\), along with its determinant. Another is to use \(f_i(\hat\eta_i)\) as the integrated mean. The local model has random mean \(f_i(\hat\eta_i)+J_i(\eta_i-\hat\eta_i)\); averaging over the population distribution, whose mean is zero, gives \(f_i(\hat\eta_i)-J_i\hat\eta_i\), not \(f_i(\hat\eta_i)\).
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- For a scalar random effect, \(J=3\), \(r(\hat\eta)=0.2\), \(R=6\), \(\Omega=2\), and \(\hat\eta=0.2\). Does the conditional-mode score condition hold?
- If \(f(\hat\eta)=5\), \(J=2\), and \(\hat\eta=0.4\), what is the FOCE mean of the integrated local model?
- For \(R=3\), \(J=2\), and \(\Omega=0.5\) in the scalar case, what is the FOCE variance?
- Why is the conditional likelihood evaluated at the mode not, by itself, the FOCE subject contribution?
Solutions. (1) \(Jr/R=3(0.2)/6=0.1\), while \(\Omega^{-1}\hat\eta=(1/2)(0.2)=0.1\), so the score condition holds. (2) The integrated mean is \(f(\hat\eta)-J\hat\eta=5-2(0.4)=4.2\). (3) \(S=R+J^2\Omega=3+2^2(0.5)=5\). (4) FOCE integrates the local Gaussian model over the random-effect distribution; evaluating the conditional likelihood at one value omits that integration and its resulting covariance and normalization terms.