9  Simple Linear Regression with Tidymodels

β€œKnowledge is a process of piling up facts; wisdom lies in their simplification.” - Martin Luther King Jr.

In earlier chapters, we fit simple linear regression models using base R functions such as lm(). In this chapter, we introduce the tidymodels framework for fitting the same kind of model in a more organized and reusable way.

tidymodels is a collection of R packages that work together to support the modeling process, including:

We will illustrate the workflow using the mtcars data first used in Example 3.1.

ExampleExample 9.1: The modeling goal

Suppose we want to model fuel economy using horsepower.

In this chapter:

  • the response variable is mpg, miles per gallon,
  • the predictor variable is hp, horsepower,
  • the statistical model is simple linear regression.

The model is \[\begin{align*} mpg = \beta_0+\beta_1 hp+\varepsilon. \end{align*}\]

The tidymodels workflow will let us fit this model, extract coefficients, make predictions, calculate metrics, and create diagnostic plots in a consistent framework.

9.1 Specifying the Model

In tidymodels, the first step is to define the type of model we want to fit. The parsnip package handles model specification.

Unlike the traditional approach where we directly call a function such as lm(), parsnip separates the model definition from the model fitting. This makes the modeling process more modular.

9.1.1 Model Specification with parsnip

The primary function for specifying a linear regression model is linear_reg(). This function says what kind of model we want, but it does not fit the model by itself.

We also specify an engine, which is the computational method used to estimate the model. For ordinary least squares regression, the engine is "lm".

ExampleExample 9.2: Specifying a linear regression model 
lm_model <- linear_reg() |>
  set_engine("lm")

lm_model
Linear Regression Model Specification (regression)

Computational engine: lm

The model specification says:

  • use linear regression as the model type,
  • use base R’s lm() function as the computational engine.

At this point, no data have been used and no coefficients have been estimated. We have only described the model we plan to fit.

NoteModel specification is not model fitting

A common beginner mistake is to think that linear_reg() fits the model. It does not.

The model specification is like saying, β€œI want to use linear regression.” The actual fitting happens later when we combine the model with a formula and data.

9.1.2 Hyperparameters and Engine Customization

Some models have hyperparameters, which are settings chosen before the model is fit. Simple linear regression with the "lm" engine has no tuning parameters, but more complex models do.

For example, ridge regression uses a penalty parameter:

ridge_model <- linear_reg(penalty = 0.1, mixture = 0) |>
  set_engine("glmnet")

Although these hyperparameters are not needed for simple linear regression, parsnip makes it easier to move from simple models to more complex models later in the course.

9.1.3 Choosing the Right Model and Engine

One strength of parsnip is that it separates the conceptual model from the computational engine.

For example:

bayesian_lm_model <- linear_reg() |>
  set_engine("stan")

The conceptual model is still linear regression, but the estimation method is Bayesian rather than ordinary least squares. This kind of modularity becomes more useful as the models become more complex.

9.2 Defining the Workflow

In the tidymodels ecosystem, a workflow is an organizing structure that ties together the pieces of the modeling process.

A workflow can contain:

  • the model specification,
  • the formula,
  • preprocessing steps, and
  • eventually the fitted model.

9.2.1 Creating a Workflow for Simple Linear Regression

For simple linear regression, our workflow consists of two main components:

  1. the model specification, and
  2. the formula that defines the response and predictor.
ExampleExample 9.3: Building a workflow

We will predict mpg using hp from the mtcars dataset.

lm_workflow <- workflow() |>
  add_model(lm_model) |>
  add_formula(mpg ~ hp)

lm_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Formula
Model: linear_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
mpg ~ hp

── Model ───────────────────────────────────────────────────────────────────────
Linear Regression Model Specification (regression)

Computational engine: lm

In this code:

  • workflow() initializes an empty workflow,
  • add_model() adds the model specification,
  • add_formula() defines the response and predictor.

The formula mpg ~ hp means that mpg is the response and hp is the predictor.

9.2.2 Formula Interface

The formula interface uses the tilde symbol, ~, to separate the response variable from the predictor variables.

For example:

mpg ~ hp

means β€œmodel mpg using hp.”

If we had multiple predictors, the formula could be written as

