16 The Linearity Assumption

“Non-linear means it’s hard to solve.” - Arthur Mattuck

16.1 Scatterplot Matrix

The first step in assessing the linearity assumption in a multiple regression model is to examine the relationships between the predictor variables and the response variable using a scatterplot matrix. This matrix provides a grid of pairwise scatterplots, allowing you to visualize potential relationships between all pairs of variables, especially the relationship between each predictor and the response variable.

In the context of multiple regression, we expect that the relationship between each predictor and the response is approximately linear. If the scatterplot for a given predictor and the response shows a straight-line trend, this suggests a linear relationship, which aligns with the assumption. However, if the plot exhibits a clear curved or nonlinear pattern, it may indicate that the linearity assumption is violated for that predictor, and you may need to transform the variable to better capture the relationship.

What does linearity mean here?

In regression, “linear” means linear in the coefficients. The predictors themselves can be transformed.

For example,

\[ E(y)=\beta_0+\beta_1x \]

is linear in the coefficients, but so is

\[ E(y)=\beta_0+\beta_1\log(x) \]

and

\[ E(y)=\beta_0+\beta_1x+\beta_2x^2. \]

The second and third models can represent curved relationships while still being linear regression models because the \(\beta\) coefficients enter the model linearly.

A scatterplot matrix is a convenient way to quickly scan for any nonlinearity before fitting the model. You can create this matrix using the GGally package in R, which extends ggplot2 to allow for a grid of scatterplots. Each scatterplot shows how one variable changes in response to another, and from this, you can judge whether a linear transformation is needed for any predictor variables before fitting the regression model.

Scatterplot matrix for checking linearity

library(GGally)
library(tidyverse)

ggpairs(mtcars)

In this example, the ggpairs function creates a scatterplot matrix for the mtcars dataset. Each plot in the matrix shows the relationship between two variables. By focusing on the plots that involve the response variable (mpg in this case), you can assess whether the relationships between the response and each predictor (like wt or hp) appear linear. If any plot shows a nonlinear trend, it suggests that a transformation might be necessary to achieve a linear relationship.

This initial diagnostic step is crucial because it allows you to anticipate issues with linearity before fitting the regression model, ensuring a better fit and more reliable interpretation of the results.

16.1.1 Discrete Predictors

When a predictor variable has only a few discrete values, such as in the case of gear or cyl in the mtcars dataset, the scatterplot matrix will show points aligned vertically or horizontally at specific values. These discrete predictors can make it harder to judge linearity directly because the relationships are less continuous. Instead of a smooth trend, look for general patterns in how the response variable changes across the levels of the predictor. If the response values vary systematically across the predictor levels (e.g., a noticeable upward or downward shift), it may still indicate a linear trend. However, if the response varies non-linearly (e.g., higher values in the middle category and lower values on the ends), it may suggest that a transformation, such as using polynomial terms, or treating the predictor as a categorical variable, could improve the model.

For example, the cyl predictor in the mtcars dataset has only three distinct values (4, 6, and 8). If the corresponding mpg values show a clear linear decrease as cyl increases, this supports linearity. However, if the trend is irregular, a transformation or alternative approach may be necessary.

16.2 Identifying Transformations

When the scatterplot matrix suggests nonlinearity between a predictor and the response, it may be necessary to transform the predictor variable to improve the linear relationship. Common transformations include logarithmic, square root, and polynomial transformations. These transformations help capture nonlinear patterns and make relationships more linear, ensuring the regression model provides accurate estimates.

In the tidyverse, transformations can be easily applied using mutate() from the dplyr package. You can then incorporate these transformed variables into your regression model within the Tidymodels framework.

Choosing a transformation

Transformations should be chosen to match the pattern in the data, not merely to improve \(R^2\).

Pattern in scatterplot	Possible next step	Interpretation caution
Curves that flatten out	Try a log transformation such as \(\log(x)\).	A one-unit change in \(\log(x)\) is not the same as a one-unit change in \(x\).
Counts or right-skewed positive predictors	Try a square-root or log transformation.	Log transformations require positive values.
U-shaped or inverted U-shaped pattern	Try a polynomial term such as \(x^2\).	Polynomial terms can behave strangely outside the observed range of the data.
Few distinct predictor values	Consider treating the predictor as categorical.	A numeric code may not represent equal spacing between levels.

After transforming, always refit the model and check residual plots again.

Example 16.1

Log transformation in mtcars

The scatterplot matrix showed a nonlinear trend for disp versus mpg. Let’s transform this variable using a log transformation.

# Add a log-transformed variable for displacement (disp)
mtcars |>
  mutate(log_disp = log(disp)) |>
  select(mpg, log_disp, disp) |>
  ggpairs()

We see the scatterplot between mpg and log_disp appears more linear than with the untransformed disp.

