Simple linear regression fit manually via matrix equations does not match lm() output

Question

I am trying to fit a linear model using matrices to my data set even though I can use OLS and do it without matrices as a simple tutorial for myself to better understand both R and matrix notation.

This is the model I am trying to fit:

$$\bf Y=X\boldsymbol\beta+\varepsilon$$

where $\bf Y$ is a $1\times n$ matrix, $\bf X$ is a $n\times k$ matrix (where $k$ is the number of $\beta$'s, which in this case is 2), $\boldsymbol \beta$ is a $k\times 1$ matrix and lastly our error term is $n\times 1$. I understand this portion.

When I simply use the lm() command to fit my data, I get the following from the summary() command:

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0503 -1.4390  0.4921  1.0589  3.9446 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9849     0.8219   3.632  0.00191 ** 
x             0.5612     0.1084   5.178 6.32e-05 ***

So the summary() is telling me that the $\beta$ matrix is a $2\times 1$ matrix, with the first number (which is $\beta_0$) as 2.0949 and the second number (which is $\beta_1)$ as 0.1084. My question is this:

We know that the matrix $\beta$ is actually:

$$\boldsymbol\beta=(\bf{X}^T\bf{X})^{-1}(\bf{X}^T\bf{Y})$$

and when I try to simply carry out this calculate by hand using R using b=(t(x)*x)^-1*t(x)*y, I get a $1\times 20$ vector (where $20$ of course is $n$, the number of observations). Why am I not getting a $2\times 1$ matrix like I should be getting?

Answer 1

You've made two mistakes in your R code for b.

1. solve is used for matrix inversion. Raising X to the $-1$ power inverts each element of X, which can occasionally be useful, but is not what we want here.
2. R uses the operator %*% for matrix multiplication. Otherwise, it does element-wise multiplication and requires your arrays to be conformable according to R's handling of vectors. Again, occasionally useful.

This is all explained in more detail in the R documentation, or the nigh-uncountable intro to R materials online, such as this one.

b<-solve(t(X)%*%X)%*%t(X)%*%y is the literal representation of the normal equations, while b<-solve(crossprod(X), crossprod(X,y)) is faster and more idiomatic, and matches the output of lm(). But you definitely don't want to explicitly compute the normal equations directly for numerical reasons (and R won't warn you until you lose all significant figures). The lm method uses QR decomposition by default, which is more numerically stable. If you apply the correct "hand cranked" computation and it still differs from the R output, it's plausibly because of loss of numerical precision due to the explicit inversion of $X^TX$.

Don't forget to add a column of $1$ for the intercept!

Answer 2

You need to include a column of ones in your matrix x, which corresponds to the intercept. The lm() function does this for you automatically, but you need to add this yourself when you calculate the answer using the normal equations.

See this previous question on this site: Using the normal equations to calculate coefficients in multiple linear regression. There is example code which you should be able to copy-and-paste directly.

This is a great introductory book that explains why you have to add the column of ones: http://www.amazon.com/Applied-Regression-Models-Edition-Student/dp/0073014664

You will also need to follow the advice given by @user777

Answer 3

This R code can be used to calculate Y (a vector of y values, the fitted values) and Beta (a vector of the coefficients) via matrix regression for a given dataset which I called insert.dataset. This should work even if you add additional numeric variables to the formula.

library(matlib) # enables function inv() to calculate a matrix's inverse

model <- lm(formula = y ~ x, data = insert.dataset) # save linear model
beta0 <- rep(1, nrow(model$model)) # column of 1s representing coefficient beta0 (intercept)
X <- as.matrix(cbind(beta0, model$model[,-1]), nrow=nrow(model$model)) # create X matrix, replacing column of outcomes with beta0

# Matrix equation to create Y (fitted values), using X and coefficients
Y <- X %*% model$coefficients
model$fitted.values # should be identical to Y

# Matrix equation to create Beta (coefficients), using X and Y
Beta <- inv(t(X) %*% X) %*% t(X) %*% model$fitted.values
model$coefficients # should be identical to Beta

It should be noted that the output from lm already has vectors representing Y (model$fitted.values) and Beta model$coefficients, but some small modification is needed to obtain X (the matrix of observed predictor values). That small modification replaces the column containing the observed y values with a column full of 1s (this is needed for multiplication with beta0) and is why you are not getting the 2 * 1 matrix that you want.