Conditional Modes in FOCE

Why an Individual Conditional Mode Matters

In a nonlinear mixed-effects model, each individual's observations provide information about that individual's random effects, while the population distribution supplies a constraint on how far those effects are expected to lie from the population mean. The individual conditional mode balances these two sources of information: it is the random-effect value that maximizes the individual's observation likelihood multiplied by the random-effect density.

This mode is central to FOCE because the individual prediction is later linearized around it. Here the focus is the mode itself: its definition, its score equation, and its dependence on the population parameters. A mode is a feature of a conditional density, not a direct observation of an individual's true random effect.

Definition: For individual $i$, let $y_i$ denote the vector of observations, $\boldsymbol{\eta}_i$ the vector of random effects, and $\boldsymbol{\theta}$ the population parameters. The individual conditional mode, also called the conditional mode of the random effects, is $$ \widehat{\boldsymbol{\eta}}_i(\boldsymbol{\theta}) = \underset{\boldsymbol{\eta}_i}{\operatorname{arg\,max}}\; p(y_i\mid\boldsymbol{\eta}_i,\boldsymbol{\theta})\, p(\boldsymbol{\eta}_i\mid\boldsymbol{\theta}). $$ Equivalently, it maximizes the log of this product. If the random effects have a proper density, the product is proportional to the conditional density of $\boldsymbol{\eta}_i$ given $y_i$ and $\boldsymbol{\theta}$; its maximizer is the conditional mode.

Model and Score Equation

Write the observation model as $y_{ij}=f_{ij}(\boldsymbol{\theta},\boldsymbol{\eta}_i)+\epsilon_{ij}$, where $f_{ij}$ is the model-predicted observation at the $j$th record for individual $i$. A common specification has $\boldsymbol{\eta}_i\sim N(\mathbf{0},\Omega)$. For the score equation below, suppose residual errors are conditionally independent normal variables with mean zero and constant variance $\sigma^2$, and are independent of the random effects. These assumptions make the individual log likelihood a sum of Gaussian contributions. Other residual models can be used, but their full log likelihood—including any parameter-dependent variance terms—must then be differentiated instead.

For a scalar random effect $\eta_i$, the individual log likelihood plus log random-effect density, after dropping constants that do not depend on $\eta_i$, is

Q_i(\eta_i) = -\frac{1}{2\sigma^2}\sum_j \left[y_{ij}-f_{ij}(\boldsymbol{\theta},\eta_i)\right]^2 -\frac{\eta_i^2}{2\omega}, \qquad \eta_i\sim N(0,\omega).

A differentiable interior mode must satisfy $Q_i'(\eta_i)=0$. Differentiating each squared residual gives the score equation below. It shows the balance explicitly: the observation contribution pulls the mode toward values that improve the predictions, while the prior-density contribution pulls it toward zero.

Q_i'(\eta_i) = \frac{1}{\sigma^2}\sum_j \left[y_{ij}-f_{ij}(\boldsymbol{\theta},\eta_i)\right] \frac{\partial f_{ij}(\boldsymbol{\theta},\eta_i)}{\partial\eta_i} -\frac{\eta_i}{\omega} =0.

For a vector of random effects with a normal distribution $N(\mathbf{0},\Omega)$, and the same constant-variance Gaussian residual assumptions, the corresponding vector equation is

\nabla_{\boldsymbol{\eta}_i}\log p(y_i\mid\boldsymbol{\eta}_i,\boldsymbol{\theta}) -\Omega^{-1}\boldsymbol{\eta}_i =\mathbf{0}.

The equation is generally nonlinear. Solving it finds stationary points, not automatically the global maximum. A candidate maximum should have a suitable objective value and negative curvature locally; optimization should also be checked for convergence and possible alternative modes.

Worked Example: A Mode From Two Observations

Consider hypothetical observations $y_{i1}=12$ and $y_{i2}=14$ from a model $f_{ij}(\theta,\eta_i)=\theta+\eta_i$, with $\theta=10$, residual variance $\sigma^2=4$, and random-effect variance $\omega=1$. Assume the independent Gaussian residual and normal random-effect conditions stated above.

The score equation is $$ \frac{(12-10-\eta_i)+(14-10-\eta_i)}{4}-\eta_i=0. $$ Substituting and collecting terms gives $$ \frac{6-2\eta_i}{4}-\eta_i =1.5-1.5\eta_i=0, $$ so the conditional mode is $\widehat{\eta}_i=1$. The second derivative is $-2/4-1=-1.5$, which is negative, so this stationary point is a local maximum. In this quadratic log objective it is also the unique maximum. The two observations jointly pull the mode above zero; the normal random-effect density penalizes larger departures.

Nonlinear Predictions and Numerical Modes

For a nonlinear prediction, the score equation usually has no closed-form solution. A numerical optimizer searches for a value that maximizes $Q_i$, using its gradient and, depending on the method, curvature information. A useful calculation checks both the score near zero and the curvature near the proposed mode. The following example uses one hypothetical observation so that the score equation can be inspected directly.

Worked Example: A Nonlinear Exponential Prediction

Suppose $y_i=12$ and $f_i(\eta_i)=10e^{\eta_i}$, with additive residual variance $\sigma^2=1$ and random-effect variance $\omega=1/16$. The random effect has distribution $N(0,1/16)$, and the residual is independent normal error with constant variance. The score equation is $$ (12-10e^{\eta_i})\,10e^{\eta_i}-16\eta_i=0. $$ At $\eta_i=0$, the score is $20$, so zero is not the mode. Solving the equation numerically gives $\widehat{\eta}_i\approx0.163624$. For example, at $\eta_i=0.1635$, $10e^{0.1635}\approx11.7763$. The left side of the rearranged equation is then $$ (12-10e^{0.1635})\,10e^{0.1635}\approx2.6349, $$ while the right side is $16(0.1635)=2.616$. The left side is still larger, so this value is slightly below the root. Near $\eta_i=0.163624$, the prediction is approximately $11.7777$, and both sides are approximately $2.618$.

