<jupyter_start><jupyter_text>Lab: Cross-Validation and the BootstrapIn this lab, we explore the resampling techniques covered in this chapter.Some of the commands in this lab may take a while to run on your computer 1 The Validation Set ApproachWe explore the use of the validation set approach in order to estimate thetest error rates that result from fitting various linear models on the Autodata set.Before we begin, we use the set.seed() function in order to set a seed forR’s random number generator, so that the reader of this book will obtainprecisely the same results as those shown below. It is generally a good ideato set a random seed when performing an analysis such as cross-validationthat contains an element of randomness, so that the results obtained canbe reproduced precisely at a later time.We begin by using the sample() function to split the set of observationsinto two halves, by selecting a random subset of 196 observations out ofthe original 392 observations. We refer to these observations as the trainingset.<jupyter_code># install.packages('ISLR', repos='http://cran.us.r-project.org')
library(ISLR)
set.seed(1)
train = sample(392,196)<jupyter_output><empty_output><jupyter_text>(Here we use a shortcut in the sample command; see ?sample for details.)We then use the subset option in lm() to fit a linear regression using onlythe observations corresponding to the training set.<jupyter_code>lm.fit =lm(mpg~horsepower ,data=Auto , subset = train )<jupyter_output><empty_output><jupyter_text>We now use the predict() function to estimate the response for all 392observations, and we use the mean() function to calculate the MSE of the196 observations in the validation set. Note that the -train index belowselects only the observations that are not in the training set<jupyter_code>attach(Auto)
mean((mpg - predict(lm.fit ,Auto))[-train]^2)<jupyter_output><empty_output><jupyter_text>Therefore, the estimated test MSE for the linear regression fit is 23.27. Wecan use the poly() function to estimate the test error for the quadraticand cubic regressions.<jupyter_code>lm.fit2=lm(mpg~poly(horsepower,2) ,data=Auto , subset =train )
mean((mpg - predict(lm.fit2 ,Auto))[-train]^2)
lm.fit3=lm(mpg~poly(horsepower,3) ,data=Auto , subset =train )
mean((mpg - predict(lm.fit3 ,Auto))[-train]^2)<jupyter_output><empty_output><jupyter_text>These error rates are 18.72 and 18.79, respectively. If we choose a differenttraining set instead, then we will obtain somewhat different errors on thevalidation set<jupyter_code>set.seed(2)
train = sample(392,196)
lm.fit =lm(mpg~horsepower, subset = train )
mean((mpg - predict(lm.fit ,Auto))[-train]^2)
mean((mpg - predict(lm.fit2 ,Auto))[-train]^2)
mean((mpg - predict(lm.fit3 ,Auto))[-train]^2)<jupyter_output><empty_output><jupyter_text>Using this split of the observations into a training set and a validationset, we find that the validation set error rates for the models with linear,quadratic, and cubic terms are 25.73, 19.95, and 19.98, respectively.These results are consistent with our previous findings: a model thatpredicts mpg using a quadratic function of horsepower performs better thana model that involves only a linear function of horsepower, and there islittle evidence in favor of a model that uses a cubic function of horsepower.  2. Leave-One-Out Cross-ValidationThe LOOCV estimate can be automatically computed for any generalizedlinear model using the glm() and cv.glm() functions. In the lab for Chapter 4, we used the glm() function to perform logistic regression by passingin the family="binomial" argument. But if we use glm() to fit a modelwithout passing in the family argument, then it performs linear regression,just like the lm() function. So for instance,<jupyter_code>glm.fit =glm(mpg~horsepower , data= Auto)
coef(glm.fit )
# and
lm.fit =lm(mpg~horsepower , data= Auto)
coef(lm.fit )<jupyter_output><empty_output><jupyter_text>yield identical linear regression models. In this lab, we will perform linearregression using the glm() function rather than the lm() function becausethe former can be used together with cv.glm(). The cv.glm() function ispart of the boot library<jupyter_code>library(boot)
glm.fit =glm(mpg~horsepower, data= Auto)
cv.err =cv.glm(Auto, glm.fit)
cv.err$delta<jupyter_output><empty_output><jupyter_text>The cv.glm() function produces a list with several components. The twonumbers in the delta vector contain the cross-validation results. In this case the numbers are identical (up to two decimal places) and correspondto the LOOCV statistic given in (5.1). Below, we discuss a situation inwhich the two numbers differ. Our cross-validation estimate for the testerror is approximately 24.23.We can repeat this procedure for increasingly complex polynomial fits.To automate the process, we use the for() function to initiate a for loopwhich iteratively fits polynomial regressions for polynomials of order i = 1 for loopto i = 5, computes the associated cross-validation error, and stores it inthe ith element of the vector cv.error. We begin by initializing the vector.This command will likely take a couple of minutes to run.<jupyter_code>cv.error =rep(0 ,5)
for(i in 1:5){
    glm.fit =glm(mpg~poly(horsepower, i),data=Auto)
    cv.error[i]=cv.glm(Auto,glm.fit)$delta[1]
}
cv.error<jupyter_output><empty_output><jupyter_text>we see a sharp drop in the estimated test MSE betweenthe linear and quadratic fits, but then no clear improvement from usinghigher-order polynomials.  3. k-Fold Cross-ValidationThe cv.glm() function can also be used to implement k-fold CV. Below weuse k = 10, a common choice for k, on the Auto data set. We once again seta random seed and initialize a vector in which we will store the CV errorscorresponding to the polynomial fits of orders one to ten.<jupyter_code>set.seed(17)
cv.error.10= rep(0 ,10)
for(i in 1:10){
    glm.fit=glm(mpg~poly(horsepower, i), data=Auto)
    cv.error.10[i] = cv.glm(Auto, glm.fit,K=10)$delta[1]
}
cv.error.10<jupyter_output><empty_output><jupyter_text>Notice that the computation time is much shorter than that of LOOCV.(In principle, the computation time for LOOCV for a least squares linearmodel should be faster than for k-fold CV, due to the availability of theformula (5.2) for LOOCV; however, unfortunately the cv.glm() functiondoes not make use of this formula.) We still see little evidence that usingcubic or higher-order polynomial terms leads to lower test error than simplyusing a quadratic fit. We saw in Section 5.3.2 that the two numbers associated with delta areessentially the same when LOOCV is performed. When we instead performk-fold CV, then the two numbers associated with delta differ slightly. The first is the standard k-fold CV estimate, as in (5.3). The second is a biascorrected version. On this data set, the two estimates are very similar to each other.  4. The BootstrapWe illustrate the use of the bootstrap in the simple example of Section 5.2,as well as on an example involving estimating the accuracy of the linearregression model on the Auto data set.  Estimating the Accuracy of a Statistic of InterestOne of the great advantages of the bootstrap approach is that it can beapplied in almost all situations. No complicated mathematical calculationsare required. Performing a bootstrap analysis in R entails only two steps.First, we must create a function that computes the statistic of interest.Second, we use the boot() function, which is part of the boot library, toperform the bootstrap by repeatedly sampling observations from the dataset with replacement.The Portfolio data set in the ISLR package is described in Section 5.2.To illustrate the use of the bootstrap on this data, we must first createa function, alpha.fn(), which takes as input the (X, Y ) data as well asa vector indicating which observations should be used to estimate α. Thefunction then outputs the estimate for α based on the selected observations.<jupyter_code>alpha.fn=function(data,index){
    X = data$X[index]
    Y = data$Y[index]
    return((var(Y)-cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y)))
}<jupyter_output><empty_output><jupyter_text>This function returns, or outputs, an estimate for α based on applying(5.7) to the observations indexed by the argument index. For instance, thefollowing command tells R to estimate α using all 100 observations.<jupyter_code>alpha.fn(Portfolio ,1:100)<jupyter_output><empty_output><jupyter_text>The next command uses the sample() function to randomly select 100 observations from the range 1 to 100, with replacement. This is equivalentto constructing a new bootstrap data set and recomputing ˆ α based on thenew data set.<jupyter_code>set.seed(1)
alpha.fn(Portfolio ,sample(100 ,100 , replace =T))<jupyter_output><empty_output><jupyter_text>We can implement a bootstrap analysis by performing this command manytimes, recording all of the corresponding estimates for α, and computing the resulting standard deviation. However, the boot() function automatesthis approach. Below we produce R = 1, 000 bootstrap estimates for α.<jupyter_code>boot(Portfolio,alpha.fn,R=1000)<jupyter_output><empty_output><jupyter_text>The final output shows that using the original data, ˆα = 0.5758, and thatthe bootstrap estimate for SE(ˆα) is 0.0937.  Estimating the Accuracy of a Linear Regression ModelThe bootstrap approach can be used to assess the variability of the coefficient estimates and predictions from a statistical learning method. Herewe use the bootstrap approach in order to assess the variability of theestimates for β0 and β1, the intercept and slope terms for the linear regression model that uses horsepower to predict mpg in the Auto data set. Wewill compare the estimates obtained using the bootstrap to those obtainedusing the formulas for SE(βˆ0) and SE(βˆ1) described in Section 3.1.2.We first create a simple function, boot.fn(), which takes in the Auto dataset as well as a set of indices for the observations, and returns the interceptand slope estimates for the linear regression model. We then apply thisfunction to the full set of 392 observations in order to compute the estimates of β0 and β1 on the entire data set using the usual linear regressioncoefficient estimate formulas from Chapter 3. Note that we do not need the{ and } at the beginning and end of the function because it is only one linelong.<jupyter_code>boot.fn= function (data ,index)
    return (coef(lm(mpg~horsepower , data=data , subset = index )))
