Instead of exposing models to a single testing target, k-fold cross-validation maps out a systematically robust evaluation architecture:
Default
DatasetWe 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