Quasi-Likelihood and Sandwich Covariance Under Misspecification

Why Misspecification Changes Likelihood Inference

A likelihood model need not describe the true data-generating distribution exactly. A fitted model can still provide a useful summary, but its target is then not necessarily a true parameter value: it is the parameter that makes the working model best approximate the data according to expected log likelihood. Moreover, the usual inverse-information covariance can be wrong. The reason is that two different features of the likelihood matter: how variable the score is under the actual distribution, and how sharply the expected log likelihood bends near its maximizer.

This tutorial derives that target and the resulting covariance. Throughout, $X_1,\ldots,X_n$ are independent and identically distributed from a true distribution $P_0$. The working model has density or mass function $q_\theta(x)$, indexed by $\theta\in\Theta\subseteq\mathbb{R}^p$, and $\ell(\theta;x)=\log q_\theta(x)$ denotes its log likelihood contribution for one observation. The observations follow $P_0$, not necessarily any member of the working model.

Definition: The pseudo-true parameter $\theta^\ast$ is a maximizer of the expected log likelihood $$ Q(\theta)=E_0[\ell(\theta;X_1)], $$ where $E_0$ denotes expectation under $P_0$. When $P_0$ has a density and the relevant expectations are finite, maximizing $Q(\theta)$ is equivalent to minimizing the Kullback–Leibler divergence from $P_0$ to the working model. This parameter is a best-fitting working-model value, not a claim that the model is true.

The sample maximum-likelihood estimator maximizes $\sum_{i=1}^n\ell(\theta;X_i)$, or equivalently its average. Under consistency and suitable differentiability, it converges to $\theta^\ast$. We focus on the regular case in which $\theta^\ast$ is a unique interior maximizer and the expected log likelihood is differentiable there.

Two Matrices with Different Roles

Write the individual score and Hessian as

s(\theta;x)=\nabla_\theta\ell(\theta;x), \qquad H(\theta;x)=\nabla_\theta^2\ell(\theta;x).

Assume differentiation can be interchanged with expectation near $\theta^\ast$; for example, an integrable envelope for the score together with the usual local differentiability and domination conditions can justify this interchange. Since $\theta^\ast$ is an interior maximizer, its expected score vanishes. Define the expected negative curvature $A$ and score variability $B$ by

E_0[s(\theta^\ast;X_1)]=0,\qquad A=-E_0[H(\theta^\ast;X_1)],\qquad B=E_0[s(\theta^\ast;X_1)s(\theta^\ast;X_1)^{\mathsf T}].

The score has mean zero, so $B$ is also its covariance matrix. The matrix $A$ measures local curvature of the expected log likelihood; $B$ measures random score fluctuation under the true distribution. Both are $p\times p$ matrices. Under a correctly specified regular model, the information identity gives $A=B$. Under misspecification, that equality generally does not hold.

Proposition: If $Q$ is differentiable at an interior maximizer $\theta^\ast$, and differentiation under expectation is valid there, then $E_0[s(\theta^\ast;X_1)]=0$.

Proof. Differentiation under expectation gives $\nabla Q(\theta^\ast)=E_0[s(\theta^\ast;X_1)]$. For any direction $v\in\mathbb{R}^p$, the scalar function $t\mapsto Q(\theta^\ast+tv)$ has a local maximum at $t=0$, because $\theta^\ast$ is interior. Its derivative at zero is therefore zero, and that derivative equals $v^{\mathsf T}\nabla Q(\theta^\ast)$. Since this holds for every $v$, $\nabla Q(\theta^\ast)=0$, proving the claim. $\square$

Worked Example: Fitting a Normal Distribution to Uniform Data

Suppose the true observations are uniform on $(-\sqrt{3},\sqrt{3})$, but the working model is $N(\mu,v)$, with $v>0$ the variance parameter. The true distribution has $E_0[X]=0$, $E_0[X^2]=1$, and $E_0[X^4]=9/5$. Determine the pseudo-true value and the sandwich covariance for the maximum-likelihood estimator of $(\mu,v)$.

Ignoring no parameter-dependent terms, the working log density is $\ell(\mu,v;x)=-\frac12\log(2\pi v)-\frac{(x-\mu)^2}{2v}$. The expected criterion is maximized at $\mu^\ast=E_0[X]=0$, and its variance component is maximized at $v^\ast=E_0[(X-\mu^\ast)^2]=1$. Thus $\theta^\ast=(0,1)$, although the true distribution is not normal.