boot.fn(Auto, 1:392)<jupyter_output><empty_output><jupyter_text>The boot.fn() function can also be used in order to create bootstrap estimates for the intercept and slope terms by randomly sampling from amongthe observations with replacement. Here we give two examples.<jupyter_code>set.seed(1)
boot.fn(Auto,sample(392,392,replace=T))
boot.fn(Auto,sample(392,392, replace=T))<jupyter_output><empty_output><jupyter_text>Next, we use the boot() function to compute the standard errors of 1,000bootstrap estimates for the intercept and slope terms.<jupyter_code>boot(Auto,boot.fn ,1000)<jupyter_output><empty_output><jupyter_text>This indicates that the bootstrap estimate for SE(βˆ0) is 0.84, and thatthe bootstrap estimate for SE(βˆ1) is 0.0073. As discussed in Section 3.1.2,standard formulas can be used to compute the standard errors for theregression coefficients in a linear model. These can be obtained using thesummary() function.<jupyter_code>summary(lm(mpg~horsepower ,data =Auto))$coef<jupyter_output><empty_output><jupyter_text>The standard error estimates for βˆ0 and βˆ1 obtained using the formulasfrom Section 3.1.2 are 0.717 for the intercept and 0.0064 for the slope.Interestingly, these are somewhat different from the estimates obtainedusing the bootstrap. Does this indicate a problem with the bootstrap? Infact, it suggests the opposite. Recall that the standard formulas given inEquation 3.8 on page 66 rely on certain assumptions. For example, theydepend on the unknown parameter σ2, the noise variance. We then estimateσ2 using the RSS. Now although the formula for the standard errors do notrely on the linear model being correct, the estimate for σ2 does. We see inFigure 3.8 on page 91 that there is a non-linear relationship in the data, andso the residuals from a linear fit will be inflated, and so will ˆ σ2. Secondly,the standard formulas assume (somewhat unrealistically) that the xi arefixed, and all the variability comes from the variation in the errors $\epsilon$i. Thebootstrap approach does not rely on any of these assumptions, and so it islikely giving a more accurate estimate of the standard errors of βˆ0 and βˆ1than is the summary() function.Below we compute the bootstrap standard error estimates and the standard linear regression estimates that result from fitting the quadratic modelto the data. Since this model provides a good fit to the data (Figure 3.8),there is now a better correspondence between the bootstrap estimates andthe standard estimates of SE(βˆ0), SE(βˆ1) and SE(βˆ2).<jupyter_code>boot.fn= function (data ,index )
    coefficients(lm(mpg~horsepower+I(horsepower^2),
                    data=data, subset=index))
set.seed(1)
boot(Auto,boot.fn,1000)
summary(lm(mpg~horsepower +I(horsepower^2) ,data= Auto))$coef<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06737000753579503, 'fraction_numerical': 0.01846269781461944, 'mean_word_length': 5.48314606741573, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '334617', 'license': 'CC0-1.0', 'n_tokens_mistral': 3835, 'n_tokens_neox': 3556, 'n_words': 1731, 'path': 'Exercise_sessions/Lab Chapter 5.ipynb', 'repo_name': 'Valkenborg/IML', 'size': '13270'}