Principle

The principle of simple linear regression is to find the line (i.e., determine its equation) which passes as close as possible to the observations, that is, the set of points formed by the pairs \((x_i, y_i)\).

In the first step, there are many potential lines. Three of them are plotted:

To find the line which passes as close as possible to all the points, we take the square of the vertical distance between each point and each potential line. Note that we take the square of the distances to make sure that a negative gap (i.e., a point below the line) is not compensated by a positive gap (i.e., a point above the line). The line which passes closest to the set of points is the one which minimizes the sum of these squared distances.

The resulting regression line is presented in blue in the following plot, and the dashed gray lines represent the vertical distance between the points and the fitted line. These vertical distances between each observed point and the fitted line determined by the least squares method are called the residuals of the linear regression model and denoted \(\epsilon\).

By definition, there is no other line with a smaller total distance between the points and the line. This method is called the least squares method, or OLS for ordinary least squares.

Equation

The regression model can be written in the form of the equation:

\[Y = \beta_0 + \beta_1 X + \epsilon\]

with:

• \(Y\) the dependent variable
• \(X\) the independent variable
• \(\beta_0\) the intercept (the mean value of \(Y\) when \(x = 0\)), also sometimes denoted \(\alpha\)
• \(\beta_1\) the slope (the expected increase in \(Y\) when \(X\) increases by one unit)
• \(\epsilon\) the residuals (the error term of mean 0 which describes the variations of \(Y\) not captured by the model, also referred as the noise)

When we determine the line which passes closest to all the points (we say that we fit a line to the observed data), we actually estimate the unknown parameters \(\beta_0\) and \(\beta_1\) based on the data at hand. Remember from your geometry classes, to draw a line you only need two parameters—the intercept and the slope.

These estimates (and thus the blue line shown in the previous scatterplot) can be computed by hand with the following formulas:

\[\begin{align} \widehat\beta_1 &= \frac{\sum^n_{i = 1} (x_i - \bar{x})(y_i - \bar{y})}{\sum^n_{i = 1}(x_i - \bar{x})^2} \\ &= \frac{\left(\sum^n_{i = 1}x_iy_i\right) - n\bar{x}\bar{y}}{\sum^n_{i = 1}(x_i - \bar{x})^2} \end{align}\]

and

\[\widehat\beta_0 = \bar{y} - \widehat\beta_1 \bar{x}\]

with \(\bar{x}\) and \(\bar{y}\) denoting the sample mean of \(x\) and \(y\), respectively.

(If you struggle to compute \(\widehat\beta_0\) and \(\widehat\beta_1\) by hand, see this Shiny app which helps you to easily find these estimates based on your data.)

Interpretations of coefficients \(\widehat\beta\)

The intercept \(\widehat\beta_0\) is the mean value of the dependent variable \(Y\) when the independent variable \(X\) takes the value 0. Its estimation has no interest in evaluating whether there is a linear relationship between two variables. It has, however, an interest if you want to know what the mean value of \(Y\) could be when \(x = 0\).⁴

The slope \(\widehat\beta_1\), on the other hand, corresponds to the expected variation of \(Y\) when \(X\) varies by one unit. It tells us two important informations:

1. The sign of the slope indicates the direction of the line—a positive slope (\(\widehat\beta_1 > 0\)) indicates that there is a positive relationship between the two variables of interest (they vary in the same direction), whereas a negative slope (\(\widehat\beta_1 < 0\)) means that there is a negative relationship between the two variables (they vary in opposite directions).
2. The value of the slope provides information on the speed of evolution of the variable \(Y\) as a function of the variable \(X\). The larger the slope in absolute value, the larger the expected variation of \(Y\) for each unit of \(X\). Note, however, that a large value does not necessarily mean that the relationship is statistically significant (more on that in the section about significance of the relationship).

This is similar to the correlation coefficient, which gives information about the direction and the strength of the relationship between two variables.

