Asymptotic Normality of GLM Maximum-Likelihood Estimators

From a Score to a Distribution

Consistency says that a maximum-likelihood estimator gets close to the true parameter as the sample grows. It does not, by itself, describe the estimator’s fluctuations around that parameter. Asymptotic normality supplies that description: after scaling the estimation error by the square root of the sample size, its distribution approaches a multivariate normal distribution.

The argument has two ingredients. The score at the true parameter is a sum of contributions from observations, so a central limit theorem can describe its fluctuations. Near the true parameter, the change in the score is governed by the likelihood’s curvature. Expanding the score and combining these facts turns the score’s limit into a limit for the estimator.

Setup: For each sample size $n$, let $Y_1,\ldots,Y_n$ be independent responses from a correctly specified GLM, with true parameter $\theta_0\in\mathbb{R}^p$, where $p$ is fixed and $\theta_0$ is an interior point of the parameter space. Write $\ell_n(\theta)$ for the log-likelihood, $U_n(\theta)=\partial\ell_n(\theta)/\partial\theta$ for its score, and $H_n(\theta)=\partial^2\ell_n(\theta)/\partial\theta\partial\theta^{\mathsf T}$ for its Hessian.

The design may be fixed or random; the conditions below are stated in terms of probability under the true model, so they can cover either case. The dimension $p$ stays fixed as $n$ grows. Assume the log-likelihood is twice continuously differentiable in a neighborhood of $\theta_0$, and that the estimator $\hat{\theta}_n$ is consistent: $\hat{\theta}_n\xrightarrow{p}\theta_0$. Assume it is an interior likelihood maximizer with probability tending to one, so $U_n(\hat{\theta}_n)=0$ with probability tending to one.

Let $J$ be a fixed, positive-definite $p\times p$ matrix. The two central limit and curvature conditions are

\frac{1}{\sqrt{n}}U_n(\theta_0)\xrightarrow{d}N(0,J), \qquad \sup_{\theta\in N} \left\| \frac{1}{n}H_n(\theta)-A(\theta) \right\|\xrightarrow{p}0, \qquad A(\theta_0)=-J,

for some open neighborhood $N$ of $\theta_0$ and a deterministic matrix-valued function $A$ continuous at $\theta_0$. Here $\|\cdot\|$ is any matrix norm. The first condition says that the score fluctuates on the order of $\sqrt n$; the second says that the normalized curvature is stable near the true parameter. These are sufficient conditions, not a claim that every GLM or every design automatically satisfies them.

The Asymptotic Normality Result

Theorem: Under the setup and conditions above, $$ \sqrt{n}(\hat{\theta}_n-\theta_0)\xrightarrow{d}N(0,J^{-1}). $$

To prove the result, use the integral form of the score expansion along the line joining $\theta_0$ to $\hat{\theta}_n$. Define

\bar{H}_n = \int_0^1 \frac{1}{n}H_n\bigl(\theta_0+t(\hat{\theta}_n-\theta_0)\bigr)\,dt.

On the event that the estimator is an interior maximizer and the segment lies in $N$, the fundamental theorem of calculus gives

0 = \frac{1}{\sqrt{n}}U_n(\hat{\theta}_n) = \frac{1}{\sqrt{n}}U_n(\theta_0) + \bar{H}_n\sqrt{n}(\hat{\theta}_n-\theta_0).

Show the average curvature converges.
Consistency puts the entire segment inside $N$ with probability tending to one. The uniform Hessian condition makes its normalized Hessian uniformly close to $A(\theta)$; continuity of $A$ at $\theta_0$ then gives $\bar{H}_n\xrightarrow{p}-J$.

Invert the curvature.
Since $J$ is positive definite, $-J$ is nonsingular. Matrix inversion is continuous at nonsingular matrices, so $\bar{H}_n^{-1}\xrightarrow{p}-J^{-1}$, on events whose probability tends to one.

Combine the limits.
The score equation gives $\sqrt{n}(\hat{\theta}_n-\theta_0)=-\bar{H}_n^{-1}U_n(\theta_0)/\sqrt{n}$. The score central limit theorem and Slutsky’s theorem yield a mean-zero normal limit with covariance $(-J^{-1})J(-J^{-1})^{\mathsf T}=J^{-1}$.

