Worked Asymptotics for Gaussian Linear Regression

Gaussian Regression and Its Large-Sample Questions

In Gaussian linear regression, the mean is described by a linear predictor and the errors have a common variance. This setting makes the connection between design and precision especially transparent: the coefficient covariance depends directly on the matrix of observed predictors. It also allows us to see, in one model, how coefficient and dispersion estimates behave jointly.

We use a fixed-design formulation. For each sample size $n$, let $x_i$ be a nonrandom $p$-vector of predictors, and let $Y_i$ be the response. The number $p$ of regression coefficients is fixed as $n$ grows. Write $X_n$ for the $n\times p$ matrix whose $i$th row is $x_i^\mathsf{T}$, and $\beta$ for the coefficient vector.

Definition: The Gaussian linear regression model is $$ Y_i=x_i^\mathsf{T}\beta+\varepsilon_i,\qquad \varepsilon_i\ \text{independent and}\ \varepsilon_i\sim N(0,\phi),\quad \phi>0. $$ Equivalently, conditional on the fixed design, $Y_n\sim N(X_n\beta,\phi I_n)$. The parameter $\phi$ is the error variance, or dispersion. The coefficient maximum-likelihood estimator is $\hat\beta=(X_n^\mathsf{T}X_n)^{-1}X_n^\mathsf{T}Y_n$ when $X_n$ has full column rank, and the dispersion maximum-likelihood estimator is $\hat\phi=\mathrm{RSS}/n$, where $\mathrm{RSS}=(Y_n-X_n\hat\beta)^\mathsf{T}(Y_n-X_n\hat\beta)$.

For a regular large-sample limit, assume $X_n$ has full column rank for all sufficiently large $n$, and that $$ Q_n=\frac{1}{n}X_n^\mathsf{T}X_n\xrightarrow{}Q, $$ where $Q$ is a finite, positive-definite $p\times p$ matrix. This condition says that the per-observation design information stabilizes and does not collapse in any coefficient direction. Since $p$ is fixed, it also implies that $n>p$ eventually.

Exact Distributions and the Information Matrix

The Gaussian likelihood gives particularly simple estimators. Up to a constant not depending on the parameters, the log-likelihood is $$ \ell(\beta,\phi)=-\frac{n}{2}\log\phi-\frac{1}{2\phi}(Y_n-X_n\beta)^\mathsf{T}(Y_n-X_n\beta). $$ For each fixed $\phi$, maximizing over $\beta$ is equivalent to minimizing the residual sum of squares, which yields the least-squares estimator above. Substitution of this estimate into the likelihood and maximization over $\phi$ gives $\hat\phi=\mathrm{RSS}/n$. The denominator is $n$ because this is the likelihood estimator; the unbiased estimator instead uses $n-p$.

Result: Under the Gaussian model and full column rank of $X_n$, $$ \hat\beta\sim N\!\left(\beta,\phi(X_n^\mathsf{T}X_n)^{-1}\right),\qquad \frac{\mathrm{RSS}}{\phi}\sim\chi^2_{n-p}. $$ Moreover, $\hat\beta$ and $\mathrm{RSS}$ are independent.

To justify these facts, let $P_X=X_n(X_n^\mathsf{T}X_n)^{-1}X_n^\mathsf{T}$, the orthogonal projection onto the column space of $X_n$. The coefficient error is $$ \hat\beta-\beta=(X_n^\mathsf{T}X_n)^{-1}X_n^\mathsf{T}\varepsilon_n, $$ so it is normal with mean zero and covariance $\phi(X_n^\mathsf{T}X_n)^{-1}$. The residual vector is $(I_n-P_X)\varepsilon_n$, and therefore $\mathrm{RSS}=\varepsilon_n^\mathsf{T}(I_n-P_X)\varepsilon_n$. The matrix $I_n-P_X$ is a symmetric projection of rank $n-p$. An orthogonal change of coordinates diagonalizes it into $n-p$ ones and $p$ zeros. The transformed error coordinates are still independent $N(0,\phi)$, giving the stated chi-squared distribution. In these coordinates, the coefficient error depends only on the $p$ coordinates in the column space, while the residual sum of squares depends only on the $n-p$ orthogonal coordinates. These coordinate sets are independent under the multivariate normal distribution, proving independence.

Worked Example: Exact Coefficient and Dispersion Estimates

Take $n=10$, with an intercept and one predictor. The predictor is $-1$ for the first five observations and $1$ for the other five. The responses in the two groups are, respectively, $1,2,1,0,1$ and $3,4,2,3,3$. Compute the estimates and their model-based covariance.

