Linear regression in R

Linear regression is an algorithm to create a statistical model that allows you to infer a relationship between a dependent variable (sometimes called a response variable) and one or more independent variables (also called explanatory variables) and their interactions. Linear regression is one of the most popular modeling techniques in machine learning and R is one of the most powerful programming languages for data science, so they’re a natural pairing.

Not every problem can be solved using linear regression because not every problem is a linear problem where the output is a linear variable. For a linear regression model to be appropriate, we need to believe that a linear relationship exists between the dependent variable and the independent variables. This means that a line can be fit between the two (or more) variables.

For example, we might want to model how much a baby should weigh by how old they are in days. The more days, (that is, the older the baby) the more kilograms (that is, the greater the weight). Not all problems have linear outputs and not all variables conform to the requirements of a linear regression.

In order to confirm that our relationship is in fact a linear one and that our data meet the requirements for a linear regression to work, we can perform some tests to confirm that we can run an accurate linear regression. Once we’ve confirmed that there are no violations of our assumptions, we can begin fitting a linear regression model to our data and evaluating how well the models fit that data and thus predict our dependent variable. To learn more about the underlying functionality of how Linear Regressions work you may want to read our tutorial on Linear Regression from Scratch.

In this tutorial, we'll use the Cars data set which is included with R and create a linear regression model that can provide estimates of the effect that independent variables have on a dependent variable. We'll also look at strategies to validate models and deal with interactions between variables and outliers in your data.

Prerequisites

Create an IBM Cloud account so that you can create a watsonx.ai project.

Steps

Step 1. Set up your environment

While you can choose from several tools, this tutorial walks you through how to set up an IBM account to use a Jupyter Notebook. Jupyter Notebooks are widely used within data science to combine code, text, images, and data visualizations to formulate a well-formed analysis.

1. Log in to watsonx.ai using your IBM Cloud account.

2. Create a watsonx.ai project.

3. Create a Jupyter Notebook using the R Kernel.

This step will open a Notebook environment where you can load your data set and copy the code from this tutorial to create linear regressions.

Step 2. Loading and evaluating data

To begin, we can load the Cars dataset that is included with R into a data frame. Once we have the Cars data loaded, we can look at the first few rows of the data:

cars_data <- mtcars
head(cars_data)

This gives us the first 5 rows of the data frame:

                    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

The meaning of each column is shown below in this table:

We’d like to explore the gas mileage of the different cars in the dataset so to begin, we’d want to check the distribution of values in the mpg. We want to do some exploratory data analysis to check whether there are extreme outliers or clusters of data or whether it’s evenly distributed:

plot(cars_data$mpg)

We can also explore the values using a histogram:

hist(cars_data$mpg)

Our data doesn’t quite follow a normal distribution, but it looks close enough. Once we begin fitting some exploratory models we will be able to see more. Let’s look at an initial question: are lighter cars more fuel efficient than heavier ones? First, we see if there’s a linear relationship between the two:

scatter.smooth(x=cars_data$mpg, y=cars_data$wt, main="Does Weight affect MPG?")  # scatterplot

There is a linear relationship, so we can fit a linear regression.

Step 3. Create a simple linear regression

We do this with a call to lm() that takes, at a minimum, two parameters: the formula and the data set to which to fit the formula.

A linear regression takes the following form: y = b0 + b1*x + e.

• y: the outcome variable
• b0: the intercept
• x: the predictor variable
• b1: the coefficient for the predictor variable
• e: the error term

A simple regression formula in R takes the following form, Y ~ X, where Y is the outcome variable and X is the predictor variable. Since we want to calculate the coefficient of wt (weight) with regard to mpg, we would use mpg ~ wt. This will calculate the regression and return a model variable that we can pass to summary() to see the model summary and information about the fitted values for each variable.

cars_lm <- lm("mpg ~ wt", data=cars_data)
summary(cars_lm)

This returns the following information about our model:

Call:
lm(formula = "mpg ~ wt", data = cars_data)

Residuals:
 Min      1Q  Median      3Q     Max
-4.5432 -2.3647 -0.1252  1.4096  6.8727

Coefficients:
         Estimate Std. Error t value             Pr(>|t|)
