An Introduction to Statistical Learning (ISLR) Solutions: Chapter 5

Problem 5

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.

Solution (a)

library(ISLR)
library(MASS)
attach(Default)
set.seed(1)
#i
lr<-glm(default~income+balance,family = binomial,data=Default)
summary(lr)

## 
## 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

pd<-predict(lr,Default$default,type="response")
pd.class<-ifelse(pd>0.5,"Yes","No")
round(mean(Default$default!=pd.class),4)

## [1] 0.0263

The logistic regression model is built using the entire data set and missclassification rate is calculated as 0.0263

(b)

Using 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.

Fit a multiple logistic regression model using only the training observations.

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.

Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.

Solution (b)

subset<-sample(nrow(Default),nrow(Default)*0.9)
default.train<-Default[subset,]
default.test<-Default[-subset,]
#iii
lr.90<-glm(default~income+balance,family = binomial,data=default.train)
predict.90<-predict(lr.90,default.test,type="response")
class.90<-ifelse(predict.90>0.5,"Yes","No")
#iv
table(default.test$default,class.90,dnn=c("Actual","Predicted"))

##       Predicted
## Actual  No Yes
##    No  965   6
##    Yes  19  10

round(mean(class.90!=default.test$default),4)

## [1] 0.025

The Data is splitted into training and validation sets in 50:10 ratio. The missclassification rate is calculated as 0.025

