(a) Explain how k-fold cross-validation is
implemented.
In this method, k-fold Cross Validation(CV) is used to randomly divide
the set of data into k roughly equal groups, or folds. The first fold is
used as a validation set, while the following k-1 folds are used to fit
the algorithm. Then, using the observations in the held-out fold, the
mean squared error, or \(MSE_1\), is
calculated. This process is repeated k times, with each time the
validation set consisting of a distinct collection of observations. K
estimations of the test error \(MSE_1\), \(MSE_2\),…, \(MSE_k\) are produced as a result of this
approach. By averaging these values, the k-fold CV estimate is
created.
(b) What are the advantages and disadvantages of k-fold cross
validation relative to:
i. The validation set approach?
Advantage of K-fold CV approach over Validation set
approach: Depending on precisely which observations are
included in the training set and which observations are included in the
validation set approach, the validation estimate of the test error rate
can vary greatly. Additionally, the test error rate for the model fit on
the full data set may be overestimated by the validation set error
rate.In comparison to validation set approach, K-fold CV have lower
variability in the test error estimates.
Disadvantage: Validation set approach is conceptually
simple and easy to implement in comparison to K-fold Cross
Validation.
ii. LOOCV?
Advantage: An exception to the general rule of k-fold
cross-validation is the LOOCV cross-validation strategy, where k=n. The
statistical learning approach must be fitted n times. The computational
cost of this could be high. Furthermore, k-fold CV frequently provides
estimates of the test error rate that are more precise than LOOCV.
Disadvantage: LOOCV should be preferred over k-fold CV
if bias reduction is the primary goal because it has a tendency to be
less biased. Therefore, the choice of k in k-fold cross-validation
involves a trade-off between bias and variance; typically, employing k=5
or k=10 results in test error rate estimates that neither have an
extremely high bias nor a very high variance.
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.library(ISLR)(a) Fit a logistic regression model that uses
income and balance to predict
default.
set.seed(1)
fit.glm = glm(default ~ income + balance, data = Default, family = "binomial")
summary(fit.glm)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4725 -0.1444 -0.0574 -0.0211 3.7245
##
## 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.
inTrain = sample(dim(Default)[1], dim(Default)[1] / 2)
default_train = Default[inTrain,]
default_test = Default[-inTrain,]ii. Fit a multiple logistic regression model using only the training observations.
glm_model = glm(default ~ income + balance, data = default_train, family = "binomial")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.
glm.probs = predict(glm_model, default_test, type="response")
dim(Default)## [1] 10000 4glm.pred=rep("No",5000)
glm.pred[glm.probs>0.5] = "Yes"iv. Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
mean(glm.pred != default_test$default)## [1] 0.0254(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.
for(i in 1:3){
inTrain = sample(dim(Default)[1], dim(Default)[1] / 2)
default_train = Default[inTrain,]
default_test = Default[-inTrain,]
glm_model = glm(default ~ income + balance, data = default_train, family = "binomial")
glm.probs = predict(glm_model, default_test, type="response")
glm.pred=rep("No",5000)
glm.pred[glm.probs>0.5] = "Yes"
print(paste("Split Number: ", i))
print(paste("Validation Error Rate", mean(glm.pred != default_test$default)))
}## [1] "Split Number: 1"
## [1] "Validation Error Rate 0.0274"
## [1] "Split Number: 2"
## [1] "Validation Error Rate 0.0244"
## [1] "Split Number: 3"
## [1] "Validation Error Rate 0.0244"The three distinct splits that are performed by iterating in a for
loop each result in a different error rate, demonstrating that the rate
varies depending on which observations are in the training and
validation sets.
(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.
inTrain = sample(dim(Default)[1], dim(Default)[1] / 2)
default_train = Default[inTrain,]
default_test = Default[-inTrain,]
glm_model2 = glm(default ~ income + balance + student, data = default_train, family = "binomial")
glm.probs = predict(glm_model2, default_test, type="response")
glm.pred=rep("No",5000)
glm.pred[glm.probs>0.5] = "Yes"
table(glm.pred, default_test$default)##
## glm.pred No Yes
## No 4808 121
## Yes 18 53 mean(glm.pred != default_test$default)## [1] 0.0278There is no apparent decrease in test error rate after using the dummy variable for student.
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 coefficients 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 coefficients
associated with income and balance in a
multiple logistic regression model that uses both
predictors.
set.seed(1)
fit.glm = glm(default ~ income + balance, data = Default, family = "binomial")
summary(fit.glm)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4725 -0.1444 -0.0574 -0.0211 3.7245
##
## 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 coefficient estimates for
income and balance in the multiple logistic
regression model.
boot.fn = function(data, index) {
return(coef(glm(default ~ balance + income, data = data, family = binomial, subset = index)))
}(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)
boot(Default, boot.fn, 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* 5.647103e-03 1.855765e-05 2.298949e-04
## t3* 2.080898e-05 1.680317e-07 4.866284e-06(d) Comment on the estimated standard errors obtained using
the glm() function and using your bootstrap
function.
While the standard errors acquired using the glm() function are
4.348e-01, 4.985e-06, and 2.274e-04, those obtained using the bootstrap
function are 4.344722e-01, 2.298949e-04, and 4.866284e-06. These
standard errors are quite similar.
Boston housing data
set, from the ISLR2 library.library(ISLR2)
attach(Boston)(a) Based on this data set, provide an estimate for the
population mean of medv. Call this estimate
ˆµ.
mean_medv = mean(medv)
mean_medv## [1] 22.53281(b) Provide 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
number of observations.
std_mean = sd(medv)/sqrt(length(medv))
std_mean## [1] 0.4088611(c) Now estimate the standard error of ˆµ using the bootstrap. How does this compare to your answer from (b)?
set.seed(1)
boot.fn = function(data, index) (return(mean(data[index])))
boot(medv, boot.fn, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.007650791 0.4106622Standard error estimate in (b) is similar to the standard error rate
computed by bootstrap function ie. 0.41.
(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 [ˆµ − 2SE(ˆµ), µˆ + 2SE(ˆµ)].
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.53281CImean_medv <- c(22.5328 - 2 * 0.4106622, 22.5328 + 2 * 0.4106622)
CImean_medv## [1] 21.71148 23.35412(e) Based on this data set, provide an estimate,
ˆµmed, for the median value of medv in the
population.
median_medv = median(medv)
median_medv## [1] 21.2(f) We 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.
boot_fn2 = function(data,index)return(median(data[index]))
boot(medv, boot_fn2, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot_fn2, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.0386 0.3770241We obtain at an estimated median value of 21.2, which is the same as
the result in (e), and a standard error of 0.38, which is reasonably
minimal when compared to the median value.
(g) Based 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.)
percent.hat = quantile(medv, c(0.1))
percent.hat## 10%
## 12.75(h) Use the bootstrap to estimate the standard error of ˆµ0.1. Comment on your findings.
boot_fn3 = function(data, index) return(quantile(data[index], c(0.1)))
boot(medv,boot_fn3, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot_fn3, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.0186 0.4925766detach(Boston)