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?

  • $\begingroup$ 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. $\endgroup$ – user603 Jun 24 '13 at 14:31
  • 1
    $\begingroup$ 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? $\endgroup$ – Riccardo Jun 24 '13 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.

  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)

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

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 

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.

  • $\begingroup$ 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? $\endgroup$ – Riccardo Jun 24 '13 at 15:19
  • $\begingroup$ by computing $\hat{e}\hat{e}'$. $\endgroup$ – user603 Jun 24 '13 at 16:39
  • $\begingroup$ 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. $\endgroup$ – Rajiv Sambasivan Jun 24 '13 at 16:48
  • $\begingroup$ 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. $\endgroup$ – Riccardo Jun 25 '13 at 18:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.