(Intercept)  37.2851     1.8776  19.858 < 0.0000000000000002 ***
wt           -5.3445     0.5591  -9.559       0.000000000129 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446
F-statistic: 91.38 on 1 and 30 DF,  p-value: 0.0000000001294

The summary shows us the regression coefficients for each independent variable, the p-value for that estimation, residual std estimated error. Below the summary for the variables is the r-squared, adjusted r-squared, degrees of freedom, and the f-statistic for the entire model.

We can see that our wt is statistically significant with an estimate of -5.3444716. This means that an increase of 1 in weight leads to a decrease of -5.3444716 in mpg. Before we can be too sure of the accuracy of this linear regression we will want to run model diagnostics. We want to determine whether we have an accurate model and if the data violates any of the linear regression assumptions. There are three types of assumption violations that we want to look for:

• Outliers or unusual observations that have too much influence on the analysis. These few observations might drastically change the outcome of the model and make it generalize less well.
• Model misspecification. We’re assuming that E[y] = Xβ but we may be wrong and we might need to include different independent variables in the model to get a satisfactory result.
• Error distribution. We have assumed that the errors will be normally distributed but autocorrelation, heteroscedasticity, and non-normality might be present. We assume that our variables are homoscedastic, are not correlated with one another, and that the errors roughly follow a normal distribution.

Many times problems with one of these violations can be addressed by transforming either the dependent or independent variables or removing outliers or using weighted values.

We can use the ggfortify library, which is built on ggplot2, to create nice looking visualizations of the residual plots of the difference between our actual data and our predicted values to confirm how well our model fits our data.

Outliers and high leverage points

An outlier is a point that has an extreme outcome variable value. The presence of outliers may affect the interpretation of the model, because it increases the residual standard error (RSE). Outliers can be identified by examining the standardized residual (or studentized residual), which is the residual divided by its estimated standard error. Standardized residuals can be interpreted as the number of standard errors away from the regression line. Observations whose standardized residuals are greater than 3 in absolute value are possible outliers.

A data point has high leverage if it has extreme predictor x values. This can be detected by examining the leverage statistic or the hat-value. A value of this statistic above 2(p + 1)/n indicates an observation with high leverage; where, p is the number of predictors and n is the number of observations.

Outliers and high leverage points can be identified by inspecting the Residuals vs Leverage plot:

library(ggfortify)
autoplot(cars_lm, label = FALSE)[1]

In our example, the residuals do roughly track our fitted values but they are generally fairly high and higher at the minimum and maximum fitted values. This can be a sign that our model isn't capturing the variance in our data especially well.

QQ Plot

The QQ plot of residuals can be used to visually check the normality assumption. The normal probability plot of residuals should approximately follow a straight line.

autoplot(cars_lm, label = FALSE)[2]

Homogeneity of variance

This assumption can be checked by examining the scale-location plot, also known as the spread-location plot.

autoplot(cars_lm, label = FALSE)[3]

Influential values

An influential value is a value which inclusion or exclusion of the value can alter the results of the regression analysis. Such a value is associated with a large residual. Not all outliers (or extreme data points) are influential in linear regression analysis.

autoplot(cars_lm, label = FALSE)[4]

Step 4. Create multiple regression

We should suspect that there are more factors at play in the mpg of each car than just the weight based on what we saw in the diagnostic plots. Multiple regression allows us to fit multiple indepedent variables to the dependent variable. We can specify which terms we want to include with the + symbol: y ~ x1 + x2.

We can also include all variables in a data set using the . operator: y ~ ..

This will include every independent variable in our multiple regression. A handy time saving technique is to find the best model using step-wise regression which chooses the best model by AIC. We will use the fit_all model and then go backwards to the best fit model.

cars_all<-lm("mpg ~ .", data=cars_data)
summary(cars_all)

This returns the following information about our model:

Call:
lm(formula = "mpg ~ .", data = cars_data)

Residuals:
 Min      1Q  Median      3Q     Max
-3.4506 -1.6044 -0.1196  1.2193  4.6271

Coefficients:
         Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.30337   18.71788   0.657   0.5181
