This tutorial uses R’s built-in mtcars data to illustrate simple and multiple regression, diagnostics, categorical predictors, interactions, and model comparison. Code and output appear together when you knit this file to HTML.

To use: Open this file in RStudio and choose Knit > Knit to HTML. Install rmarkdown if necessary with install.packages("rmarkdown"). Analysis uses base R; document rendering uses rmarkdown and knitr.

1 Prepare the data

data(mtcars)
cars <- mtcars
cars$transmission <- factor(cars$am, levels = c(0, 1),
                            labels = c("Automatic", "Manual"))
knitr::kable(head(cars[, c("mpg", "wt", "hp", "transmission")]),
             caption = "Selected observations")
Selected observations
mpg wt hp transmission
Mazda RX4 21.0 2.620 110 Manual
Mazda RX4 Wag 21.0 2.875 110 Manual
Datsun 710 22.8 2.320 93 Manual
Hornet 4 Drive 21.4 3.215 110 Automatic
Hornet Sportabout 18.7 3.440 175 Automatic
Valiant 18.1 3.460 105 Automatic 
summary(cars[, c("mpg", "wt", "hp", "transmission")])
##       mpg              wt              hp           transmission
##  Min.   :10.40   Min.   :1.513   Min.   : 52.0   Automatic:19   
##  1st Qu.:15.43   1st Qu.:2.581   1st Qu.: 96.5   Manual   :13   
##  Median :19.20   Median :3.325   Median :123.0                  
##  Mean   :20.09   Mean   :3.217   Mean   :146.7                  
##  3rd Qu.:22.80   3rd Qu.:3.610   3rd Qu.:180.0                  
##  Max.   :33.90   Max.   :5.424   Max.   :335.0
Variable Meaning Role
mpg Miles per gallon Continuous response
wt Weight in thousands of pounds Continuous predictor
hp Gross horsepower Continuous predictor
transmission Automatic or manual Categorical predictor

2 Simple linear regression

Question: Is vehicle weight associated with fuel economy?

\[ mpg_i = \beta_0 + \beta_1 wt_i + \varepsilon_i. \]

2.1 Scatterplot and correlation

plot(mpg ~ wt, data = cars, pch = 19, col = "steelblue",
     xlab = "Weight (1,000 pounds)", ylab = "Fuel economy (mpg)",
     main = "Vehicle Weight and Fuel Economy")

cor(cars$wt, cars$mpg)
## [1] -0.8676594

A downward pattern indicates a negative association. Correlation measures linear association; it does not establish causation.

2.2 Fit and interpret the model

fit_simple <- lm(mpg ~ wt, data = cars)
summary(fit_simple)
## 
## Call:
## lm(formula = mpg ~ wt, data = cars)
## 
## 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
confint(fit_simple)
##                 2.5 %    97.5 %
## (Intercept) 33.450500 41.119753
## wt          -6.486308 -4.202635

The fitted equation is

\[ \widehat{mpg} = 37.285 - 5.344 wt. \]

An additional 1,000 pounds is associated with a predicted change of -5.34 mpg. The model explains approximately 75.3% of observed mpg variation. The intercept describes a zero-weight car and has little practical meaning.

plot(mpg ~ wt, data = cars, pch = 19, col = "steelblue",
     xlab = "Weight (1,000 pounds)", ylab = "Fuel economy (mpg)",
     main = "Simple Regression Fit")
abline(fit_simple, col = "firebrick", lwd = 2)

2.3 Correlation, R-squared, and significance tests

r <- cor(cars$wt, cars$mpg)
t_slope <- summary(fit_simple)$coefficients["wt", "t value"]
f_model <- unname(summary(fit_simple)$fstatistic["value"])
c(correlation = r, correlation_squared = r^2,
  R_squared = summary(fit_simple)$r.squared)
##         correlation correlation_squared           R_squared 
##          -0.8676594           0.7528328           0.7528328
c(t_squared = t_slope^2, F_statistic = f_model)
##   t_squared F_statistic 
##    91.37533    91.37533

