We now review k-fold cross-validation.
Q3(a) Explain how k-fold cross-validation is
implemented.
A3(a) Th K-fold cross-validation is implemented by
splitting the data set into groups/folds of the same size. First part is
removed and the train is conducted on the remainder of the data set
using the first group as the test, this is repeated with each fold/group
(similarity to leave one out, but K is set to the number of observations
vs just one observation that occurs in the leave one out cross
validation. Within k-fold cross validation observations are specific to
their fold and do not repeat. Results end up being an average of the
estimates.
Q3(b)What are the advantages and disadvantages of
k-fold cross validation relative to:
Q3(b)i The validation set approach?
A3(b)i Relative to validation set approach, the K-fold
CV validation:
Advantages
- K-fold CV is less bias since it uses bigger training sets
- K-fold CV has a much lower variability with smaller datasets
- All data is used to train and test the model
Disadvantages
- K-fold CV is not as easy of a concepts to understand or explain.
- K-fold CV takes up more computation/needs more time to consume large
data sets
Q3(b)ii Relative to LOOCV, the k-fold validation:
A3(b)ii
Advantages - k-fold CV often gives more accurate
estimates of the test error rate than does LOOCV
- K-fold CV has a smaller variance than LOOCV.
- K-fold CV is less computationally intensive than LOOCV.
Disadvantages
- LOOCV is less bias than K-fold CV (when k < n).
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)
library(boot)
set.seed(1)
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 ...attach(Default)Q5(a) Fit a logistic regression model that uses income and balance to predict default.
glm.fit = glm(default~income+balance,
family = 'binomial',
data = Default)
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: 8Q5(b) Using the validation set approach, estimate
the test error of this model. In order to do this, you must perform the
following steps.
A5(b) The estimated test error of this model is
2.7%.
Q5(b)iSplit the sample set into a training set and a
validation set.
A5(b)i Accomplished below:
data = 1:10000
test = 1:3000
train = data[data != test]
test = Default[1:3000,]
test.target = default[1:3000]Q5(b)iiFit a multiple logistic regression model
using only the training observations.
A5(b)ii Accomplished below:
glm.fit2 = glm(default~income+balance,
data = Default,family='binomial',
subset=train)
summary(glm.fit2)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3532 -0.1499 -0.0614 -0.0235 3.6830
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.108e+01 4.955e-01 -22.368 <2e-16 ***
## income 1.666e-05 5.821e-06 2.862 0.0042 **
## balance 5.416e-03 2.605e-04 20.790 <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: 2037.0 on 6999 degrees of freedom
## Residual deviance: 1128.7 on 6997 degrees of freedom
## AIC: 1134.7
##
## Number of Fisher Scoring iterations: 8Q5(b)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.
A5(b)iii Accomplished below:
glm.prob=predict(glm.fit2, Default[-train, ], type = "response")
glm.pred=rep("No", length(glm.prob))
glm.pred[glm.prob >.5]="Yes"
table(glm.pred, test.target)## test.target
## glm.pred No Yes
## No 2892 74
## Yes 7 27Q5(b)iv Compute the validation set error, which is
the fraction of the observations in the validation set that are
misclassified.
A5(b)iv 2.7% is the fraction of the observation in the
validation set that are misclassified
1 - mean(glm.pred==test.target)## [1] 0.027Q5(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.
A5(c) Repeating the various splits of the observations
and running the models resulted in a finding of test error rates
consistently in a small range, ranging from 2.43-2.87%.
Split 2
data = 1:10000
test = 3001:6000
train = data[data != test]
test = Default[3001:6000,]
test.target = default[3001:6000]
glm.fit2 = glm(default~income+balance,
data = Default,
family='binomial',
subset=train)
summary(glm.fit2)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5808 -0.1379 -0.0523 -0.0183 3.7644
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.205e+01 5.475e-01 -22.017 < 2e-16 ***
## income 2.409e-05 6.078e-06 3.964 7.38e-05 ***
## balance 5.912e-03 2.866e-04 20.628 < 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: 1082.3 on 6997 degrees of freedom
## AIC: 1088.3
##
## Number of Fisher Scoring iterations: 8glm.prob=predict(glm.fit2, Default[-train, ], type = "response")
glm.pred=rep("No", length(glm.prob))
glm.pred[glm.prob >.5]="Yes"
1 - mean(glm.pred==test.target)## [1] 0.02866667Split 3
data = 1:10000
test = 6001:9000
train = data[data != test]
test = Default[6001:9000,]
test.target = default[6001:9000]
glm.fit2 = glm(default~income+balance,
data = Default,
family='binomial',
subset=train)
summary(glm.fit2)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4494 -0.1497 -0.0624 -0.0232 3.6802
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.129e+01 5.044e-01 -22.384 < 2e-16 ***
## income 2.105e-05 5.928e-06 3.551 0.000383 ***
## balance 5.512e-03 2.630e-04 20.958 < 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: 1137.2 on 6997 degrees of freedom
## AIC: 1143.2
##
## Number of Fisher Scoring iterations: 8glm.prob=predict(glm.fit2, Default[-train, ], type = "response")
glm.pred=rep("No", length(glm.prob))
glm.pred[glm.prob >.5]="Yes"
1 - mean(glm.pred==test.target)## [1] 0.02433333Split 4
data = 1:10000
test = 2001:4000
train = data[data != test]
test = Default[2001:4000,]
test.target = default[2001:4000]
glm.fit2 = glm(default~income+balance,
data = Default,
family='binomial',
subset=train)
summary(glm.fit2)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4548 -0.1458 -0.0575 -0.0213 3.7219
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.146e+01 4.838e-01 -23.690 < 2e-16 ***
## income 1.886e-05 5.539e-06 3.405 0.000662 ***
## balance 5.636e-03 2.559e-04 22.026 < 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: 2327.1 on 7999 degrees of freedom
## Residual deviance: 1274.8 on 7997 degrees of freedom
## AIC: 1280.8
##
## Number of Fisher Scoring iterations: 8glm.prob=predict(glm.fit2, Default[-train, ], type = "response")
glm.pred=rep("No", length(glm.prob))
glm.pred[glm.prob >.5]="Yes"
1 - mean(glm.pred==test.target)## [1] 0.026Q5(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.
A5(d) Introducing the dummy variable for
student did not change or reduce the error rate. The
results for the 3001-10000:1-3000 train:test split were 2.7% error rate
for both.
Default=Default
Default$studentDV = ifelse(Default$student == 'Yes', 1, 0)
data <- 1:10000
test <- 1:3000
train <- data[data != test]
test <- Default[1:3000,]
test.target <- default[1:3000]
glm.fit2 = glm(default~income+balance+studentDV,
data = Default,
family='binomial',
subset=train)
summary(glm.fit2)##
## Call:
## glm(formula = default ~ income + balance + studentDV, family = "binomial",
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3525 -0.1467 -0.0592 -0.0227 3.7131
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.032e+01 5.604e-01 -18.420 < 2e-16 ***
## income -3.825e-06 9.554e-06 -0.400 0.68887
## balance 5.528e-03 2.668e-04 20.717 < 2e-16 ***
## studentDV -7.534e-01 2.772e-01 -2.718 0.00658 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2037.0 on 6999 degrees of freedom
## Residual deviance: 1121.4 on 6996 degrees of freedom
## AIC: 1129.4
##
## Number of Fisher Scoring iterations: 8glm.prob=predict(glm.fit2, Default[-train, ], type = "response")
glm.pred=rep("No", length(glm.prob))
glm.pred[glm.prob >.5]="Yes"
1 - mean(glm.pred==test.target)## [1] 0.027We 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.
Q6(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.
A6(a) The estimated standard errors for from
glm() function are income = 4.985 and balance = 2.274.
set.seed(1)
glm.fit6 = glm(default~income+balance,
family = 'binomial',
data = Default)
summary(glm.fit6)##
## 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: 8Q6(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.
A6(b) Accomplished below:
boot.fn=function(data,index){
coef(glm(default ~ income+balance,
family = 'binomial',
data = data,
subset = index))
}Q6(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.
A6(c) Using the boot() function along with
boot.fn, the estimated standard errors are income =
4.8663e-06 and balance = 2.2989e-04.
set.seed(1)
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.945460e-02 4.344722e-01
## t2* 2.080898e-05 1.680317e-07 4.866284e-06
## t3* 5.647103e-03 1.855765e-05 2.298949e-04Q6(d) Comment on the estimated standard errors
obtained using the glm() function and using your bootstrap
function.
A6(d) The Bootstrap function has similar estimated
standard errors for both variables. The estimated standard errors for
from glm() function are income = 4.985 and balance = 2.274.
The bootstrap function had had negligible difference, the estimated
standard errors were income = 4.903 and balance = 2.279.
detach(Default)We will now consider the Boston housing data set, from
the ISLR2 library.
attach(Boston)Q9(a) Based on this data set, provide an estimate
for the population mean of medv. Call this estimate
ˆµ.
A9(a) the estimate for the population mean, ˆµ, is
22.53281.
mean.hat = mean(medv)
mean.hat## [1] 22.53281Q9(b) 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.
A9(b) The estimate of the standard error of ˆµ is
0.4088611.
estSE= sd(medv)/sqrt((length(medv)))
estSE## [1] 0.4088611Q9(c) Now estimate the standard error of ˆµ using
the bootstrap. How does this compare to your answer from (b)?
A9(c) The standard error of ˆµ using the bootstrap
(1,000 samples) resulted in a slightly lower standard error compared to
part b. The standard error of ˆµ in part b was .4088611, using bootstrap
it is 0.4069546.
set.seed(1)
boot.fn.mn=function(data,index){
mean(data[index])
}
b.mn = boot(medv, boot.fn.mn, R=10000)
b.mn##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot.fn.mn, R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 -0.002350771 0.4069546Q9(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 [ˆµ − 2SE(ˆµ), µˆ +
2SE(ˆµ)].
A9(d) When comparing the confidence intervals between
the two they were fairly close, the bootstrap method has a slightly
larger span than than the t.test function. Bootstrap 95% confidence
interval is 21.7189 to 23.34672, t.test function 95% confidence interval
is 21.72953 23.33608.
lower_ci = b.mn$t0 - 2*0.4069546
lower_ci## [1] 21.7189upper_ci = b.mn$t0 + 2*0.4069546
upper_ci## [1] 23.34672t.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.53281Q9(e) Based on this data set, provide an estimate,
ˆµmed, for the median value of medv in the
population.
A9(e) The estimate, ˆµmed, for the median value of
medv in the population is 21.2.
medvMed = median(medv)
medvMed## [1] 21.2Q9(f) 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.
A9(f) The results are similar to part G, with the exact
same median of 21.2. The std error established by the bootstrap is
.3700318.
set.seed(1)
boot.fn.mn=function(data,index){
median(data[index])
}
b.md = boot(medv, boot.fn.mn, R=10000)
b.md##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot.fn.mn, R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.014705 0.3780355Q9(g) 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.).
A9(g) The tenth percentile of medv in
Boston census tracts (ˆµ0.1) is 12.75.
medv10 = quantile(medv, probs = .1)
medv10## 10%
## 12.75Q9(h) Use the bootstrap to estimate the standard
error of ˆµ0.1. Comment on your findings.
A9(h) While the standard error is small at .4768, it is
slightly large compared to the size of ˆµ0.1 at 12.75.
set.seed(1)
b.se10 = function(data,index){
quantile(data[index], probs = .1)
}
boot(medv, b.se10, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = b.se10, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.0339 0.4767526