0:55 to 1:12

library(vembedr)
embed_url(
  "https://www.youtube.com/watch?v=TJgUiZgX5rE"
)

\[ Y=f(x) + \epsilon \]

Here Y is Dist. and X is Speed

Stop Dist. vs. Speed :

data give the speed of cars and the distances taken to stop. Note that the data were recorded in the 1920s.

Research Question :

Ie how fucked are pedestrians if ur speeding?

Modeling Philosophy

# data(cars) 
?cars # investigate data set

Sampling Framework

Test vs. Training Data

# sample!
set.seed(1)

idx <- sample(
  1:nrow(cars),
  size = 40,
  replace = FALSE
)

sample <- cars[idx, ]
test <- cars[-idx, ]

# Common axis limits for all three plots
x_limits <- range(cars$speed)
y_limits <- c(0, max(cars$dist))

par(mfrow = c(1, 3))

# Training data
plot(
  sample$speed,
  sample$dist,
  main = "Training Data",
  xlab = "Speed (mph)",
  ylab = "Stopping Distance (ft)",
  col = "forestgreen",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

# Test data
plot(
  test$speed,
  test$dist,
  main = "Test Data",
  xlab = "Speed (mph)",
  ylab = "Stopping Distance (ft)",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

# Training + test data
plot(
  sample$speed,
  sample$dist,
  main = "Train + Test",
  xlab = "Speed (mph)",
  ylab = "Stopping Distance (ft)",
  col = "forestgreen",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

# Overlay test data
points(
  test$speed,
  test$dist,
  col = "black",
  pch = 19
)

par(mfrow = c(1, 1))

Resource, Suggested Reading : Ch 5 (Resampling Methods), Sec 1 (Cross-Validation)

Modeling Framework

Candidate Models :

High Variance, Low Bias – Overfit Example

Overly Complex : High Degree Polynomial Model

\[ \hat{y} = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \cdots + \beta_p x^p \]

\[ \hat{y} = \beta_0 + \sum_{j=1}^{p} \beta_j x^j \]

mdl_poly <- lm(
  dist ~ poly(speed, 14),
  data = sample
)
summary_mdl_poly <- mdl_poly |> summary()

High Bias, Low Variance – Underfit Example

\[ \hat{y}=mx+b \]

Overly Simple : Linear Model

lm_mdl <- lm(dist~speed,data=sample)
summary_mdl_lm <- lm_mdl |> summary()

The shit – Physics Example

Literature Informed : Quadratic Model

mdl_quadratic <- lm(
  dist ~ poly(speed, 2),
  data = sample
)
summary_mdl_quadratic <- mdl_quadratic |> summary()

# Create a smooth grid of speed values
speed_grid <- seq(
  min(cars$speed),
  max(cars$speed),
  length.out = 500
)

# Predictions from each model
pred_linear <- predict(
  lm_mdl,
  newdata = data.frame(speed = speed_grid)
)

pred_quadratic <- predict(
  mdl_quadratic,
  newdata = data.frame(speed = speed_grid)
)

pred_poly <- predict(
  mdl_poly,
  newdata = data.frame(speed = speed_grid)
)

Training data performance

# Plot training data
plot(
  sample$speed,
  sample$dist,
  main = "Training",
  xlab = "Speed (mph)",
  ylab = "Stopping Distance (ft)",
  col = "forestgreen",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

# Linear model
lines(
  speed_grid,
  pred_linear,
  col = "red",
  lwd = 2
)

# Quadratic model
lines(
  speed_grid,
  pred_quadratic,
  col = "blue",
  lwd = 2
)

# High-degree polynomial
lines(
  speed_grid,
  pred_poly,
  col = "purple",
  lwd = 2
)

# Legend
legend(
  "topleft",
  legend = c(
    "Training Data",
    "Linear Model",
    "Quadratic Model",
    "Degree-14 Polynomial"
  ),
  col = c(
    "forestgreen",
    "red",
    "blue",
    "purple"
  ),
  pch = c(19, NA, NA, NA),
  lty = c(NA, 1, 1, 1),
  lwd = c(NA, 2, 2, 2),
  cex = 0.7
)

Bias & Var. Tradeoff

par(mfrow = c(1, 3))

# Linear model
plot(
  test$speed,
  test$dist,
  main = "Bias Model",
  xlab = "Speed",
  ylab = "Stopping Distance",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

lines(
  test$speed[order(test$speed)],
  predict(
    lm_mdl,
    newdata = test
  )[order(test$speed)],
  col = "red",
  lwd = 2
)

# Quadratic model
plot(
  test$speed,
  test$dist,
  main = "Physics Model",
  xlab = "Speed",
  ylab = "Stopping Distance",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

lines(
  test$speed[order(test$speed)],
  predict(
    mdl_quadratic,
    newdata = test
  )[order(test$speed)],
  col = "blue",
  lwd = 2
)

# Degree-14 polynomial
plot(
  test$speed,
  test$dist,
  main = "High Variance",
  xlab = "Speed",
  ylab = "Stopping Distance",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

lines(
  test$speed[order(test$speed)],
  predict(
    mdl_poly,
    newdata = test
  )[order(test$speed)],
  col = "purple",
  lwd = 2
)

par(mfrow = c(1, 1))

Model Adequacy

\[ \text{MSE}(\hat{\theta})=\mathbb{E}(\hat{\theta}-\theta)^2 \]

\(\theta\) = parameter – ie. what u want to est.
\(\hat{\theta}=\) estimate

Decomposition

\[ \text{MSE}(\hat{\theta})=B^2_{\hat{\theta}}+V_{\hat{\theta}} + \epsilon \]

\(\epsilon\) natural error

\[ B_{\hat{\theta}}=\mathbb{E}(\hat{\theta})-\theta \]

Distance from True Value

Interpretation

\[ \text{RMSE}(\hat{\theta})=\sqrt{\text{MSE}(\hat{\theta})} \]

Average deviation of estimate from the mean

# Sample performance
sample_performance <- data.frame(
  Model = c(
    "High Bias",
    "Physics",
    "High Variance"
  ),
  MSE = c(
    mean(residuals(lm_mdl)^2),
    mean(residuals(mdl_quadratic)^2),
    mean(residuals(mdl_poly)^2)
  ),
  RMSE = c(
    sqrt(mean(residuals(lm_mdl)^2)),
    sqrt(mean(residuals(mdl_quadratic)^2)),
    sqrt(mean(residuals(mdl_poly)^2))
  )
)

# Test performance
test_performance <- data.frame(
  Model = c(
    "High Bias",
    "Physics",
    "High Variance"
  ),
  MSE = c(
    mean((test$dist - predict(lm_mdl, newdata = test))^2),
    mean((test$dist - predict(mdl_quadratic, newdata = test))^2),
    mean((test$dist - predict(mdl_poly, newdata = test))^2)
  ),
  RMSE = c(
    sqrt(mean((test$dist - predict(lm_mdl, newdata = test))^2)),
    sqrt(mean((test$dist - predict(mdl_quadratic, newdata = test))^2)),
    sqrt(mean((test$dist - predict(mdl_poly, newdata = test))^2))
  )
)

Sample Performance

knitr::kable(
  sample_performance,
  digits = 2,
  align = c("l", "r", "r")
)

Model	MSE	RMSE
High Bias	272.43	16.51
Physics	265.13	16.28
High Variance	164.24	12.82

Test Performance

knitr::kable(
  test_performance,
  digits = 2,
  align = c("l", "r", "r")
)

Model	MSE	RMSE
High Bias	47.58	6.90
Physics	24.37	4.94
High Variance	519660641.27	22796.07

# Plot training data
plot(
  sample$speed,
  sample$dist,
  main = "Training",
  xlab = "Speed (mph)",
  ylab = "Stopping Distance (ft)",
  col = "forestgreen",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

# Linear model
lines(
  speed_grid,
  pred_linear,
  col = "red",
  lwd = 2
)

# Quadratic model
lines(
  speed_grid,
  pred_quadratic,
  col = "blue",
  lwd = 2
)

# High-degree polynomial
lines(
  speed_grid,
  pred_poly,
  col = "purple",
  lwd = 2
)

# Legend
legend(
  "topleft",
  legend = c(
    "Training Data",
    "Linear Model",
    "Quadratic Model",
    "Degree-14 Polynomial"
  ),
  col = c(
    "forestgreen",
    "red",
    "blue",
    "purple"
  ),
  pch = c(19, NA, NA, NA),
  lty = c(NA, 1, 1, 1),
  lwd = c(NA, 2, 2, 2),
  cex = 0.7
)

par(mfrow = c(1, 3))

# Linear model
plot(
  test$speed,
  test$dist,
  main = "Bias Model",
  xlab = "Speed",
  ylab = "Stopping Distance",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

lines(
  test$speed[order(test$speed)],
  predict(
    lm_mdl,
    newdata = test
  )[order(test$speed)],
  col = "red",
  lwd = 2
)

# Quadratic model
plot(
  test$speed,
  test$dist,
  main = "Physics Model",
  xlab = "Speed",
  ylab = "Stopping Distance",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

lines(
  test$speed[order(test$speed)],
  predict(
    mdl_quadratic,
    newdata = test
  )[order(test$speed)],
  col = "blue",
  lwd = 2
)

# Degree-14 polynomial
plot(
  test$speed,
  test$dist,
  main = "High Variance",
  xlab = "Speed",
  ylab = "Stopping Distance",
  col = "black",
  pch = 19,
  xlim = x_limits,
  ylim = y_limits
)

lines(
  test$speed[order(test$speed)],
  predict(
    mdl_poly,
    newdata = test
  )[order(test$speed)],
  col = "purple",
  lwd = 2
)

par(mfrow = c(1, 1))

Model Comparison

ls()

##  [1] "idx"                   "lm_mdl"                "mdl_poly"             
##  [4] "mdl_quadratic"         "pred_linear"           "pred_poly"            
##  [7] "pred_quadratic"        "sample"                "sample_performance"   
## [10] "speed_grid"            "summary_mdl_lm"        "summary_mdl_poly"     
## [13] "summary_mdl_quadratic" "test"                  "test_performance"     
## [16] "x_limits"              "y_limits"

mdl_quadratic

## 
## Call:
## lm(formula = dist ~ poly(speed, 2), data = sample)
## 
## Coefficients:
##     (Intercept)  poly(speed, 2)1  poly(speed, 2)2  
##           44.60           127.62            17.08

mdl_poly

## 
## Call:
## lm(formula = dist ~ poly(speed, 14), data = sample)
## 
## Coefficients:
##       (Intercept)   poly(speed, 14)1   poly(speed, 14)2   poly(speed, 14)3  
##           44.6000           127.6185            17.0832            19.0524  
##  poly(speed, 14)4   poly(speed, 14)5   poly(speed, 14)6   poly(speed, 14)7  
##           24.1107            23.1418            14.6597            13.3240  
##  poly(speed, 14)8   poly(speed, 14)9  poly(speed, 14)10  poly(speed, 14)11  
##           -2.0001           -19.1855           -12.9484            38.8046  
## poly(speed, 14)12  poly(speed, 14)13  poly(speed, 14)14  
##           10.4622             0.9072             2.7378

anova(lm_mdl, mdl_quadratic, mdl_poly)

## Analysis of Variance Table
## 
## Model 1: dist ~ speed
## Model 2: dist ~ poly(speed, 2)
## Model 3: dist ~ poly(speed, 14)
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)
## 1     38 10897.1                           
## 2     37 10605.3  1     291.8 1.1105 0.3020
## 3     25  6569.7 12    4035.6 1.2798 0.2895

Data Dist

cars$speed|>hist()

cars$dist|>hist()

Bias-Variance Tradeoff

Isaiah C. Mireles