For simple regression with an intercept, \(R^2=r^2\) and \(F=t^2\). The overall F-test and two-sided slope t-test give the same p-value for \(H_0:\beta_1=0\).

3 Confidence and prediction intervals

Question: What mpg does the model predict for a car weighing 3,000 pounds?

new_car <- data.frame(wt = 3)
predict(fit_simple, newdata = new_car, interval = "confidence")
##        fit      lwr      upr
## 1 21.25171 20.12444 22.37899
predict(fit_simple, newdata = new_car, interval = "prediction")
##        fit      lwr      upr
## 1 21.25171 14.92987 27.57355
weight_grid <- data.frame(wt = seq(min(cars$wt), max(cars$wt), length.out = 100))
mean_ci <- predict(fit_simple, newdata = weight_grid, interval = "confidence")
individual_pi <- predict(fit_simple, newdata = weight_grid, interval = "prediction")
plot(mpg ~ wt, data = cars, pch = 19, col = "gray40",
     ylim = range(c(cars$mpg, individual_pi)),
     xlab = "Weight (1,000 pounds)", ylab = "Fuel economy (mpg)",
     main = "Confidence and Prediction Intervals")
lines(weight_grid$wt, mean_ci[, "fit"], lwd = 2)
matlines(weight_grid$wt, mean_ci[, c("lwr", "upr")],
         col = "blue", lty = 2)
matlines(weight_grid$wt, individual_pi[, c("lwr", "upr")],
         col = "firebrick", lty = 3)
legend("topright", c("Fitted mean", "95% confidence interval", "95% prediction interval"),
       col = c("black", "blue", "firebrick"), lty = c(1, 2, 3), bty = "n")

Avoid predictions far outside the observed weight range.

4 Model diagnostics

The standard model assumes a linear conditional mean, independent errors, constant error variance, and normal errors for conventional exact inference. Predictors need not be normally distributed.

old_par <- par(mfrow = c(2, 2))
plot(fit_simple)

par(old_par)
Plot What to examine
Residuals vs Fitted Curvature or systematic patterns; residuals should scatter around zero
Normal Q-Q Points approximately following the reference line
Scale-Location Similar residual spread across fitted values
Residuals vs Leverage Observations with large residuals and high leverage

Independence must also be assessed from the study design. Investigate influential observations rather than automatically deleting them.

cook_values <- cooks.distance(fit_simple)
head(sort(cook_values, decreasing = TRUE), 5)
##   Chrysler Imperial      Toyota Corolla            Fiat 128 Lincoln Continental 
##          0.53190563          0.25987495          0.19298580          0.07192515 
##      Ford Pantera L 
##          0.03713611

5 Multiple linear regression

Question: Is horsepower associated with mpg after accounting for weight?

\[ mpg_i = \beta_0 + \beta_1 wt_i + \beta_2 hp_i + \varepsilon_i. \]

fit_multiple <- lm(mpg ~ wt + hp, data = cars)
summary(fit_multiple)
## 
## Call:
## lm(formula = mpg ~ wt + hp, data = cars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.941 -1.600 -0.182  1.050  5.854 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
## wt          -3.87783    0.63273  -6.129 1.12e-06 ***
## hp          -0.03177    0.00903  -3.519  0.00145 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8148 
## F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12
confint(fit_multiple)
##                   2.5 %      97.5 %
## (Intercept) 33.95738245 40.49715778
## wt          -5.17191604 -2.58374544
## hp          -0.05024078 -0.01330512
anova(fit_simple, fit_multiple)
## Analysis of Variance Table
## 
## Model 1: mpg ~ wt
## Model 2: mpg ~ wt + hp
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1     30 278.32                                
## 2     29 195.05  1    83.274 12.381 0.001451 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict(fit_multiple, newdata = data.frame(wt = 3, hp = 120),
        interval = "prediction")
##        fit      lwr      upr
## 1 21.78102 16.38174 27.18031