mpg ~ hp + wt + qsec

which would fit a multiple regression model.

9.2.3 Preprocessing in Workflows

Simple linear regression may not require much preprocessing, but workflows can also include preprocessing through the recipes package.

Common preprocessing steps include:

  • standardizing predictors,
  • handling missing data,
  • encoding categorical predictors,
  • transforming variables.

For example, if am were treated as a categorical predictor, we could create dummy variables using a recipe. This is a small example, but it shows the key idea: preprocessing can live inside the modeling workflow instead of being handled separately.

ExampleExample 9.4: Using a recipe to create dummy variables

First, treat am as a categorical variable rather than a numeric 0/1 variable.

mtcars_recipe_data <- mtcars |>
  mutate(
    am = factor(
      am,
      levels = c(0, 1),
      labels = c("automatic", "manual")
    )
  )

mtcars_recipe <- recipe(mpg ~ hp + am, data = mtcars_recipe_data) |>
  step_dummy(all_nominal_predictors())

mtcars_recipe |>
  prep() |>
  bake(new_data = NULL) |>
  select(mpg, hp, starts_with("am_")) |>
  head() |>
  knitr::kable(digits = 3)
mpg hp am_manual
21.0 110 1
21.0 110 1
22.8 93 1
21.4 110 0
18.7 175 0
18.1 105 0

The recipe has converted the categorical predictor am into a dummy variable that can be used by the linear regression model.

Now we can put the recipe inside a workflow and fit the model.

recipe_workflow <- workflow() |>
  add_model(lm_model) |>
  add_recipe(mtcars_recipe)

recipe_fit <- recipe_workflow |>
  fit(data = mtcars_recipe_data)

extract_fit_engine(recipe_fit) |>
  tidy() |>
  knitr::kable(digits = 4)
term estimate std.error statistic p.value
(Intercept) 26.5849 1.4251 18.6548 0
hp -0.0589 0.0079 -7.4952 0
am_manual 5.2771 1.0795 4.8883 0

This model is no longer simple linear regression because it uses both hp and am as predictors. The point of the example is not to replace our mpg ~ hp model, but to show why recipes are useful: the workflow remembers the preprocessing steps and applies them consistently.

ExampleExample 9.5: Why workflows help

Imagine fitting a model where you:

  1. replace missing values,
  2. create dummy variables,
  3. standardize predictors,
  4. fit the model,
  5. predict on new data.

If those steps are written separately, it is easy to accidentally preprocess the training data differently from the new data.

A workflow helps prevent that mistake by keeping the preprocessing and modeling steps bundled together.

9.3 Fitting the Model

Once the model and workflow have been defined, the next step is to fit the model to data.

In tidymodels, this is done with fit().

ExampleExample 9.6: Fitting the workflow 
lm_fit <- lm_workflow |>
  fit(data = mtcars)

lm_fit
══ Workflow [trained] ══════════════════════════════════════════════════════════
Preprocessor: Formula
Model: linear_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
mpg ~ hp

── Model ───────────────────────────────────────────────────────────────────────

Call:
stats::lm(formula = ..y ~ ., data = data)

Coefficients:
(Intercept)           hp  
   30.09886     -0.06823

The fitted workflow contains the model specification, formula, and estimated model.

Behind the scenes, tidymodels uses the "lm" engine, so the fitted model is ultimately based on base R’s lm() function.

9.3.1 Extracting Model Coefficients

To extract the estimated coefficients in a tidy table, we can use tidy().

For a fitted workflow, it is often useful to extract the fitted engine first and then apply tidy().

ExampleExample 9.7: Extracting and interpreting coefficients 
lm_engine <- extract_fit_engine(lm_fit)

coefs <- tidy(lm_engine)

knitr::kable(coefs, digits = 4)
term estimate std.error statistic p.value
(Intercept) 30.0989 1.6339 18.4212 0
hp -0.0682 0.0101 -6.7424 0

The estimated slope for hp is -0.0682. This means that for each additional unit of horsepower, the predicted fuel economy decreases by about 0.0682 miles per gallon, on average.

