Finding the Conditional Mode by Newton Iteration

Why Find a Conditional Mode?

In a generalized linear mixed model, integrating over random effects produces a likelihood contribution that is generally an integral rather than a closed-form expression. A useful first step in numerical integration is to locate where the integrand is largest. For a cluster with one random intercept, this means finding the value of that intercept that maximizes the log of the cluster’s likelihood contribution multiplied by its random-effects density.

This tutorial develops Newton iteration for that scalar maximization. The iteration uses the first and second derivatives of the log integrand: the score indicates the direction of local increase, while the second derivative describes local curvature. Finding a mode is a numerical optimization task; the iteration is not, by itself, an approximation to the integral.

Definition: Let \(b_j\) be the random intercept for cluster \(j\), and let \(y_j\) denote that cluster’s outcomes. The conditional log integrand is \(h_j(b_j)=\log f(y_j\mid b_j;\beta)+\log f(b_j;\sigma_b^2)\), where terms not depending on \(b_j\) may be omitted. A conditional mode is a maximizer \(\hat b_j\) of \(h_j\). In the scalar case, an interior mode satisfies \(h_j'(\hat b_j)=0\); it is a strict local maximum if also \(h_j''(\hat b_j)<0\).

Deriving the Newton Update

Suppress the cluster subscript temporarily and write \(h(b)\) for the log integrand. Near a current value \(b^{(t)}\), its second-order Taylor approximation is

h(b^{(t)}+\delta)\approx h(b^{(t)})+h'(b^{(t)})\delta +\frac{1}{2}h''(b^{(t)})\delta^2.

The maximum of this quadratic approximation occurs where its derivative with respect to \(\delta\) is zero. Solving \(h'(b^{(t)})+h''(b^{(t)})\delta=0\) gives the Newton displacement and update:

\delta^{(t)}=-\frac{h'(b^{(t)})}{h''(b^{(t)})}, \qquad b^{(t+1)}=b^{(t)}-\frac{h'(b^{(t)})}{h''(b^{(t)})}.

When \(h''(b^{(t)})<0\), the update can also be written \(b^{(t+1)}=b^{(t)}+h'(b^{(t)})/I(b^{(t)})\), where \(I(b^{(t)})=-h''(b^{(t)})>0\) is the local curvature magnitude. If the score is positive, this form moves upward; if the score is negative, it moves downward. The approximation is local, however: a large step can leave the region where the quadratic expansion is informative.

Binomial Random-Intercept Example

Suppose a hypothetical cluster has \(y=3\) successes in \(m=4\) trials. Given its random intercept \(b\), let the success probability be \(\operatorname{logit}^{-1}(-0.5+b)\), and assume \(b\sim N(0,1)\). Conditional on \(b\), the four outcomes are independent Bernoulli observations. The Gaussian density is proper, and the binomial likelihood is positive for finite \(b\), so the log integrand is defined for every finite starting value.

After dropping terms independent of \(b\), the conditional log integrand, its score, and its second derivative are

h(b)=3(-0.5+b)-4\log\{1+\exp(-0.5+b)\}-\frac{b^2}{2},

h'(b)=3-4p(b)-b,\qquad h''(b)=-4p(b)\{1-p(b)\}-1, \qquad p(b)=\frac{1}{1+\exp\{-(b-0.5)\}}.

The second derivative is negative for every finite \(b\): \(p(b)\{1-p(b)\}>0\), so \(h''(b)<-1\). Thus the log integrand is strictly concave and has at most one stationary point. It also tends to negative infinity in both tails because of the Gaussian term, so a maximum exists and is unique. These facts justify seeking the mode by solving the score equation.

Worked Example: Newton Iteration for a Binomial Cluster

Start at \(b^{(0)}=0\). The success probability is \(p(0)=1/(1+\exp(0.5))\approx0.377541\). Therefore \(h'(0)=3-4(0.377541)-0\approx1.48984\), and \(h''(0)=-4(0.377541)(0.622459)-1\approx-1.94001\).

The first Newton update is \(b^{(1)}=0-1.48984/(-1.94001)\approx0.76796\). Evaluating the same score and curvature at this value gives a second update of approximately \(b^{(2)}=0.750634\). At that point, \(b^{(2)}-0.5=0.250634\), so \(p(b^{(2)})\approx0.56233\). The score is \(h'(b^{(2)})=3-4(0.56233)-0.750634\approx3.583\times10^{-5}\).

The curvature magnitude there is approximately \(I(b^{(2)})=1+4(0.56233)(0.43767)\approx1.98445\). Hence \(b^{(3)}=b^{(2)}+h'(b^{(2)})/I(b^{(2)}) \approx0.750634+0.00003583/1.98445 \approx0.750652\). The tiny score and small change in \(b\) indicate that the unique mode is about \(0.7507\). The displayed digits are numerical approximations, not exact values.

Poisson and Gaussian Examples

The same procedure applies when the conditional likelihood is Poisson. For a Poisson log-link model with a random intercept, the conditional mean contains \(\exp(b)\); the score therefore contains the observed count minus the fitted mean. The Gaussian random-effect density contributes a quadratic penalty and a linear term to the score.

Worked Example: A Poisson Cluster with Two Observations

Consider hypothetical counts \(y_1=2\) and \(y_2=3\). Suppose the two conditional Poisson means at \(b=0\) are 1 and 3, respectively, and both are multiplied by \(\exp(b)\) when the random intercept changes. Let \(b\sim N(0,1)\). The observations are conditionally independent given \(b\). Their log integrand, omitting constants, is \(h(b)=5b-4\exp(b)-b^2/2\), so \(h'(b)=5-4\exp(b)-b\) and \(h''(b)=-4\exp(b)-1<0\). Strict concavity and the Gaussian tail penalty ensure a unique mode.

At \(b^{(0)}=0\), the score is \(5-4-0=1\), and the second derivative is \(-5\). Thus \(b^{(1)}=0-1/(-5)=0.2000\). At \(b=0.2000\), \(\exp(b)\approx1.221403\), giving score \(5-4(1.221403)-0.2000\approx-0.08561\) and curvature \(-4(1.221403)-1\approx-5.88561\). The next update is \(b^{(2)}=0.2000-(-0.08561)/(-5.88561)\approx0.18545\).

At \(b=0.18545\), \(\exp(b)\approx1.20376\), so the score is approximately \(5-4(1.20376)-0.18545\approx-0.00049\). One further update is about \(0.18545-(-0.00049)/(-5.815)\approx0.18537\). The iterates are converging to a mode near \(0.1854\); the residual score and the size of the update are both small.

Worked Example: A Gaussian Random-Intercept Model

For contrast, suppose two hypothetical responses, after subtracting their fixed-effect means, are 2 and 4. Let their independent residual errors have variance \(\tau^2=1\), and let \(b\sim N(0,2)\). The conditional log integrand, omitting constants, is \(h(b)=-\{(2-b)^2+(4-b)^2\}/2-b^2/4\). Its score is \(h'(b)=(2-b)+(4-b)-b/2=6-2.5b\), and its curvature is \(h''(b)=-2.5\). This curvature is negative everywhere, so the unique mode solves \(6-2.5b=0\), giving \(\hat b=2.4\).

From \(b^{(0)}=0\), Newton’s update gives \(b^{(1)}=0-6/(-2.5)=2.4\), the exact mode in one step. This happens because the log integrand is exactly quadratic: its second-order Taylor approximation is the function itself. The one-step result is not a general property of GLMMs with non-Gaussian outcomes.

Checking Convergence and Recognizing Problems

A numerical stopping rule should check more than whether consecutive iterates look similar. A common practice is to monitor the absolute score \(|h'(b^{(t)})|\), the absolute step \(|b^{(t+1)}-b^{(t)}|\), and the change in \(h\). A small step alone can be misleading if the curvature is extremely large; a small score alone can also be insufficient if the curvature is close to zero and the mode is poorly identified. Reported convergence tolerances depend on the scale of the problem.

Check the sign of curvature. At a candidate interior mode, \(h''(b)\) should be negative. A positive value indicates a local minimum for the quadratic approximation, not a local maximum. Near-zero curvature can make the Newton displacement unstable.
Track the log integrand. In a well-behaved concave problem, iterations should head toward its maximum. A large step that sharply lowers \(h\) is a warning that the local approximation has been used too far from its useful range.
Check arithmetic and derivatives. A sign error in \(h''\) reverses the update direction. Compare an analytic derivative with a finite-difference check when the implemented score is in doubt.
Use the model’s structure. In the binomial and Poisson examples, strict concavity rules out multiple stationary points. Without concavity, a converged iteration may find only a local mode, and different starting values may lead to different answers.

If iterations diverge, oscillate, or encounter nearly zero curvature, the computed mode should not be accepted merely because software stopped. Possible remedies include revisiting the derivatives and starting value, or using an optimization method with safeguards. Fisher scoring and damped steps provide alternative strategies for mode finding; they modify how the update is constructed or applied rather than changing the definition of the conditional mode.

How Newton Iteration Fits Among Alternatives

Newton iteration is attractive when the score and curvature are available and informative: near a well-behaved mode it can converge quickly. Its main weakness is that it relies on a local quadratic approximation. A poor starting value, nonconcavity, or weak curvature can produce an unhelpful step. Direct numerical optimization instead seeks a maximum of \(h(b)\) through an optimizer’s chosen search strategy; it can be useful when a Newton step is unreliable, but it still requires checks that a maximum was found.

For a scalar random intercept, solving \(h'(b)=0\) with a root-finding method is another option, provided the relevant root corresponds to a maximum. When the log integrand is strictly concave, its score is decreasing and any root is the unique mode. Newton’s method uses both the score and its derivative to find that root. In a multidimensional random-effects problem the analogous update involves the score vector and Hessian matrix, but the scalar examples here do not establish that a multivariate iteration will be stable.

Mode finding is also only one part of an integration calculation. The mode identifies a location; it does not supply the curvature correction or normalization needed to approximate the integral. Those are separate steps. In addition, a single mode may be inadequate if the integrand has multiple important peaks or is strongly nonquadratic near its maximum.

Common Mistake

Common pitfall: Using \(b^{(t+1)}=b^{(t)}+h'(b^{(t)})/h''(b^{(t)})\) because the score is positive, reasoning that a positive score means “move up.” This has the wrong sign when \(h''<0\). The Newton update is \(b^{(t+1)}=b^{(t)}-h'(b^{(t)})/h''(b^{(t)})\); equivalently, it is \(b^{(t)}+h'(b^{(t)})/I(b^{(t)})\) when \(I=-h''>0\). For a positive score at a concave point, the latter correctly increases \(b\).

Check Your Understanding

Try each question, then check the solutions below.

For a scalar log integrand \(h(b)\), write the Newton update for finding a stationary point. What curvature condition identifies a strict local maximum?
For \(h(b)=2b-3\exp(b)-b^2/2\), calculate the first Newton update starting at \(b^{(0)}=0\).
Why does Newton iteration find the Gaussian example’s mode in one step?
Name two quantities, in addition to the iterate itself, that should be checked when assessing convergence.
If \(h''(b)>0\) at a stationary point, what does that imply in the scalar case?

Solutions. (1) \(b^{(t+1)}=b^{(t)}-h'(b^{(t)})/h''(b^{(t)})\); a strict local maximum requires \(h''(b)<0\). (2) The score is \(h'(b)=2-3\exp(b)-b\) and \(h''(b)=-3\exp(b)-1\). At zero these are \(-1\) and \(-4\), respectively, so \(b^{(1)}=0-(-1)/(-4)=-0.2500\). (3) Its log integrand is quadratic, so the Taylor approximation used by Newton is exact and the update reaches the quadratic’s vertex. (4) Check, for example, the score magnitude and the change in log integrand, as well as the step size. (5) The point is a strict local minimum under the usual second-derivative condition, not a strict local maximum.