1. Introduction

Imagine you coach a running club and want to know why some runners are faster than others. Is it how many kilometres they run each week? How long they’ve been running? Or whether they have a coach?

Linear regression answers questions like this. It draws the best fitting straight line through your data so you can see how one measurement (the outcome, here 5 km time) changes with others (the predictors).

In other words, a regression equation looks like this:

5 km time = intercept + (slope x predictor) + error

  • The intercept is the predicted outcome when the predictor is 0.
  • The slope is how much the outcome changes for each one-unit increase in the predictor.
  • The error is the gap between the line and each runner’s actual time.

By the end of this lesson, you will be able to perform and fit a simple regression, a multiple regression and an interaction in R, and then compare them.

Set up. Install the packages below (only once in the console). Then load them each time you start R.

install.packages("tidyverse")
install.packages("dplyr")
install.packages("sjPlot")
install.packages("performance")
install.packages("GGally")
install.packages("patchwork") # needed for check model
install.packages("see") # needed for check model
library(tidyverse)
library(dplyr)
library(sjPlot)
library(performance)
library(GGally)

2. The data

We simulate data for 200 recreational runners, so we don’t need to download or import anything. (With real data, you would import your own file with read.csv).

  • time_5k: Outcome -> 5 km race time in minutes (lower = faster).
  • weekly_km: Kilometres run per week.
  • experience: Years of running experience.
  • sex: Female or Male.
  • coached: Does the runner have a coach? (Yes or No).
set.seed(2026) # This gives us the same "random" numbers / results every time

n <- 200 # number of runners

# Randomly give each runner a sex and coaching status
sex <- sample(c("Female", "Male"), n, replace = TRUE)
coached <-  sample(c("Yes", "No"), n, replace = TRUE)

# Random years of experience between 0 and 15
experience <- round(runif(n, 0, 15))

# Weekly km: more experienced runners tend to run more, plus some random variation.

# pmax makes sure nobody runs less than 5 km a week.
weekly_km <- round(pmax(5, 15 + 1.5 * experience + rnorm(n, 0, 8)))

# 5 km time (mins), built from a starting time of 30 mins:

# - males are 2.5 mins faster
# - each year of experience makes a runner 0.25 mins faster
# - each km per week above 25 makes a runner 0.08 mins faster
# - coached runners are 0.8 mins faster, and get an extra 0.12 mins of benefit per km above 25 (this is what we will look for later)
# - rnorm(n, 0, 1.5) adds random variation, as no two runners are identical

time_5k <- 30 +
  ifelse(sex == "Male", -2.5, 0) - 
  0.25 * experience - 
  0.08 * (weekly_km - 25) - 
  ifelse(coached == "Yes", 0.8 + 0.12 * (weekly_km - 25), 0) +
  rnorm(n, 0, 1.5)

# Put everything together into one table called 'runners'
runners <- data.frame(time_5k = round(time_5k, 2),
                      weekly_km, experience, sex, coached) |> dplyr::as_tibble()

3. Exploring the data

Before building any model, we explore the data. This is called exploratory data analysis (EDA). This helps us understand what the variables look like, spot anything unusual, and to get a feel for what the models are likely to find.

Always look at your data before modelling. First, the average 5 km time for coached and uncoached runners:

# This gives us our descriptive statistics of the 5 km time, so we can see a typical time and the spread of the times.
summary(runners$time_5k)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.15   24.33   26.52   26.38   28.73   32.56
runners |> 
  summarise(mean = mean(time_5k),
            sd = sd(time_5k), .by = coached)
## # A tibble: 2 × 3
##   coached  mean    sd
##   <chr>   <dbl> <dbl>
## 1 No       27.1  2.41
## 2 Yes      25.5  3.29

Boxplot. A box plot compares the 5 km times of the two groups. We use ggplot() to set up and create all types of plots and visualisations.

ggplot(runners, aes(x = coached, y = time_5k, fill = coached)) +
  geom_boxplot()

Relationships between variables. ggpairs() from the GGally package draws a grid comparing every pair of the variables listed in columns.

ggpairs(data = runners,
        columns = c("time_5k", "weekly_km", "experience"))

Look at the correlations with time_5k. Runners who run more kilometres, and runners with more experience, tend to have faster times (lower). Hence, these correlations should be negative. Now look at the correlation between weekly_km and experience: they are related to each other too, as experienced runners usually run more. This is important: when two predictors overlap, adding both to a model can change how each one looks. We will see this in Section 5.