The group means are $1$ and $3$. With an intercept and a predictor coded $-1,1$, the fitted values are the two group means. Thus $$ \hat\beta_0=\frac{1+3}{2}=2,\qquad \hat\beta_1=\frac{3-1}{2}=1. $$ The residuals are $0,1,0,-1,0$ in the first group and $0,1,-1,0,0$ in the second. Hence $\mathrm{RSS}=1+1+1+1=4$, and the maximum-likelihood dispersion estimate is $\hat\phi=4/10=0.4$.

The design matrix has $X^\mathsf{T}X=\operatorname{diag}(10,10)$, since the predictor sums to zero and its squared values sum to ten. The estimated covariance matrix is therefore $$ \widehat{\operatorname{Cov}}(\hat\beta)=\hat\phi(X^\mathsf{T}X)^{-1} =0.4\begin{pmatrix}1/10&0\\0&1/10\end{pmatrix} =\begin{pmatrix}0.04&0\\0&0.04\end{pmatrix}. $$ Both coefficient standard errors are $0.2$. The unbiased dispersion estimate would instead be $\mathrm{RSS}/(n-p)=4/8=0.5$; it is not the maximum-likelihood estimate.

Joint Large-Sample Behavior

The exact coefficient distribution already contains its asymptotic behavior. The design condition gives $X_n^\mathsf{T}X_n=nQ_n$, so $$ \sqrt n(\hat\beta-\beta)\sim N(0,\phi Q_n^{-1}). $$ Because $Q_n$ converges to the positive-definite matrix $Q$, its inverse converges to $Q^{-1}$. Thus the limiting coefficient covariance on the $\sqrt n$ scale is $\phi Q^{-1}$.

Theorem: Suppose $p$ is fixed, the observations are independent with $Y_n\sim N(X_n\beta,\phi I_n)$, $X_n$ has full column rank for all sufficiently large $n$, and $Q_n=X_n^\mathsf{T}X_n/n\to Q$ for a positive-definite $Q$. Then $$ \begin{pmatrix} \sqrt n(\hat\beta-\beta)\\ \sqrt n(\hat\phi-\phi) \end{pmatrix} \xrightarrow{d} N\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix} \phi Q^{-1}&0\\ 0&2\phi^2 \end{pmatrix} \right). $$

Proof. The exact normal distribution of the coefficient estimator and $Q_n^{-1}\to Q^{-1}$ imply $$ \sqrt n(\hat\beta-\beta)\xrightarrow{d}N(0,\phi Q^{-1}). $$ Set $\nu_n=n-p$. The exact residual result gives $\mathrm{RSS}/\phi\sim\chi^2_{\nu_n}$, and $\hat\phi=\mathrm{RSS}/n$. Therefore $$ \sqrt n(\hat\phi-\phi) =\phi\frac{\sqrt n}{n}\left(\chi^2_{\nu_n}-n\right) =\phi\frac{\sqrt n}{n}\left[ (\chi^2_{\nu_n}-\nu_n)-p \right]. $$ Since $p$ is fixed, the term $\phi p/\sqrt n$ tends to zero. The central limit theorem for sums of independent squared standard normal variables gives $$ \frac{\chi^2_{\nu_n}-\nu_n}{\sqrt{2\nu_n}}\xrightarrow{d}N(0,1). $$ Also, $\sqrt{2\nu_n}/\sqrt n\to\sqrt 2$. It follows that $\sqrt n(\hat\phi-\phi)\xrightarrow{d}N(0,2\phi^2)$. For every $n$, $\hat\beta$ and $\hat\phi$ are independent by the exact normal projection argument. Their joint characteristic function consequently factors into the product of their marginal characteristic functions. Taking limits gives the joint normal distribution with zero off-diagonal blocks displayed in the theorem. $\square$

The theorem describes simultaneous uncertainty, not just two separate standard errors. For large $n$, the coefficient covariance is approximately $\phi Q^{-1}/n$, while the variance of $\hat\phi$ is approximately $2\phi^2/n$. In practice $\phi$ is unknown and is commonly replaced by $\hat\phi$ when estimating coefficient covariance. This substitution estimates uncertainty; it does not change the theoretical limiting covariance.

Worked Example: How Predictor Location Changes Precision

Compare two designs with six observations and an intercept and slope. In the first, the predictor values are $-2,-1,0,0,1,2$. In the second, each value is increased by 10, giving $8,9,10,10,11,12$. Assume the same Gaussian error variance $\phi$ in both designs.

For the first design, the predictor sum is zero and the sum of squared values is ten, so $$ X^\mathsf{T}X=\begin{pmatrix}6&0\\0&10\end{pmatrix},\qquad \operatorname{Var}(\hat\beta)=\phi(X^\mathsf{T}X)^{-1} =\phi\begin{pmatrix}1/6&0\\0&1/10\end{pmatrix}. $$