To perform a linear regression in R, we use the lm() function (which stands for linear model). The function requires to set the dependent variable first then the independent variable, separated by a tilde (~).

Applied to our example of weight and car’s consumption, we have:

model <- lm(mpg ~ wt, data = dat)

The summary() function gives the results of the model:

summary(model)

## 
## Call:
## lm(formula = mpg ~ wt, data = dat)
## 
## 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  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## 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: 1.294e-10

In practice, we usually check the conditions of application before interpreting the coefficients (because if they are not respected, results may be biased).

In this article, however, I present the interpretations before testing the conditions because the point is to show how to interpret the results, and less about finding a valid model.

The results can be summarized as follows (see the column Estimate in the table Coefficients):

• The intercept \(\widehat\beta_0 =\) 37.29 indicates that, for a hypothetical car weighting 0 lbs, we can expect, on average, a consumption of 37.29 miles/gallon. This interpretation is shown for illustrative purposes, but as a car weighting 0 lbs is impossible, the interpretation has no meaning. In practice, we would therefore refrain from interpreting the intercept in this case. See another interpretation of the intercept when the independent variable is centered around its mean in this section.
• The slope \(\widehat\beta_1 =\) -5.34 indicates that:
  • There is a negative relationship between the weight and the distance a car can drive with a gallon (this was expected given the negative trend of the points in the scatterplot shown previously).
  • But more importantly, a slope of -5.34 means that, for an increase of one unit in the weight (that is, an increase of 1000 lbs), the number of miles per gallon decreases, on average, by 5.34 units. In other words, for an increase of 1000 lbs, the number of miles/gallon decreases, on average, by 5.34.

dat_centered <- dat

dat_centered$wt_centered <- dat$wt - mean(dat$wt)

mod_centered <- lm(mpg ~ wt_centered,
  data = dat_centered
)

summary(mod_centered)

## 
## Call:
## lm(formula = mpg ~ wt_centered, data = dat_centered)
## 
## 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)  20.0906     0.5384  37.313  < 2e-16 ***
## wt_centered  -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## 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: 1.294e-10

Significance of the relationship

As mentioned earlier, the value of the slope does not, by itself, make it possible to assess the significance of the linear relationship.

In other words, a slope different from 0 does not necessarily mean it is significantly different from 0, so it does not mean that there is a significant relationship between the two variables in the population. There could be a slope of 10 that is not significant, and a slope of 2 that is significant.

Significance of the relationship also depends on the variability of the slope, which is measured by its standard error and generally noted \(se(\widehat\beta_1)\).

Without going too much into details, to assess the significance of the linear relationship, we divide the slope by its standard error. This ratio is the test statistic and follows a Student distribution with \(n - 2\) degrees of freedom:⁵

\[T_{n - 2} = \frac{\widehat\beta_1}{se(\widehat\beta_1)}\]

For a bilateral test, the null and alternative hypotheses are:⁶

• \(H_0 : \beta_1 = 0\) (there is no (linear) relationship between the two variables)
• \(H_1 : \beta_1 \ne 0\) (there is a (linear) relationship between the two variables)

Roughly speaking, if this ratio is greater than 2 in absolute value then the slope is significantly different from 0, and therefore the relationship between the two variables is significant (and in that case it is positive or negative depending on the sign of the estimate \(\widehat\beta_1\)).

The standard error and the test statistic are shown in the column Std. Error and t value in the table Coefficients.

Fortunately, R gives a more precise and easier way to assess to the significance of the relationship. The information is provided in the column Pr(>|t|) of the Coefficients table. This is the p-value of the test. As for any statistical test, if the p-value is greater than or equal to the significance level (usually \(\alpha = 0.05\)), we do not reject the null hypothesis, and if the p-value is lower than the significance level, we reject the null hypothesis.

If we do not reject the null hypothesis, we do not reject the hypothesis of no relationship between the two variables (because we do not reject the hypothesis of a slope of 0). On the contrary, if we reject the null hypothesis of no relationship, we can conclude that there is a significant linear relationship between the two variables.

