HW6
Elaine Perera
5/11/2021
Chapter 07 (page 297): 6, 10 Question #6
library(ISLR)
library(MASS)
library(ggplot2)
library(data.table)
library(leaps)
library(glmnet)## Loading required package: Matrix## Loaded glmnet 4.1-1library(boot)
set.seed(5)
data(Wage)
#a
deltas = rep(NA, 10)
for (i in 1:10) {
glm.model = glm(wage~poly(age, i), data=Wage)
deltas[i] = cv.glm(Wage, glm.model, K=10)$delta[2]
}
data <- data.table(seq(1:10),deltas,keep.rownames = TRUE)
sd(deltas)## [1] 25.6542min(deltas)## [1] 1592.924deltas## [1] 1675.786 1600.476 1595.701 1593.974 1594.573 1593.530 1592.924 1596.653
## [9] 1593.402 1593.527ggplot(data, aes(V1,deltas))+geom_line()
fit.1 = lm(wage~poly(age, 1), data=Wage)
fit.2 = lm(wage~poly(age, 2), data=Wage)
fit.3 = lm(wage~poly(age, 3), data=Wage)
fit.4 = lm(wage~poly(age, 4), data=Wage)
fit.5 = lm(wage~poly(age, 5), data=Wage)
fit.6 = lm(wage~poly(age, 6), data=Wage)
fit.7 = lm(wage~poly(age, 7), data=Wage)
fit.8 = lm(wage~poly(age, 8), data=Wage)
fit.9 = lm(wage~poly(age, 9), data=Wage)
fit.10 = lm(wage~poly(age, 10), data=Wage)
anova(fit.1, fit.2, fit.3, fit.4, fit.5, fit.6, fit.7, fit.8, fit.9, fit.10)## Analysis of Variance Table
##
## Model 1: wage ~ poly(age, 1)
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Model 6: wage ~ poly(age, 6)
## Model 7: wage ~ poly(age, 7)
## Model 8: wage ~ poly(age, 8)
## Model 9: wage ~ poly(age, 9)
## Model 10: wage ~ poly(age, 10)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.7638 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.9005 0.001669 **
## 4 2995 4771604 1 6070 3.8143 0.050909 .
## 5 2994 4770322 1 1283 0.8059 0.369398
## 6 2993 4766389 1 3932 2.4709 0.116074
## 7 2992 4763834 1 2555 1.6057 0.205199
## 8 2991 4763707 1 127 0.0796 0.777865
## 9 2990 4756703 1 7004 4.4014 0.035994 *
## 10 2989 4756701 1 3 0.0017 0.967529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggplot(Wage, aes(age,wage)) + geom_point(color="green") +
stat_smooth(method = "lm", formula = y ~ poly(x, 3), size = 1)
#b
cvs = rep(NA, 10)
for (i in 2:10) {
Wage$cut.point = cut(Wage$age, i)
lm.fit = glm(wage~cut.point, data=Wage)
cvs[i] = cv.glm(Wage, lm.fit, K=10)$delta[2]
}
cvs.data <- data.table(c(2,3,4,5,6,7,8,9,10),cvs[-1])
ggplot(cvs.data, aes(V1,V2)) + geom_line()
fit <- glm(wage~cut(age, 8), data=Wage)
prediction <- predict(fit, data.frame(age = seq(from=18,to=80)))
pre <- data.table(seq(from=18,to=80),prediction,keep.rownames = TRUE)
ggplot(Wage, aes(age,wage)) + geom_point(color="green") + geom_line(data = pre, aes(V1,prediction), size= 1)
Question #10
#a
set.seed(1)
dt <- data.table(College)
n <- length(dt$Outstate)
train <- sample(n, n/2)
dt.train <- College[train, ]
dt.test <- College[-train, ]
reg.fit <- regsubsets(Outstate ~ ., data = dt.train, nvmax = 17, method = "forward")
reg.summary <- summary(reg.fit)
ggplot(data.frame(cp =reg.summary$cp, nrVar=1:17), aes(x=nrVar, y=cp))+xlab("Number of Variables") + ylab("Cp") + geom_line()
which.min(reg.summary$cp)## [1] 14ggplot(data.frame(bic =reg.summary$bic, nrVar=1:17), aes(x=nrVar, y=bic))+xlab("Number of Variables") + ylab("BIC") + geom_line()
which.min(reg.summary$bic)## [1] 6ggplot(data.frame(adjr2 =reg.summary$adjr2, nrVar=1:17), aes(x=nrVar, y=adjr2))+xlab("Number of Variables") + ylab("adjr2") + geom_line()
which.max(reg.summary$adjr2)## [1] 14co <- coef(reg.fit, id = 6)
names(co)## [1] "(Intercept)" "PrivateYes" "Room.Board" "Terminal" "perc.alumni"
## [6] "Expend" "Grad.Rate"#b
library(gam)## Loading required package: splines## Loading required package: foreach## Loaded gam 1.20gam.fit <- gam(Outstate ~ Private + s(Room.Board, df = 2) + s(PhD, df = 2) +
s(perc.alumni, df = 2) + s(Expend, df = 2) + s(Grad.Rate, df = 2), data = dt.train)
par(mfrow = c(2, 3))
plot(gam.fit, se = T, col = "red")
par(mfrow = c(1, 1))
#c
gam.pred <- predict(gam.fit, dt.test)
gam.err <- mean((dt.test$Outstate - gam.pred)^2)
gam.err## [1] 3438192lm.pred <- predict(lm(Outstate~Private+Room.Board+PhD+perc.alumni+Expend+Grad.Rate, data = dt.train), dt.test)
lm.err <- mean((dt.test$Outstate - lm.pred)^2)
lm.err## [1] 3841483#d
summary(gam.fit)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(PhD,
## df = 2) + s(perc.alumni, df = 2) + s(Expend, df = 2) + s(Grad.Rate,
## df = 2), data = dt.train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6500.0 -1344.2 -104.2 1350.9 8792.9
##
## (Dispersion Parameter for gaussian family taken to be 3946537)
##
## Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1483898477 on 376.0001 degrees of freedom
## AIC: 7007.986
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1843025465 1843025465 466.998 < 2.2e-16 ***
## s(Room.Board, df = 2) 1 1692725515 1692725515 428.914 < 2.2e-16 ***
## s(PhD, df = 2) 1 397159048 397159048 100.635 < 2.2e-16 ***
## s(perc.alumni, df = 2) 1 351049729 351049729 88.951 < 2.2e-16 ***
## s(Expend, df = 2) 1 419096845 419096845 106.194 < 2.2e-16 ***
## s(Grad.Rate, df = 2) 1 81869675 81869675 20.745 7.101e-06 ***
## Residuals 376 1483898477 3946537
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## Private
## s(Room.Board, df = 2) 1 2.005 0.15763
## s(PhD, df = 2) 1 3.781 0.05258 .
## s(perc.alumni, df = 2) 1 0.314 0.57572
## s(Expend, df = 2) 1 47.156 2.725e-11 ***
## s(Grad.Rate, df = 2) 1 1.057 0.30446
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1