Assignment4

Problem 3

(a) Explain how k-fold cross-validation is implemented.

• K-fold cross-validation works by randomly dividing the full data set into k roughly equal-sized folds. The model is then fit k times, each time using k-1 folds as the training data and the remaining fold as the validation set. The test error is estimated by averaging the k validation errors together.

(b) What are the advantages and disadvantages of k-fold cross-validation relative to:

i. The validation set approach?

• K-fold CV has two big advantages over the validation set approach: it makes use of more data for training (since most observations are used for fitting in each fold), and it produces a less variable test error estimate because the result is averaged over k splits rather than a single one. The downside is that it takes longer to compute since the model must be fir k times instead of once.

ii. LOOCV?

• Compared to LOOCV, k-fold CV is much faster to compute because it fits the model one k times instead of once for every observation. It may have slightly more bias than LOOKV, but it often produces a more stable estimate of the test error.

---

Problem 5

(a) Fit a logistic regression model that uses income and balance to predict default.

library(ISLR2)
set.seed(42)

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

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) Using the validation set approach, estimate the test error of this model.

set.seed(42)

# i. Split into training and validation sets
n <- nrow(Default)
train_idx <- sample(n, n / 2)
default_train <- Default[train_idx, ]
default_val   <- Default[-train_idx, ]

# ii. Fit logistic regression on training set only
glm_train <- glm(default ~ income + balance, data = default_train, family = binomial)

# iii. Predict default status on validation set
val_probs <- predict(glm_train, default_val, type = "response")
val_pred  <- ifelse(val_probs > 0.5, "Yes", "No")

# iv. Compute validation set error
mean(val_pred != default_val$default)

[1] 0.026

• The validation set error is approximately 2.6%, meaning the logistic regression model incorrectly classifies about 2.6% of the observations in the validation set.

(c) Repeat the process in (b) three times, using three different splits of the observations into a training set and a validation set. Comment on the results obtained.

set.seed(1)
n <- nrow(Default)
train_idx <- sample(n, n / 2)
default_train <- Default[train_idx, ]
default_val   <- Default[-train_idx, ]

# ii. Fit logistic regression on training set only
glm_train <- glm(default ~ income + balance, data = default_train, family = binomial)

# iii. Predict default status on validation set
val_probs <- predict(glm_train, default_val, type = "response")
val_pred  <- ifelse(val_probs > 0.5, "Yes", "No")

mean(val_pred != default_val$default)

[1] 0.0254

set.seed(2)
n <- nrow(Default)
train_idx <- sample(n, n / 2)
default_train <- Default[train_idx, ]
default_val   <- Default[-train_idx, ]

# ii. Fit logistic regression on training set only
glm_train <- glm(default ~ income + balance, data = default_train, family = binomial)

# iii. Predict default status on validation set
val_probs <- predict(glm_train, default_val, type = "response")
val_pred  <- ifelse(val_probs > 0.5, "Yes", "No")

# iv. Compute validation set error
mean(val_pred != default_val$default)

[1] 0.0238

set.seed(3)
n <- nrow(Default)
train_idx <- sample(n, n / 2)
default_train <- Default[train_idx, ]
default_val   <- Default[-train_idx, ]

# ii. Fit logistic regression on training set only
glm_train <- glm(default ~ income + balance, data = default_train, family = binomial)

# iii. Predict default status on validation set
val_probs <- predict(glm_train, default_val, type = "response")
val_pred  <- ifelse(val_probs > 0.5, "Yes", "No")

# iv. Compute validation set error
mean(val_pred != default_val$default)

[1] 0.0264

• The validation error ranged from about 2.3% to 2.8% across the three splits. This shows that the validation set approach can produce slightly different test error estimates depending on how the data are split.

(d) Now consider a logistic regression model that predicts the probability of default using income, balance, and a dummy variable for student. Estimate the test error for this model using the validation set approach. Comment on whether or not including a dummy variable for student leads to a reduction in the test error rate.

set.seed(42)
train_idx_d <- sample(n, n / 2)
glm_student <- glm(default ~ income + balance + student,
                   data = Default[train_idx_d, ], family = binomial)
probs_student <- predict(glm_student, Default[-train_idx_d, ], type = "response")
pred_student  <- ifelse(probs_student > 0.5, "Yes", "No")
mean(pred_student != Default[-train_idx_d, ]$default)

[1] 0.0258

• Adding the student dummy variable does not reduce the test error by much, it remains around 2.6%, essentially the same as the model without student. This suggests that student status does not add significant predictive power for default.

---

Problem 6

(a) Using the summary() and glm() functions, determine the estimated standard errors for the coefficients associated with income and balance in a multiple logistic regression model that uses both predictors.