5.1 Multicollinearity

\[ VIF_j = \frac{1}{1-R_j^2}, \]

where \(R_j^2\) comes from regressing predictor \(j\) on the other predictors.

cor(cars$wt, cars$hp)
## [1] 0.6587479
r2_wt <- summary(lm(wt ~ hp, data = cars))$r.squared
r2_hp <- summary(lm(hp ~ wt, data = cars))$r.squared
c(wt = 1 / (1 - r2_wt), hp = 1 / (1 - r2_hp))
##       wt       hp 
## 1.766625 1.766625

Larger VIF values indicate greater inflation of coefficient variance. With two predictors, both VIF values are equal. A threshold such as 10 is a warning indicator, not an automatic deletion rule.

6 Categorical predictor: transmission

Question: Do transmission types differ in mean mpg after accounting for weight?

Let \(Z=0\) for automatic and \(Z=1\) for manual transmission.

\[ mpg_i = \beta_0 + \beta_1 wt_i + \beta_2 Z_i + \varepsilon_i. \]

fit_group <- lm(mpg ~ wt + transmission, data = cars)
summary(fit_group)
## 
## Call:
## lm(formula = mpg ~ wt + transmission, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5295 -2.3619 -0.1317  1.4025  6.8782 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        37.32155    3.05464  12.218 5.84e-13 ***
## wt                 -5.35281    0.78824  -6.791 1.87e-07 ***
## transmissionManual -0.02362    1.54565  -0.015    0.988    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.098 on 29 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7358 
## F-statistic: 44.17 on 2 and 29 DF,  p-value: 1.579e-09
confint(fit_group)
##                        2.5 %    97.5 %
## (Intercept)        31.074114 43.568989
## wt                 -6.964951 -3.740672
## transmissionManual -3.184815  3.137584
predict(fit_group, newdata = data.frame(
  wt = c(3, 3),
  transmission = factor(c("Automatic", "Manual"),
                        levels = levels(cars$transmission))),
  interval = "confidence")
##        fit      lwr      upr
## 1 21.26312 19.35277 23.17346
## 2 21.23950 19.24208 23.23693

transmissionManual estimates manual minus automatic mean mpg at the same weight. This model assumes equal weight slopes and a constant adjusted group difference.

# Reusable function: draw lines within each group's observed weight range
plot_transmission <- function(model, title) {
  group_colors <- c("steelblue", "firebrick")
  plot(mpg ~ wt, data = cars, pch = 19,
       col = group_colors[as.integer(cars$transmission)],
       xlab = "Weight (1,000 pounds)", ylab = "Fuel economy (mpg)",
       main = title)
  for (j in seq_along(levels(cars$transmission))) {
    group_name <- levels(cars$transmission)[j]
    weights <- cars$wt[cars$transmission == group_name]
    grid <- data.frame(
      wt = seq(min(weights), max(weights), length.out = 100),
      transmission = factor(rep(group_name, 100),
                            levels = levels(cars$transmission)))
    lines(grid$wt, predict(model, newdata = grid),
          col = group_colors[j], lwd = 2)
  }
  legend("topright", levels(cars$transmission),
         col = group_colors, pch = 19, lwd = 2, bty = "n")
}
plot_transmission(fit_group, "Additive Model: Parallel Lines")

7 Interaction: different weight slopes

\[ mpg_i = \beta_0 + \beta_1 wt_i + \beta_2 Z_i + \beta_3 wt_i Z_i + \varepsilon_i. \]

