Bias in Finite-Sample Generalized Linear Model Estimates

Why Finite-Sample Bias Matters

A maximum-likelihood estimate can be consistent and asymptotically normal yet still be systematically displaced from the true parameter in a finite sample. In generalized linear models (GLMs), likelihood curvature is one source of this displacement: even when the score has mean zero at the true coefficient, solving the nonlinear score equation does not generally produce an estimator with mean equal to that coefficient.

The leading bias is typically of order \(n^{-1}\), whereas the standard error is of order \(n^{-1/2}\). Thus the bias vanishes faster than the standard error in large samples, but it can still shift Wald inference noticeably in a finite sample. The derivation below treats independent observations in a canonical-link GLM, where the score has a particularly simple form. For a noncanonical link, derivatives of the score weights also contribute, so the canonical-link formula should not be applied unchanged.

Setup and Leading-Bias Formula

Let \(Y_1,\ldots,Y_n\) be independent observations from a canonical exponential-family model with fixed dispersion \(\phi>0\). Write \(x_i\) for the fixed \(p\)-vector of covariates, \(\eta_i=x_i^\mathsf{T}\beta\), and \(\mu_i=E[Y_i]\). In the canonical parametrization, the variance is \(\operatorname{Var}(Y_i)=\phi V_i\), where \(V_i=V(\mu_i)\), and \(d\mu_i/d\eta_i=V_i\). Assume \(p\) is fixed as \(n\) grows. Define \(V_i'=dV(\mu_i(\eta_i))/d\eta_i\), with the derivative taken with respect to \(\eta_i\).

Definition: The finite-sample bias of an estimator \(\hat\beta\) for \(\beta\) is \(E_\beta[\hat\beta]-\beta\), when this expectation exists and is finite. A leading bias expansion writes this expectation as a specified order-\(n^{-1}\) term plus a smaller remainder. It is not a claim that the exact finite-sample bias equals that leading term.

For the canonical model, the score and expected information for \(\beta\) are

U(\beta)=\frac{1}{\phi}\sum_{i=1}^n x_i(Y_i-\mu_i), \qquad I(\beta)=\frac{1}{\phi}\sum_{i=1}^n V_i x_i x_i^\mathsf{T}.

The information is assumed nonsingular. Put \(q_i=x_i^\mathsf{T}I(\beta)^{-1}x_i\). Under regularity conditions, the leading bias is

E_\beta[\hat\beta-\beta] = -\frac{1}{2\phi}I(\beta)^{-1} \sum_{i=1}^n x_i V_i' q_i +o(n^{-1}).

Here the formula is an asymptotic expansion for an estimator whose expectation is well-defined. Sufficient conditions include bounded covariates and the needed derivatives of \(V\), information eigenvalues of order \(n\), enough moment and tail control for the Taylor remainder to have expectation \(o(n^{-1})\), and a consistent score root with probability tending to one. If the ordinary MLE fails to be finite on some samples, the formula does not by itself define its unconditional bias. One must instead analyze a finite estimator on all samples, or explicitly restrict attention to an interior event and justify the resulting conditional expansion.

Deriving the Expansion

The key point is that the expected score is zero at the true parameter, but the score equation is nonlinear in \(\beta\). Let \(\delta=\hat\beta-\beta\). Taylor-expand the score around the true value. In a canonical model, its first derivative is \(-I(\beta)\), and its second derivative, applied twice to a direction \(\delta\), is \(-\phi^{-1}\sum_i x_iV_i'(x_i^\mathsf{T}\delta)^2\). Thus the score equation gives

0 = U(\beta)-I(\beta)\delta -\frac{1}{2\phi}\sum_{i=1}^n x_iV_i'(x_i^\mathsf{T}\delta)^2 +R_n.

The stated regularity conditions ensure that the expectation of the remainder, after solving for \(\delta\), is \(o(n^{-1})\). To obtain the leading term, first use \(\delta_1=I^{-1}U(\beta)\) in the quadratic term. Replacing \(\delta\) by \(\delta_1\) there changes the expectation only by the higher-order remainder. Since \(E[U(\beta)]=0\) and \(\operatorname{Cov}(U(\beta))=I(\beta)\),

E[(x_i^\mathsf{T}\delta_1)^2] = x_i^\mathsf{T}I^{-1}\operatorname{Cov}(U)I^{-1}x_i = x_i^\mathsf{T}I^{-1}x_i = q_i.

Solving the Taylor equation for \(\delta\), and then taking expectations, gives the formula.

Derivation: Ignoring only the remainder whose expected contribution is \(o(n^{-1})\), rearrange the score equation as

\delta = I^{-1}U(\beta) -\frac{1}{2\phi}I^{-1}\sum_{i=1}^n x_iV_i'(x_i^\mathsf{T}\delta)^2.

The first term has expectation zero. Substituting \(\delta_1\) into the quadratic term and using the covariance calculation above yields

E[\delta] = -\frac{1}{2\phi}I^{-1}\sum_{i=1}^n x_iV_i'q_i+o(n^{-1}).

The factor \(1/\phi\) comes from the second derivative of the score and must remain in the bias formula. Since \(I^{-1}\) and each \(q_i\) both scale with \(\phi\), this explicit factor gives the leading term its stated dispersion scaling:

\boxed{ E[\hat\beta-\beta] = -\frac{1}{2\phi}I^{-1}\sum_{i=1}^n x_iV_i'q_i +o(n^{-1}). }

The boxed expression is the dispersion-explicit version used below. In particular, since \(I^{-1}\) and each \(q_i\) both scale with \(\phi\), the resulting bias scales with \(\phi\). The expected-information convention is essential when tracking this factor.

Worked Examples

Worked Example: Poisson Intercept and a Conditional Interior MLE

Let \(Y_1,\ldots,Y_n\) be independent Poisson observations with mean \(\lambda=e^\beta\), and consider the intercept-only model. The MLE is \(\hat\beta=\log(\bar Y)\) when \(T=\sum_iY_i>0\); it is not finite when \(T=0\). Therefore an unconditional finite-sample bias formula for this MLE is not valid. Instead, consider its expectation conditional on \(T>0\), where the estimate is finite.

Here \(x_i=1\), \(V_i=\lambda\), and \(V_i'=d\lambda/d\beta=\lambda\). The information is \(I=n\lambda\), so \(q_i=1/(n\lambda)\). Substitution into the leading term gives

-\frac{1}{2}I^{-1}\sum_{i=1}^n V_i'q_i = -\frac{1}{2n\lambda} \sum_{i=1}^n\frac{\lambda}{n\lambda} = -\frac{1}{2n\lambda}.

This is indeed the conditional leading bias. To check the boundary issue, \(T\sim\operatorname{Poisson}(n\lambda)\), and \(P(T=0)=e^{-n\lambda}\). Conditional on \(T>0\), \(E[T\mid T>0]=n\lambda/(1-e^{-n\lambda})\); hence the conditional mean of \((T-n\lambda)/(n\lambda)\) is exponentially small in \(n\). Its conditional second moment is \(1/(n\lambda)+o(n^{-1})\). Expanding the logarithm around \(T/(n\lambda)=1\), with the tail away from that neighborhood contributing an exponentially small amount, gives

E\!\left[\log\!\left(\frac{T}{n\lambda}\right)\middle|T>0\right] = -\frac{1}{2n\lambda}+o(n^{-1}).

This calculation is explicitly conditional. It does not repair or imply an unconditional expectation for the ordinary MLE, which equals \(-\infty\) on \(T=0\) if extended in that way.

Worked Example: Binomial Logit Intercept on Interior Samples

Suppose \(Y\sim\operatorname{Binomial}(n,p)\), with \(0<p<1\), and use a logit intercept \(\beta=\log\{p/(1-p)\}\). The estimate \(\hat\beta=\log\{Y/(n-Y)\}\) is finite only when \(0<Y<n\). We therefore discuss its conditional expectation given \(0<Y<n\), not an unconditional bias for this infinite-at-the-boundary MLE.

For one Bernoulli observation, \(V=p(1-p)\) and \(V'=p(1-p)(1-2p)\). The information for the \(n\) observations is \(I=np(1-p)\); each observation has \(q_i=1/I\). The formula gives

-\frac{1}{2}I^{-1}\sum_{i=1}^n V'q_i = -\frac{1}{2np(1-p)} \frac{n\,p(1-p)(1-2p)}{np(1-p)} = \frac{2p-1}{2np(1-p)}.

To verify the conditional interpretation, write \(\hat p=Y/n\) and expand \(\operatorname{logit}(\hat p)\) around \(p\). The first derivative is \(1/[p(1-p)]\), and the second derivative is \((2p-1)/[p^2(1-p)^2]\). Since \(E[\hat p-p]=0\) and \(E[(\hat p-p)^2]=p(1-p)/n\), the quadratic contribution is \((2p-1)/[2np(1-p)]\). Conditioning on \(0<Y<n\) changes these moments only by exponentially small terms for fixed \(p\in(0,1)\), because \(P(Y=0)=(1-p)^n\) and \(P(Y=n)=p^n\). The Taylor expansion is valid on a neighborhood of \(p\), while the probability of leaving that neighborhood is exponentially small. Thus the same leading term holds conditionally. Neither boundary event may be ignored when claiming an unconditional expectation for the unmodified MLE.

Worked Example: Gaussian Identity-Link Regression

Let \(Y_i\) be independent \(N(x_i^\mathsf{T}\beta,\sigma^2)\) observations with fixed, full-rank design matrix \(X\), and known \(\sigma^2\). This is a canonical GLM with identity link and \(V_i=1\), so \(V_i'=0\). The leading bias expression is therefore zero.

In this model the conclusion is exact, not merely asymptotic. The likelihood score equation is \(X^\mathsf{T}(y-X\beta)/\sigma^2=0\), giving \(\hat\beta=(X^\mathsf{T}X)^{-1}X^\mathsf{T}Y\). Since \(E[Y]=X\beta\),

E[\hat\beta] = (X^\mathsf{T}X)^{-1}X^\mathsf{T}X\beta = \beta.

The contrast with Poisson and binomial models is informative: nonlinear mean-variance structure can create the leading curvature term, whereas the Gaussian identity-link score equation is linear in the coefficients.

How Bias Affects Wald Inference

For a scalar contrast \(c^\mathsf{T}\beta\), Wald inference uses the estimate \(c^\mathsf{T}\hat\beta\) and its estimated standard error. If the coefficient bias is \(b_n=E[\hat\beta-\beta]\), then the contrast is shifted by \(c^\mathsf{T}b_n\). Since the standard error is usually of order \(n^{-1/2}\) and the leading bias is of order \(n^{-1}\), the bias measured in standard-error units is of order \(n^{-1/2}\). It vanishes asymptotically, but it can shift finite-sample Wald statistics and the coverage of intervals centered at the estimate.

For the Poisson intercept example with \(n=10\) and \(\lambda=1\), the conditional leading bias is \(-1/(2n\lambda)=-0.05\). The asymptotic standard error is \(1/\sqrt{n\lambda}\approx0.316\), so the bias is about \(-0.05/0.316\approx-0.158\) standard errors. This is a modest shift in this example, not a universal correction to apply to every Wald calculation. The standard error itself is also estimated, and finite-sample performance depends on the model and design.

Common pitfall: A zero-mean score does not imply an unbiased MLE. The incorrect argument is “\(E[U(\beta)]=0\), so the solution \(\hat\beta\) must have mean \(\beta\).” That would follow for a linear score equation, but generally not for a nonlinear one. The score expansion shows the missing quadratic term explicitly. A separate pitfall is to report an unconditional bias for an MLE that is infinite on samples of positive probability; the correction is to specify a finite estimator or state clearly that the calculation is conditional or formal.

Bias expansions help explain why Wald procedures can behave unevenly in small samples, especially when the likelihood is curved or the estimate lies near a boundary. They do not by themselves prescribe a bias correction: changing an estimator can also change its variability, and inference must account for both. For noncanonical links, the full score expansion must include derivatives of the observation-specific score weights.

Check Your Understanding

Use the definitions and calculations above to check your understanding.

Why does the canonical-link score have expected information \(I=\phi^{-1}\sum_iV_ix_ix_i^\mathsf{T}\)?
In the leading bias derivation, where does the factor \(q_i=x_i^\mathsf{T}I^{-1}x_i\) come from?
Why is the Poisson intercept calculation not an unconditional bias expansion for the ordinary MLE?
For the binomial logit intercept, what happens to the leading bias term when \(p=1/2\), and what is the interpretation?
Why does a bias of order \(n^{-1}\) still have the potential to affect finite-sample Wald inference?