(a) Explain how k-fold cross-validation is implemented
We have a set of data that is randomly split into k number of similarly sized groups. One of the groups is used as a validation set, and the other groups (k-1 folds) are combined to form the training set. The MSE is calculated on the validation set. This process will repeat until all of the groups have been used as the validation set at least one time. In total, this will give us k number of MSE values (one for each validation set). Finally, we can take the mean of these k number of MSE’s to estimate our test error for the whole model.
(b) What are the advantages and disadvantages of k-fold crossvalidation relative to:
The validation set approach has two drawbacks: the estimate of the test error rate can be highly variable, and the approach tends to lead to overestimates of the test error. k-fold validation is slightly more computationally intensive because we have k number of validation sets instead of just one. However, we get a more accurate estimate of our test error because we are using more data than the validation set approach. In practice, k-fold cross validation is used frequently when n is small compared to the validation set approach.
k-fold cross validation is much more efficient than LOOCV because with LOOCV you have to run the model n number of times which can be too computationally intensive if we have a lot of observations. k-fold CV also gives more accurate estimates of the test error rate than does LOOCV because LOOCV has very low bias and consequently much more variance. k-fold cross validation has moderate bias and variance (using k=5 or k=10 seem to offer the best results).
(a) Fit a logistic regression model that uses income and balance to predict default.
library(ISLR)## Warning: package 'ISLR' was built under R version 4.0.3glm.model = glm(default~income + balance, data=Default, family="binomial")
summary(glm.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:
set.seed(1)
train = sample(nrow(Default),0.8*nrow(Default))glm.model = glm(default~income + balance, data=Default, family="binomial", subset=train)
summary(glm.model)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4758 -0.1413 -0.0563 -0.0210 3.4620
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.168e+01 4.893e-01 -23.879 < 2e-16 ***
## income 2.547e-05 5.631e-06 4.523 6.1e-06 ***
## balance 5.613e-03 2.531e-04 22.176 < 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: 2313.6 on 7999 degrees of freedom
## Residual deviance: 1239.2 on 7997 degrees of freedom
## AIC: 1245.2
##
## Number of Fisher Scoring iterations: 8preds = predict(glm.model, newdata=Default[-train,], type="response")
preds.choice = ifelse(preds>0.5, "Yes","No")
caret::confusionMatrix(as.factor(preds.choice),as.factor(Default[-train,1]), positive = 'Yes')## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1928 50
## Yes 2 20
##
## Accuracy : 0.974
## 95% CI : (0.966, 0.9805)
## No Information Rate : 0.965
## P-Value [Acc > NIR] : 0.01371
##
## Kappa : 0.4252
##
## Mcnemar's Test P-Value : 7.138e-11
##
## Sensitivity : 0.2857
## Specificity : 0.9990
## Pos Pred Value : 0.9091
## Neg Pred Value : 0.9747
## Prevalence : 0.0350
## Detection Rate : 0.0100
## Detection Prevalence : 0.0110
## Balanced Accuracy : 0.6423
##
## 'Positive' Class : Yes
## mean(preds.choice != Default[-train, ]$default)## [1] 0.026The validation set approach gave us a 2.6% test error rate.
(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.
set.seed(2)
train = sample(nrow(Default),0.8*nrow(Default))
glm.model = glm(default~income + balance, data=Default, family="binomial", subset=train)
preds = predict(glm.model, newdata=Default[-train,], type="response")
preds.choice = ifelse(preds>0.5, "Yes","No")
caret::confusionMatrix(as.factor(preds.choice),as.factor(Default[-train,1]), positive = 'Yes')## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1929 33
## Yes 9 29
##
## Accuracy : 0.979
## 95% CI : (0.9717, 0.9848)
## No Information Rate : 0.969
## P-Value [Acc > NIR] : 0.0041153
##
## Kappa : 0.5699
##
## Mcnemar's Test P-Value : 0.0003867
##
## Sensitivity : 0.4677
## Specificity : 0.9954
## Pos Pred Value : 0.7632
## Neg Pred Value : 0.9832
## Prevalence : 0.0310
## Detection Rate : 0.0145
## Detection Prevalence : 0.0190
## Balanced Accuracy : 0.7315
##
## 'Positive' Class : Yes
## mean(preds.choice != Default[-train, ]$default)## [1] 0.021set.seed(3)
train = sample(nrow(Default),0.8*nrow(Default))
glm.model = glm(default~income + balance, data=Default, family="binomial", subset=train)
preds = predict(glm.model, newdata=Default[-train,], type="response")
preds.choice = ifelse(preds>0.5, "Yes","No")
caret::confusionMatrix(as.factor(preds.choice),as.factor(Default[-train,1]), positive = 'Yes')## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1941 40
## Yes 6 13
##
## Accuracy : 0.977
## 95% CI : (0.9694, 0.9831)
## No Information Rate : 0.9735
## P-Value [Acc > NIR] : 0.1837
##
## Kappa : 0.352
##
## Mcnemar's Test P-Value : 1.141e-06
##
## Sensitivity : 0.2453
## Specificity : 0.9969
## Pos Pred Value : 0.6842
## Neg Pred Value : 0.9798
## Prevalence : 0.0265
## Detection Rate : 0.0065
## Detection Prevalence : 0.0095
## Balanced Accuracy : 0.6211
##
## 'Positive' Class : Yes
## mean(preds.choice != Default[-train, ]$default)## [1] 0.023set.seed(4)
train = sample(nrow(Default),0.8*nrow(Default))
glm.model = glm(default~income + balance, data=Default, family="binomial", subset=train)
preds = predict(glm.model, newdata=Default[-train,], type="response")
preds.choice = ifelse(preds>0.5, "Yes","No")
caret::confusionMatrix(as.factor(preds.choice),as.factor(Default[-train,1]), positive = 'Yes')## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1931 35
## Yes 16 18
##
## Accuracy : 0.9745
## 95% CI : (0.9666, 0.981)
## No Information Rate : 0.9735
## P-Value [Acc > NIR] : 0.42557
##
## Kappa : 0.4014
##
## Mcnemar's Test P-Value : 0.01172
##
## Sensitivity : 0.3396
## Specificity : 0.9918
## Pos Pred Value : 0.5294
## Neg Pred Value : 0.9822
## Prevalence : 0.0265
## Detection Rate : 0.0090
## Detection Prevalence : 0.0170
## Balanced Accuracy : 0.6657
##
## 'Positive' Class : Yes
## mean(preds.choice != Default[-train, ]$default)## [1] 0.0255We get a different test error rate for each split but they are relatively close to each other.
(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.
set.seed(1)
train = sample(nrow(Default),0.8*nrow(Default))
glm.model = glm(default~income + balance + student, data=Default, family="binomial", subset=train)
preds = predict(glm.model, newdata=Default[-train,], type="response")
preds.choice = ifelse(preds>0.5, "Yes","No")
caret::confusionMatrix(as.factor(preds.choice),as.factor(Default[-train,1]), positive = 'Yes')## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1928 53
## Yes 2 17
##
## Accuracy : 0.9725
## 95% CI : (0.9644, 0.9792)
## No Information Rate : 0.965
## P-Value [Acc > NIR] : 0.03522
##
## Kappa : 0.3726
##
## Mcnemar's Test P-Value : 1.562e-11
##
## Sensitivity : 0.2429
## Specificity : 0.9990
## Pos Pred Value : 0.8947
## Neg Pred Value : 0.9732
## Prevalence : 0.0350
## Detection Rate : 0.0085
## Detection Prevalence : 0.0095
## Balanced Accuracy : 0.6209
##
## 'Positive' Class : Yes
## mean(preds.choice != Default[-train, ]$default)## [1] 0.0275If anything it looks like the inclusion of the student variable may have increased the test error.
(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.
glm.model = glm(default~income + balance, data=Default, family="binomial")
summary(glm.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) 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)))(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.890536e-02 4.347841e-01
## t2* 2.080898e-05 1.566550e-07 4.864757e-06
## t3* 5.647103e-03 1.843880e-05 2.301132e-04(d) Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.
The estimated standard error values for the coefficients are very similar for both of the functions.
(a) Based on this data set, provide an estimate for the population mean of medv. Call this estimate “mu”.
library(MASS)
mu.hat <- mean(Boston$medv)
mu.hat## [1] 22.53281(b) Provide an estimate of the standard error of ^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.
(sd(Boston$medv))/(sqrt(nrow(Boston)))## [1] 0.4088611(c) Now estimate the standard error of (mu-hat) using the bootstrap. How does this compare to your answer from (b)?
set.seed(1)
boot.fn <- function(data, index) {
mu <- mean(data[index])
return (mu)
}
boot(Boston$medv, boot.fn, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.007650791 0.4106622The standard error is very similar using both methods.
(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 [(mu-hat) ??? 2SE(mu-hat), (mu-hat) + 2SE(mu-hat)].
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.53281confint.mu.hat <- c(22.53281 - (2*0.4106622) , 22.53281 + (2*0.4106622))
confint.mu.hat## [1] 21.71149 23.35413Both confidence intervals are very similar between the two methods.
(e) Based on this data set, provide an estimate, (mu-hat)med, for the median value of medv in the population.
mu.med = median(Boston$medv)
mu.med## [1] 21.2(f) We now would like to estimate the standard error of ^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.
boot.fn <- function(data ,index)
return(median(data[index]))
set.seed(10003)
boot(Boston$medv ,boot.fn, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.0241 0.3826952The median is calculated as 21.2 which is the same as found in (e). The standard deviation of the median of medv is approximately 0.3829652 which is relatively small.
(g) Based on this data set, provide an estimate for the tenth percentile of medv in Boston suburbs. Call this quantity (mu-hat)0.1. (You can use the quantile() function.)
mu.hat1 <- quantile(Boston$medv, 0.1)
mu.hat1## 10%
## 12.75(h) Use the bootstrap to estimate the standard error of ^mu-0.1. Comment on your findings.
boot.fn <- function(data ,index)
return(quantile(data[index], .1))
set.seed(10003)
boot(Boston$medv ,boot.fn, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 -0.002 0.4998058The estimated 10th percentile value is 12.75 which is the same as the value calculated in (g). The standard error is 0.4998058 which is relatively small.