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.
a)Fit a logistic regression model that uses income and balance to predict default
library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.4.2set.seed(42)
fit <- glm(default ~ income + balance, data = Default, family = "binomial")train <- sample(nrow(Default), nrow(Default) / 2)
fit <- glm(default ~ income + balance, data = Default, family = "binomial", subset = train)
pred <- ifelse(predict(fit, newdata = Default[-train, ], type = "response") > 0.5, "Yes", "No")
table(pred, Default$default[-train])##
## pred No Yes
## No 4817 110
## Yes 20 53mean(pred != Default$default[-train])## [1] 0.026replicate(3, {
train <- sample(nrow(Default), nrow(Default) / 2)
fit <- glm(default ~ income + balance, data = Default, family = "binomial", subset = train)
pred <- ifelse(predict(fit, newdata = Default[-train, ], type = "response") > 0.5, "Yes", "No")
mean(pred != Default$default[-train])
})## [1] 0.0260 0.0294 0.0258The results obtained are variable and depend on the samples allocated to training vs. test.
replicate(3, {
train <- sample(nrow(Default), nrow(Default) / 2)
fit <- glm(default ~ income + balance + student, data = Default, family = "binomial", subset = train)
pred <- ifelse(predict(fit, newdata = Default[-train, ], type = "response") > 0.5, "Yes", "No")
mean(pred != Default$default[-train])
})## [1] 0.0278 0.0256 0.0250Including student does not seem to make a substantial improvement to the test error.
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.
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.
fit <- glm(default ~ income + balance, data = Default, family = "binomial")
summary(fit)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default)
##
## 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: 8The standard errors obtained by bootstrapping are β^1= 5.0e-6 and β^2= 2.3e-4.
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(x, i) {
fit <- glm(default ~ income + balance, data = x[i, ], family = "binomial")
coef(fit)[-1]
}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(42)
boot(Default, boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Default, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 2.080898e-05 2.737444e-08 5.073444e-06
## t2* 5.647103e-03 1.176249e-05 2.299133e-04d)Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.
The standard errors obtained by bootstrapping are similar to those estimated by glm.
In Sections 5.3.2 and 5.3.3, we saw that the cv.glm() function can be used in order to compute the LOOCV test error estimate. Alternatively, one could compute those quantities using just the glm() and predict.glm() functions, and a for loop. You will now take this approach in order to compute the LOOCV error for a simple logistic regression model on the Weekly data set. Recall that in the context of classification problems, the LOOCV error is given in (5.4).
a)Fit a logistic regression model that predicts Direction using Lag1 and Lag2.
fit <- glm(Direction ~ Lag1 + Lag2, data = Weekly, family = "binomial")b)Fit a logistic regression model that predicts Direction using Lag1 and Lag2 using all but the first observation.
fit <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-1, ], family = "binomial")c)Use the model from (b) to predict the direction of the first observation. You can do this by predicting that the first observation will go up if P(Direction=“Up” | Lag1 , Lag2 )>0.5.Was this observation correctly classified?
predict(fit, newdata = Weekly[1, , drop = FALSE], type = "response") > 0.5## 1
## TRUEYes the observation was correctly classified
i.Fit a logistic regression model using all but the ith observation to predict Direction using Lag1 and Lag2
Compute the posterior probability of the market moving up for the ith observation.
Use the posterior probability for the ith observation in order to predict whether or not the market moves up.
Determine whether or not an error was made in predicting the direction for the ith observation. If an error was made,then indicate this as a 1, and otherwise indicate it as a 0.
error <- numeric(nrow(Weekly))
for (i in 1:nrow(Weekly)) {
fit <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-i, ], family = "binomial")
p <- predict(fit, newdata = Weekly[i, , drop = FALSE], type = "response") > 0.5
error[i] <- ifelse(p, "Down", "Up") == Weekly$Direction[i]
}mean(error)## [1] 0.4499541The LOOCV test error rate is 45% which implies that our predictions are marginally more often correct than not.
We will now consider the Boston housing data set, from the ISLR2 library.
(mu <- mean(Boston$medv))## [1] 22.53281sd(Boston$medv) / sqrt(length(Boston$medv))## [1] 0.4088611set.seed(42)
(bs <- boot(Boston$medv, function(v, i) mean(v[i]), 10000))##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = function(v, i) mean(v[i]),
## R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.002175751 0.4029139The standard error using the bootstrap (0.403) is very close to that obtained from the formula above (0.409).
se <- sd(bs$t)
c(mu - 2 * se, mu + 2 * se)## [1] 21.72698 23.33863median(Boston$medv)## [1] 21.2set.seed(42)
boot(Boston$medv, function(v, i) median(v[i]), 10000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = function(v, i) median(v[i]),
## R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 -0.01331 0.3744634The estimated standard error of the median is 0.374. This is lower than the standard error of the mean.
quantile(Boston$medv, 0.1)## 10%
## 12.75set.seed(42)
boot(Boston$medv, function(v, i) quantile(v[i], 0.1), 10000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = function(v, i) quantile(v[i],
## 0.1), R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.013405 0.497298We get a standard error of ~0.5. This is higher than the standard error of the median. Nevertheless the standard error is quite small, thus we can be fairly confidence about the value of the 10th percentile.