In our example, the p-value = 1.29e-10 < 0.05 so we reject the null hypothesis at the significance level \(\alpha = 5\%\). We therefore conclude that there is a significant relationship between a car’s weight and its fuel consumption.

Tip: In order to make sure I interpret only parameters that are significant, I tend to first check the significance of the parameters thanks to the p-values, and then interpret the estimates accordingly. For completeness, note that the test is also performed on the intercept. The p-value being smaller than 0.05, we also conclude that the intercept is significantly different from 0.

Conditions of application

Unfortunately, linear regression cannot be used in all situations.

In addition to the requirement that the dependent variable must be a continuous quantitative variables, simple linear regression requires that the data satisfy the following conditions:

1. Linearity: The relationship between the two variables should be linear (at least roughly). For this reason it is always necessary to represent graphically the data with a scatterplot before performing a simple linear regression.⁷

2. Independence: Observations must be independent. It is the sampling plan and the experimental design that usually provide information on this condition. If the data come from different individuals or experimental units, they are usually independent. On the other hand, if the same individuals are measured at different periods, the data are probably not independent.
3. Normality of the residuals: For large sample sizes, confidence intervals and tests on the coefficients are (approximately) valid whether the error follows a normal distribution or not (a consequence of the central limit theorem, see more in Ernst and Albers (2017) and Lumley et al. (2002))! For small sample sizes, residuals should follow a normal distribution. This condition can be tested visually (via a QQ-plot and/or a histogram), or more formally (via the Shapiro-Wilk test for instance).
4. Homoscedasticity of the residuals: The variance of the errors should be constant. There is a lack of homoscedasticity when the dispersion of the residuals increases with the predicted values (fitted values). This condition can be tested visually (by plotting the standardized residuals vs. the fitted values) or more formally (via the Breusch-Pagan test).
5. No influential points: If the data contain outliers, it is essential to identify them so that they do not, on their own, influence the results of the regression. Note that an outlier is not an issue per se if the point is in the alignment of the regression line for example because it does not influence the regression line. It becomes a problem in the context of linear regression if it influences in a substantial manner the estimates (and in particular the slope of the regression line). This can be tackled by identifying outliers (via the Cook’s distance⁸ or the leverage index⁹ for instance), and comparing the results with and without the potential outliers. Do the results remain the same with the two approaches? If yes, outliers are not really an issue in this case. If results are much different, you can use the Theil-Sen estimator, robust regression or quantile regression which are all more robust to outliers.

Tip: I remember the first 4 conditions thanks to the acronym “LINE”, for Linearity, Independence, Normality and Equality of variance.

If any of the condition is not met, the tests and the conclusions could be erroneous so it is best to avoid using and interpreting the model. If this is the case, sometimes the conditions can be met by transforming the data (e.g., logarithmic transformation, square or square root, Box-Cox transformation, etc.) or by adding a quadratic term to the model.

If it does not help, it could be worth thinking about removing some variables or adding other variables, or even considering other types of models such as non-linear models.

Keep in mind that in practice, conditions of application should be verified before drawing any conclusion based on the model. I refrain here from testing the conditions on our data because it will be covered in details in the context of multiple linear regression (see this section).

Visualizations

If you are a frequent reader of the blog, you may know that I like to draw (simple but efficient) visualizations to illustrate my statistical analyses. Linear regression is not an exception.

There are numerous ways to visualize the relationship between the two variables of interest, but the easiest one I found so far is via the visreg() function from the package of the same name:

library(visreg)
visreg(model)

I like this approach for its simplicity—only a single line of code.

However, other elements could be displayed on the regression plot (for example the regression equation and the \(R^2\)). This can easily be done with the stat_regline_equation() and stat_cor() functions from the {ggpubr} package:

# load necessary libraries
library(ggpubr)