The estimated intercept is 30.0989. In this context, the intercept is the predicted mpg when hp = 0. Since hp = 0 is not realistic for the cars in this dataset, the intercept is needed mathematically but is not very meaningful as a practical interpretation.

ExampleExample 9.8: Comparing Tidymodels to base R

Since we used the "lm" engine, the fitted coefficients agree with base R’s lm() function.

base_lm_fit <- lm(mpg ~ hp, data = mtcars)

bind_rows(
  tidy(lm_engine) |> mutate(method = "tidymodels workflow"),
  tidy(base_lm_fit) |> mutate(method = "base lm")
) |>
  select(method, term, estimate, std.error, statistic, p.value) |>
  knitr::kable(digits = 4)
method term estimate std.error statistic p.value
tidymodels workflow (Intercept) 30.0989 1.6339 18.4212 0
tidymodels workflow hp -0.0682 0.0101 -6.7424 0
base lm (Intercept) 30.0989 1.6339 18.4212 0
base lm hp -0.0682 0.0101 -6.7424 0

The advantage of tidymodels is not that it changes ordinary least squares. The advantage is that it gives us a consistent framework that scales well as modeling tasks become more complex.

9.3.2 Base R and Tidymodels Side-by-Side

For ordinary least squares regression, tidymodels can use base R’s lm() function as the engine. That means the estimated line is the same, but the organization of the modeling process is different.

Modeling task Base R approach Tidymodels approach
Specify the model Usually done inside lm() linear_reg() and set_engine("lm")
Attach the formula lm(mpg ~ hp, data = mtcars) add_formula(mpg ~ hp) or add_recipe()
Fit the model lm() fits immediately fit(data = mtcars) fits the workflow
Extract coefficients coef(), summary(), or broom::tidy() extract_fit_engine() and tidy()
Make predictions predict(object, newdata = ...) predict(object, new_data = ...)
Add preprocessing Manually before fitting Store steps in a recipe()
Reuse the process Recreate steps by hand Reuse the workflow object

The two approaches are not enemies. Base R is direct and compact. Tidymodels is more structured and becomes especially helpful when the modeling process includes preprocessing, resampling, tuning, or repeated prediction on new data.

9.4 Making Predictions

After fitting the model, we can use predict() to generate fitted values for existing data or predictions for new data.

ExampleExample 9.9: Predicting on the original data 
mtcars_with_preds <- lm_fit |>
  predict(new_data = mtcars) |>
  bind_cols(mtcars |> select(mpg, hp))

head(mtcars_with_preds) |>
  knitr::kable(digits = 3)
.pred mpg hp
22.594 21.0 110
22.594 21.0 110
23.754 22.8 93
22.594 21.4 110
18.159 18.7 175
22.935 18.1 105

The .pred column contains the predicted value of mpg from the fitted model.

Comparing .pred to the observed mpg helps us understand how far the model’s predictions are from the actual values.

ExampleExample 9.10: Predicting for new horsepower values

Suppose we want predicted fuel economy for cars with horsepower values 100, 150, and 200.

new_cars <- tibble(hp = c(100, 150, 200))

lm_fit |>
  predict(new_data = new_cars) |>
  bind_cols(new_cars) |>
  select(hp, predicted_mpg = .pred) |>
  knitr::kable(digits = 3)
hp predicted_mpg
100 23.276
150 19.865
200 16.453

These are point predictions from the fitted line. They are not prediction intervals; they do not show uncertainty around the predictions.

9.5 Performance Metrics

Once predictions have been made, we can evaluate model performance.

Common regression metrics include:

  • RMSE: root mean squared error,
  • MAE: mean absolute error,
  • \(R^2\): proportion of variation explained.
ExampleExample 9.11: Calculating model metrics 
model_metrics <- lm_fit |>
  predict(new_data = mtcars) |>
  bind_cols(mtcars) |>
  yardstick::metrics(truth = mpg, estimate = .pred)

knitr::kable(model_metrics, digits = 4)
.metric .estimator .estimate
rmse standard 3.7403
rsq standard 0.6024
mae standard 2.9075

The RMSE is approximately 3.74, so the model’s fitted values are typically off by about that many miles per gallon.

The \(R^2\) value is approximately 0.602, so the model explains about 60.2% of the sample variation in mpg.

