Chapter 05 (page 197): 3, 5, 6, 9
3. We now review k-fold cross-validation.
(a) Explain how k-fold cross-validation is implemented.
K-fold cross-validation splits the dataset into k different parts. K can be any number, but is often 5 or 10. The model is then fit multiple times, where each time, one of the k parts is removed, the model is fit, then tested against the removed portion of the data. After the model is fit k times, an average of the mean squared error (MSE) for each of the fit models is taken and used as an estimate for the error rate.
(b) What are the advantages and disadvantages of k-fold cross-validation relative to:
i. The validation set approach?
K-fold cross-validation will result in less bias and lower MSE variance than the validation set approach. Also, when fitting a model on a smaller dataset, it may be impractical to split the data into dedicated training and testing sets as is needed in the validation set approach. K-fold is more computationally intensive than using the validaiton set approach.
ii. LOOCV?
While the LOOCV method has less bias than k-fold cross-validation, k-fold has lower variance in the MSEs. Additionally, for large datasets, k-fold will be much faster and computationally less taxing than LOOCV, since fewer models need to be fit (k fits, instead of one for each observation).
5. 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.
library(ISLR)## Warning: package 'ISLR' was built under R version 3.6.3set.seed(1)
attach(Default)(a) Fit a logistic regression model that uses income and balance to predict default.
glm.fit.log = glm(default~income+balance, data=Default, family = binomial)
summary(glm.fit.log)##
## 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.
train=sample(nrow(Default), .8*nrow(Default))
default.test=default[!train]ii. Fit a multiple logistic regression model using only the training observations.
glm.fit.log.val.set = glm(default~income+balance, data=Default, subset=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.
prob.log.val.set=predict(glm.fit.log.val.set, Default[-train,], type ='response')
pred.log.val.set=rep('No',nrow(Default[-train,]))
pred.log.val.set[prob.log.val.set>0.5] = 'Yes'
table(pred.log.val.set, Default[-train,]$default)##
## pred.log.val.set No Yes
## No 1928 50
## Yes 2 20iv. Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
mean(pred.log.val.set!=Default[-train,]$default)## [1] 0.026We 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 co-efficients 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)
glm.fit.default = glm(default~income+balance, data = Default, family = binomial)
summary(glm.fit.default)##
## 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: 8Based on the above model, the estimated standard errors are \(4.985 X 10^{-6}\) for income and \(2.274 X 10^{-4}\) for balance.
(b) Write a function, boot.fn(), that takes as input the Default dataset 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){
cofs = coef(glm(default~income+balance, data = data, subset = index, family = binomial))
return(cofs)
}
boot.fn(Default, 1:nrow(Default))## (Intercept) income balance
## -1.154047e+01 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)
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-04Using the bootstrap method, the standard errors are \(4.987 X 10^{-6}\) for income and \(2.299 X 10^{-4}\) for balance.
(d) Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.
The estimated standard errors obtained using the glm() function and the boot() function are nearly identical.
9. We will now consider the Boston housing data set, from the MASS library.
library(MASS)
attach(Boston)
set.seed(1)(a) Based on this data set, provide an estimate for the population mean of medv. Call this estimate \(\hat{\mu}\).
mean.medv = mean(medv)
mean.medv## [1] 22.53281\(\hat{\mu}\) = 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.
stderr.medv = sd(medv)/sqrt(nrow(Boston))
stderr.medv## [1] 0.4088611The standard error of \(\hat{\mu}\) = .4088611
(c) Now estimate the standard error of \(\hat{\mu}\) using the bootstrap. How does this compare to your answer from (b)?
boston.boot.fn = function(col,index){
cof = mean(col[index])
return(cof)
}
bootstrap = boot(medv, boston.boot.fn, R = 1000)
bootstrap##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boston.boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.007650791 0.4106622Using the bootstrap method, the standard error of \(\hat{\mu}\) = .3982354, which is close to the answer of .4088611 from question B.
(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})]\).
left.ci = 22.53281 - 2*.394032
right.ci = 22.53281 + 2*.394032
ci.table = matrix(c(left.ci, right.ci), nrow = 1, ncol =2)
ci.table## [,1] [,2]
## [1,] 21.74475 23.32087An approximation of the 95% confidence interval for the mean of medv is \(21.74475 \leq \hat{\mu} \leq 23.32087\)
(e) Based on this data set, provide an estimate, \(\hat{\mu}_{med}\), for the median value of medv in the population.
med.medv = median(medv)
med.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.
med.boot.fn = function(data,index){
est.med = median(data[index])
return(est.med)
}
boot(medv, med.boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = med.boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.0386 0.3770241The boostrap method produced a median of 21.2, with a small standard error of .3610336.
(g) Based on this data set, provide an estimate for the tenth percentile of medv in Boston suburbs. Call this quantity \(\hat{\mu}_{0.1}\). (You can use the quantile() function.)
percentile.ten = quantile(medv,c(.1))
percentile.ten## 10%
## 12.75\(\hat{\mu}_{0.1}\) = 12.75
(h) Use the bootstrap to estimate the standard error of \(\hat{\mu}_{0.1}\). Comment on your findings.
precentile.boot.fn = function(data,index){
precentile = quantile(data[index],c(.1))
return(precentile)
}
boot(medv, precentile.boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = precentile.boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.0186 0.4925766The estimated standard error of \(\hat{\mu}_{0.1}\) is very small (~.5) compared to the estimated value of \(\hat{\mu}_{0.1}\)