Firth Correction in Mixed-Effects Models

Why a Firth Adjustment Needs Care in a Mixed Model

In ordinary logistic regression, separation can send maximum-likelihood coefficient estimates toward infinity even though fitted probabilities remain bounded. Firth’s correction modifies the likelihood to reduce first-order estimation bias and, in common full-rank logistic models, produces finite estimates when separation prevents an ordinary finite maximum. A mixed-effects model presents an additional challenge: its likelihood for the fixed effects and variance parameters is obtained by integrating over the random effects. A correction defined for an ordinary logistic likelihood cannot simply be attached to the conditional likelihood and assumed to correct the mixed model.

The distinction is consequential. A conditional model describes outcomes given particular random-effect values; a GLMM’s marginal likelihood accounts for the distribution of those random effects. The information matrix, and therefore the Jeffreys-information adjustment, depends on which likelihood is being corrected. The derivation below first establishes the ordinary logistic adjustment, then shows what must change in a mixed-effects setting.

Definition: For a regular likelihood $L(\psi)$, let $\ell(\psi)=\log L(\psi)$ and let $I(\psi)$ be its expected Fisher information matrix. Firth’s penalized log likelihood is $$ \ell_F(\psi)=\ell(\psi)+\frac{1}{2}\log|I(\psi)|. $$ The added term is the Jeffreys-information penalty. Estimates maximize $\ell_F$, rather than the unadjusted log likelihood. The adjustment is defined only where the information determinant is positive; rank deficiency or other nonregularity requires separate treatment.

Firth’s Adjustment for Logistic Regression

For independent Bernoulli observations, suppose $Y_i$ has success probability $\mu_i$, with $\operatorname{logit}(\mu_i)=x_i^\mathsf{T}\beta$. Here $x_i$ is the predictor vector for observation $i$, and the rows of the design matrix form $X$. The ordinary score is $U(\beta)=X^\mathsf{T}(y-\mu)$, and the expected information is $I(\beta)=X^\mathsf{T}WX$, where $W$ is diagonal with entries $\mu_i(1-\mu_i)$. The Firth score is the derivative of $\ell_F$:

U_F(\beta)=U(\beta)+A(\beta), \qquad A_j(\beta)=\frac{1}{2}\operatorname{tr}\left\{I(\beta)^{-1}\frac{\partial I(\beta)}{\partial\beta_j}\right\}.

In logistic regression this can also be written as

U_F(\beta)=X^\mathsf{T}\left[y-\mu+h\mathbin{\odot}\left(\frac{1}{2}-\mu\right)\right],

