3
$\begingroup$

I'm working on a linear model like:

enter image description here

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:

enter image description here

enter image description here

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:

enter image description here

enter image description here

enter image description here

And

enter image description here

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Roland Jan 4 '18 at 12:35
  • $\begingroup$ 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? $\endgroup$ – s_t Jan 4 '18 at 12:45
  • $\begingroup$ Please add all code you used to the question. $\endgroup$ – Roland Jan 4 '18 at 12:46
  • 1
    $\begingroup$ Where did you get the formula for sigma? It's not correct. See help("summary.lm") (under "Value") for the correct definition. $\endgroup$ – Roland Jan 4 '18 at 12:51
  • 1
    $\begingroup$ Unweighted regression is the same as all weights equal to 1. $\endgroup$ – Roland Jan 4 '18 at 14:33
1
$\begingroup$

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

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

enter image description here

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
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – s_t Jan 5 '18 at 8:34
  • $\begingroup$ perfect! Maybe you should edit it again and state it as solved :). Greetings from germany! $\endgroup$ – bucky Jan 5 '18 at 8:38
0
$\begingroup$

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

enter image description here

$\endgroup$

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.