6 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.
The three ingredients of inference are:
- 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.
Standard Error
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 Distribution
The Gaussian distribution is also called the normal distribution.
Normal Distribution
A random variable Z has a normal distribution if its PDF is f_Z(z) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(z-\mu)^2}{2\sigma^2} \right), \qquad z\in\mathbb R.
We write Z\sim\mathcal N(\mu,\sigma^2), where \mu\in\mathbb R and \sigma^2>0. Its mean and variance are E[Z]=\mu, \qquad \mathrm{Var}(Z)=\sigma^2.
The mean determines the center of the distribution. The variance determines its dispersion.
The normal distribution with mean zero and variance one is the standard normal distribution: Z\sim\mathcal N(0,1).
Let \phi denote its PDF and let \Phi denote its CDF. Symmetry implies \phi(a)=\phi(-a), \qquad \Phi(-a)=1-\Phi(a). \tag{6.1}
Thus, intervals and tail probabilities can be calculated symmetrically around zero.
Every normal variable can be standardized: Z\sim\mathcal N(\mu,\sigma^2) \quad\Longrightarrow\quad \frac{Z-\mu}{\sigma}\sim\mathcal N(0,1).
For a standard normal variable, P(-1.96\leq Z\leq1.96)\approx0.95.
6.3 Central Limit Theorem
Convergence in distribution, Z_n\overset{d}{\to}Z, means that the sampling distribution of Z_n approaches the distribution of Z.
Central Limit Theorem (CLT)
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.2}
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.2 does not apply.
To see the normal approximation emerge, choose a non-normal population and increase the sample size:
An analogous result holds for averages of random vectors with finite covariance matrices. Their joint sampling distribution approaches a multivariate normal distribution.
6.4 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_i u_i. \tag{6.3}
Define \boldsymbol\Omega :=E[u_i^2\boldsymbol X_i\boldsymbol X_i'].
Because E[\boldsymbol X_i u_i]=\boldsymbol 0, \boldsymbol\Omega is the covariance matrix of the moment vector \boldsymbol X_i u_i.
The multivariate CLT applies to the sample average of \boldsymbol X_i u_i in Equation 6.3.
Asymptotic Normality of OLS
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.4} 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.4.
The approximate covariance matrix of \widehat{\boldsymbol\beta} is \boldsymbol\Sigma_\beta/n.
6.5 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'.
HC0 Covariance Matrix
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.5}
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} }.
HC0 is the natural estimator for deriving the robust covariance matrix.
HC1 Adjustment
HC1 applies a degrees-of-freedom adjustment: \widehat{\mathrm{Var}}_{HC1} (\widehat{\boldsymbol\beta}) = \frac{n}{n-k} \widehat{\mathrm{Var}}_{HC0} (\widehat{\boldsymbol\beta}), where k is the number of estimated coefficients, including the intercept.
The factor n/(n-k) compensates for estimating k coefficients before calculating the residuals.
HC0 and HC1 are asymptotically equivalent. For routine applications, we use HC1.
HC3, Leverage, and the Jackknife
HC3 uses an observation-specific correction based on leverage.
Let \boldsymbol H := \boldsymbol X (\boldsymbol X'\boldsymbol X)^{-1} \boldsymbol X' denote the hat matrix. It maps observed outcomes into fitted values: \widehat{\boldsymbol Y} =\boldsymbol H\boldsymbol Y.
The i-th diagonal entry h_{ii} =\boldsymbol X_i' (\boldsymbol X'\boldsymbol X)^{-1} \boldsymbol X_i is the leverage value or hat value of observation i.
The leverage values:
- depend only on the regressors,
- satisfy 0\leq h_{ii}\leq1,
- sum to k and therefore have sample average k/n.
The average k/n is a useful reference, not a fixed cutoff for high leverage.
In a simple regression with an intercept, h_{ii} = \frac{1}{n} + \frac{(Z_i-\overline Z)^2} {\sum_{\ell=1}^n(Z_\ell-\overline Z)^2}.
Leverage is therefore high when the regressor value lies far from the center of the sample.
In multiple regression, high leverage indicates an unusual combination of regressor values.
A high-leverage observation is not necessarily influential. Its influence also depends on its residual.
Why Residuals Shrink
Let h_{i\ell} denote entry (i,\ell) of the hat matrix.
The fitted value for observation i contains its own outcome: \widehat Y_i = h_{ii}Y_i + \sum_{\ell\neq i}h_{i\ell}Y_\ell.
A high-leverage observation therefore pulls its fitted value toward its own outcome. This tends to make its OLS residual small.
Let \widehat{\boldsymbol\beta}_{(-i)} be the OLS estimate after observation i is removed.
Provided that the regression remains estimable after this removal, \widehat{\boldsymbol\beta}_{(-i)} -\widehat{\boldsymbol\beta} = -(\boldsymbol X'\boldsymbol X)^{-1} \boldsymbol X_i \frac{\widehat u_i}{1-h_{ii}}. \tag{6.6}
The residual for observation i based on this leave-one-out fit is \widehat u_{i,(-i)} := Y_i-\boldsymbol X_i'\widehat{\boldsymbol\beta}_{(-i)} = \frac{\widehat u_i}{1-h_{ii}}. \tag{6.7}
Equivalently, \widehat u_i =(1-h_{ii})\widehat u_{i,(-i)}.
The ordinary residual is therefore most strongly reduced when leverage is high.
Under homoskedasticity, E[\widehat u_i^2\mid\boldsymbol X] =(1-h_{ii})\sigma^2.
Dividing by one factor 1-h_{ii} would correct this homoskedastic shrinkage.
HC3 instead uses the squared leave-one-out residual from Equation 6.7 and therefore divides by two factors.
HC3 Correction
HC3 Covariance Matrix
HC3 uses the squared leave-one-out residuals: \widehat{\boldsymbol\Omega}_{HC3} := \frac{1}{n}\sum_{i=1}^n \frac{\widehat u_i^2}{(1-h_{ii})^2} \boldsymbol X_i\boldsymbol X_i'.
The HC3 covariance matrix is \widehat{\mathrm{Var}}_{HC3} (\widehat{\boldsymbol\beta}) := \frac{1}{n} \widehat{\boldsymbol Q}^{-1} \widehat{\boldsymbol\Omega}_{HC3} \widehat{\boldsymbol Q}^{-1}.
When h_{ii} is high, (1-h_{ii})^2 is small. HC3 therefore multiplies the squared OLS residual by a larger adjustment.
HC3 does not downweight the observation or change the OLS estimate.
It changes only the estimated sampling uncertainty.
HC3 is also called a jackknife-type robust covariance estimator.
The jackknife removes one observation at a time and measures how much the estimate changes.
From Equation 6.6, \begin{align*} &\sum_{i=1}^n \left( \widehat{\boldsymbol\beta}_{(-i)} -\widehat{\boldsymbol\beta} \right) \left( \widehat{\boldsymbol\beta}_{(-i)} -\widehat{\boldsymbol\beta} \right)'\\ &\qquad= \widehat{\mathrm{Var}}_{HC3} (\widehat{\boldsymbol\beta}). \end{align*}
This identity explains the jackknife terminology.
HC0, HC1, and HC3 are asymptotically equivalent. Their differences are finite-sample adjustments.
HC1 remains our routine default. HC3 is especially useful in smaller samples or when leverage values are unusually large.
A noticeable difference between HC1 and HC3 shows that the estimated uncertainty is sensitive to the finite-sample leverage correction.
Neither estimator repairs a misspecified regression function.
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 the lack of a credible causal 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.
All three columns compare school districts with the same percentages of English learners and students eligible for subsidized lunch.
According to the fitted regression, those with one more student per teacher have, on average, test scores 0.998 points lower.
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.6 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).
Large-Sample Confidence Interval
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.8}
For a 95% interval, \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.
It is therefore incorrect to say that the fixed parameter has a 95% probability of lying in the particular interval we calculated.
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 mean test-score difference for one more student per teacher in this model.
The comparison is between districts with the same percentages of English learners and students eligible for subsidized lunch.
It does not establish that changing class size would cause test scores to change.
6.7 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.9}
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.9 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. A causal interpretation requires a credible causal argument for the underlying comparison.
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.10}
This result is called the delta method.
For example, the exact percentage change in the conditional geometric mean for a one-unit increase in a log-linear 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.
6.8 Classical Inference
The preceding confidence intervals use HC1 standard errors and a large-sample normal approximation.
The traditional bachelor-level approach strengthens the assumptions and obtains exact finite-sample inference.
Homoskedastic Covariance Matrix
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.
Here, k is the number of estimated coefficients, including the intercept.
The homoskedastic standard error \mathrm{se}_{hom}(\widehat\beta_j) is the square root of diagonal entry j.
Without homoskedasticity, this covariance matrix estimator is generally invalid.
Gaussian Regression Model
Let \boldsymbol u=(u_1,\ldots,u_n)' collect the population errors.
Gaussian Regression Model
Conditional on the regressor matrix \boldsymbol X, the errors are jointly normal with mean zero and the same variance: \boldsymbol u\mid\boldsymbol X \sim\mathcal N(\boldsymbol 0,\sigma^2\boldsymbol I_n).
This assumption concerns the conditional errors. It does not require the regressors to be normally distributed.
Chi-Squared Distribution
Chi-Squared Distribution
Let Z_1,\ldots,Z_\nu be independent standard normal variables.
Then V=\sum_{\ell=1}^{\nu}Z_\ell^2 \sim\chi_\nu^2.
The parameter \nu is called the degrees of freedom.
The chi-squared distribution takes only nonnegative values.
Under the Gaussian regression model, \frac{(n-k)s_{\widehat u}^2}{\sigma^2} \ \Big|\ \boldsymbol X \sim\chi_{n-k}^2.
Student t-Distribution
Student t-Distribution
Let Z\sim\mathcal N(0,1), \qquad V\sim\chi_\nu^2, and suppose that Z and V are independent.
Then T=\frac{Z}{\sqrt{V/\nu}} \sim t_\nu.
The t_\nu distribution is symmetric and has heavier tails than the standard normal distribution.
It approaches \mathcal N(0,1) as \nu increases.
Exact t-Interval
Let t_{\nu,p} denote the p-quantile of the t_\nu distribution.
Under the Gaussian regression model, Z_j := \frac{\widehat\beta_j-\beta_j} {\sigma \sqrt{[(\boldsymbol X'\boldsymbol X)^{-1}]_{jj}}} \ \Big|\ \boldsymbol X \sim\mathcal N(0,1).
Moreover, Z_j and (n-k)s_{\widehat u}^2/\sigma^2 are conditionally independent.
Therefore, the definition of the t-distribution gives the following result.
Exact t-Inference
Suppose the Gaussian regression model holds, n>k, and the regressor matrix has full column rank.
Then \frac{\widehat\beta_j-\beta_j} {\mathrm{se}_{hom}(\widehat\beta_j)} \ \Big|\ \boldsymbol X \sim t_{n-k}.
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].
This interval has conditional coverage 1-\alpha.
The two approaches can now be compared:
| Main practical route | Classical route |
|---|---|
| Heteroskedasticity allowed | Homoskedasticity required |
| HC1 standard errors | Homoskedastic standard errors |
| Large-sample normal approximation | Exact finite-sample t-distribution |
| No conditional normality required | Conditional normality required |
Under generic heteroskedasticity, using a t critical value is a finite-sample convention rather than an exact result.
Homoskedastic standard errors can be consistent without normality. Conditional normality is the additional assumption that makes the t-result exact.
6.9 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.