3(a). Explain how k-fold cross-validation is implemented. With k-fold cross-validation you first divide the data into different K parts. The first fold is then treated as a validation set with the remaining methods being fit onto the folds. MSE is then computed from the held out fold. This is repeated k number of times with a different group used as the validation set. We then average all the different MSE’s from the different folds. This results in our k-fold cross-validation.
(b). What are the advantages and disadvantages of k-fold cross-validation relative to: (i). The validation set approach? A validation set has the potential to over estiamte the test error rate while k-fold cross-validation can give us a more accurate rate, given the number of observations.
(ii). LOOCV? The advantage of k-fold cross-validation relative to LOOCV is that is a lot more computationally intuitive while also producing more accurate estimates than the LOOCV. There is little to suggest that LOOCV has a strong advantage over K-fold cross-validation.
head(Default)## default student balance income
## 1 No No 729.5265 44361.625
## 2 No Yes 817.1804 12106.135
## 3 No No 1073.5492 31767.139
## 4 No No 529.2506 35704.494
## 5 No No 785.6559 38463.496
## 6 No Yes 919.5885 7491.559str(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 ...(a).Fit a logistic regression model that uses income and balance to predict default.
set.seed(1)
deafult.glm = glm(default ~ income + balance, data = Default, family = 'binomial')
summary(deafult.glm)##
## 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).
(i). Split the sample set into a training set and a validation set.
index = sample(nrow(Default), nrow(Default)*.7)
default.train = Default[index, ]
default.test = Default[-index, ](ii).Fit a multiple logistic regression model using only the training observations
default.glm7 = glm(default ~ income + balance, data = default.train, family = 'binomial')
summary(default.glm7)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = default.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4481 -0.1402 -0.0561 -0.0211 3.3484
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.167e+01 5.214e-01 -22.379 < 2e-16 ***
## income 2.560e-05 6.012e-06 4.258 2.06e-05 ***
## balance 5.574e-03 2.678e-04 20.816 < 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: 2030.3 on 6999 degrees of freedom
## Residual deviance: 1079.6 on 6997 degrees of freedom
## AIC: 1085.6
##
## Number of Fisher Scoring iterations: 8(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.
default.glm.test = factor(ifelse(predict(default.glm7, newdata = default.test,
type = 'response')> .5, 'Yes', 'No'))
default.glm.test[1:16]## 1 10 11 12 13 17 20 23 25 28 31 33 35 37 42 43
## No No No No No No No No No No No No No No No No
## Levels: No Yes(IV).Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
round(mean(default.test$default != default.glm.test), 5)## [1] 0.02667(c).
Using 80/20 split with a valadation set error of .0265
index = sample(nrow(Default), nrow(Default)*.8)
default.train = Default[index, ]
default.test = Default[-index, ]default.glm8 = glm(default ~ income + balance, data = default.train, family = 'binomial')
summary(default.glm8)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = default.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4958 -0.1403 -0.0552 -0.0198 3.7427
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.164e+01 4.908e-01 -23.725 < 2e-16 ***
## income 2.006e-05 5.546e-06 3.616 0.000299 ***
## balance 5.731e-03 2.592e-04 22.111 < 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: 2340.6 on 7999 degrees of freedom
## Residual deviance: 1252.2 on 7997 degrees of freedom
## AIC: 1258.2
##
## Number of Fisher Scoring iterations: 8default.glm.test = factor(ifelse(predict(default.glm8, newdata = default.test,
type = 'response')> .5, 'Yes', 'No'))
default.glm.test[1:16]## 5 7 9 11 14 17 23 26 33 42 43 53 54 58 59 61
## No No No No No No No No No No No No No No No No
## Levels: No Yesround(mean(default.test$default != default.glm.test), 5)## [1] 0.0265Using 60/40 split with a valadation error rate of .027.
index = sample(nrow(Default), nrow(Default)*.6)
default.train = Default[index, ]
default.test = Default[-index, ]default.glm6 = glm(default ~ income + balance, data = default.train, family = 'binomial')
summary(default.glm6)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = default.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9851 -0.1300 -0.0479 -0.0167 3.4056
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.214e+01 5.997e-01 -20.252 < 2e-16 ***
## income 1.712e-05 6.472e-06 2.645 0.00818 **
## balance 6.110e-03 3.175e-04 19.244 < 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: 1780.59 on 5999 degrees of freedom
## Residual deviance: 910.29 on 5997 degrees of freedom
## AIC: 916.29
##
## Number of Fisher Scoring iterations: 8default.glm.test = factor(ifelse(predict(default.glm6, newdata = default.test,
type = 'response')> .5, 'Yes', 'No'))
default.glm.test[1:16]## 1 2 3 4 8 9 11 12 13 14 19 20 26 27 33 34
## No No No No No No No No No No No No No No No No
## Levels: No Yesround(mean(default.test$default != default.glm.test), 5)## [1] 0.027(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.
After running it with the variable for student using a 70/30 split we actually produced a higher error rate of .0275.
index = sample(nrow(Default), nrow(Default)*.7)
default.train = Default[index, ]
default.test = Default[-index, ]default.glm.dum = glm(default ~ income + balance + student, data = default.train, family = 'binomial')
summary(default.glm.dum)##
## Call:
## glm(formula = default ~ income + balance + student, family = "binomial",
## data = default.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1889 -0.1329 -0.0514 -0.0181 3.7918
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.111e+01 6.099e-01 -18.217 <2e-16 ***
## income 1.680e-06 1.016e-05 0.165 0.869
## balance 5.893e-03 2.863e-04 20.579 <2e-16 ***
## studentYes -6.762e-01 2.908e-01 -2.325 0.020 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2016.8 on 6999 degrees of freedom
## Residual deviance: 1049.0 on 6996 degrees of freedom
## AIC: 1057
##
## Number of Fisher Scoring iterations: 8default.glm.test = factor(ifelse(predict(default.glm.dum, newdata = default.test,
type = 'response')> .5, 'Yes', 'No'))
default.glm.test[1:16]## 3 10 13 18 19 23 28 31 32 35 40 42 45 47 53 54
## No No No No No No No No No No No No No No No No
## Levels: No Yesround(mean(default.test$default != default.glm.test), 5)## [1] 0.02833(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.default = glm(default ~ income + balance, data = Default, family = 'binomial')
summary(glm.default)##
## 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)
{coef(glm(default~income+balance, data = data, family = 'binomial',
subset = index))}
set.seed(1)
boot.fn(Default,sample(1000,500,replace = T))## (Intercept) income balance
## -1.029272e+01 1.558656e-05 5.236697e-03(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.
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.933672e-02 4.343265e-01
## t2* 2.080898e-05 1.671659e-07 4.867306e-06
## t3* 5.647103e-03 1.850346e-05 2.298227e-04(d).
The standard errors obtained using the glm() function and using the bootstap function are fairly close.
(a). Based on this data set, provide an estimate for the population mean of medv. Call this estimate ˆμ.
attach(Boston)
mean.medv = mean(medv)
mean.medv## [1] 22.53281(b). Provide an estimate of the standard error of ˆμ. Interpret this result.
stdr.er = sd(medv)/ sqrt(nrow(Boston))
stdr.er## [1] 0.4088611(c). Now estimate the standard error of ˆμ using the bootstrap. How does this compare to your answer from (b)? For the mean we return .408 while using bootstrap we return .41. The errors are very close.
boot.fn = function(vector, index) {
mean(vector[index])
}
set.seed(1)
(boot.result = boot(medv, boot.fn, 1000))##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.007650791 0.4106622(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). Again, the the results are very similar.
boot.result.95 = sd(boot.result$t)
round(c(mean(medv)-2*boot.result.95,mean(medv)+2*boot.result.95),4)## [1] 21.7115 23.3541t.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(e). Based on this data set, provide an estimate, ˆμmed, for the median value of medv in the population.
median(medv)## [1] 21.2(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.
Here we return a smaller error of .378.
boot.fn = function(vector, index) {
median(vector[index])
}
set.seed(1)
(boot.results.er = boot(medv, boot.fn, 1000))##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 0.02295 0.3778075(g). Based on this data set, provide an estimate for the tenth percentile of medv in Boston suburbs. Call this quantity ˆμ0.1. (You can use the quantile() function.)
quantile(medv, .1)## 10%
## 12.75(h). Use the bootstrap to estimate the standard error of ˆμ0.1. Comment on your findings.
Using quantile we return a higher error similar to that using the mean.
boot.fn = function(vector, index) {
quantile(vector[index], .1)
}
set.seed(1)
(boot.result.quan = 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.0339 0.4767526