Using the normal equations to calculate coefficients in multiple linear regression

Question

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

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

You forgot to include a vector of 1 representing the intercept. The rest looks correct. — Michael M, Dec 30 '13 at 16:26
I changed only the 'X' as you suggested, but then I got different coefficients again. — S12000, Dec 30 '13 at 16:36
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](http://stats.stackexchange.com/questions/26176/). — gung - Reinstate Monica, Feb 10 '14 at 04:32

ahwillia · Accepted Answer · 2013-12-30T17:06:34.503

@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

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. — S12000, Dec 30 '13 at 16:40
Those two answers are the same up to floating point roundoff error. You should be comparing the coefficients *relative to their standard errors.* — whuber, Dec 30 '13 at 16:56

Using the normal equations to calculate coefficients in multiple linear regression

Update

1 Answers1

Linked