For the shifted design, the predictor sum is $60$ and its sum of squares is $610$, so $$ X^\mathsf{T}X=\begin{pmatrix}6&60\\60&610\end{pmatrix}. $$ Its determinant is $6(610)-60^2=60$, and hence $$ (X^\mathsf{T}X)^{-1} =\frac{1}{60}\begin{pmatrix}610&-60\\-60&6\end{pmatrix} =\begin{pmatrix}61/6&-1\\-1&1/10\end{pmatrix}. $$ The slope variance remains $\phi/10$, because shifting every predictor by a constant does not change its centered sum of squares. The intercept variance increases from $\phi/6$ to $61\phi/6$, and the intercept and slope estimates are now negatively correlated. The intercept represents the fitted mean at predictor value zero, which lies far outside the shifted design's central range; its precision is therefore much poorer. This is a change in the interpretation and precision of the intercept, not a loss of information about the slope.

Using the Design Limit and Estimated Dispersion

The matrix $Q$ summarizes the asymptotic precision of all coefficient combinations. For a fixed vector $a$, the linear combination $a^\mathsf{T}\hat\beta$ has asymptotic variance $\phi a^\mathsf{T}Q^{-1}a/n$. If $Q$ has a small eigenvalue, then $Q^{-1}$ has a large eigenvalue in the corresponding direction: some combination of coefficients is weakly informed by the design. A large sample size alone does not guarantee good precision if the design fails to accumulate information in every direction.

Worked Example: Uncertainty in a Fitted Mean

Return to the balanced design with five observations at each of $x=-1$ and $x=1$, so $X^\mathsf{T}X=10I_2$. Estimate the mean response at $x=1$, using the estimates from the first worked example.

The fitted mean is $\hat\mu(1)=\hat\beta_0+\hat\beta_1=2+1=3$. It is the linear combination $a^\mathsf{T}\hat\beta$ with $a=(1,1)^\mathsf{T}$. Using $\hat\phi=0.4$, the estimated variance is $$ \hat\phi\,a^\mathsf{T}(X^\mathsf{T}X)^{-1}a =0.4\begin{pmatrix}1&1\end{pmatrix} \begin{pmatrix}0.1&0\\0&0.1\end{pmatrix} \begin{pmatrix}1\\1\end{pmatrix} =0.4(0.2)=0.08. $$ The estimated standard error is $\sqrt{0.08}\approx0.283$. This quantifies uncertainty in the mean response at $x=1$; predicting a future individual response would additionally include its error variance.

Common Mistakes and Practical Interpretation

Using the wrong denominator for dispersion. The Gaussian maximum-likelihood estimate is $\mathrm{RSS}/n$. The unbiased estimate is $\mathrm{RSS}/(n-p)$. They answer different finite-sample estimation criteria, although they have the same first-order large-sample limit when $p$ is fixed.
Ignoring the design matrix. Coefficient precision is not determined by $n$ alone. It is determined by $(X_n^\mathsf{T}X_n)^{-1}$, or asymptotically by $Q^{-1}/n$, multiplied by the dispersion.
Assuming all coefficients are independent. The joint asymptotic coefficient-dispersion covariance has zero off-diagonal blocks in this correctly specified Gaussian model. The coefficient estimates themselves need not be mutually independent: their covariance is governed by $Q^{-1}$.
Applying the limit without its conditions. The result requires fixed $p$, a stable positive-definite design limit, and the stated Gaussian independent-error model. If a design direction becomes uninformative, $Q$ is not positive definite and this limit does not apply in the displayed form.

Key takeaway: In Gaussian linear regression, the coefficient estimator is exactly normal and independent of the residual sum of squares. Under a stable, full-rank design, their joint limiting covariance is block diagonal, with coefficient precision determined by $Q^{-1}$ and dispersion uncertainty determined by $2\phi^2$.

Check Your Understanding

Use the assumptions and calculations in this tutorial to answer the following questions.

Why does the Gaussian maximum-likelihood estimate of dispersion use $n$, while the unbiased estimate uses $n-p$?
For a fixed design, what condition on $X_n^\mathsf{T}X_n/n$ ensures that all coefficient directions have a nondegenerate limiting covariance?
What is the limiting variance of $\sqrt n(\hat\phi-\phi)$, and which distributional result gives it?
Explain why shifting all predictor values by a constant can alter intercept precision without altering slope precision in a simple regression.
How does the covariance of a fitted mean at a predictor vector $x_0$ follow from the coefficient covariance?