From An Introduction to Statistical Learning with Applications in R by Drs. Gareth James, Daniela Witten, Trevor Hastie, and Rob Tibshirani:
Introduction
Statistical Learning
Linear Regression
The book is freely available at https://www.statlearning.com/
Also, check out the following videos:
Introduction to statistical learning
Supervised and unsupervised learning
Linear regression
Statistical learning and machine learning refer to a set of computational tools for modeling data and gaining insights from data.
Generally speaking, statistical learning is considered a subfield of statistics, while machine learning is a subfield of computer science.
Traditionally, statisticians focus on understanding relationships and making inference, with an emphasis on mathematical properties and interpretability. On the other hand, computer scientists, focus on prediction and performance with large data.
Regardless, statistical learning and machine learning have a significant overlap. Both requires knowledge of mathematics and programming.
In this course, we will use statistical learning and machine learning intergeably.
Main tasks in statistical learning include:
In the above, regression and classification are considered supervised learning, in which we already have some examples to train the model.
Clustering and dimension reduction are considered unsupervised learning, in which we do not have any examples to train the model.
\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \epsilon,\]
In the above, \(\beta_0\) is the intercept, \(\beta_1, \beta_2, \dots, \beta_p\) are coefficients, and \(\epsilon\) is the unexplained random error.
The unknown intercept \(\beta_0\) and coefficients \(\beta_1, \beta_2, \dots, \beta_p\) can be estimated by the method of least squares.
\[\sum_{i = 1}^n (y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2} - \dots - \beta_p x_{ip})^2\]
Advertising is available from https://www.statlearning.com/resources-second-edition. It also consists of sales of that product in 200 different markets. It also includes advertising budgets in thousand dollars for the product in each of those markets for three different media: TV, radio, and newspapers.
load(url('https://github.com/hungtong/DS-01-101/raw/refs/heads/main/dataset/Advertising.rdata'))
head(Advertising) sales TV radio newspaper
1 22.1 230.1 37.8 69.2
2 10.4 44.5 39.3 45.1
3 9.3 17.2 45.9 69.3
4 18.5 151.5 41.3 58.5
5 12.9 180.8 10.8 58.4
6 7.2 8.7 48.9 75.0mod <- lm(sales ~ TV + radio + newspaper, data = Advertising)
summary(mod)
Call:
lm(formula = sales ~ TV + radio + newspaper, data = Advertising)
Residuals:
Min 1Q Median 3Q Max
-8.8277 -0.8908 0.2418 1.1893 2.8292
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.938889 0.311908 9.422 <2e-16 ***
TV 0.045765 0.001395 32.809 <2e-16 ***
radio 0.188530 0.008611 21.893 <2e-16 ***
newspaper -0.001037 0.005871 -0.177 0.86
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16# Intercept = 2.94. If TV = 0, radio = 0, and newspaper = 0 (we spend nothing on advertising), on average, the sales are 2.94 thousand dollars.
# Coefficient TV = 0.05. For fixed values of radio and newspaper, for an additional thousand dollars increase (one-unit change) in TV, on average, the sales increase by 0.05 x 1000 = 50 dollars.
# Coefficient radio = 0.19. For fixed values of TV and newspaper, for additional thousand dollars increase in radio, on average, the sales increase by 0.19 x 1000 = 190 dollars.
# Coefficient newspaper = -0.001. For fixed values of TV and radio, for additional thousand dollars spent on newspaper, the sales decrease by -0.0001 * x 1000 = 1 dollar. (recall that the correlation between sales and newspaper is pretty low)The standard errors of the coefficients tell us the average amount their estimates differ from the actual values. Thus, smaller standard errors mean more accurate estimates.
A small p-value, such as less than a significance level \(\alpha = 0.05\), indicates that there is statistically significant evidence of a substantial association between the predictor and the response.
# The p-values for TV and radio are pretty small (2e-16 which is much smaller than 0.05), so they seem to be more useful predictors. On the other hand, the p-value for newspaper is 0.86, which is very large, so newspaper does not seem to be very useful.# In the fitted model, the residual standard error is 1.686
# The multiple R-squared is 0.8972, meaning about 89.72% of variability in the response can be explained by the multiple linear regression model with TV, radio, and newspaper.
# The multiple R-squared will always increase with more predictors. The adjusted R-squared only increase if the predictor is useful.💻 Hands-On
The mtcars was extracted from the 1974 Motor Trend US magazine, and comprises fuel consumption and 10 aspects of automobile design and performance for 32 automobiles (1973–74 models).
Fit a multiple linear regression to learn how horsepower hp, weight wt, and number of cylinders cyl affect fuel efficiency mpg. Draw some insights from the model.
wt has a very small p-value, so it appears to be a useful predictor. All the coefficient estimates are negative, meaning heavier cars with more horsepower and cyclinders tend to be less fuel efficient (lower mpg). The multiple R-squared is 0.8431, which is close to 1. About 84.31% of variability in the mpg can be explained by the multiple linear regression model with hp, wt, and cyl.
mod <- lm(mpg ~ hp + wt + cyl, data = mtcars)
summary(mod)
Call:
lm(formula = mpg ~ hp + wt + cyl, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-3.9290 -1.5598 -0.5311 1.1850 5.8986
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.75179 1.78686 21.687 < 2e-16 ***
hp -0.01804 0.01188 -1.519 0.140015
wt -3.16697 0.74058 -4.276 0.000199 ***
cyl -0.94162 0.55092 -1.709 0.098480 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.512 on 28 degrees of freedom
Multiple R-squared: 0.8431, Adjusted R-squared: 0.8263
F-statistic: 50.17 on 3 and 28 DF, p-value: 2.184e-11mod <- lm(sales ~ TV + radio + newspaper, data = Advertising)
fitted.values(mod)new_data <- data.frame(
TV = c(145.0, 25.3, 210.5, 95.6, 12.4),
radio = c(23.1, 35.2, 18.4, 12.7, 41.8),
newspaper = c(45.0, 60.3, 30.5, 20.1, 72.4)
)
predict(mod, new_data)y <- Advertising$sales
y_hat <- fitted.values(mod)
y - y_hat\[\text{RMSE} = \sqrt{\dfrac{1}{n} \sum_{i = 1}^n (y_i - \hat{y}_i)^2}\]
rmse <- sqrt( mean((y - y_hat)^2) )
rmse[1] 1.66857In practice, fitting a model to the entire data set and evaluating how well the model predicts those same observations is not very useful, because we already know the outcomes for those observations.
We are more interested in seeing how the model predicts outcomes for new observations that were not used to fit the model.
For this reason, we typically split the data into a training set and a validation (test) set. The model is fit using the training set, and predictive performance is evaluated on the test set. This gives a more realistic measure of how well the model will perform on new data.
n <- nrow(Advertising)
# 0.7 * 200 = 140
# 0.3 * 200 = 60
set.seed(101)
train_ids <- sample(1:n, size = 140)
test_ids <- setdiff(1:n, train_ids)
train_data <- Advertising[train_ids, ]
test_data <- Advertising[test_ids, ]
mod <- lm(sales ~ TV + radio + newspaper, data = train_data)
# root mean squared errors for training data (1.7153)
(train_data$sales - fitted(mod))^2 |>
mean() |>
sqrt()[1] 1.7153# root mean squared errors for test data (1.609)
pred <- predict(mod, test_data[, c('TV', 'radio', 'newspaper')])
(test_data$sales - pred)^2 |>
mean() |>
sqrt()[1] 1.60903
💻 Hands-On
Consider a multiple linear regression in which the response is fuel efficiency mpg, and the predictors are how horsepower hp, weight wt, and number of cylinders cyl. Split the data set into a training set (about 80%) and a test set (about 20%). Compute the root mean squared errors for the training and test sets.
n <- nrow(mtcars)
# train: 0.8 * 32 = 25.6 -> 26
# test: 6
set.seed(101)
train_ids <- sample(1:n, size = 26)
test_ids <- setdiff(1:n, train_ids)
train_data <- mtcars[train_ids, ]
test_data <- mtcars[test_ids, ]
mod <- lm(mpg ~ hp + wt + cyl, data = train_data)
# root mean squared errors for train data
(train_data$mpg - fitted(mod))^2 |>
mean() |>
sqrt()[1] 2.205632# root mean squared errors for test data
pred <- predict(mod, test_data[, c('hp', 'wt', 'cyl')])
(test_data$mpg - pred)^2 |>
mean() |>
sqrt()[1] 2.953884