4. Simple linear regression

Question: how is weekly distance associated with 5 km time? A simple regression has one predictor. We fit it with lm() , and the model is written outcome ~ predictor.

m1 <- lm(time_5k ~ weekly_km, data = runners)

# Estimates and p-values
summary(m1) 
## 
## Call:
## lm(formula = time_5k ~ weekly_km, data = runners)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6722 -1.4733  0.1221  1.6694  4.9451 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 31.24360    0.39145   79.81   <2e-16 ***
## weekly_km   -0.18409    0.01367  -13.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.127 on 198 degrees of freedom
## Multiple R-squared:  0.4781, Adjusted R-squared:  0.4755 
## F-statistic: 181.4 on 1 and 198 DF,  p-value: < 2.2e-16
# 95% CI
confint(m1)
##                  2.5 %     97.5 %
## (Intercept) 30.4716486 32.0155580
## weekly_km   -0.2110474 -0.1571374
# Check model/results in table
tab_model(m1)
  time 5 k
Predictors Estimates CI p
(Intercept) 31.24 30.47 – 32.02 <0.001
weekly km -0.18 -0.21 – -0.16 <0.001
Observations 200
R2 / R2 adjusted 0.478 / 0.475

Results: Each extra kilometre run per week is linked to faster 5 km times, by about 0.18 minutes. Weekly distance alone explains about 48% of the difference between runners.

Always check the model’s assumptions. check_model() draws the diagnostic plots, and each panel explains what a good plot looks like.

check_model(m1)

5. Multiple regression

Add more predictors with +. Each coefficient is then the effect of that predictor holding the others constant. This is how we perform Multiple regression.

m2 <- lm(time_5k ~ weekly_km + experience + sex + coached, data = runners)

summary(m2)
## 
## Call:
## lm(formula = time_5k ~ weekly_km + experience + sex + coached, 
##     data = runners)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6748 -1.1430  0.0010  0.9763  4.5760 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.11525    0.31835 104.023  < 2e-16 ***
## weekly_km   -0.13615    0.01232 -11.049  < 2e-16 ***
## experience  -0.21590    0.03345  -6.454 8.45e-10 ***
## sexMale     -2.09038    0.22212  -9.411  < 2e-16 ***
## coachedYes  -1.22325    0.22362  -5.470 1.37e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.548 on 195 degrees of freedom
## Multiple R-squared:  0.7276, Adjusted R-squared:  0.7221 
## F-statistic: 130.2 on 4 and 195 DF,  p-value: < 2.2e-16
tab_model(m1, m2)
  time 5 k time 5 k
Predictors Estimates CI p Estimates CI p
(Intercept) 31.24 30.47 – 32.02 <0.001 33.12 32.49 – 33.74 <0.001
weekly km -0.18 -0.21 – -0.16 <0.001 -0.14 -0.16 – -0.11 <0.001
experience -0.22 -0.28 – -0.15 <0.001
sex [Male] -2.09 -2.53 – -1.65 <0.001
coached [Yes] -1.22 -1.66 – -0.78 <0.001
Observations 200 200
R2 / R2 adjusted 0.478 / 0.475 0.728 / 0.722

tab_model(m1, m2) puts the two models side by side so we can see how the estimates changed.

Results: In the first model, each extra kilometre per week was linked to about 0.18 minutes faster. Once experience, sex, and coaching are included, this drops to about 0.14 minutes faster. It decreases as weekly kilometres were also picking up the effect of experience (experienced runner = run more and faster). Holding the other variables constant, males are about 2.09 minutes faster than females, and coached runners are about 1.22 minutes faster than uncoached runners.

6. Interactions

What if coaching changes how much each weekly kilometre helps? When the effect of one predictor depends on another, that’s an interaction.

First, plot it. facet_wrap() gives each group its own panel, and geom_smooth(method = 'lm') draws a regression line in each one. If the slopes differ between the panels, there can be an interaction.

