Asymptotic Normality of Generalized Linear Model Estimates

Why Asymptotic Normality Matters

A generalized linear model (GLM) coefficient estimate is usually obtained by maximizing a likelihood whose exact finite-sample distribution is difficult to describe. Asymptotic normality gives a tractable approximation when the sample grows: after centering at the true coefficient and scaling by the square root of sample size, the estimate approaches a multivariate normal distribution. The covariance in this limit is determined by the limiting information per observation.

The result is not automatic from taking a Newton or Fisher-scoring step. It depends on the model being identifiable, the design supplying information in every coefficient direction, a central limit theorem for the score, and a suitable local behavior of the likelihood curvature. The proof also needs care: a vector-valued score generally cannot be expanded using one common intermediate point on the line between two coefficient vectors.

Score and Information in a GLM

Consider independent responses \(Y_1,\ldots,Y_n\) with fixed covariate vectors \(x_i\in\mathbb{R}^p\), where \(p\) does not grow with \(n\). Let \(\eta_i=x_i^{\mathsf T}\beta\), and suppose the response distribution belongs to an exponential-dispersion family with known dispersion \(\phi\):

f_i(y_i;\beta) = \exp\left\{ \frac{y_i\theta_i-b(\theta_i)}{\phi} +c(y_i,\phi) \right\}, \qquad \mu_i=b'(\theta_i), \qquad \eta_i=g(\mu_i).

Here \(g\) is the link function and \(\beta\in\mathbb{R}^p\) is the coefficient vector. Assume the link and mean relationship are differentiable in the neighborhood under consideration, and write \(\dot{\mu}_i=d\mu_i/d\eta_i\). The variance function satisfies \(\operatorname{Var}(Y_i)=\phi b''(\theta_i)\).

Definition: The score is \(U_n(\beta)=\nabla\ell_n(\beta)\), where \(\ell_n\) is the log-likelihood. The observed information is \(-\nabla^2\ell_n(\beta)\), and the expected information is \(I_n(\beta)=E_\beta[-\nabla^2\ell_n(\beta)]\). Information per observation is the expected information divided by \(n\).

Differentiating the contribution of observation \(i\) gives a useful expression for the score. Since \(d\theta_i/d\mu_i=1/b''(\theta_i)\), the chain rule yields

\frac{\partial \ell_i}{\partial \beta} = \frac{Y_i-\mu_i}{\phi b''(\theta_i)} \dot{\mu}_i x_i = \frac{(Y_i-\mu_i)\dot{\mu}_i}{\operatorname{Var}(Y_i)}x_i.

The score has mean zero under the model. Independence then gives its covariance as the sum of the individual score covariances. Under the usual likelihood regularity conditions that permit the information identity, this covariance equals expected information:

I_n(\beta) = \sum_{i=1}^n \frac{\dot{\mu}_i^{\,2}}{\operatorname{Var}(Y_i)} x_i x_i^{\mathsf T}.

For the exponential-dispersion representation above, the weight is also \(\dot{\mu}_i^{\,2}/\{\phi b''(\theta_i)\}\). Thus design conditions matter through the weighted directions \(x_i\): having many observations does not help estimate a coefficient direction that the design fails to distinguish.

Conditions for the Limit

Write \(\beta_0\) for the true coefficient and \(\hat{\beta}\) for a maximum-likelihood estimator. The following formulation makes the main requirements explicit. It permits a fixed design that changes with \(n\), while keeping the coefficient dimension fixed.

Regularity conditions: Assume the responses are independent and the likelihood is correctly specified at an interior point \(\beta_0\) of an open parameter space. The model is identifiable in a neighborhood of \(\beta_0\), and a consistent MLE exists with probability tending to one; with probability tending to one it is interior and satisfies \(U_n(\hat{\beta})=0\). Assume the likelihood is twice continuously differentiable in that neighborhood and the score identities hold. The normalized expected information satisfies \(I_n(\beta_0)/n\to J\), where \(J\) is positive definite. Assume a multivariate Lindeberg condition for the independent score contributions \(s_i(\beta_0)=\nabla\ell_i(\beta_0)\): for every \(\epsilon>0\),

\(\displaystyle \frac1n\sum_{i=1}^n E_{\beta_0}\!\left[\|s_i(\beta_0)\|^2\,\mathbf{1}\{\|s_i(\beta_0)\|>\epsilon\sqrt n\}\right]\to0.\)