set.seed(42)
glm_se <- glm(default ~ income + balance, data = Default, family = binomial)
summary(glm_se)$coefficients

                 Estimate   Std. Error    z value      Pr(>|z|)
(Intercept) -1.154047e+01 4.347564e-01 -26.544680 2.958355e-155
income       2.080898e-05 4.985167e-06   4.174178  2.990638e-05
balance      5.647103e-03 2.273731e-04  24.836280 3.638120e-136

(b) Write a function, boot.fn(), that takes as input the Default data set as well as an index of the observations, and that outputs the coefficient estimates for income and balance in the multiple logistic regression model.

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

# Test on full data
boot.fn(Default, 1:nrow(Default))

      income      balance 
2.080898e-05 5.647103e-03 

(c) Use the boot() function together with your boot.fn() function to estimate the standard errors of the logistic regression coefficients for income and balance.

library(boot)
set.seed(42)
boot(Default, boot.fn, R = 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


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


Bootstrap Statistics :
        original       bias     std. error
t1* 2.080898e-05 2.737444e-08 5.073444e-06
t2* 5.647103e-03 1.176249e-05 2.299133e-04

(d) Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.

• The standard errors from the bootstrap method are very similar to those from glm(). This suggests that the estimates are consistent, and either method provides similar results for this data set.

---

Problem 9

(a) Based on this data set, provide an estimate for the population mean of medv. Call this estimate µ̂.

library(ISLR2)
mu_hat <- mean(Boston$medv)
print(mu_hat)

[1] 22.53281

• The estimated population mean of medv is approximately 22.53, meaning the average median home value across Boston is about $22,530.

(b) Provide an estimate of the standard error of µ̂. Interpret this result.

se_hat <- sd(Boston$medv) / sqrt(nrow(Boston))
print(se_hat)

[1] 0.4088611

• The estimated standard error of µ̂ is approximately 0.409. This means that the sample mean of medv would typically vary by about.409 from the true population mean if we repeatedly took random samples from the population.

(c) Now estimate the standard error of µ̂ using the bootstrap. How does this compare to your answer from (b)?

library(boot)

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

set.seed(42)
boot_result <- boot(Boston, mean.fn, R = 1000)
boot_result

ORDINARY NONPARAMETRIC BOOTSTRAP


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


Bootstrap Statistics :
    original     bias    std. error
t1* 22.53281 0.02671186   0.4009216

• The bootstrap standard error is very close to the estimate from part (b), with both values being about 0.41. This suggests that both methods give similar estimates of the standard error.

(d) Based on your bootstrap estimate from (c), provide a 95% confidence interval for the mean of medv. Compare it to the results obtained using t.test(Boston$medv).

ci_lower <- mu_hat - 2 * sd(boot_result$t)
ci_upper <- mu_hat + 2 * sd(boot_result$t)

# t.test CI
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 for the mean of medv is approximately [21.73, 23.34], which is very similar to the interval produced by t.test(). This suggests that both methods give similar estimates for the true population mean.

(e) Based on this data set, provide an estimate, µ̂_med, for the median value of medv in the population.

mu_med <- median(Boston$medv)
print(mu_med)

[1] 21.2

• The estimated population median of medv is 21.2, meaning half of Boston has a median home value below $21,200.

(f) We now would like to estimate the standard error of µ̂_med. Estimate the standard error of the median using the bootstrap. Comment on your findings.

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

set.seed(42)
boot_median <- boot(Boston, median.fn, R = 1000)
boot_median

ORDINARY NONPARAMETRIC BOOTSTRAP


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


Bootstrap Statistics :
    original  bias    std. error
t1*     21.2  0.0106   0.3661785

• The bootstrap standard error of the median is approximately 0.38, which is similar to the standard error of the mean. This suggests that both the sample mean and sample median are fairly stable estimates for this data set.

(g) Based on this data set, provide an estimate for the tenth percentile of medv in Boston census tracts. Call this quantity µ̂_0.1.

mu_10 <- quantile(Boston$medv, 0.10)
print(mu_10)

  10% 
12.75 

• The estimated tenth percentile of medv is 12.75, meaning that approx 10% of census tracts have a median home value at or below $12,750.

(h) Use the bootstrap to estimate the standard error of µ̂_0.1. Comment on your findings.

percentile10.fn <- function(data, index) {
  quantile(data$medv[index], 0.10)
}

set.seed(1)
boot_p10 <- boot(Boston, percentile10.fn, R = 1000)
boot_p10

ORDINARY NONPARAMETRIC BOOTSTRAP


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


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

• The bootstrap estimated the standard error of th e10th percentile to be approx .50. This is slightly larger than the standard errors for the mean and median, suggesting that the 10th percentile is more variable from sample to sample.