The integral expansion also explains why the sign of the Hessian matters: at a local likelihood maximum its limiting value is negative, whereas the covariance of the limiting estimator is positive definite. The proof uses consistency to place the estimator in a region where the curvature is controlled; it does not require assuming that the estimator already has a square-root-$n$ rate.

When Does the Score Have a Central Limit Theorem?

For independent observations, write the score as $U_n(\theta_0)=\sum_{i=1}^n s_{ni}$, where $s_{ni}$ is the score contribution of observation $i$ at the true parameter. A useful sufficient condition for the score limit is a multivariate Lyapunov condition. Suppose $E[s_{ni}]=0$, $n^{-1}\sum_i\operatorname{Cov}(s_{ni})\to J$, and, for some $\delta>0$ and every fixed vector $a\in\mathbb{R}^p$,

\frac{1}{n^{1+\delta/2}} \sum_{i=1}^n E\left[|a^{\mathsf T}s_{ni}|^{2+\delta}\right] \longrightarrow 0.

Result: If $J$ is positive definite and these conditions hold, then $n^{-1/2}U_n(\theta_0)\xrightarrow{d}N(0,J)$.

Proof. Fix any nonzero vector $a$. The variance of $a^{\mathsf T}U_n(\theta_0)/\sqrt{n}$ is

\frac{1}{n}\sum_{i=1}^n \operatorname{Var}(a^{\mathsf T}s_{ni}) \longrightarrow a^{\mathsf T}Ja>0.

The Lyapunov condition for the independent scalar summands $a^{\mathsf T}s_{ni}/\sqrt{n}$ holds: its numerator, before division by the total variance to the power $1+\delta/2$, is $n^{-1-\delta/2}\sum_i E|a^{\mathsf T}s_{ni}|^{2+\delta}$, which tends to zero. The scalar Lyapunov central limit theorem therefore gives $a^{\mathsf T}U_n(\theta_0)/\sqrt{n}\xrightarrow{d}N(0,a^{\mathsf T}Ja)$. Since this holds for every $a$, the Cramér–Wold theorem gives the stated multivariate limit. $\square$

In a GLM, each score contribution is built from the response and its design vector. The moment condition must therefore be checked for the actual response distribution and design sequence; independence alone does not establish it.

Worked GLM Examples

Worked Example: Bernoulli Logistic Intercept

Suppose $Y_1,\ldots,Y_n$ are independent Bernoulli responses with common success probability $p_0=0.4$, and use the logistic intercept parameter $\theta$, so $p(\theta)=e^\theta/(1+e^\theta)$. The true parameter is $\theta_0=\log(0.4/0.6)=\log(2/3)$.

The log-likelihood, score, and Hessian are

\ell_n(\theta)=\sum_{i=1}^n\{Y_i\theta-\log(1+e^\theta)\}, \qquad U_n(\theta)=\sum_{i=1}^n\{Y_i-p(\theta)\}, \qquad H_n(\theta)=-np(\theta)\{1-p(\theta)\}.

At $\theta_0$, $E[U_n(\theta_0)]=0$ and $\operatorname{Var}\{U_n(\theta_0)\}=n(0.4)(0.6)=0.24n$. Thus $J=0.24$, and the bounded Bernoulli score contributions satisfy the central limit condition. Also $H_n(\theta)/n=-p(\theta)\{1-p(\theta)\}$ is continuous and equals $-0.24$ at $\theta_0$. The theorem gives

\sqrt{n}(\hat{\theta}_n-\theta_0) \xrightarrow{d}N\left(0,\frac{1}{0.24}\right) = N\left(0,\frac{25}{6}\right).

The variance is for the log-odds parameter $\theta$, not for the success probability $p$. The likelihood has an interior finite maximizer with probability tending to one because the sample proportion converges in probability to $0.4$, which is strictly between zero and one.

Worked Example: Poisson Log-Mean Intercept

Let $Y_1,\ldots,Y_n$ be independent Poisson responses with mean $\lambda(\theta)=e^\theta$, and take the true mean to be $\lambda_0=3$, so $\theta_0=\log 3$. Ignoring terms that do not depend on $\theta$, the log-likelihood and its derivatives are

