Empirical Bayes Modes and Means in Nonlinear Models

Why Two Individual Estimates Can Differ

In a nonlinear mixed-effects model, an individual's observations update the population distribution for that individual's random effects. The resulting conditional distribution contains more information than a single empirical Bayes estimate (EBE). When one number is needed, two common summaries are the conditional mode and the conditional mean. They coincide for a symmetric Gaussian conditional distribution, but need not coincide when the likelihood is nonlinear and the conditional distribution is skewed or multimodal.

The distinction matters for individual predictions and interpretation. A mode identifies a point where the conditional density is highest; a mean averages over the whole conditional distribution. With skewness, the mean can be pulled toward a long tail. With multiple modes, it can fall between regions of high density, where the conditional density is low.

Definition: Let \(\boldsymbol{\eta}_i\) denote the random-effects vector for individual \(i\), and let \(\mathbf{y}_i\) denote that individual's observations. Conditional on fitted population parameters \(\hat{\boldsymbol{\theta}}\) and the model, the conditional mean is \(E(\boldsymbol{\eta}_i\mid\mathbf{y}_i,\hat{\boldsymbol{\theta}})\). The conditional mode is a value of \(\boldsymbol{\eta}_i\) that maximizes the conditional density \(p(\boldsymbol{\eta}_i\mid\mathbf{y}_i,\hat{\boldsymbol{\theta}})\). These are plug-in empirical Bayes summaries because the population parameters are treated as fitted values.

Conditional Distribution and Interpretation

For a model with prior random-effects density \(p(\boldsymbol{\eta}_i\mid\hat{\boldsymbol{\theta}})\) and individual likelihood \(p(\mathbf{y}_i\mid\boldsymbol{\eta}_i,\hat{\boldsymbol{\theta}})\), the conditional density is

p(\boldsymbol{\eta}_i\mid\mathbf{y}_i,\hat{\boldsymbol{\theta}}) = \frac{ p(\mathbf{y}_i\mid\boldsymbol{\eta}_i,\hat{\boldsymbol{\theta}}) p(\boldsymbol{\eta}_i\mid\hat{\boldsymbol{\theta}}) }{ \int p(\mathbf{y}_i\mid\mathbf{u},\hat{\boldsymbol{\theta}}) p(\mathbf{u}\mid\hat{\boldsymbol{\theta}})\,d\mathbf{u} }.

The denominator makes the density integrate to one over the random-effects space. In a nonlinear model, the likelihood as a function of \(\boldsymbol{\eta}_i\) can be asymmetric or have multiple high-density regions, even when the population random-effects distribution is Gaussian. Thus a Gaussian prior does not guarantee a Gaussian conditional distribution.

For a vector of random effects, the conditional mean is a vector of componentwise expectations. A joint conditional mode is a vector maximizing the full joint density. It is not generally the same as taking the mode of each separate marginal conditional density; the joint and marginal summaries answer different questions.

What Each Summary Represents

The conditional mean is the optimal point estimate under squared-error loss. To see why, let \(a\) be a proposed estimate of a scalar \(\eta_i\). Given the observations, its expected squared error is

E\{(\eta_i-a)^2\mid\mathbf{y}_i\} = \operatorname{Var}(\eta_i\mid\mathbf{y}_i) + \{E(\eta_i\mid\mathbf{y}_i)-a\}^2.

The variance term does not depend on \(a\), and the second term is minimized at \(a=E(\eta_i\mid\mathbf{y}_i)\). By contrast, a conditional mode maximizes density; it is a maximum a posteriori estimate. It is not generally the estimator that minimizes squared error.

The mode can also depend on parameterization. If a parameter is transformed nonlinearly, the density in the transformed coordinate includes a Jacobian factor, which can change the location of the mode. A mean of a transformed parameter is likewise not generally the transformation of the mean: for example, \(E\{\exp(\eta_i)\mid\mathbf{y}_i\}\) need not equal \(\exp\{E(\eta_i\mid\mathbf{y}_i)\}\). State the scale and summary being reported.

Worked Example: A Skewed Conditional Distribution After One Observation

Suppose a hypothetical individual has a scalar random effect \(\eta_i\sim N(0,1)\). Conditional on \(\eta_i\), a binary observation \(Y_i\) has probability \(P(Y_i=1\mid\eta_i)=\Phi(2\eta_i)\), where \(\Phi\) is the standard normal cumulative distribution function. The observed value is \(Y_i=1\). Find the conditional mean and mode.

The likelihood contribution for the observed outcome is \(\Phi(2\eta_i)\). By symmetry of the prior and the identity \(\Phi(x)+\Phi(-x)=1\), the marginal probability of \(Y_i=1\) is \(1/2\). The conditional density is therefore

f(\eta_i\mid Y_i=1)=2\phi(\eta_i)\Phi(2\eta_i),

where \(\phi\) is the standard normal density. This is a skewed density, not a Gaussian density. Its conditional mean is the mean of a skew-normal distribution with shape parameter 2:

E(\eta_i\mid Y_i=1) = \frac{2}{\sqrt{5}}\sqrt{\frac{2}{\pi}} \approx 0.7136.

To find the mode, differentiate the log density and set the derivative to zero:

\frac{d}{d\eta_i}\log f(\eta_i\mid Y_i=1) = -\eta_i+\frac{2\phi(2\eta_i)}{\Phi(2\eta_i)} =0.

Solving this one-dimensional equation numerically gives a mode of approximately \(0.531\). Thus the conditional mean is about \(0.714\), while the mode is about \(0.531\). The mean is farther to the right because it reflects the entire right-skewed distribution, not just its highest point. These are estimates on the \(\eta_i\) scale, conditional on this individual's observed binary outcome and the stated model.

