Assignment on Multiple Linear Regression Analysis Using mtcars
20251MBI036 Marie Francine Kwizera
2026-06-04
HOMEWORK 1: Find a data set of which you can fit multiple linear regression and interpret your results
1. LOAD AND EXPLORE THE DATASET
data(mtcars)
head(mtcars)## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1str(mtcars)## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...summary(mtcars)## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000Variables used:
mpg — miles per gallon (outcome/response) cyl — number of cylinders disp — engine displacement (cu. in.) hp — gross horsepower wt — weight (1000 lbs) gear — number of forward gears
2. EXPLORATORY DATA ANALYSIS
2a. Correlation matrix (numeric)
round(cor(mtcars[, c("mpg", "cyl", "disp", "hp", "wt", "gear")]), 2)## mpg cyl disp hp wt gear
## mpg 1.00 -0.85 -0.85 -0.78 -0.87 0.48
## cyl -0.85 1.00 0.90 0.83 0.78 -0.49
## disp -0.85 0.90 1.00 0.79 0.89 -0.56
## hp -0.78 0.83 0.79 1.00 0.66 -0.13
## wt -0.87 0.78 0.89 0.66 1.00 -0.58
## gear 0.48 -0.49 -0.56 -0.13 -0.58 1.002b. Visual correlation plot
corrplot(
cor(mtcars[, c("mpg", "cyl", "disp", "hp", "wt", "gear")]),
method = "color",
type = "upper",
addCoef.col = "black",
tl.cex = 0.9,
number.cex = 0.8,
title = "Correlation Matrix — mtcars predictors",
mar = c(0, 0, 1, 0)
)
# 2c. Scatterplot matrix (pairs plot)
pairs(
mtcars[, c("mpg", "cyl", "disp", "hp", "wt", "gear")],
main = "Scatterplot Matrix",
pch = 19,
col = "steelblue"
)
# 2d. Distribution of the outcome variable
hist(
mtcars$mpg,
breaks = 10,
main = "Distribution of MPG",
xlab = "Miles per Gallon",
col = "steelblue",
border = "white"
)
# 3. FIT THE MULTIPLE LINEAR REGRESSION MODEL
model <- lm(mpg ~ cyl + disp + hp + wt + gear, data = mtcars)
summary(model)##
## Call:
## lm(formula = mpg ~ cyl + disp + hp + wt + gear, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.223 -1.685 -0.383 1.293 5.943
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.36255 5.97248 6.256 1.28e-06 ***
## cyl -1.11864 0.71436 -1.566 0.12946
## disp 0.01380 0.01232 1.121 0.27274
## hp -0.02787 0.01660 -1.679 0.10519
## wt -3.71428 1.04818 -3.544 0.00152 **
## gear 0.67881 1.03454 0.656 0.51750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.539 on 26 degrees of freedom
## Multiple R-squared: 0.8511, Adjusted R-squared: 0.8225
## F-statistic: 29.72 on 5 and 26 DF, p-value: 5.723e-104. COEFFICIENTS INTERPRETATION
Coefficient Table
coef_table <- summary(model)$coefficients
print(round(coef_table, 4))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.3626 5.9725 6.2558 0.0000
## cyl -1.1186 0.7144 -1.5659 0.1295
## disp 0.0138 0.0123 1.1205 0.2727
## hp -0.0279 0.0166 -1.6787 0.1052
## wt -3.7143 1.0482 -3.5435 0.0015
## gear 0.6788 1.0345 0.6561 0.517595% Confidence Intervals
confint(model, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 25.08595153 49.639155453
## cyl -2.58703596 0.349746991
## disp -0.01151918 0.039127866
## hp -0.06200258 0.006256271
## wt -5.86884738 -1.559703466
## gear -1.44772056 2.8053429945. MODEL DIAGNOSTICS
Diagnostic Plots
par(mfrow = c(2, 2))
plot(model, main = "")
par(mfrow = c(1, 1)) INTERPRETATION
Residuals vs Fitted: assesses linearity. Normal Q-Q: assesses normality of residuals. Scale-Location: assesses homoscedasticity. Residuals vs Leverage: identifies influential observations.
Breusch-Pagan Test
bptest(model)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 3.1724, df = 5, p-value = 0.6734Null Hypothesis: Constant variance of residuals.
Variance Inflation Factor (VIF)
vif(model)## cyl disp hp wt gear
## 7.824309 11.207264 6.229789 5.056423 2.800679Shapiro-Wilk Test for Normality of Residuals
shapiro.test(residuals(model))##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.91262, p-value = 0.01318Null Hypothesis: Residuals are normally distributed.
Influential Observations (Cook’s distance)
cooksd <- cooks.distance(model)
plot(cooksd, type = "h",
main = "Cook's Distance — Influential Observations",
ylab = "Cook's Distance")
abline(h = 4 / nrow(mtcars), col = "red", lty = 2)
text(which(cooksd > 4 / nrow(mtcars)),
cooksd[cooksd > 4 / nrow(mtcars)],
labels = names(cooksd[cooksd > 4 / nrow(mtcars)]),
pos = 2, cex = 0.8)
# 6. MODEL PERFORMANCE METRICS
r2 <- summary(model)$r.squared
adj_r2 <- summary(model)$adj.r.squared
rmse <- sqrt(mean(residuals(model)^2))
cat("R-squared :", round(r2, 4), "\n")## R-squared : 0.8511cat("Adjusted R-squared:", round(adj_r2, 4), "\n")## Adjusted R-squared: 0.8225cat("RMSE (train) :", round(rmse, 4), "mpg\n")## RMSE (train) : 2.289 mpg7. VISUALISE ACTUAL vs PREDICTED
mtcars$predicted <- fitted(model)
ggplot(mtcars, aes(x = predicted, y = mpg)) +
geom_point(color = "steelblue", size = 2.5, alpha = 0.8) +
geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
labs(
title = "Actual vs Predicted MPG",
subtitle = "Points on the red line = perfect prediction",
x = "Predicted MPG",
y = "Actual MPG"
) +
theme_minimal()
# 8. REDUCED MODEL (drop insignificant predictors)
A reduced model using only horsepower and weight is fitted for comparison.
model_reduced <- lm(mpg ~ hp + wt, data = mtcars)
summary(model_reduced)##
## Call:
## lm(formula = mpg ~ hp + wt, data = mtcars)
##
## 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 ***
## hp -0.03177 0.00903 -3.519 0.00145 **
## wt -3.87783 0.63273 -6.129 1.12e-06 ***
## ---
## 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-12ANOVA Comparison
anova(model_reduced, model)## Analysis of Variance Table
##
## Model 1: mpg ~ hp + wt
## Model 2: mpg ~ cyl + disp + hp + wt + gear
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 29 195.05
## 2 26 167.67 3 27.38 1.4152 0.26089. MODEL SUMMARY
cat("Observations:", nrow(mtcars), "\n")## Observations: 32cat("Predictors :", 5, "\n")## Predictors : 5cat("R² :", round(summary(model)$r.squared, 3), "\n")## R² : 0.851cat("Adj. R² :", round(summary(model)$adj.r.squared, 3), "\n")## Adj. R² : 0.822cat("F-statistic :", round(summary(model)$fstatistic[1], 2),
"on", summary(model)$fstatistic[2], "and",
summary(model)$fstatistic[3], "DF\n")## F-statistic : 29.72 on 5 and 26 DFcat("F p-value :",
round(pf(summary(model)$fstatistic[1],
summary(model)$fstatistic[2],
summary(model)$fstatistic[3],
lower.tail = FALSE), 6), "\n")## F p-value : 0HOMEWORK 2: Variable Selection Methods
1. Introduction
Variable selection refers to the process of identifying the subset of
predictors that are most relevant for explaining the response variable.
Including too many irrelevant predictors leads to overfitting, inflated
variance, and poor generalisability. I will keep using the
mtcars dataset, with mpg (miles per gallon) as
the response variable to demonstrates and compares five variable
selection approaches.
1. Stepwise Selection (AIC)
Stepwise selection uses AIC (Akaike Information Criterion) to add or remove predictors. Lower AIC = better model.
full_model <- lm(mpg ~ ., data = mtcars)
step_model <- stepAIC(full_model, direction = "both", trace = FALSE)
summary(step_model)##
## Call:
## lm(formula = mpg ~ predicted, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.223 -1.685 -0.383 1.293 5.943
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00000 1.59013 0.0 1
## predicted 1.00000 0.07637 13.1 6.11e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.364 on 30 degrees of freedom
## Multiple R-squared: 0.8511, Adjusted R-squared: 0.8461
## F-statistic: 171.5 on 1 and 30 DF, p-value: 6.107e-142. Best Subset Selection
Best subset evaluates all possible combinations of predictors and ranks models of each size by Adjusted R² and BIC.
best_sub <- regsubsets(mpg ~ ., data = mtcars, nvmax = 10)## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 1 linear dependencies foundbest_sum <- summary(best_sub)
par(mfrow = c(1, 2))
plot(best_sum$adjr2, type = "b", xlab = "Predictors", ylab = "Adjusted R²",
main = "Adjusted R²", col = "steelblue", pch = 19)
points(which.max(best_sum$adjr2), max(best_sum$adjr2), col = "red", pch = 19, cex = 1.5)
plot(best_sum$bic, type = "b", xlab = "Predictors", ylab = "BIC",
main = "BIC (lower = better)", col = "steelblue", pch = 19)
points(which.min(best_sum$bic), min(best_sum$bic), col = "red", pch = 19, cex = 1.5)
par(mfrow = c(1, 1))
cat("Best number of predictors by Adjusted R²:", which.max(best_sum$adjr2), "\n")## Best number of predictors by Adjusted R²: 4cat("Best number of predictors by BIC :", which.min(best_sum$bic), "\n")## Best number of predictors by BIC : 13. Lasso Regression
Lasso applies an L1 penalty that shrinks some coefficients exactly to zero, performing automatic variable selection.
x <- model.matrix(mpg ~ ., mtcars)[, -1]
y <- mtcars$mpg
cv_lasso <- cv.glmnet(x, y, alpha = 1)
plot(cv_lasso, main = "Lasso Cross-Validation")
cat("Optimal lambda:", round(cv_lasso$lambda.min, 4), "\n\n")## Optimal lambda: 0.644cat("Lasso coefficients (zeros = excluded variables):\n")## Lasso coefficients (zeros = excluded variables):print(coef(cv_lasso, s = "lambda.min"))## 12 x 1 sparse Matrix of class "dgCMatrix"
## lambda.min
## (Intercept) 2.3642888
## cyl .
## disp .
## hp .
## drat .
## wt .
## qsec .
## vs .
## am .
## gear .
## carb .
## predicted 0.8823188Summary
| Method | How it works | Selects variables? |
|---|---|---|
| Stepwise (AIC) | Adds/removes predictors by AIC | Yes |
| Best Subset | Tests all combinations | Yes |
| Lasso (L1) | Shrinks coefficients to zero | Yes (automatically) |
Stepwise and best subset are discrete — a variable is either in or out. Lasso is continuous — it gradually shrinks coefficients and is preferred when predictors are correlated or when many candidates exist.