Modified IRLS for Firth Logistic Regression

Why the Iteration Must Be Modified

Ordinary iteratively reweighted least squares (IRLS) for logistic regression uses a working response and weights built from the usual residuals $y_i-\pi_i$. Firth logistic regression changes the score by a leverage-dependent adjustment. Therefore, reusing the ordinary working response without modification does not solve the Firth score equations. The central task here is to put that adjusted score into a weighted least-squares form.

Assume $Y_1,\ldots,Y_n$ are independent Bernoulli responses with $P(Y_i=1)=\pi_i$, and that the logistic model has linear predictor $\eta_i=x_i^\top\beta$ and probability $\pi_i=1/(1+\exp(-\eta_i))$. The $n\times p$ design matrix $X$ has rows $x_i^\top$ and full column rank. At finite $\beta$, define $w_i=\pi_i(1-\pi_i)$, $W=\operatorname{diag}(w_1,\ldots,w_n)$, and the Fisher information $I(\beta)=X^\top W X$. Full column rank and positive $w_i$ make $I(\beta)$ positive definite.

Definition: The logistic leverage at observation $i$ is $h_i=w_i x_i^\top I(\beta)^{-1}x_i$, the $i$th diagonal entry of $H=W^{1/2}X(X^\top W X)^{-1}X^\top W^{1/2}$. The Firth-adjusted score is $$ U^*(\beta)=X^\top\{y-\pi+h\odot(1/2-\pi)\}, $$ where $\pi=(\pi_1,\ldots,\pi_n)^\top$, $h=(h_1,\ldots,h_n)^\top$, and $\odot$ denotes componentwise multiplication.

This score is the gradient of the Jeffreys-penalized log-likelihood $\ell_F(\beta)=\ell(\beta)+\tfrac12\log|I(\beta)|$. In the logistic model, differentiating this penalized likelihood gives the leverage adjustment shown above. The adjusted score, rather than the ordinary score $X^\top(y-\pi)$, is the equation the iteration must target.

From the Adjusted Score to Modified IRLS

At a current finite value $\beta$, let $\eta=X\beta$ and define the adjusted residual $r_i^*=y_i-\pi_i+h_i(1/2-\pi_i)$. Since $w_i>0$, we can divide this residual by its working weight and define an adjusted working response.

Definition: The modified IRLS working response and weights at $\beta$ are $$ z_i^*=\eta_i+\frac{r_i^*}{w_i} =\eta_i+\frac{y_i-\pi_i+h_i(1/2-\pi_i)}{\pi_i(1-\pi_i)}, \qquad w_i=\pi_i(1-\pi_i). $$ The modified Fisher-scoring update is the weighted least-squares fit of $z^*$ on $X$ with weights $W$.

“Fit” here means solve the weighted least-squares normal equations using quantities evaluated at the current $\beta$. The weights and response are then recomputed at the new iterate. They are not held fixed across all iterations.

Theorem: If $X$ has full column rank and $\beta$ is finite, the modified IRLS update is well-defined and satisfies $$ \beta_{\mathrm{new}}=(X^\top W X)^{-1}X^\top W z^* =\beta+I(\beta)^{-1}U^*(\beta). $$

Proof. By the definition of $z_i^*$, its vector form is $z^*=\eta+W^{-1}r^*$. Since $\eta=X\beta$, $$ X^\top W z^*=X^\top W X\beta+X^\top r^* =I(\beta)\beta+U^*(\beta). $$ The full-rank assumption and $w_i>0$ make $I(\beta)$ invertible. Multiplying by its inverse gives the stated update. The same invertibility ensures the weighted least-squares normal equations have a unique solution. $\square$

Worked Example: Computing Leverages and an Adjusted Score

Take three independent Bernoulli observations, with covariates $x_i=(1,t_i)^\top$, $t=(0,1,2)$, and outcomes $y=(0,0,1)^\top$. Evaluate one modified IRLS step at $\beta=(0,0)^\top$.

At this iterate every $\eta_i=0$, so $\pi_i=1/2$ and $w_i=1/4$. The information matrix and its inverse are $$ I=\tfrac14 \begin{pmatrix}3&3\\3&5\end{pmatrix} = \begin{pmatrix}0.75&0.75\\0.75&1.25\end{pmatrix}, \qquad I^{-1}= \begin{pmatrix}10/3&-2\\-2&2\end{pmatrix}. $$

For $x_1=(1,0)^\top$, $h_1=\tfrac14 x_1^\top I^{-1}x_1=5/6$. For $x_2=(1,1)^\top$, $h_2=\tfrac14(10/3-4+2)=1/3$. For $x_3=(1,2)^\top$, $h_3=\tfrac14(10/3-8+8)=5/6$. The leverages sum to $2$, the number of fitted coefficients.

Because $\pi_i=1/2$, the leverage adjustments $h_i(1/2-\pi_i)$ are zero at this particular iterate. Thus $r^*=y-\pi=(-1/2,-1/2,1/2)^\top$, and $$ U^*=X^\top r^*= \begin{pmatrix}-1/2\\1/2\end{pmatrix}, \qquad \delta=I^{-1}U^*= \begin{pmatrix}-8/3\\2\end{pmatrix}. $$ The adjusted working responses are $z^*=r^*/(1/4)=(-2,-2,2)^\top$. Their weighted least-squares normal-equation right-hand side is $X^\top Wz^*=(-1/2,1/2)^\top$, so the fitted coefficient vector is $(-8/3,2)^\top$, agreeing with $\beta+\delta$. This is a large full step; the calculation establishes the update, not that taking it undamped will always improve the penalized likelihood.

What the Leverage Adjustment Does

