library(ISLR)
library(ISLR2)##
## Attaching package: 'ISLR2'## The following objects are masked from 'package:ISLR':
##
## Auto, Creditlibrary(boot)We now review k-fold cross-validation.
(a) Explain how k-fold cross-validation is implemented.
In k-fold cross-validation, the data is randomly split into k number of subsets. The model is then trained k times, each time using k - 1 subset as the training set. In other words, going down each subset and using the next one as the training set. A different subset is held out as the validation set each time. The test error is calculated on the validation subset during each iteration. Once all k iterations are complete, the k test error estimates are averaged to produce the overall cross-validation estimate of the test error. This method makes sure every observation is used for both training and validation once.
(b) What are the advantages and disadvantages of k-fold cross-validation relative to:
i. The validation set approach?
The k-fold cv approach has a more efficient use of the data by using all observations for both training and validation. The validation set approach splits itself only once. Additionally, because all observations are used and the test error is averaged over k splits, there will be a lower variance in the test error estimate. The only disadvantage of k-fold cv is greater computation from fitting it k times instead of only once. It can become very time-consuming with large datasets.
ii. LOOCV?
The k-fold cv approach is more efficient than LOOCV because it fits the model times, whereas LOOCV fits the model n times where n is the total number of observations and removes 1 each time. LOOCV also tends to have a higher variance due to the similar training sets. The models are nearly identical, making test errors fluctuate from one left-out point to another. The disadvantage of k-fold cv is that is has a slightly higher bias. This is due to to smaller training sets.
In Chapter 4, we used logistic regression to predict the probability
of default using income and balance on the Default data set.
We will now estimate the test error of this logistic regression model
using the validation set approach. Do not forget to set a random seed
before beginning your analysis.
attach(Default)(a) Fit a logistic regression model that uses
income and balance to predict
default.
set.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: 8(b) Using the validation set approach, estimate the test error of this model. In order to do this, you must perform the following steps:
i. Split the sample set into a training set and a validation set.
The observations have been split into a 70/30 for training and validation.
set.seed(1)
train <- sample(nrow(Default), nrow(Default) * 0.7)ii. Fit a multiple logistic regression model using only the training observations.
glm.fit <- glm(default ~ income + balance, data = Default, family = "binomial", subset = train)iii. Obtain a prediction of default status for each individual in the validation set by computing the posterior probability of default for that individual, and classifying the individual to the default category if the posterior probability is greater than 0.5.
The prediction of default status of each individual is then collected
to give to the mean() function for calculating the test
error.
glm.probs <- predict(glm.fit, data = Default, type = "response")
glm.pred <- ifelse(glm.probs > 0.5, "Yes", "No")iv. Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
The validation set error was 5.066667%.
mean(glm.pred[-train] != Default$default[-train])## [1] 0.05066667(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.
The validation set error for the random splits are:
The error varies slightly across the random samples. This highlights the variability in the validation set approach due to the random data partitioning.
set.seed(2)
train1 <- sample(nrow(Default), nrow(Default) * 0.7)
glm.fit1 <- glm(default ~ income + balance, data = Default, family = "binomial", subset = train1)
glm.probs1 <- predict(glm.fit1, Default, type = "response")
glm.pred1 <- ifelse(glm.probs1 > 0.5, "Yes", "No")
error1 <- mean(glm.pred1[-train1] != Default$default[-train1])
error1## [1] 0.02133333set.seed(3)
train2 <- sample(nrow(Default), nrow(Default) * 0.7)
glm.fit2 <- glm(default ~ income + balance, data = Default, family = "binomial", subset = train2)
glm.probs2 <- predict(glm.fit2, Default, type = "response")
glm.pred2 <- ifelse(glm.probs2 > 0.5, "Yes", "No")
error2 <- mean(glm.pred2[-train2] != Default$default[-train2])
error2## [1] 0.02466667set.seed(4)
train3 <- sample(nrow(Default), nrow(Default) * 0.7)
glm.fit3 <- glm(default ~ income + balance, data = Default, family = "binomial", subset = train3)
glm.probs3 <- predict(glm.fit3, Default, type = "response")
glm.pred3 <- ifelse(glm.probs3 > 0.5, "Yes", "No")
error3 <- mean(glm.pred3[-train3] != Default$default[-train3])
error3## [1] 0.02366667(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.
The validation set error for the model including student is 5%.
Compared to the the previous models, this result shows that including
student did not lead to a reduction in the rest error rate.
The previous models were around a 2% error rate.
set.seed(1)
train <- sample(nrow(Default), nrow(Default) * 0.7) glm.fit.student <- glm(default ~ income + balance + student, data = Default, family = "binomial", subset = train)
glm.probs.student <- predict(glm.fit.student, data = Default, type = "response")
glm.pred.student <- ifelse(glm.probs.student > 0.5, "Yes", "No")
error_student <- mean(glm.pred.student[-train] != Default$default[-train])
error_student## [1] 0.05We continue to consider the use of a logistic regression model to
predict the probability of default using
income and balance on the Default data set. In particular, we will now
compute estimates for the standard errors of the income and
balance logistic regression coeffcients in two different
ways: (1) using the bootstrap, and (2) using the standard formula for
computing the standard errors in the glm() function. Do not forget to set a
random seed before beginning your analysis.
(a) Using the summary() and
glm() functions, determine the
estimated standard errors for the coeffcients associated with
income and balance in a multiple logistic
regression model that uses both predictors.
set.seed(123)
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) 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 coeffcient estimates for
income and balance in the multiple logistic
regression model.
boot.fn <- function(data, index) {
fit <- glm(default ~ income + balance, data = data, subset = index, family = "binomial")
return(coef(fit)[2:3])
}(c) Use the boot() function
together with your boot.fn() function
to estimate the standard errors of the logistic regression coeffcients
for income and balance.
boot(data = Default, statistic = 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 1.582518e-07 4.729534e-06
## t2* 5.647103e-03 1.296980e-05 2.217214e-04(d) Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.
Using the glm() function, the estimated standard errors
for the logistic regression coefficients were approximately 4.985e-06
for income and 2.274e-04 for balance. Using
the bootstrap function, the error for income was 4.7295e-06
and 2.2172e-04 for balance. The bootstrap method gave a
smaller error for income and a larger error for
balance. This small difference in standard errors is likely
due to to sampling variability in the bootstrap method.
We will now consider the Boston housing data set, from the ISLR2 library.
data("Boston")
attach(Boston)(a) Based on this data set, provide an estimate for the
population mean of medv. Call this estimate \(\hat{\mu}\).
The estimate for the population mean of medv is
22.53281.
mu_hat <- mean(Boston$medv)
mu_hat## [1] 22.53281(b) Provide an estimate of the standard error of \(\hat{\mu}\). Interpret this result.
Hint: We can compute the standard error of the sample mean by dividing the sample standard deviation by the square root of the number of observations.
The standard error of \(\hat{\mu}\) is 0.4088611.
se_mu <- sd(Boston$medv) / sqrt(nrow(Boston))
se_mu## [1] 0.4088611(c) Now estimate the standard error of \(\hat{\mu}\) using the bootstrap. How does this compare to your answer from (b)?
The standard error of \(\hat{\mu}\) using the bootstrap method is 0.4106622. In part c) the error was 0.4088611. The two estimates are very similar, which suggests that the sample mean provides a stable estimate of the population mean. The bootstrap method also reported a negligible bias that further supports the accuracy of the sample mean as an estimator.
boot.fn2 <- function(data, index) {
mean(data[index])
}
set.seed(1)
boot(Boston$medv, boot.fn2, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn2, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.007650791 0.4106622(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).
Hint: You can approximate a 95 % confidence interval using the formula \([\hat{\mu} − 2SE(\hat{\mu}), \hat{\mu} + 2SE(\hat{\mu})]\).
boot_se <- 0.4106622
c(mu_hat - 2 * boot_se, mu_hat + 2 * boot_se)## [1] 21.71148 23.35413t.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.53281Using the bootstrap estimate from part c), the 95% confidence
interval for the mean of medv is approximately [21.71148,
23.35413]. Additionally, using the t.test() function gave a similar
result of [21.72953 23.33608]. The difference of hundredths from these
results suggests that the bootstrap method provides a reliable estimate
for the confidence interval of the population mean.
(e) Based on this data set, provide an estimate, \(\hat{\mu}_{med}\), for the median value of
medv in the population.
The estimated median value of medv is 21.2.
median(Boston$medv)## [1] 21.2(f) We now would like to estimate the standard error of \(\hat{\mu}_{med}\). Unfortunately, there is no simple formula for computing the standard error of the median. Instead, estimate the standard error of the median using the bootstrap. Comment on your findings.
Since there is no exact formula to calculate the standard error of the median, I used the bootstrap method to obtain an estimate of 0.3778075. This relatively low standard error suggests that the median is a stable statistic and not highly sensitive to sampling variability.
boot.fn3 <- function(data, index) {
return(median(data[index]))
}
set.seed(1)
boot(data = Boston$medv, statistic = boot.fn3, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn3, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 0.02295 0.3778075(g) Based on this data set, provide an estimate for the tenth
percentile of medv in Boston census tracts. Call this
quantity \(\hat{\mu}_{0.1}\). (You can
use the quantile()
function.)
The estimated value for the tenth percentile of medv is 12.75. In other words, 10% of the Boston census tracts have a median home value below $12.75k.
mu_hat_0.1 <- quantile(Boston$medv, 0.1)
mu_hat_0.1## 10%
## 12.75(h) Use the bootstrap to estimate the standard error of \(\hat{\mu}_{0.1}\). Comment on your findings.
Using the bootstrap method, the estimated standard error of the 10th
percentile of medv is 0.4767526. This value quantifies the
variability of the estimate \(\hat{\mu}_{0.1}
= 12.75\) across different resamples of the data. The small
standard error suggests that the 10th percentile of median home values
in the Boston data is a stable estimate, even at the lower tail of the
distribution.
boot.fn4 <- function(data, index) {
return(quantile(data[index], 0.1))
}
set.seed(1)
boot(Boston$medv, boot.fn4, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn4, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.0339 0.4767526