We now review k-fold cross-validation.
(a) Explain how k-fold cross-validation is implemented.
K-fold Cross-validation is implemented by first randomly dividing set of observations into K groups, or folds, that are all relatively equal in size. The first fold treated as validation set, and the predictive method is then fit on remaining (K-1) folds. The mean squared error is calculated on held out fold to get a sense for how well the prediction model performs.
This procedure is repeated K times, with a different group of observations held out as a validation set each time. K is most commonly set to 5 or 10, meaning that the procedure is performed 5 or 10 times.
Ultimately this process results in K estimates of test error. These estimates are averaged together to calculate the K-fold cross-validation estimate.
(b) What are the advantages and disadvantages of k-fold cross-validation relative to:
Variability of validation estimates for error rate in K-fold cross validation is a lot lower than variability seen in validation set approach. K-fold cross-validation is a much more stable approach given the lower variability in its estimates of error rate.
Another main draw back of the validation set approach is that the model is fit on only a subset of the data. The validation set approach can lead to over estimate of test error rate because statistical models tend to perform worse when trained on fewer observations.
Leave-One-Out Cross-Validation (LOOCV) is a special case of K-fold cross-validation where \(K\) is set to equal the total number of observations in the data set, \(n\). When using K-fold cross-validation you are able to set \(K\) equal to 5 or 10, which would require less computational effort than required by LOOCV. LOOCV requires fitting a learning method for every time you have an observation in the data set. With large data sets LOOCV will be much more computationally expensive.
LOOCV will provide less biased estimates of test error given that each training set contains \(n-1\) observations. Ultimately the K-fold cross-validation approach will have less variance because with LOOCV we are averaging the outputs of \(n\) fitted models trained on data sets that are highly correlated.. K-fold results in a more intermediate level of bias.
When determining \(K\) in K-fold cross-validation, there is a bias-variance trade-off. Using \(K=5\) or \(K=10\) will typically yield test error rates that have neither excessively high bias or very high variance.
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(ISLR2)
attach(Default)
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 ...(a)
Fit a logistic regression model that uses income
and balance to predict default.
Default_model<-glm(default~income+balance,data=Default,family="binomial")
summary(Default_model)##
## 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:
Split the sample set into a training set and a validation set.
set.seed(3)
Train<-sample(1:nrow(Default),0.8*nrow(Default))
default_train<-Default[Train,]
default_test<-Default[-Train,]Fit a multiple logistic regression model using only the training observations.
glm_default<-glm(default~income+balance,data=default_train,family="binomial")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_default,default_test,type="response")
glm_pred<-factor(ifelse(glm_probs>=0.50,"Yes","No"))Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
table(glm_pred,default_test$default)##
## glm_pred No Yes
## No 1941 40
## Yes 6 13mean(glm_pred!=default_test$default)## [1] 0.023(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.
In the first split, we see a slight increase in the error rate of the validation set to 0.027.
Train<-sample(1:nrow(Default),0.8*nrow(Default))
default_train<-Default[Train,]
default_test<-Default[-Train,]
glm_default<-glm(default~income+balance,data=default_train,family="binomial")
glm_probs<-predict(glm_default,default_test,type="response")
glm_pred<-factor(ifelse(glm_probs>=0.50,"Yes","No"))
table(glm_pred,default_test$default)##
## glm_pred No Yes
## No 1939 37
## Yes 6 18mean(glm_pred!=default_test$default)## [1] 0.0215The second split yields a similar error rate for the validation set of 0.026.
Train<-sample(1:nrow(Default),0.8*nrow(Default))
default_train<-Default[Train,]
default_test<-Default[-Train,]
glm_default<-glm(default~income+balance,data=default_train,family="binomial")
glm_probs<-predict(glm_default,default_test,type="response")
glm_pred<-factor(ifelse(glm_probs>=0.50,"Yes","No"))
table(glm_pred,default_test$default)##
## glm_pred No Yes
## No 1928 47
## Yes 7 18mean(glm_pred!=default_test$default)## [1] 0.027The third split yielded the highest error rate for the validation set of 0.03.
Train<-sample(1:nrow(Default),0.8*nrow(Default))
default_train<-Default[Train,]
default_test<-Default[-Train,]
glm_default<-glm(default~income+balance,data=default_train,family="binomial")
glm_probs<-predict(glm_default,default_test,type="response")
glm_pred<-factor(ifelse(glm_probs>=0.50,"Yes","No"))
table(glm_pred,default_test$default)##
## glm_pred No Yes
## No 1920 41
## Yes 11 28mean(glm_pred!=default_test$default)## [1] 0.026The validation estimate of the test error rate varies as we try different splits of the observations into training 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.
Including a dummy variable for student in our logistic
regression model yields a test error rate of .0295. Including this
variable does not appear to reduce the test error seen from previous
validation estimates in our initial model which only included
income and balance as predictors.
glm_default<-glm(default~income+balance+student,data=default_train,family="binomial")
glm_probs<-predict(glm_default,default_test,type="response")
glm_pred<-factor(ifelse(glm_probs>=0.50,"Yes","No"))
table(glm_pred,default_test$default)##
## glm_pred No Yes
## No 1919 41
## Yes 12 28mean(glm_pred!=default_test$default)## [1] 0.0265By reviewing a summary of our model, we can see that including the
dummy variable for student is considered significant given
a p-value of 0.04. It also results in the income variable
no longer being considered significant given a p-value of 0.69.
summary(glm_default)##
## Call:
## glm(formula = default ~ income + balance + student, family = "binomial",
## data = default_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5000 -0.1380 -0.0539 -0.0195 3.7715
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.108e+01 5.629e-01 -19.688 <2e-16 ***
## income 3.306e-06 9.145e-06 0.362 0.7177
## balance 5.854e-03 2.684e-04 21.813 <2e-16 ***
## studentYes -5.870e-01 2.645e-01 -2.219 0.0265 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2320.3 on 7999 degrees of freedom
## Residual deviance: 1243.9 on 7996 degrees of freedom
## AIC: 1251.9
##
## Number of Fisher Scoring iterations: 8We 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 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.
From the summary of the multiple logistic regression model produced,
we can see that the coefficient associated with income has
an estimated standard error of 4.985e-06 and the coefficient associated
with balance has an estimated standard error of
2.274e-04.
set.seed(2)
glm_fit<-glm(default~income+balance,data=Default,family="binomial")
summary(glm_fit)##
## 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){
glm_fit<-glm(default~income+balance,data=data,family="binomial",subset=index)
return(coef(glm_fit))
}(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.
Using the boot_fn, yields estimated standard errors of
the logistic regression coefficients for income and
balance of 4.924565e-06 and 2.228767e-04 respectively.
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.956457e-02 4.228903e-01
## t2* 2.080898e-05 -4.605498e-07 4.924565e-06
## t3* 5.647103e-03 3.048004e-05 2.228767e-04(d)
Comment on the estimated standard errors obtained using the
glm() function and using your bootstrap
function.
The estimated standard errors obtains using the bootstrap function
are very similar to those obtained using the glm()
function.
detach(Default)We will now consider the Boston housing data
set, from the ISLR2 library.
library(ISLR2)
attach(Boston)
str(Boston)## 'data.frame': 506 obs. of 13 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 ...
## $ 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 ...(a)
Based on this data set, provide an estimate for the
population mean of medv. Call this estimate \(\hat{\mu}\).
mu_hat<-mean(medv)
mu_hat## [1] 22.53281(b)
Provide an estimate of the standard error of \(\hat{\mu}\). Interpret this result.
stderror_hat<-sd(medv)/sqrt(dim(Boston)[1])
stderror_hat## [1] 0.4088611(c)
Now estimate the standard error of \(\hat{\mu}\) using the bootstrap. How does this compare to your answer from (b)?
Using the bootstrap function, we yield an estimated standard error of 0.407. This estimate for standard error is very similar to the estimate provided in part (b) of 0.409.
set.seed(2)
boot_fn<-function(data,index){
mu<-mean(data[index])
return(mu)
}
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.008232806 0.407299(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).
The confidence interval calculated of to be 21.72, 23.35 is very similar to the confidence interval yielded by the bootstrap function of 21.73, 23.34.
Confint_mu_hat<-c((mu_hat-2*0.407),(mu_hat+2*0.407))
Confint_mu_hat## [1] 21.71881 23.34681t.test(medv)##
## One Sample t-test
##
## data: 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(e)
Based on this data set, provide an estimate, \(\hat{\mu}med\), for the median value of
medv in the population.
med_hat<-median(medv)
med_hat## [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.
Using the bootstrap function yields an estimated standard error for the median of 0.373. The median value from this function is also shown to be 21.2.
boot_fn2<-function(data,index){
mu<-median(data[index])
return(mu)
}
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.02605 0.3726497(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.)
percentile_hat<-quantile(medv,c(0.1))
percentile_hat## 10%
## 12.75(h)
Use the bootstrap to estimate the standard error of \(\hat{\mu}0.1\). Comment on your findings.
Using the bootstrap function yields an estimated standard error for
the tenth percentile of medv of .4997. Our bootstrap
function also shows the tenth percentile of medv to be
12.75%.
boot_fn3<-function(data,index){
mu<-quantile(data[index],c(0.1))
return(mu)
}
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.00745 0.4997417detach(Boston)