First-Order Small-Sample Bias in Logistic Regression

Why Logistic Maximum-Likelihood Estimates Can Be Biased

In logistic regression, the maximum-likelihood estimate is obtained by solving a nonlinear score equation. Even when the model is correctly specified and a finite estimate exists, the nonlinearity means that the estimate need not average to the true coefficient in repeated samples. With a fixed number of coefficients and increasing sample size, the leading bias is typically of order \(n^{-1}\).

This tutorial derives that leading term. It shows how the inverse information matrix transmits the bias across coefficients, while each observation’s contribution depends on its weighted leverage and on whether its fitted probability lies below or above one half. The formula is an approximation for the ordinary MLE, not a guarantee that a finite MLE exists for every possible sample.

Model, Information, and Leverage

Consider independent responses \(Y_i\sim\operatorname{Bernoulli}(p_i)\), for \(i=1,\ldots,n\), with fixed covariate vectors \(x_i\in\mathbb{R}^p\). The logistic model at the true coefficient vector \(\beta\) is

\operatorname{logit}(p_i)=x_i^\mathsf{T}\beta,\qquad p_i=\frac{\exp(x_i^\mathsf{T}\beta)}{1+\exp(x_i^\mathsf{T}\beta)},\qquad w_i=p_i(1-p_i).

Write \(X\) for the matrix with rows \(x_i^\mathsf{T}\), and \(W=\operatorname{diag}(w_1,\ldots,w_n)\). The score at the true coefficient and the Fisher information are

U(\beta)=\sum_{i=1}^n x_i(Y_i-p_i),\qquad I=X^\mathsf{T}WX=\sum_{i=1}^n w_i x_i x_i^\mathsf{T}.

Assume \(I\) is nonsingular. Define the weighted leverage of observation \(i\) by

h_i=w_i x_i^\mathsf{T}I^{-1}x_i.

