Chris Serrano - Assignment 4

Question 3

a

k-fold cross-validation is implemented by following these steps:The original dataset is randomly divided into k roughly equal-sized, non-overlapping subsets.

1. For each fold
   • Use i as the validation (or testing) set.
   • Use the remaining folds combined as the training set.
   • Fit the statistical learning method) on the training set.
   • Compute the test error
2. The overall k fold cross-validation estimate of the test error is calculated by averaging the individual test errors computed over all k iterations:

b

i.

• Advantages:
  1. Lower Variance
  2. Reduced Pessimistic Bias
• Disadvantages:
  1. Computational Cost

ii.

• Advantages:
  1. Computational Efficiency
  2. Better Bias-Variance Trade-off
• Disadvantages:
  1. Higher Bias

---

Question 5

a

library(ISLR2)

glm.fit <- glm(default ~ income + balance, data = Default, family = binomial)
summary(glm.fit)

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

We now estimate the test error of this model using the validation set approach. We will split the dataset 50/50, fit the model on the training set, predict on the validation set using a threshold of 0.5, and calculate the validation error rate.

set.seed(1)

train <- sample(nrow(Default), nrow(Default) / 2)

glm.train_fit <- glm(default ~ income + balance, data = Default, family = binomial, subset = train)

glm.probs <- predict(glm.train_fit, Default[-train, ], type = "response")
glm.pred <- ifelse(glm.probs > 0.5, "Yes", "No")

val_error <- mean(glm.pred != Default[-train, ]$default)
cat("Validation Set Error Rate:", val_error, "\n")

Validation Set Error Rate: 0.0254 

c

for (s in 2:4) {
  set.seed(s)
  train_split <- sample(nrow(Default), nrow(Default) / 2)
  
  fit <- glm(default ~ income + balance, data = Default, family = binomial, subset = train_split)
  probs <- predict(fit, Default[-train_split, ], type = "response")
  preds <- ifelse(probs > 0.5, "Yes", "No")
  
  err <- mean(preds != Default[-train_split, ]$default)
  cat("Seed:", s, "-> Validation Error Rate:", round(err * 100, 2), "%\n")
}

Seed: 2 -> Validation Error Rate: 2.38 %
Seed: 3 -> Validation Error Rate: 2.64 %
Seed: 4 -> Validation Error Rate: 2.56 %

The validation error rates fluctuate between approximately 2.44% and 2.74%. This variability shows that the validation set approach has some variance, as the estimated test error rate depends on which observations happen to end up in the training and validation sets.

d

set.seed(1)
train_d <- sample(nrow(Default), nrow(Default) / 2)

glm.student_fit <- glm(default ~ income + balance + student, data = Default, family = binomial, subset = train_d)

student.probs <- predict(glm.student_fit, Default[-train_d, ], type = "response")
student.pred <- ifelse(student.probs > 0.5, "Yes", "No")

val_error_student <- mean(student.pred != Default[-train_d, ]$default)
cat("Validation Set Error Rate (with Student):", val_error_student, "\n")

Validation Set Error Rate (with Student): 0.026 

Including the dummy variable does not lead to a reduction in the test error rate.

---

Question 6

a

glm.full <- glm(default ~ income + balance, data = Default, family = binomial)
summary(glm.full)$coefficients[, c("Estimate", "Std. Error")]

                 Estimate   Std. Error
(Intercept) -1.154047e+01 4.347564e-01
income       2.080898e-05 4.985167e-06
balance      5.647103e-03 2.273731e-04

b

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

c

library(boot)
set.seed(1)

boot_results <- boot(Default, boot.fn, R = 1000)
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

d

The bootstrap standard errors are similar to those obtained using standard logistic regression formulas.

---

Question 9

a

mu.hat <- mean(Boston$medv)
mu.hat

[1] 22.53281

b

se.hat <- sd(Boston$medv) / sqrt(nrow(Boston))
se.hat

[1] 0.4088611

c

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

set.seed(1)
boot.mean_results <- boot(Boston, boot.fn.mean, R = 1000)
boot.mean_results

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Boston, statistic = boot.fn.mean, R = 1000)


Bootstrap Statistics :
    original      bias    std. error
t1* 22.53281 0.007650791   0.4106622

The bootstrap standard error estimate is very close to the theoretical standard error.

d

boot_se <- sd(boot.mean_results$t)

ci.boot <- c(mu.hat - 2 * boot_se, mu.hat + 2 * boot_se)
cat("Bootstrap 95% CI:", ci.boot, "\n")

Bootstrap 95% CI: 21.71148 23.35413 

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 

The bootstrap 95% confidence interval is almost identical to the interval obtained using the standard formula.

e

median.hat <- median(Boston$medv)
median.hat

[1] 21.2

f

boot.fn.median <- function(data, index) {
  return(median(data$medv[index]))
}

set.seed(1)
boot.median_results <- boot(Boston, boot.fn.median, R = 1000)
boot.median_results

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Boston, statistic = boot.fn.median, R = 1000)


Bootstrap Statistics :
    original  bias    std. error
t1*     21.2 0.02295   0.3778075

Because there is no straightforward theoretical formula for the standard error of a median, the boot strap is an incredibly useful tool here. The standard error is small, indicating that the median estimate is relatively stable.

g

q10.hat <- quantile(Boston$medv, 0.1)
q10.hat

  10% 
12.75 

h

boot.fn.q10 <- function(data, index) {
  return(quantile(data$medv[index], 0.1))
}

set.seed(1)
boot.q10_results <- boot(Boston, boot.fn.q10, R = 1000)
boot.q10_results

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Boston, statistic = boot.fn.q10, R = 1000)


Bootstrap Statistics :
    original  bias    std. error
t1*    12.75  0.0339   0.4767526

Similar to the median, there is no simple analytical formula for the standard error of percentiles. The standard error xshows that the 10th percentile is estimated with a moderate degree of precision.