library(ggplot2) # graphics library
library(ISLR) # contains code and data from the textbook
library(knitr) # contains kable() function
library(leaps) # for regsubsets() function
library(boot) # for cv.glm
library(gam)## Loading required package: splines## Loading required package: foreach## Loaded gam 1.14library(glmnet)## Loading required package: Matrix## Loaded glmnet 2.0-5options(scipen = 4) # Suppresses scientific notationThis portion of the lab gets you to carry out the Lab in §6.5.1 of ISLR (Pages 244 - 247). You will want to have the textbook Lab open in front you as you go through these exercises. The ISLR Lab provides much more context and explanation for what you’re doing.
You will need the
Hittersdata set from theISLRlibrary in order to complete this exercise.
Please run all of the code indicated in §6.5.1 of ISLR, even if I don’t explicitly ask you to do so in this document.
View() command on the Hitters data to see what the data set looks like.#View(Hitters)qplot to construct a histogram of of the Salary variable. Does Salary appear to be normally distributed, or is the distribution skewed? What units is Salary recorded in?# Edit mepanel.cor function that properly handles missing values. Use the pairs command to construct a pairs plot for the Hitters data, displaying correlations in the lower panel and plots in the upper panel. Your pairs plot should include the variables: Salary, AtBat, Hits, HmRun, CRBI, RBI, Errors. Read the ?Hitters documentation to understand what these variables mean.panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y, use = "complete.obs"))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
# Edit meNA values using the na.omit() command. Then, use the regsubsets command to perform best subset selection. Your should go up to models of size 15.# Edit me# Edit me# Edit meeval = TRUE in the header to get this code to run once you’ve constructed the regfit.full object in previous parts of this problem.)plot(regfit.full,scale="r2")
plot(regfit.full,scale="adjr2")
plot(regfit.full,scale="Cp")
plot(regfit.full,scale="bic")# Edit meThe next portion of the lab gets you to carry out the Lab in §6.5.2 of ISLR (Page 247). You will want to have the textbook Lab open in front you as you go through these exercises. The ISLR Lab provides much more context and explanation for what you’re doing.
Hitters data.# Edit me# Edit me# Edit meThe next portion of the lab gets you to carry out the Lab in §6.5.3 of ISLR (Page 248 - 250). You will want to have the textbook Lab open in front you as you go through these exercises. The ISLR Lab provides much more context and explanation for what you’re doing.
When we use the Validation Set Approach or Cross-Validation, do we first run best subset/forward/backward selection to get a sequence of models, or do we perform these steps on the training set of each validation method?
True or False: (In the validation set approach…) In the end, the variables we wind up using in our final model (the one fit to the full data) will be the same as those selected on the training set.
True or False: When we use Cross-Validation with Best subset/Forward/Backward selection, we need to apply the variable selection method on each set of training data.
# You'll need to set eval = TRUE in the code chunk header
# in order for this code to run
set.seed(1)
train=sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE)
test=(!train)
regfit.best=regsubsets(Salary~.,data=Hitters[train,],nvmax=19)
test.mat=model.matrix(Salary~.,data=Hitters[test,])
val.errors=rep(NA,19)
for(i in 1:19){
coefi=coef(regfit.best,id=i)
pred=test.mat[,names(coefi)]%*%coefi
val.errors[i]=mean((Hitters$Salary[test]-pred)^2)
}
val.errors
which.min(val.errors)
coef(regfit.best,10)
predict.regsubsets=function(object,newdata,id,...){
form=as.formula(object$call[[2]])
mat=model.matrix(form,newdata)
coefi=coef(object,id=id)
xvars=names(coefi)
mat[,xvars]%*%coefi
}
regfit.best=regsubsets(Salary~.,data=Hitters,nvmax=19)
coef(regfit.best,10)
k=10
set.seed(1)
folds=sample(1:k,nrow(Hitters),replace=TRUE)
cv.errors=matrix(NA,k,19, dimnames=list(NULL, paste(1:19)))
for(j in 1:k){
best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,],nvmax=19)
for(i in 1:19){
pred=predict(best.fit,Hitters[folds==j,],id=i)
cv.errors[j,i]=mean( (Hitters$Salary[folds==j]-pred)^2)
}
}
mean.cv.errors=apply(cv.errors,2,mean)
mean.cv.errors
par(mfrow=c(1,1))
plot(mean.cv.errors,type='b')
reg.best=regsubsets(Salary~.,data=Hitters, nvmax=19)
coef(reg.best,11)The next portion of the lab gets you to carry out the Lab in §6.6.2 of ISLR (Page 255). You will want to have the textbook Lab open in front you as you go through these exercises. The ISLR Lab provides much more context and explanation for what you’re doing.
# Define x matrix and y vector for use with glmnet
x <- model.matrix(Salary~.,Hitters)[,-1]
y <- Hitters$Salary
# Split data into test and train
set.seed(1)
train <- sample(1:nrow(x), nrow(x)/2)
test <- (-train)
y.test <- y[test]
# Predefined grid of lambda values:
grid=10^seq(10,-2, length =100)glmnet command to fit a lasso model to the train subset of the Hitters data. You’ll want to specify lambda=grid to use the predefined sequence of \(lambda\) values constructed above. Use the plot command to produce a regularization plot of your model fit.# Edit mecv.glmnet. Use the plot command to construct a CV error curve.# Edit me# Edit me# Edit me# Edit me