cyl         -0.11144    1.04502  -0.107   0.9161
disp         0.01334    0.01786   0.747   0.4635
hp          -0.02148    0.02177  -0.987   0.3350
drat         0.78711    1.63537   0.481   0.6353
wt          -3.71530    1.89441  -1.961   0.0633 .
qsec         0.82104    0.73084   1.123   0.2739
vs           0.31776    2.10451   0.151   0.8814
am           2.52023    2.05665   1.225   0.2340
gear         0.65541    1.49326   0.439   0.6652
carb        -0.19942    0.82875  -0.241   0.8122
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.65 on 21 degrees of freedom
Multiple R-squared:  0.869,  Adjusted R-squared:  0.8066
F-statistic: 13.93 on 10 and 21 DF,  p-value: 0.0000003793

This model has a high r-squared value but many of the independent variables have high standard error rates. This might be because some of the variables are masking the effect of others. Now that we have a model with all the independent variables, we can eliminate unnecessary variables one by one using the step() function. This returns the model with the highest Akaike Information Criterion (AIC) score:

cars_best <- step(cars_all, direction = "backward", trace = 0)
summary(cars_best)

This returns the following information about our model:

Call:
lm(formula = mpg ~ wt + qsec + am, data = cars_data)

Residuals:
    Min      1Q  Median      3Q     Max
-3.4811 -1.5555 -0.7257  1.4110  4.6610

Coefficients:
            Estimate Std. Error t value   Pr(>|t|)
(Intercept)   9.6178     6.9596   1.382   0.177915
wt           -3.9165     0.7112  -5.507 0.00000695 ***
qsec          1.2259     0.2887   4.247   0.000216 ***
am            2.9358     1.4109   2.081   0.046716 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.459 on 28 degrees of freedom
Multiple R-squared:  0.8497, Adjusted R-squared:  0.8336
F-statistic: 52.75 on 3 and 28 DF,  p-value: 0.0000000000121

The best fit model by AIC has just 3 predictors: mpg ~ wt + qsec + am. The p-value of the coefficients are all significant at the 0.05 level, the residual standard error is 2.459, and the adjusted r-squared value is 0.8497. The model captures most of the variablility of the data and seems to have predictive power.

Step 5. Use interaction terms

There is one feature of this model that can be improved: clarifying the relationship between the transmission types encoded in the am variable. We know that the transmission type affects the mpg but is there a difference in how the other variables affect mpg that's different for each transmission type? One approach would be to use either a t-test or anova analysis to check, but we could probably get away with a simple visual comparison. Let’s graph a boxplot to see how different the values for MPG are for the different transmission types:

ggplot(cars_data, aes(group=am, y=mpg, fill=as.factor(am))) + geom_boxplot()

It seems the type of transmission does affect mpg, so it might be worth measuring how weight affects mpg depending on the transmission type. We do this by using an interaction term.

Interactions are formed by the product of any two variables. Interactions can be between numeric variables or categorial variables. In our example the formula is essentiall: y = b + (b1 * X) + (b2 * W) + (b3 * X * W). We can interpret each coefficient like so:

• b is the intercept, or the predicted outcome when W = 0 and X = 0.
• b1 is the simple effect or slope of X, for a one unit change in the predicted change in Y at W = 0
• b2 is the simple effect or slope of W, for a one unit change in W the predicted change in Y at X = 0
• b3 is the interaction of X and W, the change in the slope of X for a one unit increase in W

In R, we indicate the interaction between two independent variables with the * symbol, such as y ~ w \* x\. For our data set, we will use the variables that we arrived at using step-wise elimination but add the interactions between transmission type and weight: mpg ~ qsec + (wt * as.factor(am)).

Below we'll run the lm() function to create our model with interactions and then print the summary of the models parameters.

cars_interactions <- lm("mpg ~ qsec + (wt * as.factor(am))", data=cars_data)
summary(cars_interactions)

This returns the following information about our model:

