Linear Regression in R

From ISLR to practical R fluency

Advanced Statistics

Today’s plan

  • Block 1 — Simple linear regression in R
  • Block 2 — Multiple regression
  • (15 min break)
  • Block 3 — Diagnostics
  • Block 4 — Broader R fluency + guided lab

Base text: ISLR2, Chapter 3. We’ll go a bit beyond it along the way.

Block 1 — Simple Linear Regression

The model (quick recap)

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

  • \(\beta_0\), \(\beta_1\): population parameters (unknown)
  • \(\hat\beta_0\), \(\hat\beta_1\): least-squares estimates
  • Goal in R: get from raw data to \(\hat\beta_0, \hat\beta_1\), and to inference about them

We’ll use the Boston housing data (ISLR2).

Explore the data

data(Boston)
str(Boston)
'data.frame':   506 obs. of  13 variables:
 $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
 $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
 $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
 $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
 $ rm     : num  6.58 6.42 7.18 7 7.15 ...
 $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
 $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
 $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
 $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
 $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
 $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
 $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Explore the data

summary(Boston)
      crim                zn             indus            chas        
 Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
 1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
 Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
 Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
 3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
 Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
      nox               rm             age              dis        
 Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
 1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
 Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
 Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
 3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
 Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
      rad              tax           ptratio          lstat      
 Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   : 1.73  
 1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.: 6.95  
 Median : 5.000   Median :330.0   Median :19.05   Median :11.36  
 Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :12.65  
 3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:16.95  
 Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :37.97  
      medv      
 Min.   : 5.00  
 1st Qu.:17.02  
 Median :21.20  
 Mean   :22.53  
 3rd Qu.:25.00  
 Max.   :50.00

Look before you fit

pairs(Boston[, c("medv", "lstat", "age", "rm")])

Fitting the model

lm.fit <- lm(medv ~ lstat, data = Boston)
lm.fit

Call:
lm(formula = medv ~ lstat, data = Boston)

Coefficients:
(Intercept)        lstat  
      34.55        -0.95

Unpacking summary()

summary(lm.fit)