(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.

Solution (c)

set.seed(1)
subset<-sample(nrow(Default),nrow(Default)*0.8)
default.train<-Default[subset,]
default.test<-Default[-subset,]
lr.80<-glm(default~income+balance,family = binomial,data=default.train)
predict.80<-predict(lr.80,default.test,type="response")
class.80<-ifelse(predict.80>0.5,"Yes","No")
#iv
mr80<-round(mean(class.80!=default.test$default),4)


subset<-sample(nrow(Default),nrow(Default)*0.7)
default.train<-Default[subset,]
default.test<-Default[-subset,]
#iii
lr.70<-glm(default~income+balance,family = binomial,data=default.train)
predict.70<-predict(lr.70,default.test,type="response")
class.70<-ifelse(predict.70>0.5,"Yes","No")
#iv
mr70<-round(mean(class.70!=default.test$default),4)

subset<-sample(nrow(Default),nrow(Default)*0.5)
default.train<-Default[subset,]
default.test<-Default[-subset,]
#iii
lr.50<-glm(default~income+balance,family = binomial,data=default.train)
predict.50<-predict(lr.50,default.test,type="response")
class.50<-ifelse(predict.50>0.5,"Yes","No")
#iv
mr50<-round(mean(class.50!=default.test$default),4)

The data set is splitted into three different ratios and the model is fitted on training set and the misclassification rate is calculated on test set.

Missclassification rate (Train:Test–80:20): 0.026

Missclassification rate (Train:Test–70:30): 0.032

Missclassification rate (Train:Test–50:50): 0.0268

It is observed that error rate is different for different samples

(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.

Solution(d)

set.seed(1)
subset<-sample(nrow(Default),nrow(Default)*0.7)
default.train<-Default[subset,]
default.test<-Default[-subset,]
#iii
lr.70<-glm(default~income+balance+student,family = binomial,data=default.train)
summary(lr.70)

## 
## Call:
## glm(formula = default ~ income + balance + student, family = binomial, 
##     data = default.train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5184  -0.1358  -0.0522  -0.0183   3.7733  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.107e+01  5.985e-01 -18.504   <2e-16 ***
## income       2.505e-06  9.968e-06   0.251   0.8016    
## balance      5.887e-03  2.844e-04  20.700   <2e-16 ***
## studentYes  -6.805e-01  2.881e-01  -2.362   0.0182 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2050.5  on 6999  degrees of freedom
## Residual deviance: 1070.7  on 6996  degrees of freedom
## AIC: 1078.7
## 
## Number of Fisher Scoring iterations: 8

predict.70<-predict(lr.70,default.test,type="response")
class.70<-ifelse(predict.70>0.7,"Yes","No")
#iv
table(default.test$default,class.70,dnn=c("Actual","Predicted"))

##       Predicted
## Actual   No  Yes
##    No  2898    3
##    Yes   84   15

round(mean(class.70!=default.test$default),4)

## [1] 0.029

The missclassification rate is 0.0444 when dummy variable for student is included. It does not lead to a reduction in the test error rate.

Problem 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. 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.

Solution (a)

set.seed(1)
lr<-glm(default~income+balance,family = binomial,data=Default)
summary(lr)

## 
## 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

The standard error of coefficients are tabulated above for the entire data set and with the help of built in R function

(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.

Solution (b)

boot.fn<-function(data,index){
  fit<-glm(default~income+balance,data=data,family="binomial",subset=index)
  return(coef(fit))
}

The function boot.fn is defind so that it takes in the data and gives the coefficient estimate for the given 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.

Solution (c)

library(boot)
boot(Default,boot.fn,100)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Default, statistic = boot.fn, R = 100)
## 
## 
## Bootstrap Statistics :
##          original        bias     std. error
## t1* -1.154047e+01  9.699111e-02 4.101121e-01
## t2*  2.080898e-05  6.715005e-08 4.127740e-06
## t3*  5.647103e-03 -5.733883e-05 2.105660e-04

The standard error of coefficient estimate is calculated using library boot and boot.fn function defined above

(d)

Comment on the estimated standard errors obtained using the glm() function and using your bootstrap function.

Solution (d)

The standard error of coefficient estimates found from two methods are pretty close

Problem 7

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 regressionmodel that predicts Direction using Lag1 and Lag2.

Solution (a)

attach(Weekly)
set.seed(1)
weekly<-glm(Direction~Lag1+Lag2,family = binomial,data=Weekly)
summary(weekly)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.623  -1.261   1.001   1.083   1.506  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.22122    0.06147   3.599 0.000319 ***
## Lag1        -0.03872    0.02622  -1.477 0.139672    
## Lag2         0.06025    0.02655   2.270 0.023232 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1496.2  on 1088  degrees of freedom
## Residual deviance: 1488.2  on 1086  degrees of freedom
## AIC: 1494.2
## 
## Number of Fisher Scoring iterations: 4

(b)

Fit a logistic regressionmodel that predicts Direction using Lag1 and Lag2 using all but the first observation.

Solution (b)

weekly1<-glm(Direction~Lag1+Lag2,family = binomial,data=Weekly[-1,])

(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?

Solution (c)

predictweekly<-predict(weekly,Weekly[1,],type = "response")
predictweekly.class<-ifelse(predictweekly>0.5,"Up","Down")
predictweekly.class

##    1 
## "Up"

The predicted direction from the model for the first observation is Up.

(d)

Write a for loop from i = 1 to i = n, where n is the number of observations in the data set, that performs each of the following steps:

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.

Solution (d)

error<-rep(0,dim(Weekly)[1])
for (i in 1:dim(Weekly)[1]){
  fit.glm<-fit.glm <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-i, ],  family = "binomial")
  pred.up <- predict.glm(fit.glm, Weekly[i, ], type = "response") > 0.5
  true.up <- Weekly[i, ]$Direction == "Up"
  if (pred.up != true.up)
    error[i] <- 1
}

The above for loop predicts and codes 1 if error is made in prediction and 0 otherwise

(e)

Take the average of the n numbers obtained in (d)iv in order to obtain the LOOCV estimate for the test error. Comment on the results.

Solution

mean(error)

## [1] 0.4499541

LOOCV estimate for the test error is 45 %

Problem 8

We will now perform cross-validation on a simulated data set.

(a)

Generate a simulated data set as follows:

set .seed (1)

y=rnorm (100)

x=rnorm (100)

y=x-2* x^2+ rnorm (100)

In this data set, what is n and what is p? Write out the model used to generate the data in equation form.

Solution (a)

set.seed(1)
y=rnorm(100)
x<-rnorm(100)
y=x-2*x^2+rnorm(100)

Here the n=100 and p=3 (including intercept).

the equation is $Y=B~0~ + B~1~*(X) + B~2~*(X)^2 + error$

(b)

Create a scatterplot of X against Y . Comment on what you find.

Solution (b)

library(ggplot2)
simulateddata<-data.frame(x,y)
ggplot(simulateddata,aes(x=x,y=y))+
  geom_point(shape=8)+
  geom_smooth()

## `geom_smooth()` using method = 'loess'

From the plot it can be inferred that there exists a quadratic relationship between x and y variables

(c)

Set a random seed, and then compute the LOOCV errors that result from fitting the following four models using least squares:

$Y = \beta0 + \beta1X + error$
$Y = \beta0 + \beta1X + \beta2X2 + error$
$Y = \beta0 + \beta1X + \beta2X2 + \beta3X3 + error$
$Y = \beta0 + \beta1X + \beta2X2 + \beta3X3 + \beta4X4 + error$

Solution (c)

#i
fit.glm1<-glm(y~x)
cv.glm(simulateddata,fit.glm1)$delta[1]

## [1] 5.890979

#ii
fit.glm2<-glm(y~poly(x,2))
cv.glm(simulateddata,fit.glm2)$delta[1]

## [1] 1.086596

#iii
fit.glm3<-glm(y~poly(x,3))
cv.glm(simulateddata,fit.glm3)$delta[1]

## [1] 1.102585

#iv
fit.glm4<-glm(y~poly(x,4))
cv.glm(simulateddata,fit.glm4)$delta[1]

## [1] 1.114772

LOOCV error for $Y = \beta0 + \beta1X + error$ is 5.8909786

LOOCV error for $Y = \beta0 + \beta1X + \beta2X2 + error$ is 1.0865956

LOOCV error for $Y = \beta0 + \beta1X + \beta2X2 + \beta3X3 + error$ is 1.1025851

LOOCV error for $Y = \beta0 + \beta1X + \beta2X2 + \beta3X3 + \beta4X4 + error$ is 1.1147723

The quadratic relationship was observed from the scatter plot between x and y variables. The cross-validation error is least for the quadratic function $Y = \beta0 + \beta1X + \beta2X2 + error$ and is obtained as 1.0865956

(d)

Repeat (c) using another random seed, and report your results. Are your results the same as what you got in (c)? Why?

Solution (d)

set.seed(2)
#i
fit.glm1.2<-glm(y~x)
cv.glm(simulateddata,fit.glm1.2)$delta[1]

## [1] 5.890979

#ii
fit.glm2.2<-glm(y~poly(x,2))
cv.glm(simulateddata,fit.glm2.2)$delta[1]

## [1] 1.086596

#iii
fit.glm3.2<-glm(y~poly(x,3))
cv.glm(simulateddata,fit.glm3.2)$delta[1]

## [1] 1.102585

#iv
fit.glm4.2<-glm(y~poly(x,4))
cv.glm(simulateddata,fit.glm4.2)$delta[1]

## [1] 1.114772

The result obtained with different seed is very close to previous result. This was expected as the model is built using the entire data set but one observation

(e)

Which of the models in (c) had the smallest LOOCV error? Is this what you expected? Explain your answer.

Solution (e)

The The cross-validation error is least for the quadratic function $Y = \beta0 + \beta1X + \beta2X2 + error$ and is obtained as 1.0865956

This was expected, as the plot showed quadratic relationship between x and y.

(f)

Comment on the statistical significance of the coefficient estimates that results from fitting each of the models in (c) using least squares. Do these results agree with the conclusions drawn based on the cross-validation results?

Solution (f)

summary(fit.glm4)

## 
## Call:
## glm(formula = y ~ poly(x, 4))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8914  -0.5244   0.0749   0.5932   2.7796  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.8277     0.1041 -17.549   <2e-16 ***
## poly(x, 4)1   2.3164     1.0415   2.224   0.0285 *  
## poly(x, 4)2 -21.0586     1.0415 -20.220   <2e-16 ***
## poly(x, 4)3  -0.3048     1.0415  -0.293   0.7704    
## poly(x, 4)4  -0.4926     1.0415  -0.473   0.6373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.084654)
## 
##     Null deviance: 552.21  on 99  degrees of freedom
## Residual deviance: 103.04  on 95  degrees of freedom
## AIC: 298.78
## 
## Number of Fisher Scoring iterations: 2

From the summary of the model it is observed that only the linear and quadratic variables are significant in the model. The results of cross validation pointed out that least error is obtained for $Y = \beta0 + \beta1X + \beta2X2 + error$

Problem 9

We will now consider the Boston housing data set, from the MASS library.

(a)

Based on this data set, provide an estimate for the population mean of medv. Call this estimate .

Solution (a)

library(MASS)
attach(Boston)
sample_mu<-mean(medv)
sample_mu

## [1] 22.53281

The sample mean is 22.5328063

(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.

Solution (b)

sample_se<-sd(medv)/sqrt(nrow(Boston))
sample_se

## [1] 0.4088611

The estimate of the standard error of is 0.4088611

The data set is the sample representing the population of Boston. The standard error gives us the accuracy of the estimate i.e. how much the mean will vary if different sample was chosen

(c)

Now estimate the standard error of using the bootstrap. How does this compare to your answer from (b)?

library(boot)
set.seed(1)
boot.fn<-function(data,index){
  mu<-mean(data[index])
  return(mu)
}
set.seed(1)
boot(medv,boot.fn,100)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = medv, statistic = boot.fn, R = 100)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 22.53281 -0.04836957   0.3815554

The standard error obtained from bootstrapping is close to the one obtained earlier

(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$ - 2SE($\mu$), $\mu$ + 2SE($\mu$)]

Solution (d)

t.test(medv)

## 
##  One Sample t-test
## 
## data:  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(sample_mu-2*0.3815554,sample_mu+2*0.3815554)
ci.mu

## [1] 21.76970 23.29592

The confidence interval obtained is very similar to the one obtained from t test

(e)

Based on this data set, provide an estimate, $\mu$ med, for the median value of medv in the population.

Solution (e)

mu_median<-median(medv)
mu_median

## [1] 21.2

The estimate for the median is 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.

Solution (f)

library(boot)
set.seed(1)
boot.fn<-function(data,index){
  mu.median<-median(data[index])
  return(mu.median)
}
set.seed(1)
boot(medv,boot.fn,100)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = medv, statistic = boot.fn, R = 100)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     21.2  -0.026   0.3683488