The term $h_i(1/2-\pi_i)$ depends both on the fitted probability and on how strongly observation $i$ influences the fitted linear predictors. It is not a constant “add one half” to each observed outcome. In particular, at $\pi_i=1/2$ the adjustment is zero, as in the example above, even when the leverages are nonzero.

Worked Example: An Intercept-Only Firth Estimate

Suppose there are $n=4$ observations in an intercept-only model, with one success. At any finite intercept $\beta$, every fitted probability is $\pi$, every weight is $w=\pi(1-\pi)$, and every leverage is $1/4$. The adjusted score is $$ U^*(\beta)=1-4\pi+4\left\{\tfrac14(\tfrac12-\pi)\right\} =\tfrac32-5\pi. $$

Setting the score to zero gives $\hat\pi_F=3/10=0.3$, and therefore $$ \hat\beta_F=\log\left(\frac{0.3}{0.7}\right)=\log(3/7)\approx-0.8473. $$ For comparison, ordinary logistic maximum likelihood would set $\pi=1/4$, corresponding to $\log(1/3)$. The different answer follows from the leverage adjustment in the score, not from changing the response data or the logistic link.

Convergence Checks and Step Damping

A modified IRLS step is a Fisher-scoring step for the penalized likelihood, but a full step need not increase that likelihood. For a trial step length $a$, calculate $\beta(a)=\beta+a\delta$, with $\delta=I^{-1}U^*$, and recompute $\ell_F$ at that trial value. If the penalized log-likelihood decreases, reduce $a$, for example by repeatedly halving it, and try again. A practical implementation should also reject trial points that produce non-finite calculations.

Theorem: Suppose $I(\beta)$ is positive definite and the penalized log-likelihood $\ell_F$ is differentiable at $\beta$, with gradient $U^*(\beta)\ne0$. Then the Fisher-scoring direction $\delta=I(\beta)^{-1}U^*(\beta)$ is a strict local ascent direction: for all sufficiently small positive $a$, $\ell_F(\beta+a\delta)>\ell_F(\beta)$.

Proof. The directional derivative at $a=0$ is $$ \left.\frac{d}{da}\ell_F(\beta+a\delta)\right|_{a=0} =U^*(\beta)^\top\delta =U^*(\beta)^\top I(\beta)^{-1}U^*(\beta)>0. $$ The last inequality holds because $I(\beta)^{-1}$ is positive definite and the score is nonzero. By the definition of a positive derivative, $\{\ell_F(\beta+a\delta)-\ell_F(\beta)\}/a$ remains positive for all sufficiently small positive $a$. Hence the penalized log-likelihood increases for those step lengths. This local statement does not guarantee that a full step increases the objective or that an iteration converges from every starting point. $\square$

Worked Example: Checking a Damped Intercept Step

Return to four observations with one success, now starting at $\beta=0$. Here $\pi=1/2$, $I=4(1/4)=1$, and the adjusted score is $U^*=1-4(1/2)=-1$. Thus $\delta=-1$. Consider the half-step $\beta(1/2)=-0.5$.

At $\beta=-0.5$, $\pi=1/(1+\exp(0.5))\approx0.3775407$, and $4\pi(1-\pi)\approx0.9400148$. The ordinary log-likelihood for one success and three failures is $$ \ell(-0.5)=\log(\pi)+3\log(1-\pi)\approx-2.3963. $$ The Jeffreys penalty, relative to the information determinant at the starting value, is $\tfrac12\log\{4\pi(1-\pi)\}\approx-0.0309$. Thus $\ell_F(-0.5)\approx-2.3963-0.0309=-2.4272$. At $\beta=0$, the ordinary log-likelihood is $-4\log 2\approx-2.7726$, and the penalty is $\tfrac12\log(1)=0$. The half-step therefore increases the penalized log-likelihood, from about $-2.7726$ to $-2.4272$. This check uses the penalized objective, not the ordinary likelihood alone.

In practice, convergence should be assessed using more than one diagnostic. Useful checks include a small change in the coefficients, a small adjusted score (or score scaled relative to the information), and a small change in $\ell_F$. Recompute $h_i$, $W$, $z^*$, and the adjusted score at the accepted iterate. A small coefficient change by itself is not enough: damping can make a step small even when the score is not yet close to zero.

Common Mistakes

Using ordinary residuals in the working response. The numerator must be $y_i-\pi_i+h_i(1/2-\pi_i)$. Omitting its leverage term generally solves a different set of equations.
Using the wrong weights. The modified response changes, but the logistic working weight remains $w_i=\pi_i(1-\pi_i)$, evaluated at the current iterate.
Treating the update as an exact maximizer. It solves the local weighted least-squares equations. It does not maximize $\ell_F$ over all coefficient values in one step.
Accepting a step because ordinary likelihood rises. Firth fitting targets $\ell_F$, so a step should be checked against the penalized log-likelihood.
Declaring convergence from a small step alone. Also check the adjusted score and objective change; a tiny damped step may simply reflect repeated step reduction.

Key takeaway: Modified IRLS replaces the ordinary logistic residual by the leverage-adjusted residual in the working response, while retaining the logistic weights. Its weighted least-squares update equals $\beta+I^{-1}U^*$; damping and convergence diagnostics help control the iteration but do not turn each full step into a guaranteed global improvement.

Check Your Understanding

Try these questions before moving on.

Write the adjusted working response $z_i^*$ and identify which quantities are recomputed at each iterate.
Why does $X^\top Wz^*=I\beta+U^*$ imply that the weighted least-squares update is $\beta+I^{-1}U^*$?
In an intercept-only model with $n$ observations, what is each leverage, and why do the leverage adjustments sum to a simple term in the adjusted score?
Explain why a positive directional derivative at step length zero guarantees some sufficiently small step increases $\ell_F$, but does not guarantee a full step will do so.
Name three diagnostics that can be monitored before declaring a modified IRLS fit converged.