Linear Regression in R
06.14.2021
Intro
Linear Regression is model that predicts a response based on one or more predictors (columns). This model is one of the most fundemental model and is often used as a baseline in machine learning. In this article, we will learn how to create to do linear regression in R.
Data
For this tutorial, we will use the Boston data set which includes
housing data with features of the houses and their prices. We would like
to predict the medv column or the medium value.
library(MASS)
data(Boston)
str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...Basic Linear Regression in R
To create a basic linear regression in r, we use the lm method. We
supply two parameters to this method. The first parameter is a formula
medv ~ . which means model the medium value parameter by all other
parameters. Then, we supply our data set, Boston.
model = lm(medv ~ ., data = Boston)
summary(model)##
## Call:
## lm(formula = medv ~ ., data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## crim -1.080e-01 3.286e-02 -3.287 0.001087 **
## zn 4.642e-02 1.373e-02 3.382 0.000778 ***
## indus 2.056e-02 6.150e-02 0.334 0.738288
## chas 2.687e+00 8.616e-01 3.118 0.001925 **
## nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
## rm 3.810e+00 4.179e-01 9.116 < 2e-16 ***
## age 6.922e-04 1.321e-02 0.052 0.958229
## dis -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
## rad 3.060e-01 6.635e-02 4.613 5.07e-06 ***
## tax -1.233e-02 3.760e-03 -3.280 0.001112 **
## ptratio -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
## black 9.312e-03 2.686e-03 3.467 0.000573 ***
## lstat -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16From the summary above, we can see the coefficients for our linear model. We can also see the performance metrics are rmse 4.745 and the r2 value is .7406.
Modling Linear Regression in R with Caret
We will now see how to model a linear regression using the Caret
package. We will use this library as it provides us with many features
for real life modeling.
To do this, we use the train method. We pass the same parameters as
above, but in addition we pass the method = 'lm' model to tell Caret
to use a linear model.
library(caret)## Loading required package: lattice
## Loading required package: ggplot2set.seed(1)
model <- train(
medv ~ .,
data = Boston,
method = 'lm'
)
model## Linear Regression
##
## 506 samples
## 13 predictor
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 506, 506, 506, 506, 506, 506, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.914881 0.71529 3.484317
##
## Tuning parameter 'intercept' was held constant at a value of TRUEWe could use summary again to get extra details. We also see that our
RMSE is 4.915 and our Rsquared is .715.
Preprocessing with Caret
One feature that we use from Caret is preprocessing. Often in real life
data science we want to run some pre processing before modeling. We will
center and scale our data by passing the following to the train method:
preProcess = c("center", "scale").
set.seed(1)
model2 <- train(
medv ~ .,
data = Boston,
method = 'lm',
preProcess = c("center", "scale")
)
model2## Linear Regression
##
## 506 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 506, 506, 506, 506, 506, 506, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.914881 0.71529 3.484317
##
## Tuning parameter 'intercept' was held constant at a value of TRUESplitting the Data Set
Often when we are modeling, we want to split our data into a train and
test set. This way, we can check for overfitting. We can use the
createDataPartition method to do this. In this example, we use the
target medv to split into an 80/20 split, p = .80.
This function will return indexes that contains 80% of the data that we should use for training. We then use the indexes to get our training data from the data set.
set.seed(1)
inTraining <- createDataPartition(Boston$medv, p = .80, list = FALSE)
training <- Boston[inTraining,]
testing <- Boston[-inTraining,]We can then fit our model again using only the training data.
set.seed(1)
model3 <- train(
medv ~ .,
data = training,
method = 'lm',
preProcess = c("center", "scale")
)
model3## Linear Regression
##
## 407 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 407, 407, 407, 407, 407, 407, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.895637 0.7167782 3.454826
##
## Tuning parameter 'intercept' was held constant at a value of TRUENow, we want to check our data on the test set. We can use the subset
method to get the features and test target. We then use the predict
method passing in our model from above and the test features.
Finally, we calculate the RMSE and r2 to compare to the model above.
test.features = subset(testing, select=-c(medv))
test.target = subset(testing, select=medv)[,1]
predictions = predict(model3, newdata = test.features)
# RMSE
sqrt(mean((test.target - predictions)^2))## [1] 5.120283# R2
cor(test.target, predictions) ^ 2## [1] 0.712638Cross Validation
In practice, we don’t normal build our data in on training set. It is
common to use a data partitioning strategy like k-fold cross-validation
that resamples and splits our data many times. We then train the model
on these samples and pick the best model. Caret makes this easy with the
trainControl method.
We will use 10-fold cross-validation in this tutorial. To do this we
need to pass three parameters method = "repeatedcv", number = 10
(for 10-fold). We store this result in a variable.
set.seed(1)
ctrl <- trainControl(
method = "cv",
number = 10,
)Now, we can retrain our model and pass the trainControl response to
the trControl parameter. Notice the our call has added
trControl = set.seed.
set.seed(1)
model4 <- train(
medv ~ .,
data = training,
method = 'lm',
preProcess = c("center", "scale"),
trControl = ctrl
)
model4## Linear Regression
##
## 407 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 367, 366, 367, 366, 365, 367, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.738682 0.7399459 3.427828
##
## Tuning parameter 'intercept' was held constant at a value of TRUEThis results seemed to have improved our accuracy for our training data. Let’s check this on the test data to see the results.
test.features = subset(testing, select=-c(medv))
test.target = subset(testing, select=medv)[,1]
predictions = predict(model4, newdata = test.features)
# RMSE
sqrt(mean((test.target - predictions)^2))## [1] 5.120283# R2
cor(test.target, predictions) ^ 2## [1] 0.712638Tuning Hyper Parameters
We can also use caret to tune hyper parameters in models. Linear Regression doesn’t ahve any, so we don’t need to tune the model.