At this value, the score is $s(\theta^\ast;X)=(X,(X^2-1)/2)^{\mathsf T}$. The Hessian has expected negative value $A=\operatorname{diag}(1,1/2)$: the mean-mean second derivative is $-1/v$, the cross derivative has expectation zero, and the variance-variance second derivative has expectation $-1/2$ at $(0,1)$. By symmetry, $E_0[X(X^2-1)]=0$, so $B=\operatorname{diag}(1,(9/5-1)/4)=\operatorname{diag}(1,1/5)$. Since $A^{-1}=\operatorname{diag}(1,2)$, the sandwich covariance is $A^{-1}BA^{-1}=\operatorname{diag}(1,4/5)$.

Consequently, the asymptotic covariance of $(\hat\mu,\hat v)$ is $\operatorname{diag}(1,4/5)/n$. The model-based inverse-curvature matrix is instead $\operatorname{diag}(1,2)$, which would incorrectly report variance $2/n$ for $\hat v$. Here the sample variance’s fluctuation depends on the true fourth moment, not the fourth moment implied by the fitted normal model.

The Sandwich Limit

The following result gives the central asymptotic calculation. It requires more than a formal Taylor expansion: the score must satisfy a central limit theorem, the estimator must be consistent and satisfy the likelihood score equation, and the empirical Hessian must converge uniformly near the target. For example, a uniform law of large numbers on a fixed neighborhood, together with continuity of the expected Hessian there, is a sufficient route to the local Hessian convergence used below.

Theorem: Suppose $X_1,\ldots,X_n$ are iid from $P_0$, $\theta^\ast$ is an interior pseudo-true parameter, and $\hat\theta$ is a consistent interior likelihood maximizer satisfying $\sum_{i=1}^n s(\hat\theta;X_i)=0$ with probability tending to one. Assume $E_0[s(\theta^\ast;X_1)]=0$, $E_0[\|s(\theta^\ast;X_1)\|^2]<\infty$, and $A=-E_0[H(\theta^\ast;X_1)]$ is positive definite. Assume also that the empirical Hessian converges uniformly in probability to its expectation on a neighborhood of $\theta^\ast$, and that the expected Hessian is continuous at $\theta^\ast$. Then $$ \sqrt n(\hat\theta-\theta^\ast) \xrightarrow{d} N\!\left(0,A^{-1}BA^{-1}\right). $$

Proof. Define the empirical Hessian $H_n(\theta)=n^{-1}\sum_{i=1}^n H(\theta;X_i)$. The vector-valued score expansion uses an integral remainder, not a scalar mean-value theorem. Applying the fundamental theorem of calculus along the line segment from $\theta^\ast$ to $\hat\theta$ gives

0= \frac{1}{\sqrt n}\sum_{i=1}^n s(\theta^\ast;X_i) + J_n\sqrt n(\hat\theta-\theta^\ast), \qquad J_n=\int_0^1 H_n\!\left(\theta^\ast+t(\hat\theta-\theta^\ast)\right)\,dt.

By consistency, the entire segment in this integral lies in any fixed neighborhood of $\theta^\ast$ with probability tending to one. Uniform convergence of $H_n$ to $E_0[H(\theta;X_1)]$ on that neighborhood, and continuity of this expectation at $\theta^\ast$, imply $J_n\xrightarrow{p}E_0[H(\theta^\ast;X_1)]=-A$. Because $A$ is positive definite, it is nonsingular; hence $J_n^{-1}\xrightarrow{p}-A^{-1}$.

The multivariate central limit theorem applies to the iid score vectors, which have mean zero and finite covariance $B$:

\frac{1}{\sqrt n}\sum_{i=1}^n s(\theta^\ast;X_i) \xrightarrow{d}N(0,B).

Rearranging the expansion yields $\sqrt n(\hat\theta-\theta^\ast)=-J_n^{-1}n^{-1/2}\sum_i s(\theta^\ast;X_i)$. Slutsky’s theorem now gives a centered normal limit whose covariance is $(-A^{-1})B(-A^{-1})^{\mathsf T}=A^{-1}BA^{-1}$, since $A$ is symmetric. This proves the result. $\square$

