Question 3. Reviewing k-fold cross-validation
a) Explaining how k-fold cross-validation is implemented.
The data is split into k groups. Within each group, the data is
further split into training data and validation date. A model is fit to
the training data set and evaluated against the validation data. Retain
the MSE and repeat the process k more times, taking out a different part
of the data each time. You then average the K number of MSEs to obtain
estimated test error rate.
b) What are the advantages and disadvantages of k-fold cross
validation in relation to:
i) the validation set approach?
Validation set is easy to use but validation MSE is highly variable.
Methods do not perform well when trained on a smaller set of
observations.
ii) LOOCV?
LOOCV has lower bias but higher variability. LOOCV is
computationally intensive.
Question 5. We will estimate test error of logistic regression model
using validation set approach.
a) Fit a logistic regression model that uses income and balance to
predict default.
library(ISLR2)
attach(Default)
set.seed(1)
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: 8
b) Using validation set approach, estimate the test error of this
model.
0.0244
train= sample(dim(Default)[1], 0.5*dim(Default)[1])
glm.fit = glm(default~income + balance, data = Default, subset=train, family=binomial)
glm.prob = predict(glm.fit, Default[-train,], type = "response")
glm.pred = rep("No", dim(Default)[1])
glm.pred[glm.prob > 0.5] = "Yes"
mean(glm.pred != Default[-train,]$default)
## [1] 0.0254
c) Repeat process (b) three times, using three different splits of
the observations into a training set and validation set. Comment on
results obtained.
There appears to be some minor variability among the validation
estimate of the test error rates.
train= sample(dim(Default)[1], 0.5*dim(Default)[1])
glm.fit = glm(default~income + balance, data = Default, subset=train, family=binomial)
glm.prob = predict(glm.fit, Default[-train,], type = "response")
glm.pred = rep("No", dim(Default)[1])
glm.pred[glm.prob > 0.5] = "Yes"
mean(glm.pred != Default[-train,]$default)
## [1] 0.0274
train= sample(dim(Default)[1], 0.5*dim(Default)[1])
glm.fit = glm(default~income + balance, data = Default, subset=train, family=binomial)
glm.prob = predict(glm.fit, Default[-train,], type = "response")
glm.pred = rep("No", dim(Default)[1])
glm.pred[glm.prob > 0.5] = "Yes"
mean(glm.pred != Default[-train,]$default)
## [1] 0.0244
train= sample(dim(Default)[1], 0.5*dim(Default)[1])
glm.fit = glm(default~income + balance, data = Default, subset=train, family=binomial)
glm.prob = predict(glm.fit, Default[-train,], type = "response")
glm.pred = rep("No", dim(Default)[1])
glm.pred[glm.prob > 0.5] = "Yes"
mean(glm.pred != Default[-train,]$default)
## [1] 0.0244
d) Consider a logistic regression model that predict 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.
The addition of the dummy variable does not significantly change the
test error rate.
train= sample(dim(Default)[1], 0.5*dim(Default)[1])
glm.fit = glm(default~income + balance + student, data = Default, subset=train, family=binomial)
glm.prob = predict(glm.fit, Default[-train,], type = "response")
glm.pred = rep("No", dim(Default)[1])
glm.pred[glm.prob > 0.5] = "Yes"
mean(glm.pred != Default[-train,]$default)
## [1] 0.0278
Question 6. 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. Set random seed before beginning analysis.
a) Using the summary() and glm() functions, determine the estimated
standard errors for the coefficients associated with income and balance
in a multiple regression model that uses both predictors.
The standard error estimates are .617, 7.02 10^(-6), and
3.1510^(-4).
attach(Default)
## The following objects are masked from Default (pos = 3):
##
## balance, default, income, student
set.seed(1)
train = sample(dim(Default)[1],0.5*dim(Default)[1])
glm.fit = glm(default ~ income + balance, data = Default, family = binomial, subset = train)
summary(glm.fit)
##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5830 -0.1428 -0.0573 -0.0213 3.3395
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.194e+01 6.178e-01 -19.333 < 2e-16 ***
## income 3.262e-05 7.024e-06 4.644 3.41e-06 ***
## balance 5.689e-03 3.158e-04 18.014 < 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: 1523.8 on 4999 degrees of freedom
## Residual deviance: 803.3 on 4997 degrees of freedom
## AIC: 809.3
##
## Number of Fisher Scoring iterations: 8
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)
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-04
d) Comment on the estimated standard errors obtained using the glm()
function and using your bootstrap function.
The standard error estimates for the linear regression model were
.617, 7.0210^(-6), and 3.1510^(-4). The standard error
estimates for the bootstrap function are .434, 4.8610^(-6), and
2.2910^(-4). There appear to be some variances across the
estimates.
Question 9. We will now be using the Boston housing data set.
a) Based on this data set, provide an estimate for the population
mean of medv.
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:ISLR2':
##
## Boston
attach(Boston)
mu.pop = mean(Boston$medv)
mu.pop
## [1] 22.53281
b) Provide an estimate of the standard error of (mu hat). Interpret
the result. (We can compute the standard error of the sample mean by
dividing the same standard deviation by the square root of the number of
observations.)
st.err = sd(Boston$medv)/sqrt(dim(Boston)[1])
st.err
## [1] 0.4088611
c) Now estimate the standard error of (mu hat) using the bootstrap.
How does this compare to part b?
The estimated standard errors appear to be similar.
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.4106622
d) Based on bootstrap estimate from part (c), provide a 95%
confidence interval for the mean of medv. Compare it to the results
obtained using t.test(Boston$medv).
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.53281
CI.mu = c(22.53-2*0.4119, 22.53+2*0.4119)
CI.mu
## [1] 21.7062 23.3538
g) Based on this data set, provide an estimate for the tenth
percentile of medv in Boston data, (use the quantile function).
quantile(Boston$medv, c(0.1))
## 10%
## 12.75
h) Use the bootstrap method to estimate the standard error of (mu
0.1). Comment on your findings.
The estimated tenth percentile value of 22.5, which does not equal
the value obtained in (g) with a standard error of .4024.
boot.fun = function(data, index) {
mu = quantile(data[index], c(0.1))
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.02542806 0.4023801