This guide walks through an example of how to conduct multiple linear regression in R, including:
- Examining the data before fitting the model
- Fitting the model
- Checking the assumptions of the model
- Interpreting the output of the model
- Assessing the goodness of fit of the model
- Using the model to make predictions
Let’s jump in!
Setup
For this example we will use the built-in R dataset mtcars, which contains information about various attributes for 32 different cars:
#view first six lines of mtcars
head(mtcars)
# mpg cyl disp hp drat wt qsec vs am gear carb
#Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
#Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
#Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
#Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
#Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
#Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
In this example we will build a multiple linear regression model that uses mpg as the response variable and disp, hp, and drat as the predictor variables.
#create new data frame that contains only the variables we would like to use to data #view first six rows of new data frame head(data) # mpg disp hp drat #Mazda RX4 21.0 160 110 3.90 #Mazda RX4 Wag 21.0 160 110 3.90 #Datsun 710 22.8 108 93 3.85 #Hornet 4 Drive 21.4 258 110 3.08 #Hornet Sportabout 18.7 360 175 3.15 #Valiant 18.1 225 105 2.76
Examining the Data
Before we fit the model, we can examine the data to gain a better understanding of it and also visually assess whether or not multiple linear regression could be a good model to fit to this data.
In particular, we need to check if the predictor variables have a linear association with the response variable, which would indicate that a multiple linear regression model may be suitable.
To do so, we can use the pairs() function to create a scatterplot of every possible pair of variables:
pairs(data, pch = 18, col = "steelblue")
From this pairs plot we can see the following:
- mpg and disp appear to have a strong negative linear correlation
- mpg and hp appear to have a strong positive linear correlation
- mpg and drat appear to have a modest negative linear correlation
Note that we could also use the ggpairs() function from the GGally library to create a similar plot that contains the actual linear correlation coefficients for each pair of variables:
#install and load the GGally library install.packages("GGally") library(GGally) #generate the pairs plot ggpairs(data)
Each of the predictor variables appears to have a noticeable linear correlation with the response variable mpg, so we’ll proceed to fit the linear regression model to the data.
Fitting the Model
The basic syntax to fit a multiple linear regression model in R is as follows:
lm(response_variable ~ predictor_variable1 + predictor_variable2 + ..., data = data)
Using our data, we can fit the model using the following code:
model
Checking Assumptions of the Model
Before we proceed to check the output of the model, we need to first check that the model assumptions are met. Namely, we need to verify the following:
1. The distribution of model residuals should be approximately normal.
We can check if this assumption is met by creating a simple histogram of residuals:
hist(residuals(model), col = "steelblue")
Although the distribution is slightly right skewed, it isn’t abnormal enough to cause any major concerns.
2. The variance of the residuals should be consistent for all observations.
This preferred condition is known as homoskedasticity. Violation of this assumption is known as heteroskedasticity.
To check if this assumption is met we can create a fitted value vs. residual plot:
#create fitted value vs residual plot plot(fitted(model), residuals(model)) #add horizontal line at 0 abline(h = 0, lty = 2)
Ideally we would like the residuals to be equally scattered at every fitted value. We can see from the plot that the scatter tends to become a bit larger for larger fitted values, but this pattern isn’t extreme enough to cause too much concern.
Interpreting the Output of the Model
Once we’ve verified that the model assumptions are sufficiently met, we can look at the output of the model using the summary() function:
summary(model) #Call: #lm(formula = mpg ~ disp + hp + drat, data = data) # #Residuals: # Min 1Q Median 3Q Max #-5.1225 -1.8454 -0.4456 1.1342 6.4958 # #Coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) 19.344293 6.370882 3.036 0.00513 ** #disp -0.019232 0.009371 -2.052 0.04960 * #hp -0.031229 0.013345 -2.340 0.02663 * #drat 2.714975 1.487366 1.825 0.07863 . #--- #Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # #Residual standard error: 3.008 on 28 degrees of freedom #Multiple R-squared: 0.775, Adjusted R-squared: 0.7509 #F-statistic: 32.15 on 3 and 28 DF, p-value: 3.28e-09
From the output we can see the following:
- The overall F-statistic of the model is 32.15 and the corresponding p-value is 3.28e-09. This indicates that the overall model is statistically significant. In other words, the regression model as a whole is useful.
- disp is statistically significant at the 0.10 significance level. In particular, the coefficient from the model output tells is that a one unit increase in disp is associated with a -0.019 unit decrease, on average, in mpg, assuming hp and drat are held constant.
- hp is statistically significant at the 0.10 significance level. In particular, the coefficient from the model output tells is that a one unit increase in hp is associated with a -0.031 unit decrease, on average, in mpg, assuming disp and drat are held constant.
- drat is statistically significant at the 0.10 significance level. In particular, the coefficient from the model output tells is that a one unit increase in drat is associated with a 2.715 unit increase, on average, in mpg, assuming disp and hp are held constant.
Assessing the Goodness of Fit of the Model
To assess how “good” the regression model fits the data, we can look at a couple different metrics:
1. Multiple R-Squared
This measures the strength of the linear relationship between the predictor variables and the response variable. A multiple R-squared of 1 indicates a perfect linear relationship while a multiple R-squared of 0 indicates no linear relationship whatsoever.
Multiple R is also the square root of R-squared, which is the proportion of the variance in the response variable that can be explained by the predictor variables. In this example, the multiple R-squared is 0.775. Thus, the R-squared is 0.7752 = 0.601. This indicates that 60.1% of the variance in mpg can be explained by the predictors in the model.
Related: What is a Good R-squared Value?
2. Residual Standard Error
This measures the average distance that the observed values fall from the regression line. In this example, the observed values fall an average of 3.008 units from the regression line.
Related: Understanding the Standard Error of the Regression
Using the Model to Make Predictions
From the output of the model we know that the fitted multiple linear regression equation is as follows:
mpghat = -19.343 – 0.019*disp – 0.031*hp + 2.715*drat
We can use this equation to make predictions about what mpg will be for new observations. For example, we can find the predicted value of mpg for a car that has the following attributes:
- disp = 220
- hp = 150
- drat = 3
#define the coefficients from the model output
intercept #use the model coefficients to predict the value for mpg
intercept + disp*220 + hp*150 + drat*3
#[1] 18.57373
For a car with disp = 220, hp = 150, and drat = 3, the model predicts that the car would have a mpg of 18.57373.
You can find the complete R code used in this tutorial here.
Additional Resources
The following tutorials explain how to fit other types of regression models in R:
How to Perform Quadratic Regression in R
How to Perform Polynomial Regression in R
How to Perform Exponential Regression in R