The following code will attempt to replicate the results of the
lm() function in R. For this exercise, we will be using a cross sectional data set provided by R called “women”, that has height and weight data for 15 individuals.
The OLS regression equation:
where a white noise error term. For this example weight, and height. the marginal impact a one unit change in height has on weight.
## This is the OLS regression we will manually calculate:
reg = lm(weight ~ height, data=women)
Recall that the following matrix equation is used to calculate the vector of estimated coefficients of an OLS regression:
where the matrix of regressor data (the first column is all 1’s for the intercept), and the vector of the dependent variable data.
Matrix operators in R
as.matrix()coerces an object into the matrix class.
t()transposes a matrix.
%*%is the operator for matrix multiplication.
solve()takes the inverse of a matrix. Note, the matrix must be invertible.
For a more complete introduction to doing matrix operations in R, check out this page.
Back to OLS
The following code calculates the 2 x 1 matrix of coefficients, :
## Create X and Y matrices for this specific regression
X = as.matrix(cbind(1,women$height))
Y = as.matrix(women$weight)
## Choose beta-hat to minimize the sum of squared residuals
## resulting in matrix of estimated coefficients:
bh = round(solve(t(X)%*%X)%*%t(X)%*%Y, digits=2)
## Label and organize results into a data frame
beta.hat = as.data.frame(cbind(c("Intercept","Height"),bh))
names(beta.hat) = c("Coeff.","Est")
Calculating Standard Errors
To calculate the standard errors, you must first calculate the variance-covariance (VCV) matrix, as follows:
The VCV matrix will be a square k x k matrix. Standard errors for the estimated coefficients are found by taking the square root of the diagonal elements of the VCV matrix.
## Calculate vector of residuals
res = as.matrix(women$weight-bh-bh*women$height)
## Define n and k parameters
n = nrow(women)
k = ncol(X)
## Calculate Variance-Covariance Matrix
VCV = 1/(n-k) * as.numeric(t(res)%*%res) * solve(t(X)%*%X)
## Standard errors of the estimated coefficients
StdErr = sqrt(diag(VCV))
## Calculate p-value for a t-test of coefficient significance
P.Value = rbind(2*pt(abs(bh/StdErr), df=n-k,lower.tail= FALSE),
2*pt(abs(bh/StdErr), df=n-k,lower.tail= FALSE))
## concatenate into a single data.frame
beta.hat = cbind(beta.hat,StdErr,P.Value)
A Scatterplot with OLS line
## Plot results
plot(women$height,women$weight, xlab = "Height", ylab = "Weight",
main = "OLS: Height and Weight")
abline(a = bh, b = bh, col = ‘red’, lwd = 2, lty="dashed")
Now you can check the results above using the canned
summary(lm(weight ~ height, data = women))