Model Selection
library(ISLR)
library(leaps)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(readr)Data Import
HWTests <- read_table("/Volumes/BOOTCAMP/r documents/HWTest(1).tsv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## male = col_double(),
## age = col_double(),
## education = col_double(),
## currentSmoker = col_double(),
## cigsPerDay = col_double(),
## BPMeds = col_double(),
## prevalentStroke = col_double(),
## prevalentHyp = col_double(),
## diabetes = col_double(),
## totChol = col_double(),
## sysBP = col_double(),
## diaBP = col_double(),
## BMI = col_double(),
## heartRate = col_double(),
## TenYearCHD = col_double(),
## glucose = col_double()
## )HWTrains <- read_table("/Volumes/BOOTCAMP/r documents/HWTrain(1).tsv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## male = col_double(),
## age = col_double(),
## education = col_double(),
## currentSmoker = col_double(),
## cigsPerDay = col_double(),
## BPMeds = col_double(),
## prevalentStroke = col_double(),
## prevalentHyp = col_double(),
## diabetes = col_double(),
## totChol = col_double(),
## sysBP = col_double(),
## diaBP = col_double(),
## BMI = col_double(),
## heartRate = col_double(),
## TenYearCHD = col_double(),
## glucose = col_double()
## )attach(HWTrains)
trainData<- HWTrains
attach(HWTests)## The following objects are masked from HWTrains:
##
## age, BMI, BPMeds, cigsPerDay, currentSmoker, diabetes, diaBP,
## education, glucose, heartRate, male, prevalentHyp, prevalentStroke,
## sysBP, TenYearCHD, totCholperform model selection using backward elimination on the HWTrains dataset
model1= regsubsets(glucose~.,data = HWTrains,method='backward' )
res=summary(model1)
res$adjr2## [1] 0.4268744 0.4336558 0.4362070 0.4380689 0.4396029 0.4408325 0.4408644
## [8] 0.4407784coef(model1, 8)## (Intercept) male cigsPerDay diabetes sysBP diaBP
## 65.59682727 0.81988369 -0.09364035 95.89338838 0.14271552 -0.18197945
## BMI heartRate TenYearCHD
## 0.07881455 0.10529847 3.97170688res.df1 = data.frame(adj_r2 = res$adjr2, num = 1:8)
res.df1## adj_r2 num
## 1 0.4268744 1
## 2 0.4336558 2
## 3 0.4362070 3
## 4 0.4380689 4
## 5 0.4396029 5
## 6 0.4408325 6
## 7 0.4408644 7
## 8 0.4407784 8performs model selection using the forward selection method on the HWTrains dataset to fit regression models predicting glucose:
model2= regsubsets(glucose~.,data = HWTrains, method='forward')
res2=summary(model2)
res2$adjr2## [1] 0.4268744 0.4336558 0.4362070 0.4380689 0.4396029 0.4408325 0.4408644
## [8] 0.4407784coef(model2, 8)## (Intercept) male cigsPerDay diabetes sysBP diaBP
## 65.59682727 0.81988369 -0.09364035 95.89338838 0.14271552 -0.18197945
## BMI heartRate TenYearCHD
## 0.07881455 0.10529847 3.97170688res.df2 = data.frame(adj_r2 = res2$adjr2, num = 1:8)
res.df2## adj_r2 num
## 1 0.4268744 1
## 2 0.4336558 2
## 3 0.4362070 3
## 4 0.4380689 4
## 5 0.4396029 5
## 6 0.4408325 6
## 7 0.4408644 7
## 8 0.4407784 8performs model selection using the exhaustive search method on the HWTrains dataset to fit regression models predicting glucose:
model3= regsubsets(glucose~.,data = HWTrains)
res3=summary(model3)
res3$adjr2## [1] 0.4268744 0.4336558 0.4362070 0.4380689 0.4396029 0.4408325 0.4408644
## [8] 0.4407784coef(model3, 8)## (Intercept) male cigsPerDay diabetes sysBP diaBP
## 65.59682727 0.81988369 -0.09364035 95.89338838 0.14271552 -0.18197945
## BMI heartRate TenYearCHD
## 0.07881455 0.10529847 3.97170688res.df3 = data.frame(adj_r2 = res3$adjr2, num = 1:8)
res.df3## adj_r2 num
## 1 0.4268744 1
## 2 0.4336558 2
## 3 0.4362070 3
## 4 0.4380689 4
## 5 0.4396029 5
## 6 0.4408325 6
## 7 0.4408644 7
## 8 0.4407784 8#RIDGE
xtrain=model.matrix(glucose~., data=HWTrains[,-1])
ytrain=HWTrains$glucose
xtest=model.matrix(glucose~., data=HWTests[,-1])
ytest=HWTests$glucose
set.seed(1)
library(glmnet)## Loading required package: Matrix## Loaded glmnet 4.1-8ridge.fit=cv.glmnet(xtrain, ytrain, alpha=0)
plot(ridge.fit)
ridge.lambda=ridge.fit$lambda.min
ridge.lambda## [1] 1.666902ridge.pred=predict(ridge.fit, s=ridge.lambda, newx = xtest)
ridge.err=mean((ridge.pred-ytest)^2)
ridge.err## [1] 315.3806mod1 = lm(glucose ~ ., HWTrains[,c(names(which(res2$outmat[8,] == '*')),'glucose')])
summary(mod1)##
## Call:
## lm(formula = glucose ~ ., data = HWTrains[, c(names(which(res2$outmat[8,
## ] == "*")), "glucose")])
##
## Residuals:
## Min 1Q Median 3Q Max
## -118.741 -8.524 -1.284 6.539 215.414
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.59683 3.67138 17.867 < 2e-16 ***
## male 0.81988 0.82660 0.992 0.321351
## cigsPerDay -0.09364 0.03473 -2.696 0.007058 **
## diabetes 95.89339 2.27017 42.241 < 2e-16 ***
## sysBP 0.14272 0.02858 4.993 6.34e-07 ***
## diaBP -0.18198 0.05293 -3.438 0.000595 ***
## BMI 0.07881 0.10118 0.779 0.436061
## heartRate 0.10530 0.03319 3.173 0.001527 **
## TenYearCHD 3.97171 1.07876 3.682 0.000236 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.08 on 2557 degrees of freedom
## Multiple R-squared: 0.4425, Adjusted R-squared: 0.4408
## F-statistic: 253.7 on 8 and 2557 DF, p-value: < 2.2e-16preds = predict(mod1, HWTests)
MSE1 <- mean((HWTests$glucose - preds)^2)
MSE1## [1] 322.6679set.seed(1)
cvfit1 = cv.glmnet(data.matrix(HWTrains),HWTrains$glucose,
type.measure = 'mse',
nfolds=5,
alpha=0)
ridge.lambda=cvfit1$lambda.min
ridge.lambda## [1] 2.550623preds.ridge1 = predict(cvfit1, s=ridge.lambda,data.matrix(HWTests))
MSE2 <- mean((preds.ridge1-HWTests$glucose )^2)
MSE2## [1] 6.445823redge.coeff=predict(cvfit1, type="coefficients", s=ridge.lambda)[1:15,]
redge.coeff## (Intercept) male age education currentSmoker
## 9.1029411848 0.0847451330 0.0075015845 0.0089032239 0.0338940197
## cigsPerDay BPMeds prevalentStroke prevalentHyp diabetes
## -0.0126838305 0.0162192361 0.3462622352 -0.0044079767 12.5476653626
## totChol sysBP diaBP BMI heartRate
## -0.0005694605 0.0151047922 -0.0166236216 0.0137557002 0.0145627497#Lasso
set.seed(1)
lasso.fit=cv.glmnet(xtrain, ytrain, alpha=1)
plot(lasso.fit)
lasso.lambda=lasso.fit$lambda.min
lasso.lambda## [1] 0.1746273lasso.coeff=predict(lasso.fit, type="coefficients", s=lasso.lambda)
lasso.coeff## 16 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 65.80014600
## (Intercept) .
## age 0.02203803
## education .
## currentSmoker .
## cigsPerDay -0.06551281
## BPMeds .
## prevalentStroke 0.64000496
## prevalentHyp .
## diabetes 95.43093594
## totChol .
## sysBP 0.10237663
## diaBP -0.09656369
## BMI 0.02930911
## heartRate 0.08441828
## TenYearCHD 3.61028868library(pls)##
## Attaching package: 'pls'## The following object is masked from 'package:stats':
##
## loadingspcr=pcr(glucose~., data=HWTrains, scale=TRUE, validation="CV")
summary(pcr)## Data: X dimension: 2566 15
## Y dimension: 2566 1
## Fit method: svdpc
## Number of components considered: 15
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 25.52 25.1 25.12 24.69 24.49 24.08 23.97
## adjCV 25.52 25.1 25.11 24.68 24.49 24.15 24.28
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 20.07 19.83 19.79 19.56 19.54 19.55 19.56
## adjCV 19.96 19.77 19.76 19.54 19.52 19.52 19.53
## 14 comps 15 comps
## CV 19.55 19.53
## adjCV 19.53 19.50
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 21.460 34.393 42.425 49.782 56.65 63.47 70.02 75.94
## glucose 3.866 3.867 7.953 9.331 12.66 12.85 41.59 42.22
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 81.69 86.63 91.09 94.97 97.47 98.84 100.00
## glucose 42.25 43.83 43.94 43.99 43.99 44.05 44.28