ggplot(runners, aes(weekly_km, time_5k)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap("coached")

To note, a * b fits both predictors and their interaction:

m3 <- lm(time_5k ~ weekly_km * coached + experience + sex, data = runners)

summary(m3)
## 
## Call:
## lm(formula = time_5k ~ weekly_km * coached + experience + sex, 
##     data = runners)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8678 -0.9812 -0.0217  0.8996  4.1293 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          31.80423    0.33841  93.982  < 2e-16 ***
## weekly_km            -0.08235    0.01334  -6.175 3.80e-09 ***
## coachedYes            2.25055    0.52647   4.275 3.00e-05 ***
## experience           -0.21999    0.02986  -7.367 4.87e-12 ***
## sexMale              -2.28079    0.20001 -11.403  < 2e-16 ***
## weekly_km:coachedYes -0.13009    0.01824  -7.130 1.93e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.382 on 194 degrees of freedom
## Multiple R-squared:  0.7842, Adjusted R-squared:  0.7786 
## F-statistic:   141 on 5 and 194 DF,  p-value: < 2.2e-16

Results: For uncoached runners, each extra kilometre per week is linked to about 0.08 minutes faster. For coached runners it’s about 0.21 minutes faster (-0.08 plus -0.13). The p-value for the interaction is below .001, so the benefit of additional kilometres really is different for coached runners.

(In this model, the coached estimate of 2.25 is the difference at 0 km per week, which isn’t a realistic runner, so focus on the two slopes).

7. Comparing models

Which model is best? tab_model() with show.aicc = TRUE shows each model’s AICc, and lower is better.

tab_model(m1, m2, m3,
          show.aicc = TRUE)
  time 5 k time 5 k time 5 k
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 31.24 30.47 – 32.02 <0.001 33.12 32.49 – 33.74 <0.001 31.80 31.14 – 32.47 <0.001
weekly km -0.18 -0.21 – -0.16 <0.001 -0.14 -0.16 – -0.11 <0.001 -0.08 -0.11 – -0.06 <0.001
experience -0.22 -0.28 – -0.15 <0.001 -0.22 -0.28 – -0.16 <0.001
sex [Male] -2.09 -2.53 – -1.65 <0.001 -2.28 -2.68 – -1.89 <0.001
coached [Yes] -1.22 -1.66 – -0.78 <0.001 2.25 1.21 – 3.29 <0.001
weekly km × coached [Yes] -0.13 -0.17 – -0.09 <0.001
Observations 200 200 200
R2 / R2 adjusted 0.478 / 0.475 0.728 / 0.722 0.784 / 0.779
AICc 873.612 749.854 705.450

compare_performance() compares AIC, BIC, adjusted R-squared and RMSE (lower AIC, BIC and RMSE are better, and a higher R-squared is better), and test_performance() tests whether the differences are meaningful:

compare_performance(m1, m2, m3)
## # Comparison of Model Performance Indices
## 
## Name | Model | AIC (weights) | AICc (weights) | BIC (weights) |    R2
## ---------------------------------------------------------------------
## m1   |    lm | 873.5 (<.001) |  873.6 (<.001) | 883.4 (<.001) | 0.478
## m2   |    lm | 749.4 (<.001) |  749.9 (<.001) | 769.2 (<.001) | 0.728
## m3   |    lm | 704.9 (>.999) |  705.4 (>.999) | 728.0 (>.999) | 0.784
## 
## Name | R2 (adj.) |  RMSE | Sigma
## --------------------------------
## m1   |     0.475 | 2.117 | 2.127
## m2   |     0.722 | 1.529 | 1.548
## m3   |     0.779 | 1.361 | 1.382
test_performance(m1, m2, m3)
## Name | Model |     BF | df | df_diff | Criterion |   Chi2 |      p
## ------------------------------------------------------------------
## m1   |    lm |        |  3 |         |    867.49 |        |       
## m2   |    lm | > 1000 |  6 |       3 |    737.42 | 130.07 | < .001
## m3   |    lm | > 1000 |  7 |       1 |    690.87 |  46.55 | < .001
## Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Model 3, the interaction model, is the best: it has the lowest AICc (705.45) and the highest adjusted R-squared (0.779), as ww built that interaction into the data.

8. Key points

  • A regression coefficient is the change in the outcome for a one-unit change in a predictor, holding the others constant.
  • Adding predictors can change the other coefficients.
  • An interaction means the effect of one predictor depends on another
  • Compare modesl using AICc, BIC, adjusted R-squared and RMSE.
  • Regression shows association, not causation.