Multimodality: A Mean Between Plausible Values

Nonlinear observation models can map quite different random-effects values to similar predicted observations. The conditional density may then have multiple peaks. In that setting, reporting only a mode can hide another plausible region, while the mean can lie in a low-density gap between peaks.

Worked Example: A Multimodal Conditional Density

Consider a hypothetical scalar conditional density represented by a mixture

f(\eta_i\mid\mathbf{y}_i) = 0.7\,N(\eta_i;-2,0.25)+0.3\,N(\eta_i;2,0.25),

where each normal component has variance \(0.25\) (standard deviation \(0.5\)). This form can approximate a conditional distribution when two separated random-effects values explain the observations well. The conditional mean follows directly from the mixture weights:

E(\eta_i\mid\mathbf{y}_i) = 0.7(-2)+0.3(2) = -1.4+0.6 = -0.800.

The component peaks are separated by 4 units, or 8 component standard deviations. The higher, left peak dominates the density near \(-2\), so the global mode is essentially \(-2\). More precisely, the small right component shifts the left mode only about \(2.2\times10^{-14}\) units to the right; its location rounds to \(-2.000\). The right component also produces a local mode near \(2\), but its peak is lower. The mean, \(-0.800\), lies between the two peaks and is not itself representative of either high-density region. In this case, show the conditional distribution or both modes rather than presenting the mean as if it described a typical individual value.

A Gaussian Benchmark and Computation

A Gaussian conditional distribution provides a useful benchmark. Its mean and mode coincide, so a substantial difference between computed summaries should prompt checks of the conditional density, numerical integration, and optimization.

Worked Example: When Mean and Mode Agree

For comparison, suppose a hypothetical individual has prior \(\eta_i\sim N(0,1)\), and an observation \(Z_i=0.60\) follows \(Z_i\mid\eta_i\sim N(\eta_i,0.25)\). The prior and likelihood are Gaussian, so the conditional distribution is Gaussian. Its variance and mean are

V_i=\left(1+\frac{1}{0.25}\right)^{-1}=0.200, \qquad m_i=V_i\left(\frac{0}{1}+\frac{0.60}{0.25}\right) =0.200(2.40)=0.480.

A Gaussian density is symmetric and maximized at its mean. Therefore both the conditional mean and conditional mode are \(0.480\) on the \(\eta_i\) scale. This agreement follows from the Gaussian form and should not be assumed in a nonlinear model.

In practice, the conditional mean is computed by integrating \(\boldsymbol{\eta}_i\) against its conditional density. Numerical quadrature or Monte Carlo integration can approximate this expectation. A conditional mode is obtained by maximizing the log conditional density, typically with numerical optimization. These computations require attention to the whole conditional distribution: a local optimizer may find one mode but miss another, and a numerical integration method may poorly capture a narrow peak or a long tail.

Useful checks include plotting one-dimensional conditional densities where possible, using multiple optimization starting points, checking whether the reported mode has the largest density among candidate modes, and verifying that integration is stable to more nodes or draws. In higher dimensions, inspect marginal summaries and relevant joint structure; a single joint-mode calculation does not reveal all marginal behavior.

Choosing and Reporting an Estimate

Choose the summary to match the task. If the goal is a point prediction judged by squared error on the random-effects scale, the conditional mean is the direct choice. If the goal is the most probable local value under the fitted conditional density, the mode is relevant, but multimodality means there may be more than one local mode. For predictions on the response scale, calculate or approximate the conditional expectation of the response itself when that is the target; transforming an EBE may not give that expectation.

Neither summary removes uncertainty. Both are conditional on the individual's observations and on fitted population parameters; they do not, by themselves, report the spread or shape of the conditional distribution. Sparse or weakly informative individual data can leave the conditional distribution close to the population distribution, while informative data can make it more concentrated or reshape it. Do not interpret one EBE as a directly observed individual parameter.

Common Mistakes and Limitations

Calling the mean and mode interchangeable. They coincide for a Gaussian conditional density, but skewness and multiple peaks can separate them substantially.
Calling the mode the most likely exact value. For a continuous variable, the probability of any exact point is zero. The mode is the location of maximum density, not a positive probability assigned to that point.
Treating a multimodal mean as a typical value. A mean between separated peaks can have low conditional density. Inspect the distribution before assigning it a typical-individual interpretation.
Assuming an optimizer found the global mode. Nonlinear log densities can have multiple local maxima. Compare starts and evaluate the density at candidate solutions.
Ignoring the parameter scale. Means and modes need not transform in the same way. State whether the estimate is for \(\eta_i\), a transformed individual parameter, or a response prediction.
Reporting a point estimate as complete individual information. The conditional distribution can be broad, skewed, or multimodal; summarize its uncertainty or shape when relevant.

Key takeaway: The conditional mean averages over the full individual-specific distribution and minimizes conditional squared-error loss. The conditional mode identifies a density maximum. In nonlinear models, compare their values and inspect distribution shape before interpreting either as a representative individual effect.

Check Your Understanding

Answer each question, then compare with the brief solutions.

What feature of a conditional distribution guarantees that its mean and mode coincide?
For the mixture in the worked example, why is the conditional mean not a good description of either high-density component?
What point estimate minimizes conditional expected squared error, and why?
A nonlinear optimizer reports one conditional mode. Name one check that can help detect a missed mode.

Solutions. (1) A Gaussian conditional density is symmetric and has its mean at its maximum. (2) Its component peaks are near \(-2\) and \(2\), while the mean is \(-0.800\), between the peaks in a low-density region. (3) The conditional mean, because expected squared error equals the conditional variance plus the squared distance from the proposed estimate to that mean. (4) Repeat the optimization from multiple starting values and compare the log densities at the candidate solutions.