Polynomial terms are still linear regression

A model with a squared predictor can represent a curved relationship while still being a linear regression model.

quadratic_fit <- lm(mpg ~ wt + I(wt^2), data = mtcars)

broom::tidy(quadratic_fit) |>
  knitr::kable(digits = 4)

term	estimate	std.error	statistic	p.value
(Intercept)	49.9308	4.2113	11.8564	0.0000
wt	-13.3803	2.5140	-5.3223	0.0000
I(wt^2)	1.1711	0.3594	3.2580	0.0029

The term I(wt^2) creates a new predictor by squaring wt. The model is curved as a function of wt, but it is still linear in the coefficients.

16.3 Fitting the Model with Transformed Variables

After identifying and applying transformations, the next step is to fit a multiple regression model using the transformed predictors. This ensures that the model aligns with the linearity assumption, yielding more reliable predictions and inferences. Below is an example of how to fit such a model.

Example 16.2

Comparing models with and without transformed predictors

We will set up two models: one with transformed variables and one without transformed variables. We will then compare the results using \(R^2\). We will transform disp, hp, and wt since these three variables appear to be nonlinear in the scatterplot matrix. We will not include the discrete variables in this example.

library(tidymodels)

# Untransformed variables
untransformed_recipe <- recipe(
  mpg ~ disp + hp + drat + wt + qsec,
  data = mtcars
)

lm_model <- linear_reg() |>
  set_engine("lm")

untransformed_workflow <- workflow() |>
  add_recipe(untransformed_recipe) |>
  add_model(lm_model)

fit_untransformed <- untransformed_workflow |>
  fit(data = mtcars)

untransformed_metrics <- fit_untransformed |>
  glance() |>
  mutate(model = "Untransformed predictors")

untransformed_metrics |>
  select(model, r.squared, adj.r.squared, sigma, statistic, p.value) |>
  knitr::kable(digits = 4)

model	r.squared	adj.r.squared	sigma	statistic	p.value
Untransformed predictors	0.8489	0.8199	2.558	29.2177	0

Now we will fit the model with the transformed variables.

# Transformed variables
transformed_recipe <- recipe(
  mpg ~ disp + hp + drat + wt + qsec,
  data = mtcars
) |>
  step_mutate(
    log_disp = log(disp),
    log_hp = log(hp),
    log_wt = log(wt)
  ) |>
  step_rm(disp, hp, wt)

transformed_workflow <- workflow() |>
  add_recipe(transformed_recipe) |>
  add_model(lm_model)

fit_transformed <- transformed_workflow |>
  fit(data = mtcars)

transformed_metrics <- fit_transformed |>
  glance() |>
  mutate(model = "Transformed predictors")

bind_rows(untransformed_metrics, transformed_metrics) |>
  select(model, r.squared, adj.r.squared, sigma, statistic, p.value) |>
  knitr::kable(digits = 4)

model	r.squared	adj.r.squared	sigma	statistic	p.value
Untransformed predictors	0.8489	0.8199	2.5580	29.2177	0
Transformed predictors	0.8922	0.8715	2.1605	43.0487	0

Note the improvement in the coefficient of determination. It has increased with the transformed variables. The residual standard error, sigma, also decreases, which suggests that the transformed model has smaller typical residuals.

Interpreting coefficients after log transformations

Once a predictor is log-transformed, its coefficient is no longer the change in \(y\) for a one-unit increase in the original predictor.

For example, in a model with

\[ E(y)=\beta_0+\beta_1\log(x), \]

\(\beta_1\) is the expected change in \(y\) for a one-unit increase in \(\log(x)\), not a one-unit increase in \(x\). A more useful interpretation is often based on percentage changes:

A 1% increase in \(x\) is associated with an approximate change of \(0.01\beta_1\) units in the mean response.
Doubling \(x\) is associated with a change of \(\beta_1\log(2)\) units in the mean response.

So transformed predictors may improve model fit, but they also change the way coefficients must be explained.

Do not rely on R-squared alone

A larger \(R^2\) can suggest that a transformed model fits better, but it does not prove that the linearity assumption has been fixed. A model can have a larger \(R^2\) and still show a curved residual pattern.

Use numerical summaries to compare fit, but use residual plots to check whether the linearity assumption looks reasonable.

16.4 Transforming the Response Variable

So far, we have transformed predictor variables. Sometimes it may also be useful to transform the response variable. This is most common when the response is positive and the relationship appears multiplicative rather than additive, or when the spread of the residuals increases as the fitted values increase.

For example, instead of modeling mpg directly, we could model log(mpg). This changes the response scale, so the residuals, fitted values, and coefficient interpretations are all on the log scale.

Transforming the response variable

The following model uses the same predictors as before, but transforms the response variable.

log_response_fit <- lm(
  log(mpg) ~ disp + hp + drat + wt + qsec,
  data = mtcars
)