Call:
lm(formula = medv ~ lstat, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.168  -3.990  -1.318   2.034  24.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.55384    0.56263   61.41   <2e-16 ***
lstat       -0.95005    0.03873  -24.53   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.216 on 504 degrees of freedom
Multiple R-squared:  0.5441,    Adjusted R-squared:  0.5432 
F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

Pulling out individual pieces

coef(lm.fit)
(Intercept)       lstat 
 34.5538409  -0.9500494 
summary(lm.fit)$coefficients
              Estimate Std. Error   t value      Pr(>|t|)
(Intercept) 34.5538409 0.56262735  61.41515 3.743081e-236
lstat       -0.9500494 0.03873342 -24.52790  5.081103e-88
summary(lm.fit)$r.squared
[1] 0.5441463
summary(lm.fit)$sigma      # residual standard error
[1] 6.21576

Confidence intervals

confint(lm.fit)
                2.5 %     97.5 %
(Intercept) 33.448457 35.6592247
lstat       -1.026148 -0.8739505

Prediction: confidence vs. prediction interval

predict(lm.fit, newdata = data.frame(lstat = c(5, 10, 15)),
        interval = "confidence")
       fit      lwr      upr
1 29.80359 29.00741 30.59978
2 25.05335 24.47413 25.63256
3 20.30310 19.73159 20.87461
predict(lm.fit, newdata = data.frame(lstat = c(5, 10, 15)),
        interval = "prediction")
       fit       lwr      upr
1 29.80359 17.565675 42.04151
2 25.05335 12.827626 37.27907
3 20.30310  8.077742 32.52846

Plotting the fit — base R

plot(Boston$lstat, Boston$medv, pch = 20, col = "steelblue",
     xlab = "lstat", ylab = "medv")
abline(lm.fit, lwd = 2, col = "darkred")

Plotting the fit — ggplot2

ggplot(Boston, aes(x = lstat, y = medv)) +
  geom_point(alpha = 0.5, color = "steelblue") +
  geom_smooth(method = "lm", se = TRUE, color = "darkred") +
  theme_minimal()

Block 2 — Multiple Regression

Adding a predictor

lm.fit2 <- lm(medv ~ lstat + age, data = Boston)
summary(lm.fit2)

Call:
lm(formula = medv ~ lstat + age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.981  -3.978  -1.283   1.968  23.158 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 33.22276    0.73085  45.458  < 2e-16 ***
lstat       -1.03207    0.04819 -21.416  < 2e-16 ***
age          0.03454    0.01223   2.826  0.00491 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.173 on 503 degrees of freedom
Multiple R-squared:  0.5513,    Adjusted R-squared:  0.5495 
F-statistic:   309 on 2 and 503 DF,  p-value: < 2.2e-16

Regressing on everything

lm.fit.all <- lm(medv ~ ., data = Boston)
summary(lm.fit.all)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.1304  -2.7673  -0.5814   1.9414  26.2526 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  41.617270   4.936039   8.431 3.79e-16 ***
crim         -0.121389   0.033000  -3.678 0.000261 ***
zn            0.046963   0.013879   3.384 0.000772 ***
indus         0.013468   0.062145   0.217 0.828520    
chas          2.839993   0.870007   3.264 0.001173 ** 
nox         -18.758022   3.851355  -4.870 1.50e-06 ***
rm            3.658119   0.420246   8.705  < 2e-16 ***
age           0.003611   0.013329   0.271 0.786595    
dis          -1.490754   0.201623  -7.394 6.17e-13 ***
rad           0.289405   0.066908   4.325 1.84e-05 ***
tax          -0.012682   0.003801  -3.337 0.000912 ***
ptratio      -0.937533   0.132206  -7.091 4.63e-12 ***
lstat        -0.552019   0.050659 -10.897  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.798 on 493 degrees of freedom
Multiple R-squared:  0.7343,    Adjusted R-squared:  0.7278 
F-statistic: 113.5 on 12 and 493 DF,  p-value: < 2.2e-16

Marginal vs. joint effect

coef(lm(medv ~ age, data = Boston))["age"]
       age 
-0.1231627 
coef(lm.fit2)["age"]
       age 
0.03454434

Controlling for lstat changes what age “means” in the model.

Interaction terms

lm.fit.int <- lm(medv ~ lstat * age, data = Boston)
summary(lm.fit.int)

Call:
lm(formula = medv ~ lstat * age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.806  -4.045  -1.333   2.085  27.552 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 36.0885359  1.4698355  24.553  < 2e-16 ***
lstat       -1.3921168  0.1674555  -8.313 8.78e-16 ***
age         -0.0007209  0.0198792  -0.036   0.9711    
lstat:age    0.0041560  0.0018518   2.244   0.0252 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.149 on 502 degrees of freedom
Multiple R-squared:  0.5557,    Adjusted R-squared:  0.5531 
F-statistic: 209.3 on 3 and 502 DF,  p-value: < 2.2e-16

lstat * age expands to lstat + age + lstat:age automatically.

Non-linear terms: two ways

lm.fit.sq   <- lm(medv ~ lstat + I(lstat^2), data = Boston)
lm.fit.poly <- lm(medv ~ poly(lstat, 2), data = Boston)
cor(Boston$lstat, Boston$lstat^2)
[1] 0.9605726

I(lstat^2) = raw polynomial (correlated with lstat). poly() = orthogonal polynomial by default (uncorrelated terms).

Comparing coefficients

coef(lm.fit.sq)
(Intercept)       lstat  I(lstat^2) 
42.86200733 -2.33282110  0.04354689 
coef(lm.fit.poly)
    (Intercept) poly(lstat, 2)1 poly(lstat, 2)2 
       22.53281      -152.45955        64.22724

Same fit, different parameterization — coefficients aren’t directly comparable across the two.

Is the quadratic term worth it?

anova(lm.fit, lm.fit.sq)
Analysis of Variance Table

Model 1: medv ~ lstat
Model 2: medv ~ lstat + I(lstat^2)
  Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
1    504 19472                                 
2    503 15347  1    4125.1 135.2 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Qualitative predictors

data(Carseats)
contrasts(Carseats$ShelveLoc)
       Good Medium
Bad       0      0
Good      1      0
Medium    0      1
lm.fit.cs <- lm(Sales ~ . + Income:Advertising + Price:Age,
                data = Carseats)

Qualitative predictors — output

summary(lm.fit.cs)

Call:
lm(formula = Sales ~ . + Income:Advertising + Price:Age, data = Carseats)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9208 -0.7503  0.0177  0.6754  3.3413 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         6.5755654  1.0087470   6.519 2.22e-10 ***
CompPrice           0.0929371  0.0041183  22.567  < 2e-16 ***
Income              0.0108940  0.0026044   4.183 3.57e-05 ***
Advertising         0.0702462  0.0226091   3.107 0.002030 ** 
Population          0.0001592  0.0003679   0.433 0.665330    
Price              -0.1008064  0.0074399 -13.549  < 2e-16 ***
ShelveLocGood       4.8486762  0.1528378  31.724  < 2e-16 ***
ShelveLocMedium     1.9532620  0.1257682  15.531  < 2e-16 ***
Age                -0.0579466  0.0159506  -3.633 0.000318 ***
Education          -0.0208525  0.0196131  -1.063 0.288361    
UrbanYes            0.1401597  0.1124019   1.247 0.213171    
USYes              -0.1575571  0.1489234  -1.058 0.290729    
Income:Advertising  0.0007510  0.0002784   2.698 0.007290 ** 
Price:Age           0.0001068  0.0001333   0.801 0.423812    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.011 on 386 degrees of freedom
Multiple R-squared:  0.8761,    Adjusted R-squared:  0.8719 
F-statistic:   210 on 13 and 386 DF,  p-value: < 2.2e-16

Block 3 — Diagnostics

The four base-R panels

par(mfrow = c(2, 2))
plot(lm.fit2)

Reading the panels

  • Residuals vs Fitted → non-linearity
  • Q-Q → normality of residuals
  • Scale-Location → heteroscedasticity
  • Residuals vs Leverage → influential points (Cook’s distance)

Collinearity: VIF

vif(lm.fit.all)
    crim       zn    indus     chas      nox       rm      age      dis 
1.767486 2.298459 3.987181 1.071168 4.369093 1.912532 3.088232 3.954037 
     rad      tax  ptratio    lstat 
7.445301 9.002158 1.797060 2.870777

Rule of thumb: VIF > 5 (or 10) flags problematic collinearity.

Tidy diagnostics with broom

tidy(lm.fit2)
# A tibble: 3 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  33.2       0.731      45.5  2.94e-180
2 lstat        -1.03      0.0482    -21.4  8.42e- 73
3 age           0.0345    0.0122      2.83 4.91e-  3
glance(lm.fit2)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.551         0.549  6.17      309. 2.98e-88     2 -1638. 3283. 3300.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Per-observation diagnostics

aug <- augment(lm.fit2)
head(aug)
# A tibble: 6 × 9
   medv lstat   age .fitted .resid    .hat .sigma   .cooksd .std.resid
  <dbl> <dbl> <dbl>   <dbl>  <dbl>   <dbl>  <dbl>     <dbl>      <dbl>
1  24    4.98  65.2    30.3  -6.34 0.00513   6.17 0.00182       -1.03 
2  21.6  9.14  78.9    26.5  -4.92 0.00382   6.18 0.000814      -0.798
3  34.7  4.03  61.1    31.2   3.53 0.00553   6.18 0.000608       0.573
4  33.4  2.94  45.8    31.8   1.63 0.00564   6.18 0.000132       0.265
5  36.2  5.33  54.2    29.6   6.61 0.00409   6.17 0.00158        1.07 
6  28.7  5.21  58.7    29.9  -1.17 0.00437   6.18 0.0000530     -0.191

Residuals vs Fitted, tidy version

ggplot(aug, aes(.fitted, .resid)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_minimal()

Alternative view: car::residualPlots()

residualPlots(lm.fit2)
           Test stat Pr(>|Test stat|)    
lstat        13.3146           <2e-16 ***
age           1.5804           0.1146    
Tukey test   13.9913           <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Block 4 — Broader R Fluency

Formula interface tricks

lm(medv ~ . - age, data = Boston)          # drop a predictor
lm(medv ~ lstat + log(rm), data = Boston)  # transform inside formula
lm(medv ~ lstat : age, data = Boston)      # interaction only

Comparing nested models

m1 <- lm(medv ~ lstat, data = Boston)
m2 <- lm(medv ~ lstat + age + tax, data = Boston)
anova(m1, m2)
Analysis of Variance Table

Model 1: medv ~ lstat
Model 2: medv ~ lstat + age + tax
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    504 19472                                  
2    502 18683  2    789.14 10.602 3.093e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A glimpse of tidymodels

library(tidymodels)
lm_spec <- linear_reg() |> set_engine("lm")
lm_fit  <- lm_spec |> fit(medv ~ lstat + age, data = Boston)
tidy(lm_fit)

Same result as lm() — a different, increasingly common interface. Good to recognize, not required today.

Extracting results programmatically

coef(lm.fit2)
(Intercept)       lstat         age 
33.22276053 -1.03206856  0.03454434 
fitted(lm.fit2)[1:5]
       1        2        3        4        5 
30.33535 26.51520 31.17418 31.77061 29.59414 
residuals(lm.fit2)[1:5]
        1         2         3         4         5 
-6.335350 -4.915202  3.525817  1.629390  6.605862

Useful for reproducible reporting — e.g. feeding into modelsummary for publication tables.

A note looking ahead

lm() assumes independent errors.

Time series regression breaks that assumption — different diagnostics, different standard errors, different models are needed. Keep this in mind if you go further with time-dependent data.

Guided Lab — Auto dataset

Your turn

data(Auto)
str(Auto)
'data.frame':   392 obs. of  9 variables:
 $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
 $ cylinders   : int  8 8 8 8 8 8 8 8 8 8 ...
 $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
 $ horsepower  : int  130 165 150 150 140 198 220 215 225 190 ...
 $ weight      : int  3504 3693 3436 3433 3449 4341 4354 4312 4425 3850 ...
 $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
 $ year        : int  70 70 70 70 70 70 70 70 70 70 ...
 $ origin      : int  1 1 1 1 1 1 1 1 1 1 ...
 $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
 - attr(*, "na.action")= 'omit' Named int [1:5] 33 127 331 337 355
  ..- attr(*, "names")= chr [1:5] "33" "127" "331" "337" ...
  1. Fit mpg ~ horsepower. Interpret the coefficient and \(R^2\).
  2. Produce the four diagnostic plots. Evidence of non-linearity?
  3. Try mpg ~ poly(horsepower, 2). Does it improve fit (anova())?
  4. Add year and origin (as factor). Interpret origin’s coefficients.
  5. Check VIF across displacement, horsepower, weight.
  6. ggplot2 figure of mpg vs horsepower with the fitted curve overlaid.

Wrap-up

  • lm() + summary() is the core workflow — everything else builds on it
  • broom gives you tidy output for downstream work
  • car gives you diagnostic tools beyond base R
  • tidymodels exists as a parallel ecosystem — know it’s there
  • Next: how these ideas connect to time series regression

Questions?