6 Regression Inference
Section 5 interpreted estimated regression coefficients and effects.
A different sample from the same population would generally give different estimates.
Statistical inference quantifies this sampling uncertainty.
- A point estimate is a single estimate of the unknown population parameter.
- A standard error measures the sampling uncertainty of the estimate.
- A critical value is chosen from the approximating distribution and determines the interval’s width.
This section uses these ingredients to construct confidence intervals.
Hypothesis tests follow in Section 7.
6.1 Estimation Uncertainty
Let \widehat\theta_n estimate a population parameter \theta from an i.i.d. sample of size n.
The sampling distribution of \widehat\theta_n describes the estimates obtained from all possible samples of this size.
Its standard deviation \mathrm{sd}(\widehat\theta_n) measures how much the estimator varies across samples.
This quantity is generally unknown. A standard error estimates it using the observed sample.
A valid standard error satisfies \frac{\mathrm{se}(\widehat\theta_n)} {\mathrm{sd}(\widehat\theta_n)} \overset{p}{\to}1.
Thus, the standard error becomes a reliable measure of sampling uncertainty as the sample size increases.
A small standard error means that the estimator would vary relatively little across repeated samples.
It does not imply that the model or its causal interpretation is correct.
6.2 Normal Approximation and the CLT
A normally distributed random variable with mean \mu and variance \sigma^2 is written as Z\sim\mathcal N(\mu,\sigma^2).
The standard normal distribution \mathcal N(0,1) is symmetric around zero.
About 95% of its probability lies between -1.96 and 1.96.
Convergence in distribution, Z_n\overset{d}{\to}Z, means that the sampling distribution of Z_n approaches the distribution of Z.
Let H_1,H_2,\ldots be i.i.d. with E[H_i]=\theta, \qquad 0<\mathrm{Var}(H_i)=\tau^2<\infty.
Then the sample mean \overline H_n satisfies \sqrt n(\overline H_n-\theta) \overset{d}{\to}\mathcal N(0,\tau^2). \tag{6.1}
Equivalently, \frac{\overline H_n-\theta}{\tau/\sqrt n} \overset{d}{\to}\mathcal N(0,1).
The LLN from Section 3 shows that \overline H_n approaches \theta.
The CLT additionally describes the rate and the approximate shape of its sampling distribution.
Very skewed or heavy-tailed distributions may require a large sample before the normal approximation is accurate.
If \tau^2 is infinite, Equation 6.1 does not apply.
The same theorem applies to averages of random vectors. Their joint sampling distribution approaches a multivariate normal distribution.
6.3 Asymptotic Normality of OLS
Consider the exact linear CEF from Section 4: Y_i=\boldsymbol X_i'\boldsymbol\beta+U_i, \qquad E[U_i\mid\boldsymbol X_i]=0.
Recall the population and sample second-moment matrices \boldsymbol Q=E[\boldsymbol X_i\boldsymbol X_i'], \qquad \widehat{\boldsymbol Q} =\frac{1}{n}\sum_{i=1}^n \boldsymbol X_i\boldsymbol X_i'.
Substituting the regression equation into the OLS estimator gives \sqrt n(\widehat{\boldsymbol\beta}-\boldsymbol\beta) =\widehat{\boldsymbol Q}^{-1} \frac{1}{\sqrt n}\sum_{i=1}^n\boldsymbol X_iU_i. \tag{6.2}
Define \boldsymbol\Omega :=E[U_i^2\boldsymbol X_i\boldsymbol X_i'].
Because E[\boldsymbol X_iU_i]=\boldsymbol 0, \boldsymbol\Omega is the covariance matrix of the moment vector \boldsymbol X_iU_i.
The multivariate CLT applies to the sample average of \boldsymbol X_iU_i in Equation 6.2.
Suppose:
- (Y_i,\boldsymbol X_i')' are i.i.d.,
- the number of regressors k is fixed,
- the exact linear CEF is correctly specified,
- \boldsymbol Q and \boldsymbol\Omega have finite entries,
- and \boldsymbol Q is invertible.
Then \sqrt n(\widehat{\boldsymbol\beta}-\boldsymbol\beta) \overset{d}{\to} \mathcal N(\boldsymbol 0,\boldsymbol\Sigma_\beta), \tag{6.3} where \boldsymbol\Sigma_\beta =\boldsymbol Q^{-1}\boldsymbol\Omega\boldsymbol Q^{-1}.
Thus, for a large sample, \widehat{\boldsymbol\beta} \overset{\mathrm{approx.}}{\sim} \mathcal N\left( \boldsymbol\beta, \frac{1}{n}\boldsymbol\Sigma_\beta \right).
The matrix \boldsymbol\Sigma_\beta is the covariance matrix of the limiting distribution in Equation 6.3.
The approximate covariance matrix of \widehat{\boldsymbol\beta} is \boldsymbol\Sigma_\beta/n.
6.4 Robust Standard Errors
Because E[U_i\mid\boldsymbol X_i]=0, the conditional error variance is \sigma^2(\boldsymbol x) :=\mathrm{Var}(U_i\mid\boldsymbol X_i=\boldsymbol x) =E[U_i^2\mid\boldsymbol X_i=\boldsymbol x].
The errors are homoskedastic if this variance is constant.
They are heteroskedastic if it changes with \boldsymbol x.
In the general case, \boldsymbol\Omega =E[\sigma^2(\boldsymbol X_i) \boldsymbol X_i\boldsymbol X_i'].
We do not need to assume homoskedasticity to estimate \boldsymbol\Omega.
HC0 Covariance Matrix
Let \widehat u_i =Y_i-\boldsymbol X_i'\widehat{\boldsymbol\beta} denote the OLS residual.
Replacing the unknown U_i^2 by \widehat u_i^2 gives \widehat{\boldsymbol\Omega}_{HC0} :=\frac{1}{n}\sum_{i=1}^n \widehat u_i^2\boldsymbol X_i\boldsymbol X_i'.
The resulting heteroskedasticity-consistent covariance matrix estimator is \widehat{\mathrm{Var}}_{HC0} (\widehat{\boldsymbol\beta}) :=\frac{1}{n} \widehat{\boldsymbol Q}^{-1} \widehat{\boldsymbol\Omega}_{HC0} \widehat{\boldsymbol Q}^{-1}. \tag{6.4}
The HC0 standard error of \widehat\beta_j is the square root of diagonal entry j: \mathrm{se}_{HC0}(\widehat\beta_j) :=\sqrt{ \left[ \widehat{\mathrm{Var}}_{HC0} (\widehat{\boldsymbol\beta}) \right]_{jj} }.
Finite-Sample Adjustments
HC0, HC1, and HC3 use the same large-sample argument.
They differ in how they correct the squared residuals in finite samples:
- HC0 uses \widehat u_i^2 directly.
- HC1 multiplies the HC0 covariance matrix by n/(n-k).
- HC3 makes a larger adjustment for observations with high leverage.
All three estimators are asymptotically equivalent.
HC0 is the natural estimator for deriving the robust covariance matrix.
For routine applications, we use HC1.
The leverage of observation i is h_{ii} =\boldsymbol X_i' \left(\sum_{\ell=1}^n \boldsymbol X_\ell\boldsymbol X_\ell'\right)^{-1} \boldsymbol X_i.
It measures how unusual the observation’s regressor values are relative to the sample.
The leverage values lie between zero and one and have sample average k/n.
An observation with leverage well above k/n may have a strong effect on the fitted regression.
HC3 modifies \widehat{\boldsymbol\Omega}_{HC0} by replacing \widehat u_i^2 with \frac{\widehat u_i^2}{(1-h_{ii})^2}.
This corrects more strongly for the tendency of OLS residuals to be small at high-leverage observations.
HC3 can therefore be useful when high leverage is a concern, especially in smaller samples.
Robust standard errors allow the conditional error variance to change with the regressors.
They do not repair an incorrect CEF, omitted-variable bias, bad controls, or a missing causal identification argument.
Under misspecification, sandwich inference may still apply to the best linear approximation from Section 4. Its coefficients need not be marginal effects of the actual CEF.
HC standard errors also require independent observations. Clustered and time-series dependence require different covariance estimators, discussed in Section 8.
Class Size Example
We return to Model (2) from the class-size application in Section 5:
The table reports the same coefficient estimates with three standard-error estimators.
Standard errors appear in parentheses.
| Homoskedastic-only | HC1 | HC3 | |
|---|---|---|---|
| Student-teacher ratio | -0.998 | -0.998 | -0.998 |
| (0.239) | (0.270) | (0.274) | |
| Num.Obs. | 420 | 420 | 420 |
The first column assumes a constant conditional error variance.
The point estimate is -0.998 in every column.
Only the estimated sampling uncertainty changes: the standard errors are 0.239, 0.270, and 0.274.
A robust standard error can be smaller or larger than its homoskedastic-only counterpart.
For the remaining examples, we use HC1.
6.5 Confidence Intervals
Let z_p denote the p-quantile of the standard normal distribution: P(Z\leq z_p)=p, \qquad Z\sim\mathcal N(0,1).
Under the conditions above, a valid robust standard error gives \frac{\widehat\beta_j-\beta_j} {\mathrm{se}_{HC1}(\widehat\beta_j)} \overset{d}{\to}\mathcal N(0,1).
An approximate (1-\alpha) confidence interval for \beta_j is therefore I_{1-\alpha} =\left[ \widehat\beta_j -z_{1-\alpha/2}\mathrm{se}_{HC1}(\widehat\beta_j), \ \widehat\beta_j +z_{1-\alpha/2}\mathrm{se}_{HC1}(\widehat\beta_j) \right]. \tag{6.5}
For a 95% interval, \boxed{ \widehat\beta_j \ \pm\ 1.96\,\mathrm{se}_{HC1}(\widehat\beta_j) }.
A confidence interval is random because its endpoints depend on the sample.
The population parameter is fixed but unknown.
A valid interval procedure has the coverage property \lim_{n\to\infty} P(\beta_j\in I_{1-\alpha}) =1-\alpha.
Thus, about 95% of the intervals constructed from repeated large samples cover the true parameter.
The particular interval calculated from the observed sample either covers the parameter or it does not.
Class Size Confidence Interval
Code
class_size_estimate = coef(class_size_fit)[["STR"]]
class_size_se = sqrt(diag(
sandwich::vcovHC(class_size_fit, type = "HC1")
))[["STR"]]
class_size_ci = class_size_estimate +
qnorm(0.975) * c(-1, 1) * class_size_se
round(c(
estimate = class_size_estimate,
standard_error = class_size_se,
lower = class_size_ci[1],
upper = class_size_ci[2]
), 3) estimate standard_error lower upper
-0.998 0.270 -1.528 -0.469 The approximate 95% HC1 confidence interval is [-1.528,-0.469].
It describes uncertainty about the population coefficient in this specified conditional-mean model.
It does not establish that changing class size would cause test scores to change.
6.6 Confidence Intervals for Effects
Many effects from Section 5 depend on several regression coefficients.
Suppose the effect of interest is the linear combination \theta=\boldsymbol r'\boldsymbol\beta.
Its estimator and robust standard error are \widehat\theta =\boldsymbol r'\widehat{\boldsymbol\beta}, \qquad \mathrm{se}(\widehat\theta) =\sqrt{ \boldsymbol r' \widehat{\mathrm{Var}}_{HC1} (\widehat{\boldsymbol\beta}) \boldsymbol r }. \tag{6.6}
The covariance terms matter because the coefficient estimates generally move together across samples.
For the quadratic experience specification, ME_A(a)=\gamma_2+2a\gamma_3.
Using \boldsymbol\gamma=(\gamma_1,\gamma_2,\gamma_3)', set \boldsymbol r(a)=(0,1,2a)'.
Then Equation 6.6 gives the uncertainty of the estimated marginal effect at experience level a.
a = 20
r = c(0, 1, 2 * a)
effect_estimate = sum(r * coef(quadratic_fit))
effect_se = sqrt(drop(t(r) %*% quadratic_vcov %*% r))
effect_ci = effect_estimate + qnorm(0.975) * c(-1, 1) * effect_se
round(c(
estimate = effect_estimate,
standard_error = effect_se,
lower = effect_ci[1],
upper = effect_ci[2]
), 3) estimate standard_error lower upper
0.223 0.007 0.209 0.237 At 20 years of experience, the estimated marginal effect is USD 0.223 in hourly wage per year of experience.
Its approximate 95% HC1 confidence interval is [0.209,0.237].
Each vertical slice of the shaded area is a confidence interval at one fixed experience level.
It is not a simultaneous confidence band for the entire curve.
As in Section 5, these estimates describe conditional associations unless the causal identification conditions hold.
Suppose an effect is a differentiable nonlinear function \theta=g(\boldsymbol\beta).
A first-order approximation gives \mathrm{se}(g(\widehat{\boldsymbol\beta})) \approx \sqrt{ \nabla g(\widehat{\boldsymbol\beta})' \widehat{\mathrm{Var}}(\widehat{\boldsymbol\beta}) \nabla g(\widehat{\boldsymbol\beta}) }. \tag{6.7}
This result is called the delta method.
For example, the exact percentage effect in a log-level model is g(\beta_2)=100(e^{\beta_2}-1).
Inference for a population AME may additionally need to account for estimating the population distribution of the regressors.
If the errors are homoskedastic, then \boldsymbol\Omega=\sigma^2\boldsymbol Q.
The classical homoskedastic covariance matrix estimator is \widehat{\mathrm{Var}}_{hom} (\widehat{\boldsymbol\beta}) =s_{\widehat u}^2 (\boldsymbol X'\boldsymbol X)^{-1}, where s_{\widehat u}^2 =\frac{1}{n-k}\sum_{i=1}^n\widehat u_i^2.
Without homoskedasticity, this covariance matrix estimator is generally invalid.
Suppose, in addition, that \boldsymbol U\mid\boldsymbol X \sim\mathcal N(\boldsymbol 0,\sigma^2\boldsymbol I_n).
This Gaussian regression model imposes conditional normality and homoskedasticity.
Then \frac{\widehat\beta_j-\beta_j} {\mathrm{se}_{hom}(\widehat\beta_j)} \ \Big|\ \boldsymbol X \sim t_{n-k}.
The t_{n-k} distribution is symmetric and approaches \mathcal N(0,1) as n-k increases.
Using the homoskedastic standard error gives the exact conditional interval \left[ \widehat\beta_j -t_{n-k,1-\alpha/2}\mathrm{se}_{hom}(\widehat\beta_j), \ \widehat\beta_j +t_{n-k,1-\alpha/2}\mathrm{se}_{hom}(\widehat\beta_j) \right].
Under generic heteroskedasticity, using a t critical value is a finite-sample convention rather than an exact result.
6.7 Inference Workflow
For a coefficient or effect:
- Calculate the point estimate.
- Select a covariance estimator that matches the sampling structure.
- Obtain the standard error.
- Combine the estimate, standard error, and critical value.
- Interpret the interval for the population parameter actually estimated.
Section 7 uses the same ingredients for hypothesis testing.