Variance-covariance matrix of the errors in linear regression

Question

How is the var/cov error matrix calculated by statistical analysis packages in practice?

This idea is clear to me in theory. But not in practice. I mean, if I have a vector of random variables $\textbf{X}=(X_{1}, X_{2}, \ldots, X_{n})^\top$, I understand that the variance/covariance matrix $\Sigma$ will be given the external product of the deviance-from-the-mean vectors: $\Sigma=\mathrm{E}\left[(\textbf{X}-\mathrm{E}(\textbf{X}))(\textbf{X}-\mathrm{E}(\textbf{X}))^\top\right]$.

But when I have a sample, the errors of my observations are not random variables. Or better, they are, but only if I take a number of identical samples from the same population. Otherwise, they're given. So, again my question is: how can a statistical package produce a var/cov matrix starting from a list of observations (i.e. a sample) supplied by the researcher?

Answer 1

The covariance matrix for a model of the type $y = X\beta + \epsilon$ is usually computed as $$(X^t X)^{-1}\frac{\sigma^2}{d}$$ where $\sigma^2$ is the residual sum of squares, $\sigma^2=\sum_i (y_i - X_i\hat\beta)^2$ and $d$ is the degrees of freedom (typically the number of observations minus the number of parameters).

For robust and or clustered standard errors, the product $X^t X$ is modified slightly. There may also be other ways to calculate the covariance matrix, e.g. as suggested by the expectation of outer products.

Answer 2

1. OLS estimation of the error variance, $\sigma^2$:

$$s^2=\frac{\hat \varepsilon^\top\hat \varepsilon}{n-p}$$

This is included in Practical Regression and Anova using R by Julian J. Faraway, page 21 .

Example of its calculation in R, based on linear model of miles-per-gallon regressed on multiple car model specs included in the mtcars database: ols = lm(mpg ~ disp + drat + wt, mtcars). These are the manual calculations and the output of the lm() function:

> rdf = nrow(X) - ncol(X)                    # Residual degrees of freedom
> s.sq = as.vector((t(ols$residuals) %*% ols$residuals) / rdf) 
>                                            # s square (OLS estimate of sigma square)
> (sigma = sqrt(s.sq))                       # Residual standar error
[1] 2.950507
> summary(ols)

Call:
lm(formula = mpg ~ disp + drat + wt, data = mtcars)
...
Residual standard error: 2.951 on 28 degrees of freedom

2. Variance - Covariance matrix of the estimated coefficients, $\hat \beta$:

$$\mathrm{Var}\left[\hat \beta \mid X \right] =\sigma^2 \left(X^\top X\right)^{-1}$$

estimated as in page 8 of this online document as

$$\hat{\mathrm{Var}}\left[\hat \beta \mid X \right] =s^2 \left(X^\top X\right)^{-1}$$

> X = model.matrix(ols)                             # Model matrix X
> XtX = t(X) %*% X                                  # X transpose X
> Sigma = solve(XtX) * s.sq                         # Variance - covariance matrix
> all.equal(Sigma, vcov(ols))                       # Same as built-in formula
[1] TRUE
> sqrt(diag(Sigma))                                 # Calculated Std. Errors of coef's
(Intercept)        disp        drat          wt 
7.099791769 0.009578313 1.455050731 1.217156605 
> summary(ols)[[4]][,2]                             # Output of lm() function
(Intercept)        disp        drat          wt 
7.099791769 0.009578313 1.455050731 1.217156605 

Answer 3

With linear regression we are fitting a model $Y = \beta*X +\varepsilon$. $Y$ is the dependent variable, the $X$'s are the predictor (explanatory) variables. We use the data provided to us (the training set or the sample) to estimate the population $\beta$'s. The $X$'s are not considered random variables. The $Y$'s are random because of the error component.