Variance of the estimators in Linear Model with manual formulae in R

I'm working on a linear model like:

I'm applying the formulae to the data to calculate the estimators of a and b using R. First I got some data like:

x <- c(10,11,23,24)
y <- c(4,2,7,8)


Then I applied the lm() to have the results:

model <- lm(y~x)          # the model

coefficients(model)       # the coefficients
(Intercept)           x
-0.5500000   0.3411765

sqrt(diag(vcov(model)))   # the s.e.
(Intercept)           x
1.69932945  0.09333313


So I decided to calculate manually the coefficients and the s.e. I used this formulae for the coefficients:

And in R I did something like:

bHat <- cov(x,y)/var(x)
aHat <- mean(y)-(cov(x,y)/var(x))*mean(x)

bHat
[1] 0.3411765
aHat
[1] -0.55


Then I calculate the s.e.. I used those formulae:

And

So I tried those formulae in R.

se.bHat <- sqrt(sigma/sum((x-(mean(x)))^2))
se.bHat
[1] 0.211205

se.aHat <- sqrt((sigma/4)*  (1+4*mean(x)^2)/sum((x-(mean(x)))^2))
se.aHat
[1] 3.590497


But the result is not equal to the lm() output. Am I using the wrong formulae in the theory or am I applying them wrongly? Thanks in advance.

EDIT:

I've tried as suggested this as the formula σ^2 = 1/(n-p) Sum(w[i] R[i]^2) with help(summary.lm):

sigma <- (1/(4-2))*sum(residuals(model)^2)
se.bHat <- sqrt(sigma/sum((x-(mean(x)))^2))
se.bHat
[1] 0.09333313


And the result is not ok. I'm going to try also with the intercept.

• How do you calculate sigma? Note that it is squared in the formula. Also, your code for se.aHat doesn't match the formula you provide. Commented Jan 4, 2018 at 12:35
• The formula of the sigma is after the var(b) formula, however is var(y) (but I can be wrong). Edited the formula. The theory formula are correct in your opinion or wrong?
– s__
Commented Jan 4, 2018 at 12:45
• Please add all code you used to the question. Commented Jan 4, 2018 at 12:46
• Where did you get the formula for sigma? It's not correct. See help("summary.lm") (under "Value") for the correct definition. Commented Jan 4, 2018 at 12:51
• Unweighted regression is the same as all weights equal to 1. Commented Jan 4, 2018 at 14:33

due to the fact, that i can´t give a comment:

i know a different formula for the s.e of a:

I have it from the book: Applied Regression Analysis, Linear Models, and Related Methods by John Fox (1997).

Using it in R: (maybe a too many brackets :D)

> (se.aHat <- ((sigma*sum(x^2))/(n*(sum((x-mean(x))^2)))))
[1] 2.887721
> sqrt(se.aHat)
[1] 1.699329

• Thanks! The problem was that I used a wrong formula for sigma square, a bit different from mine. Using the command help("summary.lm") as suggest by Roland, I've found that the right formula is σ^2 = 1/(n-p) Sum(w[i] R[i]^2). Your formula is however correct, it gives the same result of mine.
– s__
Commented Jan 5, 2018 at 8:34
• perfect! Maybe you should edit it again and state it as solved :). Greetings from germany! Commented Jan 5, 2018 at 8:38

The problem is that I used the wrong formula for sigma square. using:

help("summary.lm")


As suggested by Roland (see comments below), you can find that the right formula is as stated by R guide:

the square root of the estimated variance of the random error

σ^2 = 1/(n-p) Sum(w[i] R[i]^2)

where R[i] is the i-th residual, residuals[i].

So