#The observations are randomly divided into k groups (folds) of roughly equal size. One fold is held out as the validation set, and the model is fit on the remaining k−1 folds. The mean squared error (or other error measure) is computed on the held-out fold. This process is repeated k times, with each fold serving as the validation set exactly once. The k resulting error estimates are then averaged to produce the k-fold CV estimate: CV(k)=(1/k)*∑MSE wher i=1 to k . #advantages of k-fold CV: Lower variance in the test error estimate, since it doesn’t depend heavily on a single random train/validation split. Also less bias, because more data is used for training in each fit (only 1/k of the data is held out at a time, vs. often ~half in the validation set approach) so the fitted models more closely resemble one fit on the full dataset. #Disadvantage: More computationally expensive, since the model must be fit k times instead of once.
##ii. LOOCV? #Advantages of k-fold CV (with k < n): Computationally cheaper: fits the model k times forexample 5 or 10 instead of n times. Lower variance: LOOCV’s n fitted models are trained on nearly identical data (each missing just one point), so their outputs are highly correlated, and averaging highly correlated quantities doesn’t reduce variance much. k-fold CV’s training sets overlap less, giving less correlated error estimates and thus lower variance in the final CV estimate.
#Disadvantage of k-fold CV: Slightly higher bias — LOOCV uses n−1 observations for training each time (almost the full dataset), giving nearly unbiased estimates of test error. k-fold CV uses less data per fit (n×(k−1)/k), so it has slightly more bias. This reflects a classic bias-variance tradeoff: LOOCV has lower bias but higher variance; k-fold CV (typically k=5 or 10) has slightly more bias but substantially lower variance, and is generally preferred in practice for this reason plus computational cost.
library(ISLR2)
data(Default)
head(Default)## default student balance income
## 1 No No 729.5265 44361.625
## 2 No Yes 817.1804 12106.135
## 3 No No 1073.5492 31767.139
## 4 No No 529.2506 35704.494
## 5 No No 785.6559 38463.496
## 6 No Yes 919.5885 7491.559dim(Default)## [1] 10000 4summary(Default)## default student balance income
## No :9667 No :7056 Min. : 0.0 Min. : 772
## Yes: 333 Yes:2944 1st Qu.: 481.7 1st Qu.:21340
## Median : 823.6 Median :34553
## Mean : 835.4 Mean :33517
## 3rd Qu.:1166.3 3rd Qu.:43808
## Max. :2654.3 Max. :73554set.seed(1)
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: 8set.seed(1)
train <- sample(nrow(Default), nrow(Default) / 2)glm.fit <- glm(default ~ income + balance, data = Default,
family = "binomial", subset = train)
summary(glm.fit)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.194e+01 6.178e-01 -19.333 < 2e-16 ***
## income 3.262e-05 7.024e-06 4.644 3.41e-06 ***
## balance 5.689e-03 3.158e-04 18.014 < 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: 1523.8 on 4999 degrees of freedom
## Residual deviance: 803.3 on 4997 degrees of freedom
## AIC: 809.3
##
## Number of Fisher Scoring iterations: 8glm.probs <- predict(glm.fit, newdata = Default[-train, ], type = "response")
glm.pred <- rep("No", nrow(Default) - length(train))
glm.pred[glm.probs > 0.5] <- "Yes"
table(glm.pred, Default[-train, ]$default)##
## glm.pred No Yes
## No 4824 108
## Yes 19 49mean(glm.pred != Default[-train, ]$default)## [1] 0.0254#validation set error rate is 2.54% (127 misclassified out of 5000) with this split
# Split 1
set.seed(2)
train <- sample(nrow(Default), nrow(Default) / 2)
glm.fit <- glm(default ~ income + balance, data = Default,
family = "binomial", subset = train)
glm.probs <- predict(glm.fit, newdata = Default[-train, ], type = "response")
glm.pred <- rep("No", nrow(Default) - length(train))
glm.pred[glm.probs > 0.5] <- "Yes"
mean(glm.pred != Default[-train, ]$default)## [1] 0.0238# Split 2
set.seed(3)
train <- sample(nrow(Default), nrow(Default) / 2)
glm.fit <- glm(default ~ income + balance, data = Default,
family = "binomial", subset = train)
glm.probs <- predict(glm.fit, newdata = Default[-train, ], type = "response")
glm.pred <- rep("No", nrow(Default) - length(train))
glm.pred[glm.probs > 0.5] <- "Yes"
mean(glm.pred != Default[-train, ]$default)## [1] 0.0264# Split 3
set.seed(4)
train <- sample(nrow(Default), nrow(Default) / 2)
glm.fit <- glm(default ~ income + balance, data = Default,
family = "binomial", subset = train)
glm.probs <- predict(glm.fit, newdata = Default[-train, ], type = "response")
glm.pred <- rep("No", nrow(Default) - length(train))
glm.pred[glm.probs > 0.5] <- "Yes"
mean(glm.pred != Default[-train, ]$default)## [1] 0.0256#The validation set error rate varies somewhat depending on which observations are included in the training vs. validation set — ranging from about 2.4% to 2.6%. This illustrates a key drawback of the validation set approach: the test error estimate has non-trivial variability because it depends heavily on the particular random split. There’s no single “correct” validation error — it’s an estimate with some noise, and this variability is exactly the kind of thing k-fold CV is designed to reduce (as we discussed in Question 3), since it averages over multiple splits rather than relying on just one
set.seed(1)
train <- sample(nrow(Default), nrow(Default) / 2)
glm.fit <- glm(default ~ income + balance + student, data = Default,
family = "binomial", subset = train)
summary(glm.fit)##
## Call:
## glm(formula = default ~ income + balance + student, family = "binomial",
## data = Default, subset = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.134e+01 6.937e-01 -16.346 <2e-16 ***
## income 1.686e-05 1.122e-05 1.502 0.1331
## balance 5.767e-03 3.213e-04 17.947 <2e-16 ***
## studentYes -5.992e-01 3.324e-01 -1.803 0.0715 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1523.77 on 4999 degrees of freedom
## Residual deviance: 800.07 on 4996 degrees of freedom
## AIC: 808.07
##
## Number of Fisher Scoring iterations: 8glm.probs <- predict(glm.fit, newdata = Default[-train, ], type = "response")
glm.pred <- rep("No", nrow(Default) - length(train))
glm.pred[glm.probs > 0.5] <- "Yes"
table(glm.pred, Default[-train, ]$default)##
## glm.pred No Yes
## No 4825 112
## Yes 18 45mean(glm.pred != Default[-train, ]$default)## [1] 0.026#Adding the student dummy variable does not reduce the test error rate — if anything, it’s very slightly higher (2.60% vs. 2.54%). This is consistent with what we saw in Chapter 4: student is a useful predictor on its own (since students tend to carry higher balances and thus higher default risk), but once balance is already in the model, student adds little additional predictive power — in fact its coefficient here isn’t even significant at the 5% level (p = 0.0715). This matches the textbook’s known finding: including student doesn’t meaningfully improve the model’s ability to classify default status once income and balance are accounted for.
#QUESTION 6:
set.seed(1)
glm.fit <- glm(default ~ income + balance, data = Default, family = "binomial")
summary(glm.fit)$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-136library(boot)
boot.fn <- function(data, index) {
fit <- glm(default ~ income + balance, data = data,
family = "binomial", subset = index)
return(coef(fit))
}
# test it on the full sample
boot.fn(Default, 1:nrow(Default))## (Intercept) income balance
## -1.154047e+01 2.080898e-05 5.647103e-03set.seed(1)
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* -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#The bootstrap standard errors are remarkably close to the standard formula-based (asymptotic) standard errors from glm() differing only in the third or fourth significant digit for all three coefficients.This close agreement makes sense: the formula-based SEs from glm() rely on asymptotic theory (maximum likelihood estimation), which assumes the model is correctly specified and relies on large-sample approximations (e.g., the Hessian/Fisher information matrix) to estimate variability. With n = 10,000, we’re well within the regime where those asymptotic approximations should be quite accurate which is exactly what we observe. The bootstrap, by contrast, makes no distributional assumptions it estimates variability directly by repeatedly resampling the data and refitting the model, capturing the empirical sampling distribution of the coefficients. The fact that the two approaches agree closely here gives us added confidence that the standard logistic regression assumptions (and the asymptotic SE formula) are reasonable for this dataset and model.
##question 9
library(ISLR2)
data(Boston)
mu.hat <- mean(Boston$medv)
mu.hat## [1] 22.53281#This is the sample mean of medv (median home value, in $1000s) across the 506 census tracts in the Boston dataset.
se.mu.hat <- sd(Boston$medv) / sqrt(nrow(Boston))
se.mu.hat## [1] 0.4088611library(boot)
set.seed(1)
boot.fn <- function(data, index) {
mean(data[index])
}
boot(Boston$medv, boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.007650791 0.4106622##Bootstrap SE(μ̂) = 0.4106622. Comparison to (b): This is extremely close to the formula-based SE of 0.4088611 — differing by less than 0.002. This is a great illustration of the bootstrap “working” as expected: it doesn’t rely on any parametric formula, yet it recovers almost exactly the same standard error as the classical formula for the standard error of the mean. This makes sense because the formula SE/√n is itself an asymptotically valid estimator, and with n=506 we’re in a regime where it’s quite accurate so an assumption-free resampling method arrives at essentially the same answer.
CI.lower <- mu.hat - 2 * 0.4106622
CI.upper <- mu.hat + 2 * 0.4106622
c(CI.lower, CI.upper)## [1] 21.71148 23.35413# Compare to R's built-in 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 two confidence intervals are nearly identical — both centered at 22.53 with a half-width of roughly 0.8. The bootstrap-based interval is very slightly wider than the t-test interval, which makes sense given the bootstrap SE (0.4107) was marginally larger than the formula SE (0.4089). This again confirms that the bootstrap, despite making no distributional assumptions, produces essentially the same inferential conclusions as the classical parametric approach herereassuring us that the normal-theory approximation is quite reasonable for this dataset
med.hat <- median(Boston$medv)
med.hat## [1] 21.2set.seed(1)
boot.fn <- function(data, index) {
median(data[index])
}
boot(Boston$medv, boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 0.02295 0.3778075#Since there’s no simple closed-form formula for the standard error of the median (unlike the mean), the bootstrap is especially valuable here — it gives us an estimate we couldn’t easily get otherwise. The SE of the median (0.378) is somewhat smaller than the SE of the mean (0.409/0.411), suggesting the median is estimated slightly more precisely than the mean in this sample which is plausible given the right-skew in medv, where the mean is more sensitive to extreme high values while the median is more robust to them
mu.hat.0.1 <- quantile(Boston$medv, 0.1)
mu.hat.0.1## 10%
## 12.75#μ̂₀.₁ = 12.75 estimates that 10% of Boston census tracts in this dataset have a median home value at or below $12,750.
set.seed(1)
boot.fn <- function(data, index) {
quantile(data[index], 0.1)
}
boot(Boston$medv, boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.0339 0.4767526