The derivative of the score is $10e^{\eta_i}(12-20e^{\eta_i})-16$, which is approximately $-152$ at the numerical root. The log objective therefore curves downward there, supporting a local maximum. For an overall check, the objective with constants omitted is $-\tfrac12(12-10e^{\eta_i})^2-8\eta_i^2$. It is approximately $-0.239$ at the candidate mode and $-2$ at zero. Both tails tend to negative infinity, since the prior penalty does so as $\eta_i\to-\infty$ and the squared residual does so as $\eta_i\to\infty$. In practical optimization, convergence and alternative stationary points still need checking rather than assuming a root is the global maximum.

Dependence on Population Parameters

The notation $\widehat{\boldsymbol{\eta}}_i(\boldsymbol{\theta})$ emphasizes that the mode depends on the population parameters. The fixed effects determine the predictions, and the random-effect covariance determines the strength and shape of the density penalty. Even with the same observations, changing either can change the maximizing random effect. The mode is therefore conditional on the specified population model and parameter values.

Worked Example: Changing the Fixed Effect or Random-Effect Variance

Use one hypothetical observation $y_i=12$ and the linear prediction $f_i(\theta,\eta_i)=\theta+\eta_i$. Assume normal random effects with variance $\omega$ and independent normal residual error with variance $\sigma^2=4$. The score equation is $$ \frac{12-\theta-\eta_i}{4}-\frac{\eta_i}{\omega}=0, \qquad \widehat{\eta}_i=\frac{\omega(12-\theta)}{4+\omega}. $$ With $\theta=10$ and $\omega=1$, substitution gives $\widehat{\eta}_i=1(12-10)/(4+1)=0.4$. Holding the data and variances fixed but changing the fixed effect to $\theta=11$ gives $\widehat{\eta}_i=1(12-11)/5=0.2$. Returning to $\theta=10$ but increasing the random-effect variance to $\omega=4$ gives $\widehat{\eta}_i=4(12-10)/(4+4)=1$. Each result is a maximum: the second derivative is $-1/4-1/\omega<0$. A larger random-effect variance weakens the penalty for moving away from zero; a changed fixed effect changes how much of the observed difference the random effect must account for.

Optimization Checks, Comparisons, and Limitations

A reliable calculation checks that the optimizer has actually found a mode of the intended conditional objective. Useful checks include:

Score: at an interior mode, the gradient should be close to zero. A large gradient suggests incomplete convergence or a coding error.
Curvature: the Hessian of the log objective should be negative definite at a strict local maximum. A singular or nearly singular Hessian can indicate weak information or numerical instability.
Objective and starts: compare objective values from multiple starting values. Distinct solutions can signal multiple stationary points; report or investigate the highest valid maximum rather than relying on the first root found.
Model specification: confirm that the likelihood includes all relevant residual-model terms. If residual variance depends on $\eta_i$, omitting its log-variance contribution gives the wrong score.

The conditional mode differs from both the marginal likelihood and a first-order expansion point. The marginal likelihood integrates the observation likelihood over the random-effect distribution; maximizing the product at one $\boldsymbol{\eta}_i$ does not perform that integration. In a linear-Gaussian model, the conditional distribution is Gaussian, so its mode and mean coincide. In a nonlinear model, the conditional distribution need not be symmetric, and its mode need not equal its mean. The first-order approximation expands predictions at the population mean random effect, usually zero; the conditional mode instead uses information from the individual's observations. FOCE uses this mode as its expansion point, but finding the mode alone is not the full FOCE likelihood approximation.

Conditional modes also have limitations. With sparse or uninformative observations, the mode can be strongly pulled toward zero by the random-effect density. It can be sensitive to the assumed covariance and residual model, and a mode is not a measure of uncertainty by itself. These limitations matter when interpreting individual estimates: a value such as $\widehat{\eta}_i$ is a conditional summary under fitted population parameters, not proof that the individual's unobserved effect equals that value.

Key takeaway: The individual conditional mode maximizes the product of an individual's observation likelihood and random-effect density. Its score balances data fit against the population penalty, so it depends on both the individual's observations and the population parameters.

Check Your Understanding

Answer each question, then compare your reasoning with the short solutions.

State the definition of an individual conditional mode and identify the two factors being maximized.
For scalar normal random effects with variance $\omega$ and independent, constant-variance Gaussian residuals, what term in the score represents the random-effect penalty?
In the exponential example, why is $\eta_i=0$ not the mode?
For the linear model in the final worked example, explain why increasing $\omega$ from 1 to 4 changes the mode from 0.4 to 1 when $\theta=10$.
Give two checks that help distinguish a valid local maximum from an optimizer that has not converged properly.

Solutions. (1) It is the value of $\boldsymbol{\eta}_i$ maximizing $p(y_i\mid\boldsymbol{\eta}_i,\boldsymbol{\theta})p(\boldsymbol{\eta}_i\mid\boldsymbol{\theta})$: the observation likelihood and random-effect density. (2) The score includes $-\eta_i/\omega$. (3) The score at zero is $20$, not zero, so zero is not even a stationary point. (4) A larger $\omega$ weakens the density penalty, allowing the observation to pull the mode farther from zero; substitution gives $4(12-10)/(4+4)=1$. (5) Check that the gradient is near zero and that the log-objective Hessian has negative curvature; also compare objective values from multiple starting values.

What You'll Learn