# create plot with regression line, regression equation and R^2
ggplot(dat, aes(x = wt, y = mpg)) +
  geom_smooth(method = "lm") +
  geom_point() +
  stat_regline_equation(label.x = 3, label.y = 32) + # for regression equation
  stat_cor(aes(label = after_stat(rr.label)), label.x = 3, label.y = 30) + # for R^2
  theme_minimal()

Principle

Multiple linear regression is a generalization of simple linear regression, in the sense that this approach makes it possible to relate one variable with several variables through a linear function in its parameters.

Multiple linear regression is used to assess the relationship between two variables while taking into account the effect of other variables. By taking into account the effect of other variables, we cancel out the effect of these other variables in order to isolate and measure the relationship between the two variables of interest. This point is the main difference with simple linear regression.

To illustrate how to perform a multiple linear regression in R, we use the same dataset than the one used for simple linear regression (mtcars). Below a short preview:

head(dat)

##                    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

We have seen that there is a significant and negative linear relationship between the distance a car can drive with a gallon and its weight (\(\widehat\beta_1 =\) -5.34, \(p\)-value < 0.001).

However, one may wonder whether there are not in reality other factors that could explain a car’s fuel consumption.

To explore this, we can visualize the relationship between a car’s fuel consumption (mpg) together with its weight (wt), horsepower (hp) and displacement (disp) (engine displacement is the combined swept (or displaced) volume of air resulting from the up-and-down movement of pistons in the cylinders, usually the higher the more powerful the car):

ggplot(dat) +
  aes(x = wt, y = mpg, colour = hp, size = disp) +
  geom_point() +
  scale_color_gradient() +
  labs(
    y = "Miles per gallon",
    x = "Weight (1000 lbs)",
    color = "Horsepower",
    size = "Displacement"
  ) +
  theme_minimal()

It seems that, in addition to the negative relationship between miles per gallon and weight, there is also:

• a negative relationship between miles/gallon and horsepower (lighter points, indicating more horsepower, tend to be more present in low levels of miles per gallon)
• a negative relationship between miles/gallon and displacement (bigger points, indicating larger values of displacement, tend to be more present in low levels of miles per gallon).

Therefore, we would like to evaluate the relation between the fuel consumption and the weight, but this time by adding information on the horsepower and displacement. By adding this additional information, we are able to capture only the direct relationship between miles/gallon and weight (the indirect effect due to horsepower and displacement is canceled out).

This is the whole point of multiple linear regression! In fact, in multiple linear regression, the estimated relationship between the dependent variable and an explanatory variable is an adjusted relationship, that is, free of the linear effects of the other explanatory variables.

Let’s illustrate this notion of adjustment by adding both horsepower and displacement in our linear regression model:

model2 <- lm(mpg ~ wt + hp + disp,
  data = dat
)

summary(model2)

## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268,	Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

We can see that now, the relationship between miles/gallon and weight is weaker in terms of slope (\(\widehat\beta_1 =\) -3.8 now, against \(\widehat\beta_1 =\) -5.34 when only the weight was considered).

The effect of weight on fuel consumption was adjusted according to the effect of horsepower and displacement. This is the remaining effect between miles/gallon and weight after the effects of horsepower and displacement have been taken into account. More detailed interpretations in this section.

Equation

Multiple linear regression models are defined by the equation

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \epsilon\]

It is similar than the equation of simple linear regression, except that there is more than one independent variables (\(X_1, X_2, \dots, X_p\)).

Estimation of the parameters \(\beta_0, \dots, \beta_p\) by the method of least squares is based on the same principle as that of simple linear regression, but applied to \(p\) dimensions. It is thus no longer a question of finding the best line (the one which passes closest to the pairs of points (\(y_i, x_i\))), but finding the \(p\)-dimensional plane which passes closest to the coordinate points (\(y_i, x_{i1}, \dots, x_{ip}\)).

This is done by minimizing the sum of the squares of the deviations of the points on the plane:

Source: Thaddeussegura

