library(tidyverse)
library(openintro)
library(ISLR)
library(ISLR2)
library(boot)
library(knitr)
set.seed(4143)We now review k-fold cross validation.
Explain how k-fold cross-validation is implemented.
k-fold cross-validation (CV) divides the data set into k groups of observations, then performs leave-one-out cross-validation according to groups, rather than individual observations. This means that one kth of the data set is reserved to be used as the validation set, while a model is trained on the remaining k - 1 groups. Then, the mean squared error is computed on the validation set. This process is repeated iteratively until each group has been reserved and used as a validation set and a corresponding MSE computed. Then, the average of the k MSEs is computed to arrive at the cross-validation estimate.
What are the advantages and disadvantages of k-fold cross-validation relative to:
The validation set approach?
The validation set approach is simple to execute, but has high MSE variability, so it is an unstable method of validation. k-fold CV is much more stable because the model is able to train on the entire data set, reducing the bias of the model.
LOOCV?
Leave-One-Out Cross-Validation (LOOCV) is the maximum- k example of k-fold CV. As such, it has minimal bias as each iteration trains on n-1 observations, or approximately the entire data set. However, the k-fold CV test error estimate actually has lower variance than the LOOCV test error estimate, due to resulting from the mean of less-overlapping (less-correlated) models.
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.
Fit a logistic regression model that uses income
and balance to predict default.
GLM.default <- glm(default ~ income + balance, data = Default, family = binomial)
summary(GLM.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: 8Using 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.
train <- sample(10000, 5000)Fit a multiple logistic regression model using only the training observations.
GLM.def.train <- glm(default ~ income + balance,
data = Default, family = binomial, subset = train)
summary(GLM.def.train)##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5058 -0.1426 -0.0581 -0.0201 3.7308
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.172e+01 6.156e-01 -19.036 < 2e-16 ***
## income 2.536e-05 7.137e-06 3.553 0.000381 ***
## balance 5.662e-03 3.153e-04 17.961 < 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: 1510.51 on 4999 degrees of freedom
## Residual deviance: 796.78 on 4997 degrees of freedom
## AIC: 802.78
##
## Number of Fisher Scoring iterations: 8Obtain 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.
predi <- predict(GLM.def.train, Default, type = "response")[-train]
predictions <- ifelse(predi < 0.5, "No", "Yes")Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
table(predictions, Default$default[-train])##
## predictions No Yes
## No 4824 115
## Yes 17 44mean(predictions != Default$default[-train])## [1] 0.0264Validation set error, 50% train and 50% test: 2.64%
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.
# Train on 70% of data
train1 <- sample(10000, 7000)
GLM.def.train1 <- glm(default ~ income + balance,
data = Default, family = "binomial", subset = train1)
summary(GLM.def.train1)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2059 -0.1465 -0.0585 -0.0213 3.7259
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.159e+01 5.215e-01 -22.23 < 2e-16 ***
## income 2.225e-05 5.982e-06 3.72 0.000199 ***
## balance 5.651e-03 2.708e-04 20.86 < 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: 2057.2 on 6999 degrees of freedom
## Residual deviance: 1110.1 on 6997 degrees of freedom
## AIC: 1116.1
##
## Number of Fisher Scoring iterations: 8pred1 <- predict(GLM.def.train1, Default, type = "response")[-train1]
result1 <- ifelse(pred1 < 0.5, "No", "Yes")
table(result1, Default$default[-train1])##
## result1 No Yes
## No 2891 66
## Yes 11 32mean(result1 != Default$default[-train1])## [1] 0.02566667Validation set error, 70% train and 30% test: 2.566667%
# Train on 80% of data
train2 <- sample(10000, 8000)
GLM.def.train2 <- glm(default ~ income + balance,
data = Default, family = "binomial", subset = train2)
summary(GLM.def.train2)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4199 -0.1506 -0.0613 -0.0234 3.6809
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.124e+01 4.716e-01 -23.838 < 2e-16 ***
## income 2.016e-05 5.511e-06 3.658 0.000254 ***
## balance 5.479e-03 2.452e-04 22.340 < 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: 2354.0 on 7999 degrees of freedom
## Residual deviance: 1304.6 on 7997 degrees of freedom
## AIC: 1310.6
##
## Number of Fisher Scoring iterations: 8pred2 <- predict(GLM.def.train2, Default, type = "response")[-train2]
result2 <- ifelse(pred2 < 0.5, "No", "Yes")
table(result2, Default$default[-train2])##
## result2 No Yes
## No 1932 41
## Yes 4 23mean(result2 != Default$default[-train2])## [1] 0.0225Validation set error, 80% train and 20% test: 2.25%
# Train on 90% of data
train3 <- sample(10000, 9000)
GLM.def.train3 <- glm(default ~ income + balance,
data = Default, family = "binomial", subset = train3)
summary(GLM.def.train3)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train3)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4638 -0.1470 -0.0588 -0.0217 3.7091
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.143e+01 4.512e-01 -25.339 < 2e-16 ***
## income 2.013e-05 5.253e-06 3.833 0.000127 ***
## balance 5.607e-03 2.359e-04 23.770 < 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: 2657.5 on 8999 degrees of freedom
## Residual deviance: 1435.5 on 8997 degrees of freedom
## AIC: 1441.5
##
## Number of Fisher Scoring iterations: 8pred3 <- predict(GLM.def.train3, Default, type = "response")[-train3]
result3 <- ifelse(pred3 < 0.5, "No", "Yes")
table(result3, Default$default[-train3])##
## result3 No Yes
## No 968 23
## Yes 3 6mean(result3 != Default$default[-train3])## [1] 0.026Validation set error, 90% train and 10% test: 2.6%
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.
#New logistic regression model with `student`
GLM.def.student <- glm(default ~ income + balance + as.numeric(student),
data = Default, family = "binomial")
summary(GLM.def.student)##
## Call:
## glm(formula = default ~ income + balance + as.numeric(student),
## family = "binomial", data = Default)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4691 -0.1418 -0.0557 -0.0203 3.7383
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.022e+01 6.360e-01 -16.072 < 2e-16 ***
## income 3.033e-06 8.203e-06 0.370 0.71152
## balance 5.737e-03 2.319e-04 24.738 < 2e-16 ***
## as.numeric(student) -6.468e-01 2.363e-01 -2.738 0.00619 **
## ---
## 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: 1571.5 on 9996 degrees of freedom
## AIC: 1579.5
##
## Number of Fisher Scoring iterations: 8GLM.def.Strain <- glm(default ~ income + balance + as.numeric(student),
data = Default, family = "binomial", subset = train1)
summary(GLM.def.Strain)##
## Call:
## glm(formula = default ~ income + balance + as.numeric(student),
## family = "binomial", data = Default, subset = train1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1427 -0.1451 -0.0572 -0.0209 3.6895
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.052e+01 7.570e-01 -13.896 <2e-16 ***
## income 7.618e-06 9.762e-06 0.780 0.435
## balance 5.733e-03 2.768e-04 20.715 <2e-16 ***
## as.numeric(student) -5.334e-01 2.802e-01 -1.903 0.057 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2057.2 on 6999 degrees of freedom
## Residual deviance: 1106.5 on 6996 degrees of freedom
## AIC: 1114.5
##
## Number of Fisher Scoring iterations: 8pred1 <- predict(GLM.def.Strain, Default, type = "response")[-train1]
result1 <- ifelse(pred1 < 0.5, "No", "Yes")
table(result1, Default$default[-train1])##
## result1 No Yes
## No 2887 67
## Yes 15 31mean(result1 != Default$default[-train1])## [1] 0.02733333Validation set error, 70% train and 30% test: 2.733333%. The model
including a dummy variable for student does not show a
reduction in test error compared to previous models that did not include
this variable.
We 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.
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.
GLM.default <- glm(default ~ income + balance,
data = Default, family = "binomial")
summary(GLM.default)$coef## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.154047e+01 4.347564e-01 -26.544680 2.958355e-155
## income 2.080898e-05 4.985167e-06 4.174178 2.990638e-05
## balance 5.647103e-03 2.273731e-04 24.836280 3.638120e-136Standard error for income coefficient: 4.985e-06
Standard error for balance coefficient: 2.274e-04
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 ~ income + balance, data = Default, family = "binomial", subset = index)))
}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.fn(Default, sample(10000, 10000, replace = TRUE))## (Intercept) income balance
## -1.121457e+01 2.351554e-05 5.422823e-03boot(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 -2.518061e-02 4.263577e-01
## t2* 2.080898e-05 7.748263e-08 4.683468e-06
## t3* 5.647103e-03 1.195997e-05 2.289205e-04Comment on the estimated standard errors obtained using the
glm() function and using your bootstrap
function.
Standard error for income glm() coefficient:
4.985e-06
Standard error for income boot() coefficient:
4.877957e-06
Standard error for balance glm() coefficient:
2.274e-04
Standard error for balance boot() coefficient:
2.282513e-04
The corresponding values appear very close.
We will now consider the Boston housing data set,
from the ISLR2 library.
set.seed(4341)Based on this data set, provide an estimate for the
population mean of medv. Call this estimate
ˆμ.
mu.hat <- mean(Boston$medv)ˆμ = 22.533
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.
stdE <- function(data, index){
X <- data$medv[index]
return((sd(X))/sqrt(506))
}
stdE(Boston, 1:506)## [1] 0.4088611The estimated standard error of ˆμ is 0.4088611, which means that the average median value of owner-occupied homes in this data set is $22,532.81 plus-or-minus $408.86.
Now estimate the standard error of ˆμ using the bootstrap. How does this compare to your answer from (b)?
set.seed(4341)
stdE(Boston, sample(506, 506, replace = TRUE))## [1] 0.3802379SE <- boot(Boston, stdE, R=1000)
SE##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston, statistic = stdE, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.4088611 -0.001462205 0.01649827The estimated standard error of ˆμ using bootstrap is exactly the
same as the standard error given by the formula
standard deviation/√(n), which is
0.4088611.
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(ˆμ)].
CI.boot.low <- mu.hat - 2*SE$t0
CI.boot.high <- mu.hat + 2*SE$t0
ttest <- t.test(Boston$medv)The 95% confidence interval for the mean of medv based
on bootstrapping is [21.72, 23.35]. The confidence interval obtained
from the t-test is [21.73, 23.34], which is 0.01 units narrower on each
side than the bootstrap confidence interval.
Based 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.2The estimated population median of medv using
bootstrapping is 21.2 (or $21,200).
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.
med <- function(data, index){
X <- data$medv[index]
return(median(X))
}
med(Boston, 1:506)## [1] 21.2set.seed(4341)
med(Boston, sample(506, 506, replace = TRUE))## [1] 21.4boot(Boston, med, R=1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston, statistic = med, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.01495 0.3740841The standard error of ˆμmed (mu.hat.med) is
0.3740841.
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.)
mu.hat.01 <- quantile(Boston$medv, probs = 0.1)
mu.hat.01## 10%
## 12.75ˆμ0.1 = 12.8
Use the bootstrap to estimate the standard error of ˆμ0.1. Comment on your findings.
#bootstrapped
quantile.boot <- function(data, index){
X <- data$medv[index]
return(quantile(X, probs = 0.1))
}
set.seed(4143)
quantile.boot(Boston, sample(506, 506, replace = TRUE))## 10%
## 12.6boot(Boston, quantile.boot, R=1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston, statistic = quantile.boot, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 -0.0118 0.4993948The standard error of ˆμ0.1 is 0.4993948.