where $h_i$ is the $i$th diagonal entry of $W^{1/2}X(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W^{1/2}$, and $\odot$ denotes elementwise multiplication. Thus the adjustment depends on the current fitted probabilities and the design’s leverage. It is not generally equivalent to adding a fixed half-success to every observation.

Worked Example: The Intercept-Only Firth Estimate

Suppose a hypothetical independent sample contains $n=9$ binary outcomes, of which $y=2$ are successes. Assume all observations share the same probability $p$, and the model contains only an intercept. The ordinary estimate is $2/9$, but the Firth estimate follows directly from the penalized likelihood.

The likelihood is proportional to $p^y(1-p)^{n-y}$, and the expected information for the intercept is $np(1-p)$. Omitting constants that do not depend on $p$, the penalized log likelihood is

\ell_F(p)=y\log(p)+(n-y)\log(1-p)+\frac{1}{2}\log\{np(1-p)\}.

Differentiating with respect to the intercept, whose derivative with respect to $p$ is $1/\{p(1-p)\}$, gives the score equation

y-np+\frac{1}{2}(1-2p)=0 \quad\Longrightarrow\quad p_F=\frac{y+0.5}{n+1}.

Substituting the counts gives $p_F=(2+0.5)/(9+1)=2.5/10=0.25$. Its fitted log odds are $\log\{0.25/(1-0.25)\}=\log(1/3)=-1.0986$. The ordinary estimate is $2/9=0.2222$, with log odds $\log\{(2/9)/(7/9)\}=\log(2/7)=-1.2528$. Both estimates are finite; this example illustrates the adjustment, not separation. The adjusted estimate is a model-based estimate of the common success probability, not a statement that the observed proportion must be changed.

A Finite Estimate Under Complete Separation

When the predictors perfectly separate successes from failures, the ordinary logistic likelihood can keep increasing as a coefficient grows without bound. The Firth penalty changes the optimization problem. The following grouped-data example shows the adjustment explicitly in a full-rank logistic model with a binary predictor.

Worked Example: Two Separated Predictor Groups

Consider hypothetical independent observations in two groups. At $x=0$, there are four failures and no successes; at $x=1$, there are four successes and no failures. Fit $\operatorname{logit}(p)=\beta_0+\beta_1x$. The observed outcomes are completely separated by $x$.

For this saturated two-group model, the Firth estimate of each group probability is its group success count plus $0.5$, divided by its group size plus $1$. This follows from the intercept-only score calculation applied to each group; the two fitted group logits are equivalent to the intercept and the sum of the intercept and slope.

For $x=0$, $p_{0,F}=(0+0.5)/(4+1)=0.1$. For $x=1$, $p_{1,F}=(4+0.5)/(4+1)=0.9$. Therefore,

\widehat{\beta}_{0,F}=\log\left(\frac{0.1}{0.9}\right)=\log(1/9)=-2.1972, \qquad \widehat{\beta}_{1,F}=\log\left(\frac{0.9}{0.1}\right)-\widehat{\beta}_{0,F} =2.1972-(-2.1972)=4.3944.

The adjusted odds ratio comparing $x=1$ with $x=0$ is $\exp(4.3944)=81$, also obtained from $(0.9/0.1)/(0.1/0.9)=9/0.1111\approx81$. The finite values result from maximizing the penalized likelihood, not from observing any failures in the $x=1$ group or successes in the $x=0$ group. The large odds ratio still reflects strong separation and should not be mistaken for high precision.

What Changes for a GLMM?

For a logistic random-intercept model, a conditional specification is $\operatorname{logit}\{P(Y_{ij}=1\mid b_j)\}=x_{ij}^\mathsf{T}\beta+b_j$, with random intercepts $b_j$ drawn from a distribution such as $N(0,\tau^2)$. Given $b_j$, outcomes in cluster $j$ are modeled by their conditional likelihood. The likelihood used for marginal maximum likelihood, however, integrates that conditional likelihood against the random-effect distribution. Writing $\theta$ for variance parameters, its form is

L_m(\beta,\theta)=\prod_j\int \left\{\prod_i P(Y_{ij}=y_{ij}\mid b_j;\beta)\right\} f(b_j;\theta)\,db_j, \qquad \ell_m(\beta,\theta)=\log L_m(\beta,\theta).

A Firth-style correction for this model must be derived from the relevant integrated likelihood and its information, or from a specifically justified bias-reducing procedure for the mixed-model parameters of interest. In a regular parameterization, the direct Jeffreys-information form would involve the information for the integrated likelihood, for example $\ell_m(\psi)+\tfrac12\log|I_m(\psi)|$ for a parameter vector $\psi$ being corrected. In practice, nuisance variance parameters, integration, approximation, and possible boundary estimates make the choice and computation of such a correction method-specific. There is no general license to substitute the ordinary fixed-effects logistic formula.

Worked Example: Why Conditioning on a Cluster Is Not Integration

Consider a hypothetical random-intercept model with one cluster containing two successes, fixed intercept $\beta_0=0$, and $b\sim N(0,1)$. Conditional on $b=0$, each observation has success probability $\operatorname{logit}^{-1}(0)=0.5$, so the conditional probability of two successes is $0.5^2=0.25$.

The integrated probability is instead

P(Y_{1}=1,Y_{2}=1) =\int \{\operatorname{logit}^{-1}(b)\}^2\phi(b;0,1)\,db =E\{p(b)^2\}.

By symmetry of the normal distribution and the identity $\operatorname{logit}^{-1}(-b)=1-\operatorname{logit}^{-1}(b)$, $E\{p(b)\}=0.5$. Hence $E\{p(b)^2\}=[E\{p(b)\}]^2+\operatorname{Var}\{p(b)\}=0.25+\operatorname{Var}\{p(b)\}>0.25$, because $p(b)$ varies when $b$ varies. The conditional contribution at $b=0$ is therefore not the integrated contribution. A penalty based on conditional information at a selected or estimated $b$ is not, merely by being a Firth penalty, the Jeffreys-information correction for the integrated likelihood.

Derivation, Assumptions, and Practical Checks

The intercept-only result also shows why the penalty changes the score. For $n$ independent Bernoulli observations, the log likelihood is $y\log p+(n-y)\log(1-p)$. The expected information for the intercept is $np(1-p)$, so the penalty contributes $\tfrac12\log p+\tfrac12\log(1-p)$, apart from a constant. Differentiating with respect to the intercept adds $\tfrac12(1-2p)$ to the ordinary score $y-np$. Solving the resulting equation gives $p_F=(y+0.5)/(n+1)$. Thus the finite estimate in the one-group case follows from the penalty’s score contribution, rather than from an arbitrary replacement of the observed data.

The ordinary logistic derivation assumes the stated independent Bernoulli likelihood, a correctly specified design, and an information matrix that is nonsingular in the parameter region used. The familiar finiteness result applies to standard full-rank logistic designs under separation; it does not repair duplicated columns, unidentified parameters, or a misspecified sampling model. In a GLMM, the random-effects distribution and conditional model must also be specified, and the marginal likelihood or its approximation must be handled by the selected method. A result for an independent logistic model does not establish that a mixed-model fit is bias-reduced.

Useful checks include whether the unadjusted fit shows separation indicators such as diverging coefficients or near-zero information in a direction, whether the corrected fit is stable to numerical tolerances and integration accuracy, and whether variance estimates lie at boundaries. Compare fitted values and likelihood calculations under the intended integrated model, not only conditional predictions at estimated random effects. A correction cannot create information absent from the design, and convergence of an optimizer alone does not verify that the intended likelihood was optimized.

Choosing and Interpreting an Approach

Firth logistic regression is the direct option for a genuinely independent binary-outcome model: it uses the ordinary logistic likelihood and its information. For clustered binary data, the target is typically a GLMM whose likelihood integrates over random effects. A mixed-model Firth-style method must address that integrated likelihood and state which parameters receive the adjustment. A proper-prior analysis is a different approach: it uses a specified prior distribution and yields inference conditional on that prior and likelihood, rather than the Jeffreys-information penalized estimate. Ridge penalization is another distinct method, not a synonym for Firth correction.

Interpret adjusted coefficients on the scale of the model. In the separated two-group example, the slope corresponds to a fitted conditional odds ratio of 81 for the specified logistic model. In a random-effects model, a conditional odds ratio and a marginal population-level odds ratio are generally different. Neither the existence of finite adjusted estimates nor an optimizer’s convergence implies that uncertainty is small. Use uncertainty procedures justified for the selected method and model.

Common pitfall: “Firth logistic regression returns finite coefficients, so I can apply the same correction to the conditional GLMM fit.” This reasoning uses the wrong likelihood for the marginal mixed model. The correction’s information term must correspond to the likelihood being corrected. Use a method explicitly developed for the integrated mixed-model likelihood, and verify its assumptions, parameterization, and numerical integration.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

In an intercept-only sample with $n=7$ and $y=0$, what is the Firth estimate of $p$, and what is its logit?
For $n=10$ and $y=6$, calculate the intercept-only Firth estimate and compare it with the ordinary sample proportion.
Why does the probability of two successes in the random-intercept example involve $E\{p(b)^2\}$, not just the squared probability at $b=0$?
State one reason a finite Firth-adjusted estimate does not demonstrate that a separated GLMM has been correctly fit.

Solutions. (1) $p_F=(0+0.5)/(7+1)=0.0625$; its logit is $\log(0.0625/0.9375)=\log(1/15)=-2.7081$. (2) $p_F=(6+0.5)/(10+1)=6.5/11=0.5909$, compared with the ordinary proportion $6/10=0.6$; the adjusted estimate is slightly closer to $0.5$. (3) The same random intercept affects both outcomes in the cluster, so their conditional joint probability is $p(b)^2$, which must be averaged over the random-effect distribution. (4) The correction may have been applied to a conditional rather than integrated likelihood, or the design may remain unidentified; finiteness alone checks neither the target likelihood nor model adequacy.