Interpretations of coefficients \(\widehat\beta\)

The least squares method results in an adjusted estimate of the coefficients. The term adjusted means after taking into account the linear effects of the other independent variables on the dependent variable, but also on the predictor variable.

In other words, the coefficient \(\beta_1\) corresponds to the slope of the relationship between \(Y\) and \(X_1\) when the linear effects of the other explanatory variables (\(X_2, \dots, X_p\)) have been removed, both at the level of the dependent variable \(Y\) but also at the level of \(X_1\).

Applied to our model with weight, horsepower and displacement as independent variables, we have:

summary(model2)

## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268,	Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

The table Coefficients gives the estimate for each parameter (column Estimate), together with the \(p\)-value of the nullity of the parameter (column Pr(>|t|)).

The hypotheses are the same as for simple linear regression, that is:

• \(H_0 : \beta_j = 0\)
• \(H_1 : \beta_j \ne 0\)

The test of \(\beta_j = 0\) is equivalent to testing the hypothesis: is the dependent variable associated with the independent variable studied, all other things being equal, that is to say, at constant level of the other independent variables.

In other words:

• the test of \(\beta_1 = 0\) corresponds to testing the hypothesis: is fuel consumption associated with a car’s weight, at a constant level of horsepower and displacement
• the test of \(\beta_2 = 0\) corresponds to testing the hypothesis: is fuel consumption associated with horsepower, at a constant level of weight and displacement
• the test of \(\beta_3 = 0\) corresponds to testing the hypothesis: is fuel consumption associated with displacement, at a constant level of weight and displacement
• (for the sake of completeness: the test of \(\beta_0 = 0\) corresponds to testing the hypothesis: is miles/gallon different from 0 when weight, horsepower and displacement are equal to 0)

In practice, we usually check the conditions of application before interpreting the coefficients (because if they are not respected, results may be biased). In this article, however, I present the interpretations before testing the conditions because the point is to show how to interpret the results, and less about finding a valid model.

Based on the output of our model, we conclude that:

• There is a significant and negative relationship between miles/gallon and weight, all else being equal. So for an increase of one unit in the weight (that is, an increase of 1000 lbs), the number of miles/gallon decreases, on average, by 3.8, for a constant level of horsepower and displacement (\(p\)-value = 0.001).
• There is a significant and negative relationship between miles/gallon and horsepower, all else being equal. So for an increase of one unit of horsepower, the distance traveled with a gallon decreases, on average, by 0.03 mile, for a constant level of weight and displacement (\(p\)-value = 0.011).
• We do not reject the hypothesis of no relationship between miles/gallon and displacement when weight and horsepower stay constant (because \(p\)-value = 0.929 > 0.05).
• (For completeness but it should be interpreted only when it makes sense: for a weight, horsepower and displacement = 0, we can expect that a car has, on average, a fuel consumption of 37.11 miles/gallon (\(p\)-value < 0.001). See a more useful interpretation of the intercept when the independent variables are centered in this section.)

This is how to interpret quantitative independent variables. Interpreting qualitative independent variables is slightly different in the sense that it quantifies the effect of a level in comparison with the reference level, sill all else being equal.

So it compares the different groups (formed by the different levels of the categorical variable) in terms of the dependent variable (this is why linear regression can be seen as an extension to the t-test and ANOVA).

For the illustration, we model the fuel consumption (mpg) on the weight (wt) and the shape of the engine (vs). The variable vs has two levels: V-shaped (the reference level) and straight engine.¹⁰

## Recoding dat$vs
library(forcats)
dat$vs <- as.character(dat$vs)
dat$vs <- fct_recode(dat$vs,
  "V-shaped" = "0",
  "Straight" = "1"
)

model3 <- lm(mpg ~ wt + vs,
  data = dat
)

summary(model3)