fit_interaction <- lm(mpg ~ wt * transmission, data = cars)
summary(fit_interaction)
## 
## Call:
## lm(formula = mpg ~ wt * transmission, data = cars)
## 
## 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 ***
## transmissionManual     14.8784     4.2640   3.489  0.00162 ** 
## wt:transmissionManual  -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
confint(fit_interaction)
##                           2.5 %    97.5 %
## (Intercept)           25.229642 37.602469
## wt                    -5.395234 -2.176581
## transmissionManual     6.143928 23.612917
## wt:transmissionManual -8.257693 -2.339028
anova(fit_group, fit_interaction)
## Analysis of Variance Table
## 
## Model 1: mpg ~ wt + transmission
## Model 2: mpg ~ wt * transmission
##   Res.Df    RSS Df Sum of Sq     F   Pr(>F)   
## 1     29 278.32                               
## 2     28 188.01  1    90.312 13.45 0.001017 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b <- coef(fit_interaction)
c(Automatic_slope = unname(b["wt"]),
  Manual_slope = unname(b["wt"] + b["wt:transmissionManual"]))
## Automatic_slope    Manual_slope 
##       -3.785908       -9.084268
plot_transmission(fit_interaction, "Interaction Model: Different Slopes")

The transmission main effect in this model is the difference at zero weight. Centering weight at 3 makes it the difference at 3,000 pounds:

cars$wt_centered <- cars$wt - 3
fit_centered <- lm(mpg ~ wt_centered * transmission, data = cars)
summary(fit_centered)
## 
## Call:
## lm(formula = mpg ~ wt_centered * transmission, data = cars)
## 
## 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)                     20.0583     0.8475  23.667  < 2e-16 ***
## wt_centered                     -3.7859     0.7856  -4.819 4.55e-05 ***
## transmissionManual              -1.0167     1.3209  -0.770  0.44794    
## wt_centered:transmissionManual  -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

Centering changes the interpretation of the intercept and transmission main effect while preserving fitted values and group slopes.

8 Compare candidate models

models <- list(
  Weight_only = fit_simple,
  Weight_and_HP = fit_multiple,
  Weight_and_transmission = fit_group,
  Weight_by_transmission = fit_interaction
)
comparison <- data.frame(
  Model = names(models),
  R_squared = vapply(models, function(m) summary(m)$r.squared, numeric(1)),
  Adjusted_R_squared = vapply(models, function(m) summary(m)$adj.r.squared, numeric(1)),
  AIC = vapply(models, AIC, numeric(1)),
  BIC = vapply(models, BIC, numeric(1)),
  row.names = NULL
)
knitr::kable(comparison, digits = 4, caption = "Candidate model comparison")
Candidate model comparison
Model R_squared Adjusted_R_squared AIC BIC
Weight_only 0.7528 0.7446 166.0294 170.4266
Weight_and_HP 0.8268 0.8148 156.6523 162.5153
Weight_and_transmission 0.7528 0.7358 168.0292 173.8921
Weight_by_transmission 0.8330 0.8151 157.4760 164.8046

Higher adjusted \(R^2\) and lower AIC/BIC are preferred under their respective criteria. These models use the same response and observations. Partial F-tests require nested models; AIC/BIC can also compare non-nested candidates. In-sample fit does not establish prediction accuracy on new data.

9 Student practice

  1. Fit mpg ~ hp and interpret its slope and 95% confidence interval.
  2. Compare that model with mpg ~ wt + hp. Does adding weight improve fit?
  3. Explain why a prediction interval is wider than a confidence interval for the mean.
  4. Explain why the transmission difference depends on weight in the interaction model.
  5. Examine diagnostics for your preferred model and discuss concerns.

Interpretation reminder: These observational vehicle data describe associations. Regression coefficients alone do not establish causal effects.

10 Reproducibility information

sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/Chicago
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37     R6_2.5.1          fastmap_1.2.0     xfun_0.49        
##  [5] cachem_1.1.0      knitr_1.49        htmltools_0.5.8.1 rmarkdown_2.29   
##  [9] lifecycle_1.0.5   cli_3.6.3         sass_0.4.9        jquerylib_0.1.4  
## [13] compiler_4.4.2    rstudioapi_0.19.0 tools_4.4.2       evaluate_1.0.1   
## [17] bslib_0.8.0       yaml_2.3.10       rlang_1.1.4       jsonlite_1.8.9