Call:
lm(formula = "mpg ~ qsec + (wt * as.factor(am))", data = cars_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5076 -1.3801 -0.5588  1.0630  4.3684 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)          9.723      5.899   1.648 0.110893    
qsec                 1.017      0.252   4.035 0.000403 ***
wt                  -2.937      0.666  -4.409 0.000149 ***
as.factor(am)1      14.079      3.435   4.099 0.000341 ***
wt:as.factor(am)1   -4.141      1.197  -3.460 0.001809 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.084 on 27 degrees of freedom
Multiple R-squared:  0.8959, Adjusted R-squared:  0.8804 
F-statistic: 58.06 on 4 and 27 DF,  p-value: 0.0000000000007168

The Adjusted R-squared for this model with interactions is 0.8804, whereas our model without the interactions had an Adjusted R-Squared of 0.8336. This is a good indication that our interaction model may be a better fit for our data. Additionally, our residual standard error is 2.084 for the model with interactions, whereas the model without had a residual standard error of 2.459.

Another way to compare the models is to look at the slopes of the regression lines for the two different transmission types using the emmip method from the emmeans library to graph the predictions of the model at different weights for each our transmission types.

install.packages("emmeans")
library(emmeans)

(mylist <- list(wt=seq(1.5,5.5,by=0.25),am=unique(cars_data$am)))
emmip(cars_interactions, am~wt, at=mylist, CIs=TRUE)

We can see that the slopes of the regression line for MPG for manual transmission cars is different than for automatic cars with regard to the weight. Thus, it makes sense for us to use an interaction term that captures these two. One thing we might want to account for is that our model is only seeing heavier cars which are manual transmissions and thus might be biased in its estimates of the MPG of those cars.

Step 6. Create a robust linear regression

We can attempt to correct for these outliers by using robust regression. Robust regression is a commonly used alternative to least squares regression when data has outliers or influential observations. It can also be used for the purpose of detecting influential observations or simply for correcting for outliers without needing to remove them from the data.

The rlm() command in the MASS package command implements several versions of robust regression. The formula for the linear regression is passed in the same way.

require(MASS)

weighted_model <- rlm(mpg ~ qsec + (wt * as.factor(am)), data=cars_data)
summary(weighted_model)

## 
## Call: rlm(formula = mpg ~ qsec + (wt * as.factor(am)), data = cars_data)
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1900 -1.1106 -0.4092  1.2357  4.7696 
## 
## Coefficients:
##                   Value   Std. Error t value
## (Intercept)       11.7019  5.7018     2.0523
## qsec               0.9265  0.2436     3.8034
## wt                -3.0813  0.6438    -4.7864
## as.factor(am)1    13.3573  3.3204     4.0228
## wt:as.factor(am)1 -3.9493  1.1568    -3.4140
## 
## Residual standard error: 1.696 on 27 degrees of freedom

We can see that rlm does not return r-squared or p-values for the variables. Robust regression uses iteratively reweighted least squares (IRLS) for maximum likelihood estimation (MLE), whereas linear regression uses ordinary lLeast squares (OLS) so the calculations are quite different. However, looking at the residual standard error, we can see that the robust regression has a smaller standard error overall as well as smaller standard errors for each of the independent variables as well. This should give us confidence that our robust model does in fact fit our data better.

Summary

Linear regression is one of the oldest and most commonly used techniques for model building in data science. Once you’ve confirmed that your data contains linear relationships, you can build linear regression models in R using the lm command.

Plotting the regression diagnostics will help observe how well your model fits your data. In the case of multiple regressions, you may find your model needs to account for interactions between variables. These can be included in your formula and evaluated like any other variable.

Finally, if you suspect that outliers are influencing your estimations and model fit, you can use techniques like robust linear regression to reweight those outliers and correct how much influence they have in your model.

Try watsonx for free

Build an artificial intelligence (AI) strategy for your business on one collaborative AI and data platform called IBM watsonx, which combines new generative AI capabilities powered by foundation models and traditional machine learning into a powerful platform spanning the AI lifecycle. With watsonx.ai, you can train, validate, tune, and deploy models with ease and build AI applications in a fraction of the time with a fraction of the data.

Try watsonx.ai, the next-generation studio for AI builders.

Next steps

Explore more articles and tutorials about watsonx on IBM Developer.

Also, to learn more about other supervised learning algorithms that you can apply to classification and regression problems, see this tutorial in the Getting started with machine learning learning path: Tutorial: Learn regression algorithms using Python and scikit-learn