## 
## Call:
## lm(formula = mpg ~ wt + vs, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7071 -2.4415 -0.3129  1.4319  6.0156 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  33.0042     2.3554  14.012 1.92e-14 ***
## wt           -4.4428     0.6134  -7.243 5.63e-08 ***
## vsStraight    3.1544     1.1907   2.649   0.0129 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.78 on 29 degrees of freedom
## Multiple R-squared:  0.801,	Adjusted R-squared:  0.7873 
## F-statistic: 58.36 on 2 and 29 DF,  p-value: 6.818e-11

Based on the output of our model, we conclude that:

• For a V-shaped engine and for an increase of one unit in the weight (that is, an increase of 1000 lbs), the number of miles/gallon decreases, on average, by 4.44 (\(p\)-value < 0.001).
• The distance traveled with a gallon of fuel increases by, on average, 3.15 miles when the engine is straight compared to a V-shaped engine, for a constant weight (\(p\)-value = 0.013).
• (For completeness but it should be interpreted only when it makes sense: for a weight = 0 and a V-shaped engine, we can expect that the car has, on average, a fuel consumption of 33 miles/gallon (\(p\)-value < 0.001). See a more useful interpretation of the intercept when the independent variables are centered in this section.)

Conditions of application

As for simple linear regression, multiple linear regression requires some conditions of application for the model to be usable and the results to be interpretable. Conditions for simple linear regression also apply to multiple linear regression, that is:

1. Linearity of the relationships between the dependent and independent variables¹¹
2. Independence of the observations
3. Normality of the residuals
4. Homoscedasticity of the residuals
5. No influential points (outliers)

But there is one more condition for multiple linear regression:

6. No multicollinearity: Multicollinearity arises when there is a strong linear correlation between the independent variables, conditional on the other variables in the model. It is important to check it because it may lead to an imprecision or an instability of the estimated parameters when a variable changes. It can be assessed by studying the correlation between each pair of independent variables, or even better, by computing the variance inflation factor (VIF). The VIF measures how much the variance of an estimated regression coefficient increases, relative to a situation in which the explanatory variables are strictly independent. A high value of VIF is a sign of multicollinearity (the threshold is generally admitted at 5 or 10 depending on the domain). The easiest way to reduce the VIF is to remove some correlated independent variables, or eventually to standardize the data.

You will often see that these conditions are verified by running plot(model, which = 1:6) and it is totally correct. However, I recently discovered the check_model() function from the {performance} package which tests these conditions all at the same time (and let’s be honest, in a more elegant way).¹²

Applied on our model2 with miles/gallon as dependent variable, and weight, horsepower and displacement as independent variables, we have:

# install.packages("performance")
# install.packages("see")
library(performance)

check_model(model2)

In addition to testing all conditions at the same time, it also gives insight on how to interpret the different diagnostic plots and what you should expect (see in the subtitles of each plot).

Based on these diagnostic plots, we see that:

• Homogeneity of variance (middle left plot) is respected
• Multicollinearity (bottom left plot) is not an issue (I tend to use the threshold of 10 for VIF, and all of them are below 10)¹³
• There is no influential points (middle right plot)
• Normality of the residuals (bottom right plot) is also not perfect due to 3 points deviating from the reference line but it still seems acceptable to me. In any case, the number of observations is large enough given the number of parameters¹⁴ and given the small deviation from normality so tests on the coefficients are (approximately) valid whether the error follows a normal distribution or not
• Linearity (top right plot) is not perfect so let’s check each independent variable separately:

# weight
ggplot(dat, aes(x = wt, y = mpg)) +
  geom_point() +
  theme_minimal()

# horsepower
ggplot(dat, aes(x = hp, y = mpg)) +
  geom_point() +
  theme_minimal()

# displacement
ggplot(dat, aes(x = disp, y = mpg)) +
  geom_point() +
  theme_minimal()

It seems that the relationship between miles/gallon and horsepower is not linear, which could be the main component of the slight linearity defect of the model.

