#############################################################################
# 16. Linear Regression
#############################################################################
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(modelr) # supports basic modeling (for more advanced: https://www.tidymodels.org/ )
# we will install and load two more packages later on
# Simple Linear Regression ------------------------------------------------
# Visualizing the linear
mtcars %>%
ggplot(aes(x = wt, y = mpg)) +
geom_point() +
labs(
y = "Miles per gallon",
x = "Car's weight (1000 lbs)"
) +
theme_minimal()
# linear regression in R
model <- mtcars %>%
lm(mpg ~ wt, data = .)
summary(model)
##
## Call:
## lm(formula = mpg ~ wt, data = .)
##
## 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
# statistical significant?
# add predicted values and residuals to the original data
mtcars %>%
add_predictions(model) %>%
add_residuals(model)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
## Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
## Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
## Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2
## pred resid
## Mazda RX4 23.282611 -2.2826106
## Mazda RX4 Wag 21.919770 -0.9197704
## Datsun 710 24.885952 -2.0859521
## Hornet 4 Drive 20.102650 1.2973499
## Hornet Sportabout 18.900144 -0.2001440
## Valiant 18.793255 -0.6932545
## Duster 360 18.205363 -3.9053627
## Merc 240D 20.236262 4.1637381
## Merc 230 20.450041 2.3499593
## Merc 280 18.900144 0.2998560
## Merc 280C 18.900144 -1.1001440
## Merc 450SE 15.533127 0.8668731
## Merc 450SL 17.350247 -0.0502472
## Merc 450SLC 17.083024 -1.8830236
## Cadillac Fleetwood 9.226650 1.1733496
## Lincoln Continental 8.296712 2.1032876
## Chrysler Imperial 8.718926 5.9810744
## Fiat 128 25.527289 6.8727113
## Honda Civic 28.653805 1.7461954
## Toyota Corolla 27.478021 6.4219792
## Toyota Corona 24.111004 -2.6110037
## Dodge Challenger 18.472586 -2.9725862
## AMC Javelin 18.926866 -3.7268663
## Camaro Z28 16.762355 -3.4623553
## Pontiac Firebird 16.735633 2.4643670
## Fiat X1-9 26.943574 0.3564263
## Porsche 914-2 25.847957 0.1520430
## Lotus Europa 29.198941 1.2010593
## Ford Pantera L 20.343151 -4.5431513
## Ferrari Dino 22.480940 -2.7809399
## Maserati Bora 18.205363 -3.2053627
## Volvo 142E 22.427495 -1.0274952
# visualize
mtcars %>%
add_predictions(model) %>%
add_residuals(model) %>%
ggplot(aes(wt, mpg)) +
geom_point() +
geom_line(aes(y = pred), color = "blue", linewidth = 1)
# do you see how these relate????
mtcars %>%
ggplot(aes(wt, mpg)) +
geom_point() +
geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'
# Connect the actual data points with their corresponding predicted values
# using geom_segment()
mtcars %>%
add_predictions(model) %>%
add_residuals(model) %>%
ggplot(aes(wt, mpg)) +
geom_segment(aes(xend = wt, yend = pred)) + # connect
geom_point() +
geom_line(aes(y = pred), color = "blue", linewidth = 1) #add predicted values
# you can access the fitted values and residuals
residuals(model)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## -2.2826106 -0.9197704 -2.0859521 1.2973499
## Hornet Sportabout Valiant Duster 360 Merc 240D
## -0.2001440 -0.6932545 -3.9053627 4.1637381
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 2.3499593 0.2998560 -1.1001440 0.8668731
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## -0.0502472 -1.8830236 1.1733496 2.1032876
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 5.9810744 6.8727113 1.7461954 6.4219792
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## -2.6110037 -2.9725862 -3.7268663 -3.4623553
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 2.4643670 0.3564263 0.1520430 1.2010593
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## -4.5431513 -2.7809399 -3.2053627 -1.0274952
fitted.values(model)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 23.282611 21.919770 24.885952 20.102650
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 18.900144 18.793255 18.205363 20.236262
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 20.450041 18.900144 18.900144 15.533127
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 17.350247 17.083024 9.226650 8.296712
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 8.718926 25.527289 28.653805 27.478021
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 24.111004 18.472586 18.926866 16.762355
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 16.735633 26.943574 25.847957 29.198941
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 20.343151 22.480940 18.205363 22.427495
# Predictions
# confidence interval for new data
predict(model,
new = data.frame(wt = 3),
interval = "confidence",
level = .95
)
## fit lwr upr
## 1 21.25171 20.12444 22.37899
# Multiple linear regression ----------------------------------------------
# considering other variable
mtcars %>%
ggplot() +
aes(x = wt, y = mpg, colour = hp, size = disp) +
geom_point() +
scale_color_gradient() +
labs(
y = "Miles per gallon",
x = "Weight (1000 lbs)",
color = "Horsepower",
size = "Displacement"
) +
theme_minimal()
# Multiple linear regression in R
model2 <- mtcars %>%
lm(mpg ~ wt + hp + disp,
data = .
)
summary(model2)
##
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.891 -1.640 -0.172 1.061 5.861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.105505 2.110815 17.579 < 2e-16 ***
## wt -3.800891 1.066191 -3.565 0.00133 **
## hp -0.031157 0.011436 -2.724 0.01097 *
## disp -0.000937 0.010350 -0.091 0.92851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8083
## F-statistic: 44.57 on 3 and 28 DF, p-value: 8.65e-11
# interpretation of coefficients
# Categorical variable
model3 <- mtcars %>%
lm(mpg ~ wt + as.factor(cyl),
data = .
)
summary(model3)
##
## Call:
## lm(formula = mpg ~ wt + as.factor(cyl), data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5890 -1.2357 -0.5159 1.3845 5.7915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.9908 1.8878 18.006 < 2e-16 ***
## wt -3.2056 0.7539 -4.252 0.000213 ***
## as.factor(cyl)6 -4.2556 1.3861 -3.070 0.004718 **
## as.factor(cyl)8 -6.0709 1.6523 -3.674 0.000999 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.557 on 28 degrees of freedom
## Multiple R-squared: 0.8374, Adjusted R-squared: 0.82
## F-statistic: 48.08 on 3 and 28 DF, p-value: 3.594e-11
# interpreting the results
# Interaction term
model4 <- mtcars %>%
lm(mpg ~ wt * am, # or, you can do mpg ~ wt + am + wt:am
data = .
)
summary(model4)
##
## Call:
## lm(formula = mpg ~ wt * am, data = .)
##
## 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 ***
## am 14.8784 4.2640 3.489 0.00162 **
## wt:am -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
# interpreting results
# assumptions
library(performance)
##
## Attaching package: 'performance'
##
## The following objects are masked from 'package:modelr':
##
## mae, mse, rmse
library(see)
check_model(model2)
# Linearity - check each independent variable separately
# weight
ggplot(mtcars, aes(x = wt, y = mpg)) +
geom_point() +
theme_minimal()
# horsepower
ggplot(mtcars, aes(x = hp, y = mpg)) +
geom_point() +
theme_minimal()
# displacement
ggplot(mtcars, aes(x = disp, y = mpg)) +
geom_point() +
theme_minimal()
# The three most common tools to select a good linear model are according to:
## the p-value associated to the model,
## the coefficient of determination R2 and
## the Akaike Information Criterion