Question 3: \(k\)-Fold Cross-Validation

(a) Implementation of \(k\)-Fold Cross-Validation

Instead of exposing models to a single testing target, k-fold cross-validation maps out a systematically robust evaluation architecture:

  1. Partitioning: The complete dataset is randomly divided into k equal, mutually exclusive, and exhaustive subsets (folds).
  2. Iterative Estimation: The model is iteratively estimated k times. In each iteration i (where i = 1, 2), the i-th fold is strictly sequestered as the validation set, while the remaining k-1 folds are amalgamated to serve as the training space.
  3. Error Aggregation: An out-of-sample performance metric (such as Mean Squared Error for regression or Misclassification Rate for classification), denoted as ERR, is computed on the validation fold.
  4. Final Estimation: The overall k-fold cross-validation estimate of the test error is calculated by averaging these distinct out-of-sample error estimates

(b) Methodological Trade-offs

  1. Relative to the Validation Set Approach
  1. Relative to Leave-One-Out Cross-Validation (LOOCV)

Question 5: Validation Set Approach on the Default Dataset

(a) Baseline Logistic Regression Model

We initialize our standard parametric model on the full space using `glm()` with a binomial link function:
library(ISLR2)
set.seed(42) # Guarantee deterministic results across environments

fit_logistic <- glm(default ~ income + balance, data = Default, family = binomial)
summary(fit_logistic)
## 
## Call:
## glm(formula = default ~ income + balance, family = binomial, 
##     data = Default)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
## income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
## balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1579.0  on 9997  degrees of freedom
## AIC: 1585
## 
## Number of Fisher Scoring iterations: 8
### ### (b) Validation Set Assessment Workflow
 
 library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.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
train_idx <- sample(nrow(Default), nrow(Default) * 0.5)
train_set <- Default[train_idx, ]
val_set   <- Default[-train_idx, ]

fit_val_split <- glm(default ~ income + balance, data = train_set, family = binomial)

post_probs <- predict(fit_val_split, newdata = val_set, type = "response")
predictions <- ifelse(post_probs > 0.5, "Yes", "No")

val_error_rate <- mean(predictions != val_set$default)
print(paste("Validation Set Error Rate:", round(val_error_rate * 100, 2), "%"))
## [1] "Validation Set Error Rate: 2.6 %"
### ### (c) Evaluation of Split Stochasticity
set.seed(101)
errors <- numeric(3)

for(i in 1:3) {
  train_idx <- sample(nrow(Default), nrow(Default) * 0.5)
  t_set <- Default[train_idx, ]
  v_set <- Default[-train_idx, ]
  
  fit <- glm(default ~ income + balance, data = t_set, family = binomial)
  probs <- predict(fit, newdata = v_set, type = "response")
  preds <- ifelse(probs > 0.5, "Yes", "No")
  errors[i] <- mean(preds != v_set$default)
}
print(errors)
## [1] 0.0250 0.0268 0.0280
## ## Question 6: Parameter Error Quantification via Bootstrap

### ### (a) Asymptotic Standard Errors via glm()

fit_full <- glm(default ~ income + balance, data = Default, family = binomial)
summary(fit_full)$coefficients[, 2]
##  (Intercept)       income      balance 
## 4.347564e-01 4.985167e-06 2.273731e-04
### ### (b) Defining the Bootstrapping Objective Vector

boot.fn <- function(data, index) {
  fit <- glm(default ~ income + balance, data = data, subset = index, family = binomial)
  return(coef(fit))
}

### ### (c) Executing the Empirical Resampling Array

library(boot)
set.seed(1)
boot_results <- boot(data = Default, statistic = boot.fn, R = 1000)
print(boot_results)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Default, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##          original        bias     std. error
## t1* -1.154047e+01 -3.945460e-02 4.344722e-01
## t2*  2.080898e-05  1.680317e-07 4.866284e-06
## t3*  5.647103e-03  1.855765e-05 2.298949e-04
## ## Question 9: Non-Parametric Resampling on the Boston Housing Market

### ### (a) Population Mean Estimation

mu_hat <- mean(Boston$medv)
print(mu_hat)
## [1] 22.53281
### ### (b) Classical Asymptotic Standard Error Estimation

se_analytical <- sd(Boston$medv) / sqrt(nrow(Boston))
print(se_analytical)
## [1] 0.4088611
### ### (c) Bootstrap Variance Assessment

boot_mean_fn <- function(data, index) {
  return(mean(data$medv[index]))
}

set.seed(1)
boot_mean_res <- boot(data = Boston, statistic = boot_mean_fn, R = 1000)
print(boot_mean_res)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston, statistic = boot_mean_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 22.53281 0.007650791   0.4106622
### (d) Inferential Confidence Intervals

se_boot <- 0.4106
ci_lower <- mu_hat - 2 * se_boot
ci_upper <- mu_hat + 2 * se_boot
print(paste("Bootstrap 95% CI: [", round(ci_lower, 4), ",", round(ci_upper, 4), "]"))
## [1] "Bootstrap 95% CI: [ 21.7116 , 23.354 ]"
# Student's t-test baseline for comparison
t.test(Boston$medv)
## 
##  One Sample t-test
## 
## data:  Boston$medv
## t = 55.111, df = 505, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  21.72953 23.33608
## sample estimates:
## mean of x 
##  22.53281
### (e) Non-Parametric Central Tendency 

median_hat <- median(Boston$medv)
print(median_hat)
## [1] 21.2
### (f) Standard Error of the Median via Bootstrap
boot_median_fn <- function(data, index) {
  return(median(data$medv[index]))
}

set.seed(1)
boot_median_res <- boot(data = Boston, statistic = boot_median_fn, R = 1000)
print(boot_median_res)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston, statistic = boot_median_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     21.2 0.02295   0.3778075
### (g) Tenth Percentile Estimation
quantile_tenth <- quantile(Boston$medv, 0.1)
print(quantile_tenth)
##   10% 
## 12.75
### (h) Standard Error of the Tenth Percentile via Bootstrap
boot_quantile_fn <- function(data, index) {
  return(quantile(data$medv[index], 0.1))
}

set.seed(1)
boot_quantile_res <- boot(data = Boston, statistic = boot_quantile_fn, R = 1000)
print(boot_quantile_res)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston, statistic = boot_quantile_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    12.75  0.0339   0.4767526