\ell_n(\theta)=\sum_{i=1}^n\{Y_i\theta-e^\theta\}, \qquad U_n(\theta)=\sum_{i=1}^n(Y_i-e^\theta), \qquad H_n(\theta)=-ne^\theta.

At the true value, each score contribution has mean zero and variance $\operatorname{Var}(Y_i)=3$. Hence $J=3$. The Poisson distribution has finite moments of every positive integer order, so, for example, the Lyapunov condition holds with $\delta=1$ for this fixed-mean model. The normalized Hessian is $-e^\theta$, continuous and equal to $-3$ at $\theta_0$. Therefore

\sqrt{n}(\hat{\theta}_n-\log 3) \xrightarrow{d}N\left(0,\frac{1}{3}\right).

The score equation gives $e^{\hat{\theta}_n}=\bar{Y}$ when $\bar{Y}>0$. Since $P(\bar{Y}=0)=e^{-3n}\to0$, the finite interior maximizer exists with probability tending to one, as required.

Worked Example: Gaussian Regression With Alternating Design

Consider the Gaussian GLM $Y_i=\theta_1+\theta_2 x_i+\varepsilon_i$, with known variance $\sigma^2=4$, independent $\varepsilon_i\sim N(0,4)$, and fixed $x_i=1$ for even $i$, $x_i=-1$ for odd $i$. Let $n$ be even and take $\theta_0=(1,0.5)^{\mathsf T}$. The design vector is $z_i=(1,x_i)^{\mathsf T}$; the score contribution is $s_{ni}=z_i\varepsilon_i/4$.

There are $n/2$ design vectors of each type. Thus $\sum_i z_i z_i^{\mathsf T}=nI_2$, and

\operatorname{Var}\{U_n(\theta_0)\} = \sum_{i=1}^n \frac{z_i z_i^{\mathsf T}}{4} = \frac{n}{4}I_2.

The normalized Hessian is exactly $-I_2/4$, so $J=I_2/4$. The score is exactly multivariate normal for every $n$; the design has full rank, and the maximum-likelihood estimator is interior. Consequently,

\sqrt{n}(\hat{\theta}_n-\theta_0) \xrightarrow{d}N(0,4I_2).

In this example the estimator is already exactly normal, with covariance $4I_2/n$ at each even sample size. The asymptotic theorem recovers that distributional scaling.

Interpretation and a Common Pitfall

The matrix $J$ describes the limiting information per observation in the chosen parameterization. More curvature means that a given random score fluctuation moves the maximizer less, which is reflected by the inverse matrix $J^{-1}$ in the limit. This is a statement about the estimator’s sampling distribution after rescaling; it is not an exact finite-sample normality claim.

Pitfall: treating consistency as a central limit theorem. The incorrect reasoning is “$\hat{\theta}_n\xrightarrow{p}\theta_0$, therefore $\sqrt{n}(\hat{\theta}_n-\theta_0)$ is asymptotically normal.” Consistency gives no rate or limiting distribution. The correction is to establish a score central limit theorem and stable, nonsingular local curvature, then use the score expansion.
Pitfall: overlooking a singular information limit. If $J$ is singular, the inverse in the theorem does not exist. The stated normal limit has not been established; the parameterization or design may fail to identify all directions at the square-root-$n$ scale.
Pitfall: applying the score equation at a boundary maximum. The equation $U_n(\hat{\theta}_n)=0$ is justified here by an interior maximizer. If a maximum occurs on the boundary, that equation need not hold, and this proof does not apply as written.

The central argument is therefore local but not automatic: consistency brings the estimator into the region where the likelihood is smooth, the score’s central limit theorem supplies its random displacement, and the nonsingular Hessian translates that displacement into estimation error.

Check Your Understanding

Use the score expansion to explain each answer, rather than relying only on the final normal-limit formula.

Which two limiting conditions combine to yield asymptotic normality of the estimator?
Why does positive definiteness of $J$ matter when solving the expanded score equation?
For the Poisson intercept example with true mean $\lambda_0$, what is the limiting variance for $\sqrt{n}(\hat{\theta}_n-\theta_0)$, where $\theta_0=\log\lambda_0$?
Why does consistency alone not establish a square-root-$n$ normal limit?
What feature of the Gaussian design example ensures that both regression parameters have nonsingular limiting information?