To improve linearity, the variable could be removed, a transformation could be applied (logarithmic and/or squared for instance) or a quadratic term could be added to the model.¹⁵ If this does not fix the issue of linearity, other types of models could be considered.

If you want to read more about these conditions of applications and how to deal with them, here is a very complete chapter on diagnostics for linear models written by Prof. Dustin Fife.

For the sake of easiness and for illustrative purposes, we assume linearity for the rest of the article.

When the conditions of application are met, we usually say that the model is valid. But not all valid models are good models. The next section deals with model selection.

How to choose a good linear model?

A model which satisfies the conditions of application is the minimum requirement, but you will likely find several models that meet this criteria. So one may wonder how to choose between different models that are all valid?

The three most common tools to select a good linear model are according to:

1. the \(p\)-value associated to the model,
2. the coefficient of determination \(R^2\) and
3. the Akaike Information Criterion

The approaches are detailed in the next sections. Note that the first two are applicable to simple and multiple linear regression, whereas the third is only applicable to multiple linear regression.

summary(model2)

## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268,	Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

summary(model2)

## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268,	Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

## vs has already been transformed into factor
## so only am is transformed here

## Recoding dat$vs
library(forcats)
dat$am <- as.character(dat$am)
dat$am <- fct_recode(dat$am,
  "Automatic" = "0",
  "Manual" = "1"
)

model4 <- lm(mpg ~ .,
  data = dat
)

model4 <- step(model4, trace = FALSE)

summary(model4)

## 
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = dat)
## 
## 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 6.95e-06 ***
## qsec          1.2259     0.2887   4.247 0.000216 ***
## amManual      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: 1.21e-11

Visualizations

There are many ways to visualize results of a linear regression. The easiest ones I am aware of are:

1. visreg() illustrates the relationships between the dependent and independent variables in different plots (one for each independent variable unless you specify which relationship you want to illustrate):

library(visreg)

visreg(model4)

2. ggcoefstats() illustrates the results in one single plot, with many statistical details:

library(ggstatsplot)

ggcoefstats(model4)

In this plot:

• when the solid line does not cross the vertical dashed line, the estimates is significantly different from 0 at the 5% significance level (i.e., \(p\)-value < 0.05)
• furthermore, a point to the right (left) of the vertical dashed line means that there is a positive (negative) relationship between the two variables
• the more extreme the point, the stronger the relationship

3. plot_summs() which also illustrates the results but in a more concise way:

library(jtools)
library(ggstance)

plot_summs(model4,
  omit.coefs = NULL
)

The advantage of this approach is that it is possible to compare coefficients of multiple models simultaneously (particularly interesting when the models are nested):

model4bis <- lm(mpg ~ wt + qsec + am + hp,
  data = dat
)

plot_summs(model4,
  model4bis,
  omit.coefs = NULL
)

To go further

Below some more advanced topics related to linear regression. Feel free to comment at the end of the article if you believe I missed an important one.

library(parameters)

model_parameters(model4, summary = TRUE)

## Parameter   | Coefficient |   SE |         95% CI | t(28) |      p
## ------------------------------------------------------------------
## (Intercept) |        9.62 | 6.96 | [-4.64, 23.87] |  1.38 | 0.178 
## wt          |       -3.92 | 0.71 | [-5.37, -2.46] | -5.51 | < .001
## qsec        |        1.23 | 0.29 | [ 0.63,  1.82] |  4.25 | < .001
## am [Manual] |        2.94 | 1.41 | [ 0.05,  5.83] |  2.08 | 0.047 
## 
## Model: mpg ~ wt + qsec + am (32 Observations)
## Residual standard deviation: 2.459 (df = 28)
## R2: 0.850; adjusted R2: 0.834

library(gt)

print_html(model_parameters(model4, summary = TRUE))

library(report)

report(model4)[1]

