I wanted to know how I can manually compute what is returned when I call the vcov() function in R on a lm object, i.e. the variance-covariance matrix.
The diagonal entries are very straightforward, one just needs to take the st. error of every parameters and square it. For the estimates under for example, entry (z,z) of the variance-covariance matrix is simply computed by taking z's estimated coefficient's standard error (0.1096658) and squaring it. But what about the non-diagonal entries? Could someone provide the code to manually compute these?
# Minimal working example code
library(data.table)
df <- data.table(y = runif(100, 0, 100),
x = runif(100, 0, 100),
z = runif(100, 0, 100))
coef(lm(y ~ x + z, df)) # coefficients
coef(summary(lm(y ~ x + z, df)))[,2] # standard errors
vcov(lm(y ~ x + z, df)) # variance-covariance matrix
# Example result
coef(lm(y ~ x + z, df))
(Intercept) x z
54.44460066 0.03052479 -0.11197309
coef(summary(lm(y ~ x + z, df)))[,2]
(Intercept) x z
8.4668880 0.1077858 0.1096658
vcov(lm(y ~ x + z, df))
(Intercept) x z
(Intercept) 71.6881931 -0.576781121 -0.585380616
x -0.5767811 0.011617776 -0.001059506
z -0.5853806 -0.001059506 0.012026594
For a "manual computation" you can check this
library(data.table)
df = data.table(y = runif(100, 0, 100),
x = runif(100, 0, 100),
z = runif(100, 0, 100))
mod = lm(y ~ ., df)
X = cbind(1, as.matrix(df[, -1]))
invXtX = solve(crossprod(X))
coef = invXtX %*% t(X) %*% df$y
resid = df$y - X %*% coef
df.res = nrow(X) - ncol(X)
manual = invXtX * sum(resid**2)/(df.res)
funct = vcov(mod)
all.equal(manual, funct, check.attributes = F)# should be TRUE