These are the diagonal entries of the weighted hat matrix \(W^{1/2}X(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W^{1/2}\). They measure how strongly a row of the design contributes in the geometry determined by the model’s information.

Worked Example: Intercept-Only Logistic Regression

Suppose \(n=20\), the model has only an intercept, and the true success probability is \(p=0.2\). Then every \(x_i=1\), \(w_i=p(1-p)=0.16\), and the information is \(I=20(0.16)=3.2\).

For each observation, \(h_i=w_i/I=0.16/3.2=0.05=1/20\). The leading bias formula derived below gives

\operatorname{Bias}(\hat\beta) \approx I^{-1}\sum_{i=1}^{20}h_i(p-\tfrac12)x_i =\frac{1}{3.2}\,20\left(\frac{1}{20}\right)(-0.3) =-0.09375.

The true intercept is \(\operatorname{logit}(0.2)=\log(0.2/0.8)\), and the leading bias is negative. This is a first-order approximation, not an exact finite-sample bias.

Leading Bias: Statement and Derivation

An expectation statement for the MLE needs care: a finite MLE may fail to exist on some samples. The theorem below makes the expectation explicit by conditioning on existence and states the approximation and moment conditions needed for the expansion. These conditions are part of the result; design regularity alone should not be taken to imply a finite MLE for every sample or sufficient control of its conditional tails.

Theorem: Consider a sequence of fixed-design logistic models with fixed coefficient dimension \(p\). For each \(n\), let \(E_n\) be the event that a finite logistic MLE \(\hat\beta\) exists, and let \(I_n\), \(U_n\), \(h_{ni}\), and \(p_{ni}\) be defined as above at the true coefficient. Suppose \(I_n\) is nonsingular and the following conditional expansion and moment conditions hold:

\(P(E_n)>0\) and \(P(E_n^c)=o(n^{-1})\).
\(E[U_n\mid E_n]=o(1)\) and \(E[U_nU_n^\mathsf{T}\mid E_n]=I_n+o(n)\), entrywise.
On \(E_n\), writing \(\delta_n=\hat\beta-\beta\), the score expansion can be written as \(\delta_n=I_n^{-1}U_n+\frac12 I_n^{-1}T_n[I_n^{-1}U_n,I_n^{-1}U_n]+r_n\), where \(E[\|r_n\|\mid E_n]=o(n^{-1})\). Here \(T_n\) is the second derivative of the score at the true coefficient, and the bracket denotes contraction of its last two indices with the displayed vectors.
The design and probabilities are regular in the sense that \(I_n\) has eigenvalues of order \(n\), the covariate vectors are bounded, and the probabilities stay bounded away from zero and one. These conditions ensure the displayed leading term is of order \(n^{-1}\); the conditional expansion and moment conditions above still need to hold.

Then the conditional first-order bias is

E[\hat\beta-\beta\mid E_n] = I_n^{-1}\sum_{i=1}^n h_{ni}\left(p_{ni}-\frac12\right)x_i +o(n^{-1}).

To derive the expression, differentiate the score. Its first derivative at \(\beta\) is \(-I_n\), and its second derivative is the three-index array

T_n=-\sum_{i=1}^n w_{ni}(1-2p_{ni})\,x_i\otimes x_i\otimes x_i.

Here \(x_i\otimes x_i\otimes x_i\) has entries \(x_{ia}x_{ib}x_{ic}\). Apply the stated expansion and take conditional expectations. The linear term contributes \(I_n^{-1}E[U_n\mid E_n]=o(n^{-1})\), because \(I_n^{-1}\) is of order \(n^{-1}\). In the quadratic term, the conditional second-moment assumption allows \(U_nU_n^\mathsf{T}\) to be replaced by \(I_n\) at leading order. In coordinates, its \(j\)th component is therefore

\begin{aligned} \frac12 \sum_{a,b,c}(I_n^{-1})_{ja}(T_n)_{abc}(I_n^{-1})_{bc} &=-\frac12\sum_{a,b,c}(I_n^{-1})_{ja} \sum_i w_{ni}(1-2p_{ni})x_{ia}x_{ib}x_{ic}(I_n^{-1})_{bc}\\ &=\sum_i w_{ni}(p_{ni}-\tfrac12)(I_n^{-1}x_i)_j\, x_i^\mathsf{T}I_n^{-1}x_i\\ &=\sum_i h_{ni}(p_{ni}-\tfrac12)(I_n^{-1}x_i)_j. \end{aligned}

In vector form, this is \(I_n^{-1}\sum_i h_{ni}(p_{ni}-\tfrac12)x_i\), the theorem’s displayed formula. The remainder has conditional expectation \(o(n^{-1})\) by assumption, so it does not change the leading term. This proves the result. \(\square\)

The formula separates the roles of the three quantities. The vector \(I_n^{-1}\) sets the scale and mixes contributions among coefficients. The leverage \(h_{ni}\) weights observations that are influential in the information geometry. The factor \(p_{ni}-\frac12\) changes sign at one half: observations with fitted probabilities above and below one half contribute in opposite directions before the information matrix is applied.

Worked Example: A Two-Point Covariate Design

Take \(m\) observations at each of \(x=-1,0,1\), with an intercept and slope, so each row is \(x_i=(1,x)^\mathsf{T}\). Let the true coefficients be \(\beta_0=0\) and \(\beta_1=\log 2\). The probabilities at \(x=-1,0,1\) are \(1/3,1/2,2/3\), respectively; the corresponding weights are \(2/9,1/4,2/9\).

Symmetry makes the off-diagonal entry of the information matrix zero. Its diagonal entries are

I_{00}=m\left(\frac{2}{9}+\frac14+\frac{2}{9}\right)=\frac{25m}{36}, \qquad I_{11}=m\left(\frac{2}{9}+\frac{2}{9}\right)=\frac{4m}{9}.

For either endpoint observation, the leverage is

h_{\pm1} =\frac{2}{9}\left(\frac{36}{25m}+\frac{9}{4m}\right) =\frac{41}{50m}.

At \(x=0\), it is \(h_0=(1/4)(36/(25m))=9/(25m)\). In the bias numerator, the intercept contributions from the two endpoints cancel because \(p-\frac12\) has opposite signs there, and the middle contribution is zero. For the slope, each endpoint observation contributes \((41/(50m))(1/6)=41/(300m)\): at \(x=-1\), both \(p-\frac12\) and \(x\) are negative, so their product is positive. Counting \(m\) observations at each of the two endpoints gives a total slope numerator of \(41/150\). Applying \(I_{11}^{-1}=9/(4m)\) gives

\operatorname{Bias}(\hat\beta_0)\approx 0,\qquad \operatorname{Bias}(\hat\beta_1)\approx \frac{9}{4m}\left(\frac{41}{150}\right) =\frac{123}{200m}.

Thus the leading slope bias is positive in this design, even though the intercept bias cancels by symmetry.

A Check on the Leverage Weights

Result: If \(I=X^\mathsf{T}WX\) is nonsingular, the weighted leverages satisfy \(\sum_{i=1}^n h_i=p\), where \(p\) is the number of coefficients.

To prove this, let \(H=W^{1/2}X(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W^{1/2}\). Its \(i\)th diagonal entry is \(h_i\), so the sum of the leverages is its trace. By cyclic invariance of the trace,

\sum_i h_i=\operatorname{tr}(H) =\operatorname{tr}\left((X^\mathsf{T}WX)^{-1}X^\mathsf{T}WX\right) =\operatorname{tr}(I_p)=p.

This identity is a useful arithmetic check, but it does not mean the bias formula is determined by the number of parameters alone: the leverage must still be paired with the probabilities and covariate vectors.

Worked Example: Two Binomial Groups

Consider a saturated logistic model for two groups, with \(m=10\) independent observations in each group. The control probability is \(p_C=0.2\), and the treatment probability is \(p_T=0.6\). Use an intercept and a treatment indicator. Each group’s fitted logit has the same leading bias as an intercept-only model with \(m\) observations, namely \((p-\frac12)/(mp(1-p))\).

For control, this is \((-0.3)/(10(0.2)(0.8))=-0.1875\). For treatment, it is \(0.1/(10(0.6)(0.4))=0.0416667\). The intercept is the control logit, so its leading bias is \(-0.1875\). The treatment coefficient is the treatment logit minus the control logit, giving leading bias

$$ 0.0416667-(-0.1875)=0.2291667. $$

The two coefficient biases differ because they represent different contrasts of the group logits. This example also shows why reporting a single direction of bias for all logistic coefficients would be misleading.

Interpretation and a Common Pitfall

The expansion explains why small-sample bias is not simply “toward zero” for every coefficient. Each observation first contributes in the direction of \((p_i-\frac12)x_i\), weighted by its leverage; the inverse information then maps that aggregate contribution into coefficient space. Changes to the covariate coding can change the coefficient components while leaving the underlying fitted probabilities unchanged.

A common mistake is to state the displayed expression as an unconditional expectation without addressing samples on which no finite MLE exists. Even a simple intercept-only model has positive probability of observing all zeros or all ones, and on either event the ordinary logistic likelihood has no finite maximizer. The correction is to specify the expectation being approximated and the conditions that justify the expansion. Here the result concerns \(E[\hat\beta-\beta\mid E_n]\), with explicit conditional moment and remainder assumptions; it does not silently assign a finite value to the MLE on nonexistence samples.

The leading term is useful for understanding the source and scale of MLE bias, and for checking calculations in small or imbalanced designs. It is an asymptotic approximation, not an exact correction or a substitute for checking whether the MLE exists. In particular, the leverage identity checks the weights’ total, but cannot by itself establish that the conditional expansion is valid.

Check Your Understanding

Use the definitions and derivation above to answer the following questions.

In the bias expression, what role does the inverse information matrix play, and what role does the factor \(p_i-\frac12\) play?
Why is it not generally valid to write an unconditional expectation of the ordinary logistic MLE without discussing existence?
For an intercept-only model, show that \(h_i=1/n\), and simplify the leading bias expression.
Why do the two endpoint contributions to the intercept bias cancel in the symmetric worked example?
Prove that the weighted leverages sum to the number of coefficients when the information matrix is nonsingular.