broom::tidy(log_response_fit) |>
  knitr::kable(digits = 4)

term	estimate	std.error	statistic	p.value
(Intercept)	3.1089	0.4823	6.4458	0.0000
disp	0.0001	0.0005	0.1215	0.9042
hp	-0.0009	0.0007	-1.3732	0.1814
drat	0.0542	0.0576	0.9402	0.3557
wt	-0.2069	0.0547	-3.7830	0.0008
qsec	0.0247	0.0202	1.2232	0.2322

broom::glance(log_response_fit) |>
  select(r.squared, adj.r.squared, sigma, statistic, p.value) |>
  knitr::kable(digits = 4)

r.squared	adj.r.squared	sigma	statistic	p.value
0.8802	0.8572	0.1125	38.2165	0

The fitted values from this model estimate \(\log(mpg)\), not mpg. If we want predictions on the original mpg scale, we need to back-transform.

log_response_predictions <- tibble(
  observed_mpg = mtcars$mpg,
  fitted_log_mpg = fitted(log_response_fit),
  fitted_mpg = exp(fitted(log_response_fit)),
  residual_log_scale = residuals(log_response_fit)
)

log_response_predictions |>
  head() |>
  knitr::kable(digits = 4)

observed_mpg	fitted_log_mpg	fitted_mpg	residual_log_scale
21.0	3.0928	22.0398	-0.0483
21.0	3.0539	21.1984	-0.0094
22.8	3.2179	24.9767	-0.0912
21.4	3.0048	20.1828	0.0586
18.7	2.8483	17.2592	0.0802
18.1	2.9587	19.2734	-0.0628

Back-transforming with exp() gives predictions in the original units, but inference is still being performed on the log scale. This is one reason response transformations should be used thoughtfully and explained carefully.

Interpreting coefficients after transforming the response

When the response is log-transformed, coefficient interpretations change again. In a model such as

\[ \log(y)=\beta_0+\beta_1x, \]

a one-unit increase in \(x\) is associated with multiplying the typical value of \(y\) by \(e^{\beta_1}\). For small values of \(\beta_1\), this is approximately a \(100\beta_1\%\) change in \(y\).

The phrase “typical value” is doing some work here. Because the model is fit on the log scale, simply exponentiating fitted values does not automatically give the arithmetic mean of \(y\) on the original scale.

16.5 Checking Linearity After Fitting the Model

After fitting the model, it is essential to validate that the transformations improved the linear relationship between the predictors and the response. This can be achieved by examining the residuals.

If the linearity assumption holds, the residuals should appear randomly scattered around zero in the residual plot, with no obvious patterns. Systematic patterns, such as curves, indicate remaining nonlinearity, suggesting that further transformations or a different model might be necessary.

Example 16.3

Checking residual plots after transforming predictors

Next, let’s examine the residuals against each predictor in the untransformed model.

# Obtain the residuals and fitted values from the fit
predictions_untransformed <- extract_fit_engine(fit_untransformed) |>
  augment()

predictions_untransformed |>
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Fitted Values", 
       y = "Residuals", 
       title = "Residuals vs. Fitted Values") +
  theme_minimal()

We see from this residual plot that there is a nonlinear pattern. That is, the residuals tend to be above 0 at the low end, then they tend to be below 0 in the middle, then they tend to be back above 0 at the high end. If there are not that many predictors, you can plot the residuals against each predictor and determine which one needs to be transformed.

Let’s first examine the fit with the untransformed variables.

library(gridExtra)

