# Variance-covariance matrix of the errors in linear regression

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?

• The errors of your observations are function of a random variables (the y's) and are therefore themselves random. Conditional on X alone, they are not given. Jun 24, 2013 at 14:31
• Yes, I fully agree on that. But what you say works in theory. If I draw, say, 100 random samples of identical size from the same population, each observation error will be a random variable with (0, sigma^2). What if, instead, I only draw one sample? In that case, the mean of the error of each observation is the error itself. Is it clear what I am saying? So, what I am trying to understand is, how does a package like Stata calculate the variance-covariance matrix using only one sample drawn from the population? Jun 24, 2013 at 15:02

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.

• Thanks for the clear explanation! Any chance you (or someone else) could expand on your last paragraph? I am currently attempting to calculate the standard errors of a larger-than-memory regression problem, and having to take the matrix inverse to calculate the covariance matrix is making the computation very memory intensive. Jul 29, 2021 at 14:27
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
 2.950507
> summary(ols)

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

1. 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}$$

$$\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
 TRUE
> sqrt(diag(Sigma))                                 # Calculated Std. Errors of coef's
(Intercept)        disp        drat          wt
7.099791769 0.009578313 1.455050731 1.217156605
> summary(ols)[][,2]                             # Output of lm() function
(Intercept)        disp        drat          wt
7.099791769 0.009578313 1.455050731 1.217156605


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.

• Hi Rajiv, thank you for the correction. So, can you explain how Stata (or any other stats package), starting from Y (and epsilon), manages to derive the variance-covariance matrix Sigma? Jun 24, 2013 at 15:19
• by computing $\hat{e}\hat{e}'$. Jun 24, 2013 at 16:39
• Agree with user603. Please check page 21 of cran.r-project.org/doc/contrib/Faraway-PRA.pdf. This is based on R but includes a good discussion of the theory behind linear regression. Jun 24, 2013 at 16:48
• Hi both, thank you, first of all. I also agree with you, user603, and I was expecting this answer. But if the var/cov matrix it's calculated by computing the external product of the error vectors, this means that the cov among the error components in most cases won't be zero as the hypothesis of independence would imply. Right? This is what my doubt revolves around. Rajiv, I looked into the good guide you suggested but couldn't find an answer. Thank you in advance for any future reply. Jun 25, 2013 at 18:03