Finally, in a fixed neighborhood of \(\beta_0\), suppose the normalized negative Hessian converges uniformly in probability to a deterministic matrix function \(K(\beta)\) that is continuous at \(\beta_0\), with \(K(\beta_0)=J\). These conditions may be verified from specific likelihood and design assumptions; they are stated here as direct conditions on the score and curvature.

Main Result and Proof

Theorem: Under the stated conditions,

\(\displaystyle \sqrt n(\hat{\beta}-\beta_0)\xrightarrow{d}N(0,J^{-1}).\)

Equivalently, the large-sample covariance scale of \(\hat{\beta}\) is \(J^{-1}/n\). This is an asymptotic approximation, not a claim that the finite-sample distribution is exactly normal.

Proof. First, the score contributions are independent, have mean zero, and their covariance matrices sum to \(I_n(\beta_0)\). The information condition gives

\operatorname{Var}_{\beta_0}\!\left(\frac{U_n(\beta_0)}{\sqrt n}\right) = \frac{I_n(\beta_0)}{n} \longrightarrow J.

The multivariate Lindeberg–Feller central limit theorem therefore gives

\frac{U_n(\beta_0)}{\sqrt n}\xrightarrow{d}N(0,J).

To handle the likelihood curvature, let \(\Delta_n=\hat{\beta}-\beta_0\). Apply the fundamental theorem of calculus to the score along the line segment from \(\beta_0\) to \(\hat{\beta}\). On the event that the MLE is interior and its score is zero,

0 = U_n(\hat{\beta}) = U_n(\beta_0) + \left[ \int_0^1 \nabla^2\ell_n(\beta_0+t\Delta_n)\,dt \right]\Delta_n.

Define the integrated normalized negative Hessian

A_n = -\frac1n\int_0^1 \nabla^2\ell_n(\beta_0+t\Delta_n)\,dt.

Consistency gives \(\|\Delta_n\|\xrightarrow{p}0\), so every point \(\beta_0+t\Delta_n\) on the segment approaches \(\beta_0\) uniformly over \(0\leq t\leq1\), in probability. By the assumed uniform convergence of the normalized negative Hessian to \(K(\beta)\), and continuity of \(K\) at \(\beta_0\), it follows that \(A_n\xrightarrow{p}J\). Positive definiteness of \(J\) implies that \(A_n\) is invertible with probability tending to one and \(A_n^{-1}\xrightarrow{p}J^{-1}\).

Rearranging the integral expansion gives

\sqrt n(\hat{\beta}-\beta_0) = A_n^{-1}\frac{U_n(\beta_0)}{\sqrt n}.

Slutsky’s theorem now yields a limiting normal distribution with covariance \(J^{-1}J(J^{-1})^{\mathsf T}=J^{-1}\), since \(J\) is symmetric and positive definite. This proves the theorem. \(\square\)

The integrated Hessian is essential. For a vector-valued score, the scalar mean-value theorem does not generally provide one point on the segment at which the entire score difference equals the Hessian at that point times \(\Delta_n\). The integral identity is valid componentwise and supplies exactly the matrix needed in the argument.

Worked Examples

Worked Example: Balanced Logistic Regression

Let \(Y_i\) be independent Bernoulli responses with \(\operatorname{logit}(p_i)=\beta_0+\beta_1 z_i\). Take an even sample size with half the covariates \(z_i=1\) and half \(z_i=-1\), and let the true coefficient be \((\beta_0,\beta_1)=(0,\log 3)\). The two success probabilities are

p_{+}=\frac{3}{4}, \qquad p_{-}=\frac{1}{4}.

For the logistic link, \(\dot{\mu}_i=p_i(1-p_i)\), and the Bernoulli variance is also \(p_i(1-p_i)\). Thus each information weight is \(3/16\). With \(x_i=(1,z_i)^{\mathsf T}\), the average of \(x_ix_i^{\mathsf T}\) over the balanced design is the identity matrix: the off-diagonal entries average to zero and both diagonal entries average to one. Consequently,

J = \frac{3}{16} \begin{pmatrix}1&0\\0&1\end{pmatrix}, \qquad J^{-1} = \frac{16}{3} \begin{pmatrix}1&0\\0&1\end{pmatrix}.