Here we computed the standard error of median using bootstrapping. This showcases how easily we can implement bootstrapping concept to find standard error when statistical softwares do no directly give us the value

(g)

Based on this data set, provide an estimate for the tenth percentile of medv in Boston suburbs. Call this quantity $\mu$ 0.1. (You can use the quantile() function.)

Solution (g)

mu.quantile<-quantile(medv,c(0.1))
mu.quantile

##   10% 
## 12.75

The estimate for the tenth percentile of medv in Boston suburbs is obtaines as 12.75

(h)

Use the bootstrap to estimate the standard error of 0.1. Comment on your findings.

Solution (h)

boot.fn <- function(data, index) {
  mu.quantile <- quantile(data[index], c(0.1))
  return (mu.quantile)
}
boot(medv, boot.fn, 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    12.75 0.01705   0.5032246

The standard error is obtained using bootstrapping and the value is small compared to the tenth percentile which means the estimate is representing the population quite accurately.

An Introduction to Statistical Learning (ISLR) Solutions: Chapter 5

Swapnil Sharma

July 22, 2017

Chapter 5 Resampling:

Cross Validation & Bootstrapping

Problem 5

(a)

Solution (a)

(b)

Solution (b)

(c)

Solution (c)

(d)

Solution(d)

Problem 6

(a)

Solution (a)

(b)

Solution (b)

(c)

Solution (c)

(d)

Solution (d)

Problem 7

(a)

Solution (a)

(b)

Solution (b)

(c)

Solution (c)

(d)

Solution (d)

(e)

Solution

Problem 8

(a)

Solution (a)

(b)

Solution (b)

(c)

Solution (c)

(d)

Solution (d)

(e)

Solution (e)

(f)

Solution (f)

Problem 9

(a)

Solution (a)

(b)

Solution (b)

(c)

(d)

Solution (d)

(e)

Solution (e)

(f)

Solution (f)

(g)

Solution (g)

(h)

Solution (h)