Worked Calculations in Scalar Models

Worked Example: A Poisson Working Model for Overdispersed Counts

Let the true count $Y$ equal zero with probability $0.6$ and five with probability $0.4$. Fit a Poisson working model with mean $\lambda>0$. The true mean is $E_0[Y]=0.4(5)=2$, and $E_0[Y^2]=0.4(25)=10$, so $\operatorname{Var}_0(Y)=10-2^2=6$.

The working log likelihood contribution is $Y\log\lambda-\lambda-\log(Y!)$, with score $s(\lambda;Y)=Y/\lambda-1$ and Hessian $-Y/\lambda^2$. Maximizing the expected log likelihood gives $\lambda^\ast=E_0[Y]=2$. Thus $A=E_0[Y]/(\lambda^\ast)^2=2/4=1/2$. The score has mean zero, and $B=\operatorname{Var}_0(Y)/(\lambda^\ast)^2=6/4=3/2$.

The sandwich variance is $A^{-1}BA^{-1}=2(3/2)2=6$. Indeed, the fitted Poisson MLE is the sample mean, whose asymptotic variance is the true variance $6$. In contrast, the inverse-curvature variance is $A^{-1}=2$, the value the Poisson model would imply if it were correct. The two answers differ because the true counts are more variable than the working Poisson model allows.

Worked Example: An Exponential Working Model for Gamma Data

Suppose $Y$ has a Gamma distribution with shape $2$ and rate $1$, so $E_0[Y]=2$ and $\operatorname{Var}_0(Y)=2$. Fit an exponential working model with rate $\lambda>0$, whose log likelihood contribution is $\ell(\lambda;y)=\log\lambda-\lambda y$.

The expected log likelihood is $\log\lambda-2\lambda$, with derivative $1/\lambda-2$. Its unique maximizer is $\lambda^\ast=1/2$. The score is $1/\lambda-Y$; the negative expected Hessian at the target is $A=1/(\lambda^\ast)^2=4$. The score variance is $B=\operatorname{Var}_0(Y)=2$. Therefore the sandwich variance is $A^{-1}BA^{-1}=(1/4)(2)(1/4)=1/8$.

The MLE is $\hat\lambda=1/\bar Y$. Applying the scalar delta method to $\bar Y$, whose asymptotic variance is $2/n$, gives asymptotic variance $(1/2^2)^2(2/n)=1/(8n)$, agreeing with the sandwich result. The working exponential model would instead give inverse-information variance $A^{-1}/n=1/(4n)$. Its implied variance for $Y$ is $4$, whereas the true variance is $2$.

Using and Interpreting the Sandwich

The covariance formula separates two sources of behavior. The middle matrix $B$ records how much the observed score varies under the actual sampling distribution. The outer factors $A^{-1}$ translate a score perturbation into a parameter perturbation according to the local curvature of the expected objective. This explains the “sandwich” terminology. In applied work, the empirical counterparts are the negative average Hessian and the average outer product of individual scores, evaluated at $\hat\theta$; their combination estimates the asymptotic covariance of $\hat\theta$, which is then divided by $n$ for its approximate finite-sample covariance.

The formula does not make a misspecified model correct, nor does it remove bias from a poor target choice. It quantifies uncertainty around the pseudo-true value under the stated iid and regularity conditions. If observations are dependent, the score covariance generally needs to account for dependence; if the estimator is on a boundary or the curvature is singular, this regular sandwich theorem may not apply.

Common pitfall: Using $A^{-1}$ as the covariance merely because it is the inverse expected negative Hessian silently assumes the information identity $A=B$. Under misspecification, compute both matrices under the true distribution and use $A^{-1}BA^{-1}$. Also, for a vector score expansion, do not assert that a single intermediate parameter gives the matrix mean-value formula; use the integral Hessian along the segment, as in the proof.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

What distinguishes the pseudo-true parameter from a parameter value in a correctly specified model?
In the sandwich covariance, which matrix describes score variability, and which describes expected curvature?
For the Poisson example, verify directly that the pseudo-true score has mean zero.
Why does the normal working model in the uniform-data example give a different variance estimate for $\hat v$ when one uses inverse curvature alone?
Why is an integral Hessian expansion appropriate for a vector score equation?