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.
sigma
? Note that it is squared in the formula. Also, your code forse.aHat
doesn't match the formula you provide. $\endgroup$ – Roland Jan 4 '18 at 12:35help("summary.lm")
(under "Value") for the correct definition. $\endgroup$ – Roland Jan 4 '18 at 12:51