WarningTraining metrics are optimistic

The metrics above are calculated on the same data used to fit the model. These are training-set metrics.

Training-set performance can make a model look better than it will perform on new data. Later, we will use tools such as data splitting, cross-validation, and resampling to estimate performance more honestly.

9.6 Model Diagnostics

After fitting a linear regression model, we still need to check whether the assumptions appear reasonable.

The main assumptions are:

  • linearity,
  • constant variance,
  • independence,
  • normality of errors.

9.6.1 Residuals and Fitted Values

We can generate residuals and fitted values using augment() from the broom package.

To use augment() with the fitted lm object, we first extract the fitted engine from the workflow.

ExampleExample 9.12: Creating diagnostic data with `augment()` 
diagnostic_data <- lm_engine |>
  augment()

head(diagnostic_data) |>
  knitr::kable(digits = 3)
..y hp .fitted .resid .hat .sigma .cooksd .std.resid
21.0 110 22.594 -1.594 0.040 3.917 0.004 -0.421
21.0 110 22.594 -1.594 0.040 3.917 0.004 -0.421
22.8 93 23.754 -0.954 0.051 3.925 0.002 -0.253
21.4 110 22.594 -1.194 0.040 3.922 0.002 -0.315
18.7 175 18.159 0.541 0.037 3.928 0.000 0.143
18.1 105 22.935 -4.835 0.043 3.820 0.037 -1.280

The augmented dataset includes columns such as:

  • .fitted: fitted values,
  • .resid: residuals,
  • .std.resid: standardized residuals.

These columns are useful for diagnostic plots.

9.6.2 Residuals versus Fitted Values

A residuals-versus-fitted-values plot helps assess linearity and constant variance.

ExampleExample 9.13: Residuals versus fitted values 
ggplot(diagnostic_data, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    x = "Fitted values",
    y = "Residuals",
    title = "Residuals vs. fitted values"
  ) +
  theme_minimal()

If the model is appropriate, we hope to see residuals scattered randomly around 0 with no strong curve or funnel shape.

9.6.3 Normality of Residuals

A Q-Q plot helps us check whether the residuals are approximately normal.

ExampleExample 9.14: Q-Q plot of residuals 
ggplot(diagnostic_data, aes(sample = .resid)) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(
    title = "Q-Q plot of residuals",
    x = "Theoretical quantiles",
    y = "Sample residual quantiles"
  ) +
  theme_minimal()

If the residuals are approximately normal, the points should roughly follow the reference line.

9.6.4 Independence of Errors

If the data are ordered in time, space, or another meaningful sequence, we should check for autocorrelation in the residuals.

ExampleExample 9.15: ACF plot of residuals 
acf(diagnostic_data$.resid, main = "ACF of residuals")

For mtcars, the row order is not a meaningful time order, so this plot is not as important as it would be for time series data. The main lesson is procedural: when observations have a meaningful order, check residual dependence.

9.6.5 Detecting Potential Outliers

Standardized residuals help identify observations that are unusually far from the fitted line.

ExampleExample 9.16: Standardized residuals 
ggplot(diagnostic_data, aes(x = .fitted, y = .std.resid)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "gray40") +
  geom_hline(yintercept = c(-2, 2), linetype = "dashed", color = "orange") +
  geom_hline(yintercept = c(-3, 3), linetype = "dotted", color = "red") +
  labs(
    x = "Fitted values",
    y = "Standardized residuals",
    title = "Standardized residuals vs. fitted values"
  ) +
  theme_minimal()

Observations outside the dashed or dotted lines should be investigated, especially if they also have high leverage or strong influence.

WarningCommon Tidymodels troubleshooting

When Tidymodels code fails, the issue is often about the object being used rather than the regression idea itself.

Symptom Likely issue What to check
predict() gives an error about new_data The argument name is wrong or the object is not fitted Use new_data = ... and make sure the workflow has been fit
New predictions fail after using a recipe The new data do not contain the variables required by the recipe Supply the original predictor columns expected by the recipe, not manually created dummy columns
Coefficients are not available directly from the workflow The fitted engine is nested inside the workflow Use extract_fit_engine() before tidy()
The model uses the wrong response or predictor The formula sides were reversed Remember: response on the left, predictors on the right
Results differ from what you expected A preprocessing step changed the variables used by the model Inspect the recipe with prep() and bake()