The design has full rank, and the probabilities stay strictly between zero and one. The score contributions are bounded, so the Lindeberg condition holds. Under the stated consistency and local regularity conditions, \(\sqrt n(\hat{\beta}-\beta)\) converges to a bivariate normal distribution with covariance \((16/3)I_2\). The corresponding asymptotic covariance scale for \(\hat{\beta}\) is \((16/3n)I_2\).

Worked Example: Poisson Regression With Two Covariate Groups

Suppose \(Y_i\) are independent Poisson responses with log link, and the covariate vector is \(x_i=(1,z_i)^{\mathsf T}\), where half the observations have \(z_i=0\) and half have \(z_i=1\). Let the true coefficient be \((\log 2,\log(3/2))^{\mathsf T}\). The two means are \(\mu_0=2\) and \(\mu_1=3\). For the Poisson log-link model, \(\dot{\mu}_i=\mu_i\) and \(\operatorname{Var}(Y_i)=\mu_i\), so the information weight is \(\mu_i\). Averaging the two group contributions gives

J = \frac12\left[ 2\begin{pmatrix}1&0\\0&0\end{pmatrix} + 3\begin{pmatrix}1&1\\1&1\end{pmatrix} \right] = \begin{pmatrix}5/2&3/2\\3/2&3/2\end{pmatrix}.

Its determinant is \((5/2)(3/2)-(3/2)^2=15/4-9/4=3/2>0\), and direct inversion gives

J^{-1} = \begin{pmatrix}1&-1\\-1&5/3\end{pmatrix}.

The distinct covariate groups identify both coefficients, and the means remain finite and positive. Poisson score contributions have finite second moments; with these two fixed covariate values and bounded means, the Lindeberg condition holds. Therefore the limiting covariance of \(\sqrt n(\hat{\beta}-\beta)\) is the displayed matrix. Its negative off-diagonal entry reflects the dependence between the estimated intercept and slope, not a negative variance.

Worked Example: Gaussian Regression as an Exact Check

Let \(Y_i\sim N(\beta_0+\beta_1z_i,\sigma^2)\) independently, with known \(\sigma^2\), and use a balanced design with half \(z_i=1\) and half \(z_i=-1\). This is a GLM with identity link. The mean derivative is one and each response variance is \(\sigma^2\), so

J = \frac{1}{\sigma^2} \left(\frac1n\sum_{i=1}^n \begin{pmatrix}1&z_i\\z_i&z_i^2\end{pmatrix}\right) = \frac{1}{\sigma^2} \begin{pmatrix}1&0\\0&1\end{pmatrix}.

Hence the theorem gives limiting covariance \(J^{-1}=\sigma^2 I_2\) for \(\sqrt n(\hat{\beta}-\beta)\). In this particular model, the least-squares coefficient estimate is exactly normal for every \(n\) with this balanced design: its covariance is \(\sigma^2(X^{\mathsf T}X)^{-1}=(\sigma^2/n)I_2\). Thus the asymptotic result agrees with an exact finite-sample calculation here. That exact normality is special to the Gaussian model and should not be assumed for other GLMs.

Interpreting the Result and Avoiding Pitfalls

The limiting matrix \(J\) summarizes information per observation. If it is nearly singular, some linear combinations of coefficients are weakly informed, and the corresponding entries of \(J^{-1}\) can be large. If it is singular, the theorem as stated does not apply: at least one coefficient direction lacks positive information, so a full-dimensional normal limit with covariance \(J^{-1}\) is unavailable.

A common mistake is to use the observed information at the estimate as though its inverse were automatically the exact covariance of the estimator. Under the theorem’s conditions, the normalized observed information converges to \(J\), so its inverse provides an asymptotically justified covariance estimate on the \(1/n\) scale. It does not describe the exact finite-sample distribution, and it does not repair failure of identifiability, nonexistence of a finite MLE, or insufficient design information.

Key takeaway: A GLM MLE has the information-based normal limit when its score obeys a central limit theorem, the normalized local curvature converges to a positive-definite matrix, and a consistent interior MLE exists. The proof uses the integral of the Hessian along the coefficient segment, not a single intermediate point.

Check Your Understanding

Use the stated conditions and examples to check each conclusion.

Why does positive definiteness of \(J\) matter when solving the integrated score equation?
For a GLM, what role do the covariate vectors and the information weights play in determining \(J\)?
What condition on the score contributions supports the multivariate central limit theorem in the result?
Why is an integrated Hessian used instead of asserting one common mean-value point for the score vector?
In the balanced logistic example, explain why the limiting off-diagonal covariance is zero.