library(ISLR)## Warning: package 'ISLR' was built under R version 4.4.3library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.4.3##
## Attaching package: 'ISLR2'## The following objects are masked from 'package:ISLR':
##
## Auto, Creditlibrary(MASS)## Warning: package 'MASS' was built under R version 4.4.3##
## Attaching package: 'MASS'## The following object is masked from 'package:ISLR2':
##
## Bostonk-fold cross-validation is a way to check how well a model works by splitting the data into k equal parts, or “folds.” For each round, one fold is used to test the model, and the other k − 1 folds are used to train it. This process repeats k times, so each fold is used once as the test set. After all k rounds, the test errors are averaged to give a final estimate of how well the model is expected to perform on new data.
i.The validation set approach?
Compared to the validation set approach, k-fold cross-validation has some clear advantages. It gives more accurate and stable results because every data point gets used for both training and testing. In contrast, the validation set method only uses part of the data for training, which can make the model less effective and lead to more variable results. The downside is that k-fold takes longer to run because the model is trained multiple times.
Compared to leave-one-out cross-validation (LOOCV), k-fold cross-validation is faster and often gives more consistent results. LOOCV trains the model n times (where n is the number of data points), so it takes a lot of time. While LOOCV uses nearly all the data for training each time and has very low bias, it can have high variability. k-fold is a good balance between speed and accuracy, which is why people often use 5-fold or 10-fold CV in practice.
str(Default)## 'data.frame': 10000 obs. of 4 variables:
## $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ student: Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 2 1 1 ...
## $ balance: num 730 817 1074 529 786 ...
## $ income : num 44362 12106 31767 35704 38463 ...set.seed(1)
model1 <- glm(default ~ income + balance, data = Default, family = "binomial")
summary(model1)##
## 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: 8b. Using the validation set approach, estimate the test error of this model. In order to do this, you must perform the following steps:
Split the sample set into a training set and a validation set.
Fit a multiple logistic regression model using only the training observations.
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.
Compute the validation set error, which is the fraction of the observations in the validation set that are misclassiffied.
# Split the data (50% training, 50% validation)
set.seed(1)
train_index <- sample(1:nrow(Default), nrow(Default) / 2)
train_data <- Default[train_index, ]
valid_data <- Default[-train_index, ]
# Fit logistic regression on training set
model2 <- glm(default ~ income + balance, data = train_data, family = "binomial")
# Predict on validation set
probs <- predict(model2, newdata = valid_data, type = "response")
preds <- ifelse(probs > 0.5, "Yes", "No")
# Compute validation set error
actual <- valid_data$default
mean(preds != actual)## [1] 0.0254Repeat 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 test error rates from the three different validation splits are fairly close to each other, ranging from 2.5% to 2.8%. This indicates that the model performs consistently across different random training/validation splits. However, there’s still a small amount of variation depending on how the data is divided. This variability is a known limitation of the validation set approach, since the outcome depends on the specific split. For more stable results, methods like k-fold cross-validation are preferred.
run_validation <- function(seed_value) {
set.seed(seed_value)
train_index <- sample(1:nrow(Default), nrow(Default) / 2)
train_data <- Default[train_index, ]
valid_data <- Default[-train_index, ]
model <- glm(default ~ income + balance, data = train_data, family = "binomial")
probs <- predict(model, newdata = valid_data, type = "response")
preds <- ifelse(probs > 0.5, "Yes", "No")
mean(preds != valid_data$default)
}
error1 <- run_validation(1)
error2 <- run_validation(2)
error3 <- run_validation(3)
c(error1, error2, error3)## [1] 0.0254 0.0238 0.0264Now 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,
A logistic regression model was fitted using income, balance, and a dummy variable for student to predict the probability of default. Using the validation set approach with a random 50% split, the test error rate for this model was approximately 0.026, or 2.6%. This error rate is nearly identical to the rates obtained when only income and balance were used as predictors. As a result, including the student variable does not appear to meaningfully reduce the test error. This suggests that student status does not provide substantial additional predictive power for default when income and balance are already included in the model.
set.seed(4)
train_index <- sample(1:nrow(Default), nrow(Default) / 2)
train_data <- Default[train_index, ]
valid_data <- Default[-train_index, ]
model3 <- glm(default ~ income + balance + student, data = train_data, family = "binomial")
probs3 <- predict(model3, newdata = valid_data, type = "response")
preds3 <- ifelse(probs3 > 0.5, "Yes", "No")
error_with_student <- mean(preds3 != valid_data$default)
error_with_student## [1] 0.0262Using 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.
library(boot)
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: 8Write a function, boot.fn(), that takes as input the Default data set as well as an index of the observations, and that outputs the coefficients estimates for income and balance in the multiple logistic regression model.
boot.fn <- function(data, index) {
glm_fit <- glm(default ~ income + balance, data = data[index, ], family = "binomial")
return(coef(glm_fit)[2:3])
}Use the boot() function together with your boot.fn() function to estimate the standard errors of the logistic regression coefficients for income and balance.
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* 2.080898e-05 -1.698372e-07 4.718297e-06
## t2* 5.647103e-03 5.198245e-06 2.313598e-04Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.
Using the glm() function, the standard error for the
income coefficient was about 0.000008, and for balance it was about
0.0002. When we used the bootstrap method, the standard error for income
was around 0.0000047, and for balance it was about 0.00023. These
numbers are fairly close to the ones from glm(). The
bootstrap gave a slightly smaller value for income and a slightly bigger
one for balance. Overall, both methods give similar results, showing
that the standard errors from glm() are reliable in this
case.
Based on this data set, provide an estimate for the population mean of medv. Call this estimate µˆ.
library(boot)
data("Boston")
str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...mu_hat <- mean(Boston$medv)
mu_hat## [1] 22.53281Provide an estimate of the standard error of µˆ. 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 umber of observations.
se_mu_hat <- sd(Boston$medv) / sqrt(nrow(Boston))
se_mu_hat## [1] 0.4088611Now estimate the standard error of µˆ using the bootstrap. How does this compare to your answer from (b)?
The standard error of the mean was about 0.4089 using the formula and about 0.4107 using the bootstrap. The two numbers are very close, which means the formula gives a good estimate. The bootstrap just confirms it by using re-sampling instead of a formula.
boot_mean <- function(data, index) {
return(mean(data[index]))
}
set.seed(1)
boot_out <- boot(Boston$medv, boot_mean, R = 1000)
boot_se <- sd(boot_out$t)
boot_se## [1] 0.4106622Based 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 [ˆµ − 2SE(ˆµ), µˆ + 2SE(ˆµ)].
The 95% confidence interval using the bootstrap was [21.7115, 23.3541], while the interval from the t.test() was [21.7295, 23.3361]. Both intervals are very similar, which shows that the bootstrap gives results close to the traditional t-test method. This means either method gives a reliable estimate of the range for the true mean of medv.
ci_boot <- c(mu_hat - 2 * boot_se, mu_hat + 2 * boot_se)
ci_boot## [1] 21.71148 23.35413t.test(Boston$medv)$conf.int## [1] 21.72953 23.33608
## attr(,"conf.level")
## [1] 0.95Based on this data set, provide an estimate, µˆmed, for the median value of medv in the population.
mu_hat_med <- median(Boston$medv)
mu_hat_med## [1] 21.2We now would like to estimate the standard error of µˆ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.
The median value of medv is 21.2. Since there’s no exact
formula to calculate the standard error of the median, we used the
bootstrap method. The bootstrap estimated the standard error to be about
0.3778. This means that the median value is expected to vary by about
0.38 if we repeatedly took new random samples from the population. The
bootstrap is helpful here because it gives a good estimate without
needing a formula.
boot_median <- function(data, index) {
return(median(data[index]))
}
set.seed(1)
boot_med <- boot(Boston$medv, boot_median, R = 1000)
se_med <- sd(boot_med$t)
se_med## [1] 0.3778075Based on this data set, provide an estimate for the tenth percentile of medv in Boston census tracts. Call this quantity µˆ0.1. (You can use the quantile() function.)
mu_hat_0.1 <- quantile(Boston$medv, 0.1)
mu_hat_0.1## 10%
## 12.75Use the bootstrap to estimate the standard error of µˆ0.1. Comment on your fIndings.
The 10th percentile of medv was estimated to be 12.75,
meaning 10% of homes are below this value. Since there’s no simple
formula to get the standard error, the bootstrap method was used. It
estimated the standard error to be about 0.48. This means the 10th
percentile might change by around 0.48 in different samples. The
bootstrap gives a good way to measure this variation.
boot_p10 <- function(data, index) {
return(quantile(data[index], 0.1))
}
set.seed(1)
boot_10 <- boot(Boston$medv, boot_p10, R = 1000)
se_p10 <- sd(boot_10$t)
se_p10## [1] 0.4767526