## [1] "We fitted a linear model (estimated using OLS) to predict mpg with wt, qsec and am (formula: mpg ~ wt + qsec + am). The model explains a statistically significant and substantial proportion of variance (R2 = 0.85, F(3, 28) = 52.75, p < .001, adj. R2 = 0.83). The model's intercept, corresponding to wt = 0, qsec = 0 and am = Automatic, is at 9.62 (95% CI [-4.64, 23.87], t(28) = 1.38, p = 0.178). Within this model:\n\n  - The effect of wt is statistically significant and negative (beta = -3.92, 95% CI [-5.37, -2.46], t(28) = -5.51, p < .001; Std. beta = -0.64, 95% CI [-0.87, -0.40])\n  - The effect of qsec is statistically significant and positive (beta = 1.23, 95% CI [0.63, 1.82], t(28) = 4.25, p < .001; Std. beta = 0.36, 95% CI [0.19, 0.54])\n  - The effect of am [Manual] is statistically significant and positive (beta = 2.94, 95% CI [0.05, 5.83], t(28) = 2.08, p = 0.047; Std. beta = 0.49, 95% CI [7.59e-03, 0.97])\n\nStandardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation."

# confidence interval for new data
predict(model4,
  new = data.frame(wt = 3, qsec = 18, am = "Manual"),
  interval = "confidence",
  level = .95
)

##        fit      lwr    upr
## 1 22.87005 21.09811 24.642

# prediction interval for new data
predict(model4,
  new = data.frame(wt = 3, qsec = 18, am = "Manual"),
  interval = "prediction",
  level = .95
)

##        fit      lwr      upr
## 1 22.87005 17.53074 28.20937

library(car)
linearHypothesis(model4, c("wt = 0", "qsec = 0"))

## Linear hypothesis test
## 
## Hypothesis:
## wt = 0
## qsec = 0
## 
## Model 1: restricted model
## Model 2: mpg ~ wt + qsec + am
## 
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
## 1     30 720.90                                 
## 2     28 169.29  2    551.61 45.618 1.55e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model5 <- lm(mpg ~ vs + am + as.factor(cyl),
  data = dat
)

summary(model5)

## 
## Call:
## lm(formula = mpg ~ vs + am + as.factor(cyl), data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2821 -1.4402  0.0391  1.8845  6.2179 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.809      2.928   7.789 2.24e-08 ***
## vsStraight         1.708      2.235   0.764  0.45135    
## amManual           3.165      1.528   2.071  0.04805 *  
## as.factor(cyl)6   -5.399      1.837  -2.938  0.00668 ** 
## as.factor(cyl)8   -8.161      2.892  -2.822  0.00884 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.097 on 27 degrees of freedom
## Multiple R-squared:  0.7701,	Adjusted R-squared:  0.736 
## F-statistic: 22.61 on 4 and 27 DF,  p-value: 2.741e-08

library(car)
Anova(model5)

## Anova Table (Type II tests)
## 
## Response: mpg
##                 Sum Sq Df F value  Pr(>F)  
## vs               5.601  1  0.5841 0.45135  
## am              41.122  1  4.2886 0.04805 *
## as.factor(cyl)  94.591  2  4.9324 0.01493 *
## Residuals      258.895 27                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model6 <- lm(mpg ~ wt + am + wt:am,
  data = dat
)

# Or in a shorter way:
model6 <- lm(mpg ~ wt * am,
  data = dat
)

summary(model6)

## 
## Call:
## lm(formula = mpg ~ wt * am, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6004 -1.5446 -0.5325  0.9012  6.0909 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  31.4161     3.0201  10.402 4.00e-11 ***
## wt           -3.7859     0.7856  -4.819 4.55e-05 ***
## amManual     14.8784     4.2640   3.489  0.00162 ** 
## wt:amManual  -5.2984     1.4447  -3.667  0.00102 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.591 on 28 degrees of freedom
## Multiple R-squared:  0.833,	Adjusted R-squared:  0.8151 
## F-statistic: 46.57 on 3 and 28 DF,  p-value: 5.209e-11

visreg(model6, "wt", by = "am")