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

0:55 to 1:12

1 Real World Data

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

1.1 Research Question :

Ie how fucked are pedestrians if ur speeding?

1.2 Modeling Philosophy

# data(cars) 
?cars # investigate data set

2 Sampling Framework

2.1 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)

3 Modeling Framework

3.1 Candidate Models :

3.1.1 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()

3.1.2 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()

3.1.3 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)
)

4 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
)

5 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))

6 Model Adequacy

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

  • \(\theta\) = parameter – ie. what u want to est.

  • \(\hat{\theta}=\) estimate

6.1 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

7 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))
  )
)

7.0.1 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

7.0.2 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))