Why Binary GLMM Shrinkage Is Nonlinear

Why Binary Group Effects Do Not Shrink by a Fixed Weight

In a Gaussian random-intercept model, a group effect can often be written as a fixed fraction of a group-specific estimate. That simple linear rule does not carry over to a binary GLMM. The likelihood for binary outcomes is curved on the log-odds scale, and its slope depends on the fitted probabilities. Consequently, the amount by which a group effect is pulled toward zero depends on the group’s observed outcomes and its covariates, not just on group size and the random-effect variance.

The central calculation is the score equation for the conditional mode. It balances the group’s evidence in the binary likelihood against the normal random-effects penalty. Deriving that equation also gives a direct way to check a proposed prediction and to understand why groups with identical sizes can have different amounts of shrinkage.

Model and Conditional-Mode Equation

For group \(j\), let \(y_{ij}\in\{0,1\}\) be the binary outcome for observation \(i\), where \(i=1,\ldots,n_j\). Suppose outcomes are conditionally independent given the group effect and covariates, with \(y_{ij}\mid b_j\sim\operatorname{Bernoulli}(p_{ij})\), and \(\operatorname{logit}(p_{ij})=x_{ij}^{T}\beta+b_j\). Here \(x_{ij}^{T}\beta\) is the fixed-effects linear predictor and \(b_j\) is the group-specific random intercept. The random-effects assumption is \(b_j\sim N(0,\tau^2)\), independently across groups. For the derivation below, treat \(\beta\) and \(\tau^2>0\) as fixed and known; in applications they are commonly replaced by fitted parameter estimates.

Write \(\eta_{ij}=x_{ij}^{T}\beta\) and \(p_{ij}(b)=\operatorname{logit}^{-1}(\eta_{ij}+b)\), where \(\operatorname{logit}^{-1}(z)=1/(1+\exp(-z))\). Omitting terms that do not depend on \(b\), the conditional log density for group \(j\) is

Differentiating each binary log-likelihood contribution gives \(y_{ij}-p_{ij}(b)\). Differentiating the normal penalty gives \(-b/\tau^2\). Thus the score is

The conditional mode solves \(U_j(\hat b_j)=0\). Equivalently, the sum of the binary residuals at the candidate group effect must equal the penalty contribution \(\hat b_j/\tau^2\). If all observations have the same fixed-effects predictor \(\eta\), and \(k_j=\sum_i y_{ij}\), this simplifies to

Why the Equation Has One Solution

The score changes monotonically with \(b\). Since the derivative of the logistic probability is \(p_{ij}(b)\{1-p_{ij}(b)\}\), differentiating the score gives

The inequality holds because every \(p_{ij}(b)\{1-p_{ij}(b)\}\) is nonnegative and \(1/\tau^2\) is strictly positive. Therefore the conditional log density is strictly concave. Also, as \(b\) tends to \(+\infty\), the term \(-b/\tau^2\) makes the score tend to \(-\infty\); as \(b\) tends to \(-\infty\), it makes the score tend to \(+\infty\). The score is continuous, so it crosses zero, and strict decrease means it crosses exactly once. This establishes a unique finite conditional mode under the stated model.

The equation is generally solved numerically. A Newton step from a current value \(b\) is \(b_{\mathrm{new}}=b-U_j(b)/U_j'(b)\). The derivative is always negative, so a numerical check should confirm that a proposed solution has a score near zero and negative curvature. Simply equating the observed group proportion to a fitted probability ignores the penalty term and is not the conditional-mode equation.

Outcomes and Fitted Probabilities Determine the Pull

At zero, the score is \(U_j(0)=\sum_i\{y_{ij}-p_{ij}(0)\}\), because the penalty contribution is zero there. This is the group’s total residual relative to the fixed-effects probabilities. A positive value pushes the mode above zero; a negative value pushes it below zero. But the score at zero alone does not determine the final estimate: as \(b\) changes, all fitted probabilities change, and the residuals are recalculated. Meanwhile, the penalty grows linearly in \(b\).

This dependence is more than a matter of the group’s number of successes. If observations have different fixed-effects predictors, two groups with the same outcomes and group size can have different residual sums and different conditional modes. Even with a common predictor, the fitted probability at the mode depends on the balance of outcomes and on the random-effect variance.

Local Curvature and the Meaning of Nonlinear Shrinkage

The score equation shows why there is no universal shrinkage fraction. In the Gaussian random-intercept case, the likelihood score is linear in the group effect, so combining it with a linear normal penalty leads to a linear weighted average. Here, changing \(b\) changes each \(p_{ij}(b)\), and therefore changes the likelihood slope. The relevant likelihood curvature is \(\sum_i p_{ij}(b)\{1-p_{ij}(b)\}\), which depends on the fitted probabilities and is largest near probability \(0.5\), decreasing near 0 and 1.

A local approximation makes the role of curvature explicit. Taylor-expand the score around zero: \(U_j(b)\approx U_j(0)+U_j'(0)b\). Solving this linear approximation for its zero gives

This is only a local approximation, not a general closed-form solution. Its numerator is the residual evidence at zero; its denominator combines the likelihood information at zero with the prior precision. Since the probabilities in the denominator change as the group effect moves away from zero, the exact balance is nonlinear.

Numerical Checks, Alternatives, and Limitations

When computing a conditional mode, evaluate the score at the reported value using unrounded probabilities where possible. A score close to zero checks the first-order condition; a negative second derivative checks that the point is a maximum. Also check that the values of \(\eta_{ij}\), \(\tau^2\), and the outcomes used in the score match the stated model. With extreme fitted probabilities or very sparse groups, rounding can make a small score appear larger or smaller than it is.

The conditional mode is not the unpenalized group-specific logistic estimate: it includes the normal penalty and is pulled toward zero. Nor is it generally a fixed fraction of that estimate. A Gaussian random-intercept predictor provides a useful comparison because its linear likelihood yields a fixed shrinkage weight under its assumptions; that weight should not be transferred to a binary GLMM. The conditional mode is also distinct from the conditional posterior mean. They are different summaries of the conditional distribution, and in a binary model there is no general reason they must coincide.

The derivation assumes the logistic link, conditionally independent Bernoulli observations, a normal random intercept with positive variance, and fixed covariates and parameters during the calculation. If the assumed mean structure or conditional independence is inappropriate, a numerically accurate root does not repair the model. In practice, the fitted fixed effects and variance are themselves estimated; treating them as fixed here isolates the nonlinear score mechanism rather than representing all sources of uncertainty in a prediction.