5
$\begingroup$

I am trying to understand how to get the coefficient of a multiple linear regression.

The formula is:

$b = (X'X)^{-1}(X')Y$

I try to calculate $b$ without package and with the lm package inside R.

Doing so, I got different results.

I want to know why. Did I made a mistake? Or does the lm package calculate differently because of the intercept?

> y <-  c(1,2,3,4,5)
> x1 <- c(1,2,3,4,5)
> x2 <- c(1,4,5,7,9)
> Y <- as.matrix(y)
> X <- as.matrix(cbind(x1,x2))
> beta = solve(t(X) %*% X) %*% (t(X) %*% Y) ; beta
            [,1]
x1  1.000000e+00
x2 -1.421085e-14
> model <- lm(y~x1+x2) ; model$coefficients
 (Intercept)           x1           x2 
1.191616e-15 1.000000e+00 1.192934e-15 

3d

Update

As Alex and the other told me, it was a question of roundoff error. Therefore, I decided to take another data from the book "Essential Statistics for business and economics" by Anderson and all. In this case, the coefficients are the same in both lm function and in my own matrix.

> y <- c(9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6, 7.6, 6.1)
> x0 <- c(1,1,1,1,1,1,1,1,1,1) 
> x1 <-  c(100,50,100,100,50,80,75,65,90,90)
> x2 <- c(4,3,4,2,2,2,3,4,3,2)
> Y <- as.matrix(y)
> X <- as.matrix(cbind(x0,x1,x2))

> beta = solve(t(X) %*% X) %*% (t(X) %*% Y);beta
         [,1]
x0 -0.8687015
x1  0.0611346
x2  0.9234254
> model <- lm(y~+x1+x2) ; model$coefficients
(Intercept)          x1          x2 
 -0.8687015   0.0611346   0.9234254 
$\endgroup$
3
  • 7
    $\begingroup$ You forgot to include a vector of 1 representing the intercept. The rest looks correct. $\endgroup$
    – Michael M
    Dec 30, 2013 at 16:26
  • $\begingroup$ I changed only the 'X' as you suggested, but then I got different coefficients again. $\endgroup$
    – S12000
    Dec 30, 2013 at 16:36
  • $\begingroup$ For general information about the effects of suppressing (ie, not fitting) the intercept in a regression model, it may be useful to read this excellent CV thread: Removal of statistically significant intercept term boosts $R^2$ in linear model. $\endgroup$ Feb 10, 2014 at 4:32

1 Answer 1

8
$\begingroup$

@MichaelMayer has it right. Try the following:

> y <-  c(1,2,3,4,5)
> x0 <- c(1,1,1,1,1)   # vector of ones representing the intercept
> x1 <- c(1,2,3,4,5)
> x2 <- c(1,4,5,7,9)
> Y <- as.matrix(y)
> X <- as.matrix(cbind(x0,x1,x2))
> beta = solve(t(X) %*% X) %*% (t(X) %*% Y) ;

Update: Even with the above changes you will still get slightly different results when you use lm due to roundoff error. If the estimated intercept and coefficient for x2 are nonzero, this is no longer visible.

x1 <- c(1,2,3,4,5)
x2 <- c(1,4,5,7,9)
y <-  x1 + x2 + rnorm(5,mean=0,sd=0.3);
Y <- as.matrix(y);
X <- as.matrix(cbind(1,x1,x2));
beta = solve(t(X) %*% X) %*% (t(X) %*% Y) ; beta
model <- lm(y~1+x1+x2) ; model$coefficients

Output is:

         [,1]
   -0.2948504
x1  0.8081534
x2  1.1741777
(Intercept)          x1          x2 
 -0.2948504   0.8081534   1.1741777 

For more info, see: https://en.wikipedia.org/wiki/Machine_epsilon

$\endgroup$
2
  • $\begingroup$ When I do it I get different coefficient between 'lm' function and the manual way. lm gives 1.191616e-15 1.000000e+00 1.192934e-15 and the manual way gives 4.263256e-14 1.000000e+00 0.000000e+00. Where the first number is the intercept then x1 and x2. $\endgroup$
    – S12000
    Dec 30, 2013 at 16:40
  • 7
    $\begingroup$ Those two answers are the same up to floating point roundoff error. You should be comparing the coefficients relative to their standard errors. $\endgroup$
    – whuber
    Dec 30, 2013 at 16:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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