The useful habit is to ask, β€œWhich object do I have right now: a model specification, a recipe, an unfitted workflow, or a fitted workflow?” Many Tidymodels errors become easier to diagnose once that distinction is clear.

9.7 Recap

This chapter introduced a tidymodels workflow for simple linear regression.

Idea Meaning
linear_reg() Specifies a linear regression model.
Engine The computational method used to fit the model, such as "lm".
Workflow An object that bundles the model specification with a formula or recipe.
recipe() Stores preprocessing steps such as dummy-variable creation, transformations, or standardization.
add_recipe() Adds a recipe to a workflow instead of using a plain formula.
fit() Fits the workflow to data.
tidy() Extracts model coefficients and inferential summaries in a clean table.
predict() Generates predictions for existing or new data.
yardstick::metrics() Calculates model performance metrics such as RMSE and \(R^2\).
augment() Adds fitted values, residuals, and diagnostic quantities to the data.
Training-set metrics Performance metrics calculated on the same data used to fit the model; often optimistic.
Base R vs. Tidymodels The fitted least-squares line can be identical, but Tidymodels organizes the modeling process more explicitly.
Diagnostics Residual-based checks for model assumptions.

9.8 Check your understanding

NoteProblems

  1. What is the difference between specifying a model with linear_reg() and fitting a model with fit()?

  2. What role does the engine play in a parsnip model specification?

  3. Why might a workflow be useful even for simple linear regression?

  4. What problem does a recipe() solve in the modeling workflow?

  5. If Tidymodels uses the "lm" engine, why should the coefficient estimates agree with base R’s lm()?

  6. In the formula mpg ~ hp, which variable is the response and which is the predictor?

  7. Why is the intercept in the mpg ~ hp model not very meaningful in practice?

  8. What does predict() return in the examples from this chapter?

  9. Why are training-set metrics often optimistic?

  10. What is the difference between tidy() and augment()?

  11. If prediction fails after you use a recipe, what is one of the first things you should check?

  12. Why is the ACF plot not especially meaningful for mtcars, but important for time-ordered data?

  1. Specification versus fitting. linear_reg() describes the type of model to use. fit() estimates the model parameters from data.

  2. The engine performs the computation. The engine specifies which underlying function or package will fit the model. For ordinary least squares regression, the "lm" engine uses base R’s lm() function.

  3. Workflows organize the process. A workflow keeps the model specification, formula, and preprocessing steps together. This makes the modeling process more reproducible and easier to modify later.

  4. Recipes store preprocessing. A recipe records steps such as dummy-variable creation, transformations, or standardization. This helps ensure that the same preprocessing is applied when the model is fit and when predictions are made on new data.

  5. Same engine, same least-squares fit. When Tidymodels uses the "lm" engine, it is using base R’s ordinary least squares machinery behind the scenes. The organization is different, but the estimated line should be the same.

  6. Left side is the response. In mpg ~ hp, mpg is the response variable and hp is the predictor variable.

  7. Zero horsepower is not realistic. The intercept is the predicted mpg when hp = 0. Since cars in this dataset do not have 0 horsepower, the intercept is needed mathematically but is not practically meaningful.

  8. Predicted values. predict() returns a tibble with a .pred column containing predicted values from the fitted model.

  9. The model has already seen the data. Training-set metrics evaluate the model on the same data used to fit it, so they usually make the model look better than it may perform on new data.

  10. Model summary versus observation-level diagnostics. tidy() summarizes model coefficients, standard errors, test statistics, and p-values. augment() gives observation-level information such as fitted values and residuals.

  11. Check the required predictor columns. New data should contain the original predictor variables expected by the workflow or recipe. If the recipe creates dummy variables, do not usually supply the dummy variables yourself; supply the original categorical predictor and let the recipe handle the encoding.

  12. Order matters for autocorrelation. The ACF plot is meaningful when observations have a meaningful order, such as time. The row order of mtcars is not a natural time sequence, so autocorrelation is not a central concern for that dataset.