p1 <- predictions_untransformed |>
  ggplot(aes(x = disp, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "disp", y = "Residuals")

p2 <- predictions_untransformed |>
  ggplot(aes(x = hp, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "hp", y = "Residuals")

p3 <- predictions_untransformed |>
  ggplot(aes(x = drat, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "drat", y = "Residuals")

p4 <- predictions_untransformed |>
  ggplot(aes(x = wt, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "wt", y = "Residuals")

p5 <- predictions_untransformed |>
  ggplot(aes(x = qsec, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "qsec", y = "Residuals")

grid.arrange(p1, p2, p3, p4, p5, nrow = 3)

Combining these residual plots with the scatterplot matrix, it appears that disp, hp, and wt are clearly nonlinear.

Let’s now see the residual plot for the transformed variables.

# Obtain the residuals and fitted values from the fit
predictions_transformed <- extract_fit_engine(fit_transformed) |>
  augment()

predictions_transformed |>
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Fitted Values", 
       y = "Residuals", 
       title = "Residuals vs. Fitted Values") +
  theme_minimal()

This residual plot shows no obvious nonlinear pattern. Thus, the transformations helped.

The two residual plots can also be shown side by side. This makes it easier to see whether the transformation reduced the curved pattern.

residual_comparison <- bind_rows(
  predictions_untransformed |>
    transmute(
      model = "Untransformed predictors",
      fitted = .fitted,
      residual = .resid
    ),
  predictions_transformed |>
    transmute(
      model = "Transformed predictors",
      fitted = .fitted,
      residual = .resid
    )
)

residual_comparison |>
  ggplot(aes(x = fitted, y = residual)) +
  geom_point() +
  geom_smooth(se = FALSE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  facet_wrap(~model) +
  labs(
    x = "Fitted Values",
    y = "Residuals",
    title = "Residual plots before and after transformation"
  ) +
  theme_minimal()

The side-by-side display emphasizes the purpose of transforming variables: we want residuals that look more like random scatter around 0 and less like a systematic curve.

16.6 Recap

In this chapter, we focused on how to assess and improve the linearity assumption in multiple regression.

Idea	Meaning
Linearity assumption	The expected response should be approximately linear in the predictors included in the model.
Linear in the coefficients	A regression model can include transformed predictors such as \(\log(x)\) or \(x^2\) and still be a linear regression model.
Scatterplot matrix	A grid of pairwise plots used to look for curved relationships before fitting the model.
Discrete predictor	A predictor with only a few possible values; it may need to be treated as categorical rather than quantitative.
Transformation	A change to a predictor, such as using \(\log(x)\), \(\sqrt{x}\), or \(x^2\), to better capture a nonlinear relationship.
Log-transformed predictor	A predictor entered as \(\log(x)\); the coefficient is interpreted in terms of changes in \(\log(x)\) or approximate percentage changes in \(x\).
Transformed response	A model that uses a response such as \(\log(y)\); fitted values and residuals are on the transformed scale.
Residual plot	A plot used after fitting the model to check whether nonlinear patterns remain.
Residuals vs. fitted values	A general diagnostic plot for checking whether residuals are randomly scattered around 0.
Residuals vs. predictors	Diagnostic plots that can help identify which predictor may need a transformation.
Side-by-side residual plots	A useful way to compare model diagnostics before and after transforming variables.
\(R^2\) caution	A larger \(R^2\) may indicate better fit, but residual plots are needed to assess whether the linearity assumption is reasonable.

16.7 Check your understanding

Problems

What does the linearity assumption mean in a multiple regression model?
Why can a model with \(\log(x)\) or \(x^2\) still be called a linear regression model?
Why is a scatterplot matrix useful before fitting a multiple regression model?
Why can predictors with only a few distinct values be hard to assess using scatterplots?
Give one situation where a log transformation might be useful.
Why should we check residual plots after transforming predictors?
What pattern in a residual plot suggests that the linearity assumption may still be violated?
Why should we avoid choosing transformations based only on \(R^2\)?
In a model with \(E(y)=\beta_0+\beta_1\log(x)\), why should we avoid saying that \(\beta_1\) is the effect of a one-unit increase in \(x\)?
What changes when we model \(\log(y)\) instead of \(y\)?
Why is a side-by-side residual plot useful after transforming predictors?

Solutions

The mean response should follow the form specified by the model. In a model using untransformed predictors, this means the relationship between each predictor and the expected response should be approximately straight-line, after accounting for the other predictors.
Linearity refers to the coefficients. A model such as \(E(y)=\beta_0+\beta_1\log(x)\) or \(E(y)=\beta_0+\beta_1x+\beta_2x^2\) is still linear in the \(\beta\) coefficients, even though it may be curved as a function of \(x\).
It helps identify potential nonlinear relationships early. Before fitting the model, a scatterplot matrix lets us scan the relationship between the response and each predictor and decide whether transformations may be needed.
There are too few x-values to see a smooth pattern. With only a few distinct predictor values, points stack vertically. It may be better to think about group differences or treat the predictor as categorical.
When the relationship bends and flattens out. A log transformation is often useful when increases in a positive predictor have a large effect at first but a smaller effect at larger values.
Transformations do not guarantee success. After refitting the model, residual plots show whether the curved pattern has actually been reduced.
A systematic curve. If residuals tend to be positive in one region, negative in another, and positive again elsewhere, the model is missing curvature.
\(R^2\) measures fit, not assumptions. A transformed model may have a larger \(R^2\) but still show a nonlinear residual pattern. Residual diagnostics are needed to assess the model form.
The predictor scale has changed. The coefficient \(\beta_1\) describes a one-unit increase in \(\log(x)\), not a one-unit increase in \(x\). It is often more useful to describe the effect using percentage changes in \(x\).
The model is now on the log-response scale. The fitted values, residuals, and coefficient interpretations refer to \(\log(y)\). To discuss predictions in the original units of \(y\), we need to back-transform carefully.
It makes improvement easier to see. A side-by-side plot lets us compare whether the transformed model has less systematic curvature in its residuals than the untransformed model.