2.2 The lm function
The function lm is the workshorse for fitting linear models. It takes as input a formula: suppose you have a data frame containing columns x (a regressor) and y (the regressand); you can then call lm(y ~ x) to fit the linear model \(y = \beta_0 + \beta_1x + \varepsilon\). The explanatory variable y is on the left hand side,
while the right hand side should contain the predictors, separated by a + sign if there are more than one.
If you provide the data frame name using data, then the shorthand y ~ . fits all the columns of the data frame (but y) as regressors.
To fit higher order polynomials or transformations, use the I function to tell R to interpret the input “as is”.
Thus, lm(y~x+I(x^2)), would fit a linear model with design matrix \((\boldsymbol{1}_n, \mathbf{x}^\top, \mathbf{x}^2)^\top\). A constant is automatically included in the regression, but can be removed by writing -1 or +0 on the right hand side of the formula.
# The function lm and its output
fit <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
fit_summary <- summary(fit)The lm output will display OLS estimates along with standard errors, \(t\) values for the Wald test of the hypothesis \(\mathrm{H}_0: \beta_i=0\) and the associated \(P\)-values. Other statistics and information about the sample size, the degrees of freedom, etc., are given at the bottom of the table.
Many methods allow you to extract specific objects. For example, the functions coef, resid, fitted, model.matrix will return \(\hat{\boldsymbol{\beta}}\), \(\boldsymbol{e}\), \(\hat{\boldsymbol{y}}\) and \(\mathbf{X}\), respectively.
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"The following simply illustrates what has been derived in Exercise series 2. R has devoted functions that are coded more efficiently.
2.2.1 Singular value decomposition
The SVD decomposition in R returns a list with elements u, d and v. u is the orthonormal \(n \times p\) matrix, d is a vector containing the diagonal elements of \(\mathbf{D}\) and v is the \(p \times p\) orthogonal matrix. Recall that the decomposition is
\[\mathbf{X} = \mathbf{UDV}^\top\]
and that \(\mathbf{VV}^\top= \mathbf{V}^\top\mathbf{V}=\mathbf{U}^\top\mathbf{U}=\mathbf{I}_p\). The matrix \(\mathbf{D}\) contains the singular values of \(\mathbf{X}\), and the diagonal elements \(\mathrm{d}_{ii}^2\) corresponds to the (ordered) eigenvalues of \(\mathbf{X}^\top\mathbf{X}\).
The following shows how to use the SVD decomposition in R. This material is optional and provided for reference only.
svdX <- svd(X)
# Projection matrix
Hx <- tcrossprod(svdX$u)
# t(U) %*% U gives p by p identity matrix
all.equal(crossprod(svdX$u), diag(p))## [1] TRUE## [1] TRUE## [1] TRUE## [1] TRUE# OLS coefficient from SVD
beta_hat_svd <- c(svdX$v %*% diag(1/svdX$d) %*% t(svdX$u) %*% y)
all.equal(beta_hat_svd, beta_hat)## [1] TRUE2.2.2 QR decomposition
R uses a QR-decomposition to calculate the OLS estimates in the function lm. There are specific functions to return coefficients, fitted values and residuals. One can also obtain the \(n \times p\) matrix \(\mathbf{Q}_1\) and the upper triangular \(p \times p\) matrix \(\mathbf{R}\) from the thinned QR decomposition,
\[\mathbf{X} = \mathbf{Q}_1\mathbf{R}.\]
The following shows how to use the QR decomposition in R. This material is optional and provided for reference only.
Hx <- X %*% solve(crossprod(X)) %*% t(X)
qrX <- qr(X)
Q1 <- qr.Q(qrX)
R <- qr.R(qrX)
# Compute beta_hat from QR
beta_hat_qr1 <- qr.coef(qrX, y) #using built-in function
beta_hat_qr2 <- c(backsolve(R, t(Q1) %*% y)) #manually
all.equal(beta_hat, beta_hat_qr1, check.attributes = FALSE)## [1] TRUE## [1] TRUE# Compute residuals
qre <- qr.resid(qrX, y)
all.equal(qre, c(y - X %*% beta_hat), check.attributes = FALSE)## [1] TRUE# Compute fitted values
qryhat <- qr.fitted(qrX, y)
all.equal(qryhat, c(X %*% beta_hat), check.attributes = FALSE)## [1] TRUE## [1] TRUE