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 a set of observations into K groups, or folds, that are all relatively equal in size. The first fold is treated as the validation set, and the predictive method is then fit on remaining \((K-1)\) folds. The mean squared error is calculated on the held out fold to get a sense for how well the predictive 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:
(i) The validation set approach?
Variability of validation estimates for error rate in K-fold cross validation is a lot lower than variability seen in estimates for test error produced by the validation set approach. K-fold cross-validation is a much more stable approach given the lower variability in its estimates of test 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 an over estimate of the test error rate because statistical models tend to perform worse when trained on fewer observations.
(ii) LOOCV?
Leave-One-Out Cross-Validation (LOOCV) method 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 larger data sets, LOOCV will be much more computationally expensive.
LOOCV will provide a less biased estimate of test error given that each training set contains \(n-1\) observations. Ultimately the K-fold cross-validation approach will have less variance because the LOOCV method averages the outputs of \(n\) fitted models trained on data sets that are all highly correlated. K-fold cross-validation 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\) has been shown to yield test error rates that have both intermediate bias and 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.
From the logistic regression model produced, it appears that both the
income and balance variables are significant
predictors given that they both have p-values less than 0.05.
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:
(i)
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,](ii)
Fit a multiple logistic regression model using only the training observations.
glm_default<-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_default,default_test,type="response")
glm_pred<-factor(ifelse(glm_probs>=0.50,"Yes","No"))(iv)
Compute the validation set error, which is the fraction of the observations in the validation set that are missclassified.
The multiple logistic regression model yields a validation set error of 0.023
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.
set.seed(6)
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 1930 44
## Yes 9 17mean(glm_pred!=default_test$default)## [1] 0.0265The second split yields a slightly lower error rate for the validation set of 0.024.
set.seed(7)
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 1927 38
## Yes 9 26mean(glm_pred!=default_test$default)## [1] 0.0235The third split also yields an error rate for the validation set of 0.027.
set.seed(10)
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 1930 44
## Yes 9 17mean(glm_pred!=default_test$default)## [1] 0.0265The 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 .0255. Including this
variable does not appear to produce a significant change in the test
error seen from previous validation estimates in our initial model which
only included income and balance as
predictors.
set.seed(9)
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 1932 44
## Yes 7 17mean(glm_pred!=default_test$default)## [1] 0.0255By reviewing a summary of our model, we can see that including the
dummy variable for student as a predictor is not considered
significant given a p-value of 0.131. It also results in the
income variable no longer being considered significant
given a p-value of 0.244.
summary(glm_default)##
## Call:
## glm(formula = default ~ income + balance + student, family = "binomial",
## data = default_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5129 -0.1389 -0.0534 -0.0189 3.7583
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.135e+01 5.671e-01 -20.020 <2e-16 ***
## income 1.062e-05 9.108e-06 1.166 0.244
## balance 5.827e-03 2.631e-04 22.145 <2e-16 ***
## studentYes -3.963e-01 2.626e-01 -1.509 0.131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2374.1 on 7999 degrees of freedom
## Residual deviance: 1257.3 on 7996 degrees of freedom
## AIC: 1265.3
##
## 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))
}
coef(glm_fit)## (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.
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 obtained 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}\).
The estimated mean for the medv variable is 22.53.
mu_hat<-mean(medv)
mu_hat## [1] 22.53281(b)
Provide an estimate of the standard error of \(\hat{\mu}\). Interpret this result.
Our calculated standard error for \(\hat{\mu}\) is .409. Given a 95% confidence interval, we could expect the real population mean to be between 21.698 and 23.302. The lower bound for this range is calculated by \(22.5-(1.96*.409)=21.698\) and the upper bound is calculated by \(22.5+(1.96*.409)=23.302\).
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 using our bootstrap estimate of
0.407 from (c) is shown to be 21.735 to 23.331. This is very similar to
the confidence interval yielded by using
t.test(Boston$medv) of 21.730 to 23.336.
Confint_mu_hat<-c((mu_hat-1.96*0.407),(mu_hat+1.96*0.407))
Confint_mu_hat## [1] 21.73509 23.33053t.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.388. Given a 95% confidence interval, we would expect the median to be between \(21.2-(1.96*.388)=20.440\) and \(21.2+(1.96*.388)=21.961\). The median value from this function is also shown to be 21.2.
set.seed(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.02465 0.3879954(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 .487. Given a 95%
confidence interval, we would expect the tenth percentile of
medv to be between \(12.75-(1.96*.487)=11.796\) and \(12.75+(1.96*.487)=13.705\). Our bootstrap
function also shows the tenth percentile of medv to be
12.75.
set.seed(3)
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.0331 0.4873263detach(Boston)