STA-6543: Assignment 6
Austin Thrash - FFK221
2024-03-21
Table of Contents:
ISLR Chapter 7
Question 6
In this exercise, you will further analyze the Wage data set considered throughout this chapter.
head(Wage)## year age maritl race education region
## 231655 2006 18 1. Never Married 1. White 1. < HS Grad 2. Middle Atlantic
## 86582 2004 24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic
## 161300 2003 45 2. Married 1. White 3. Some College 2. Middle Atlantic
## 155159 2003 43 2. Married 3. Asian 4. College Grad 2. Middle Atlantic
## 11443 2005 50 4. Divorced 1. White 2. HS Grad 2. Middle Atlantic
## 376662 2008 54 2. Married 1. White 4. College Grad 2. Middle Atlantic
## jobclass health health_ins logwage wage
## 231655 1. Industrial 1. <=Good 2. No 4.318063 75.04315
## 86582 2. Information 2. >=Very Good 2. No 4.255273 70.47602
## 161300 1. Industrial 1. <=Good 1. Yes 4.875061 130.98218
## 155159 2. Information 2. >=Very Good 1. Yes 5.041393 154.68529
## 11443 2. Information 1. <=Good 1. Yes 4.318063 75.04315
## 376662 2. Information 2. >=Very Good 1. Yes 4.845098 127.115746.(a) Age to predict Wage: Polynomial Regression
Perform polynomial regression to predict wage using age. Use cross-validation to select the optimal degree d for the polynomial. What degree was chosen, and how does this compare to the results of hypothesis testing using ANOVA? Make a plot of the resulting polynomial fit to the data.
# Hypothesis Testing Using ANOVA
ploymodel1 <- lm(wage ~ poly(age, 1), data=Wage)
ploymodel2 <- lm(wage ~ poly(age, 2), data=Wage)
ploymodel3 <- lm(wage ~ poly(age, 3), data=Wage)
ploymodel4 <- lm(wage ~ poly(age, 4), data=Wage)
ploymodel5 <- lm(wage ~ poly(age, 5), data=Wage)
anova_results <- anova(ploymodel1, ploymodel2, ploymodel3, ploymodel4, ploymodel5)
print(anova_results)## 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)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.5931 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8888 0.001679 **
## 4 2995 4771604 1 6070 3.8098 0.051046 .
## 5 2994 4770322 1 1283 0.8050 0.369682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Polynomial Regression Using Cross-Validation
set.seed(4)
cv.error.10 <- rep(0, 10)
for (i in 1:10) {
glm.fit <- glm(wage ~ poly(age, i), data = Wage)
cv.error.10[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
# Getting model number with lowest MSE
optimal.degree <- which.min(cv.error.10)
# Create Formula
formula_str <- sprintf("y ~ poly(x, %d, raw = TRUE)", optimal.degree)
# Convert to formula object
formula_obj <- as.formula(formula_str)
# Plotting
ggplot(Wage, aes(x = age, y = wage)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", formula = formula_obj, color = "blue") +
labs(title = paste("Degree-", optimal.degree, "Polynomial Fit"))
When performing Polynomial Regression using cross-validation (k-fold validation), the results suggest using a degree of 4 and it seems to be a pretty good fit. However when using hypothesis testing using ANOVA, the optimal degree is 3, as going any higher than that does not grant significant improvement to the model and will add further complexity.
6.(b) Age to predict Wage: Step Function
Fit a step function to predict wage using age, and perform cross-validation to choose the optimal number of cuts. Make a plot of the fit obtained.
set.seed(4)
cuts <- numeric(19)
for (i in 2:20) {
Wage$age_cut <- cut(Wage$age, i, include.lowest = TRUE)
lm.fit <- glm(wage ~ age_cut, data = Wage)
cuts[i - 1] <- cv.glm(Wage, lm.fit, K = 10)$delta[1]
}
bins <- 2:20
plot(bins, cuts, type = "b", pch = 19, col = "blue", xlab = "Number of Bins", ylab = "CV Error", main = "CV Error vs. Number of Bins for Age Binning")
# Getting model number with lowest MSE
optimal.cut <- which.min(cuts) + 1
# Create Formula
formula_str <- sprintf("y ~ cut(x, %d, raw = TRUE)", optimal.cut)
# Convert to formula object
formula_obj <- as.formula(formula_str)
# Plotting
ggplot(Wage, aes(x = age, y = wage)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", formula = formula_obj, color = "blue") +
labs(title = paste("Degree -", optimal.cut, ": Step Function"))
When using the step function to predict Wage using Age, the results suggest using a step function with 16 intervals. However I would most likely use an interval of 7 or 8 as the “1 Standard Deviation” rule applies and they would make the model less complex.
Question 10
This question relates to the College data set.
head(College)## Private Apps Accept Enroll Top10perc Top25perc
## Abilene Christian University Yes 1660 1232 721 23 52
## Adelphi University Yes 2186 1924 512 16 29
## Adrian College Yes 1428 1097 336 22 50
## Agnes Scott College Yes 417 349 137 60 89
## Alaska Pacific University Yes 193 146 55 16 44
## Albertson College Yes 587 479 158 38 62
## F.Undergrad P.Undergrad Outstate Room.Board Books
## Abilene Christian University 2885 537 7440 3300 450
## Adelphi University 2683 1227 12280 6450 750
## Adrian College 1036 99 11250 3750 400
## Agnes Scott College 510 63 12960 5450 450
## Alaska Pacific University 249 869 7560 4120 800
## Albertson College 678 41 13500 3335 500
## Personal PhD Terminal S.F.Ratio perc.alumni Expend
## Abilene Christian University 2200 70 78 18.1 12 7041
## Adelphi University 1500 29 30 12.2 16 10527
## Adrian College 1165 53 66 12.9 30 8735
## Agnes Scott College 875 92 97 7.7 37 19016
## Alaska Pacific University 1500 76 72 11.9 2 10922
## Albertson College 675 67 73 9.4 11 9727
## Grad.Rate
## Abilene Christian University 60
## Adelphi University 56
## Adrian College 54
## Agnes Scott College 59
## Alaska Pacific University 15
## Albertson College 5510.(a)
Split the data into a training set and a test set. Using out-of-state tuition as the response and the other variables as the predictors, perform forward stepwise selection on the training set in order to identify a satisfactory model that uses just a subset of the predictors.
# Performing 70/30 train test split
set.seed(4)
train_index <- sample(1:nrow(College), round(nrow(College) * 0.7))
train <- College[train_index, ]
test <- College[-train_index, ]
nv_max <- ncol(College) - 2
fwd <- regsubsets(Outstate ~ ., data = train, method = "forward", nvmax = nv_max)
summary1 <- summary(fwd)
summary1## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = train, method = "forward",
## nvmax = nv_max)
## 17 Variables (and intercept)
## Forced in Forced out
## PrivateYes FALSE FALSE
## Apps FALSE FALSE
## Accept FALSE FALSE
## Enroll FALSE FALSE
## Top10perc FALSE FALSE
## Top25perc FALSE FALSE
## F.Undergrad FALSE FALSE
## P.Undergrad FALSE FALSE
## Room.Board FALSE FALSE
## Books FALSE FALSE
## Personal FALSE FALSE
## PhD FALSE FALSE
## Terminal FALSE FALSE
## S.F.Ratio FALSE FALSE
## perc.alumni FALSE FALSE
## Expend FALSE FALSE
## Grad.Rate FALSE FALSE
## 1 subsets of each size up to 16
## Selection Algorithm: forward
## PrivateYes Apps Accept Enroll Top10perc Top25perc F.Undergrad
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " " " " " " "
## 4 ( 1 ) "*" " " " " " " " " " " " "
## 5 ( 1 ) "*" " " " " " " " " " " " "
## 6 ( 1 ) "*" " " " " " " " " " " " "
## 7 ( 1 ) "*" " " " " " " " " " " " "
## 8 ( 1 ) "*" " " "*" " " " " " " " "
## 9 ( 1 ) "*" " " "*" "*" " " " " " "
## 10 ( 1 ) "*" "*" "*" "*" " " " " " "
## 11 ( 1 ) "*" "*" "*" "*" "*" " " " "
## 12 ( 1 ) "*" "*" "*" "*" "*" " " " "
## 13 ( 1 ) "*" "*" "*" "*" "*" " " " "
## 14 ( 1 ) "*" "*" "*" "*" "*" " " " "
## 15 ( 1 ) "*" "*" "*" "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## P.Undergrad Room.Board Books Personal PhD Terminal S.F.Ratio
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " "
## 5 ( 1 ) " " "*" " " " " " " " " " "
## 6 ( 1 ) " " "*" " " " " " " "*" " "
## 7 ( 1 ) " " "*" " " "*" " " "*" " "
## 8 ( 1 ) " " "*" " " "*" " " "*" " "
## 9 ( 1 ) " " "*" " " "*" " " "*" " "
## 10 ( 1 ) " " "*" " " "*" " " "*" " "
## 11 ( 1 ) " " "*" " " "*" " " "*" " "
## 12 ( 1 ) " " "*" " " "*" " " "*" "*"
## 13 ( 1 ) " " "*" " " "*" "*" "*" "*"
## 14 ( 1 ) " " "*" "*" "*" "*" "*" "*"
## 15 ( 1 ) " " "*" "*" "*" "*" "*" "*"
## 16 ( 1 ) " " "*" "*" "*" "*" "*" "*"
## perc.alumni Expend Grad.Rate
## 1 ( 1 ) " " "*" " "
## 2 ( 1 ) " " "*" " "
## 3 ( 1 ) " " "*" " "
## 4 ( 1 ) "*" "*" " "
## 5 ( 1 ) "*" "*" "*"
## 6 ( 1 ) "*" "*" "*"
## 7 ( 1 ) "*" "*" "*"
## 8 ( 1 ) "*" "*" "*"
## 9 ( 1 ) "*" "*" "*"
## 10 ( 1 ) "*" "*" "*"
## 11 ( 1 ) "*" "*" "*"
## 12 ( 1 ) "*" "*" "*"
## 13 ( 1 ) "*" "*" "*"
## 14 ( 1 ) "*" "*" "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"ggplot() +
geom_path(aes(x = seq_along(summary1$bic), y = summary1$bic)) +
geom_point(aes(x = which.min(summary1$bic), y = min(summary1$bic)), size = 3, colour = "red") +
xlab("Number of predictors") +
ylab("BIC") +
theme_minimal()
When performing forward stepwise selection and using the out-of-state tutition as the response variable, the results suggest using 9 predictors which are:
Private, Accept, Enroll, Room.Board, Personal, Terminal, perc.alumni, Expend, and Grad.Rate
However if I were picking the number of predictors I would most likely chose to use 5 as it has a BIC value closer the BIC value when using 9 predictors. This would also reduce the complexity of the model and might provide more accurate results.
10.(b)
Fit a GAM on the training data, using out-of-state tuition as the response and the features selected in the previous step as the predictors. Plot the results, and explain your findings.
fits <- gam(Outstate ~ Private + s(Accept) + s(Enroll) + s(Room.Board) + s(Personal) + s(Terminal) + s(perc.alumni) + s(Expend) + s(Grad.Rate), data = train)
par(mfrow = c(3,4))
plot(fits, se = TRUE)
From the results, we can see that ‘Private’, ‘Accept’, ‘Room.Board’, ‘perc.alumni’, ‘Expend’, and ‘Grad.Rate’ all seem to exhibit a positive relationship with the response variable; For each increase in these variables we can expect to see an increase in the response on average. While ‘Enroll’ and ‘Personal’ have a negative relationship with the response; On average we expect to see a decrease in the response for each increase of these variables.
10.(c)
Evaluate the model obtained on the test set, and explain the results obtained.
preds_nonlinear <- predict(fits, newdata = test)
mse_nonlinear <- mean((preds_nonlinear - test$Outstate)^2)
fits_lm <- lm(Outstate ~ Private + Accept + Enroll + Room.Board + Personal + Terminal + perc.alumni + Expend + Grad.Rate, data = train)
preds_linear <- predict(fits_lm, newdata = test)
mse_linear <- mean((preds_linear - test$Outstate)^2)
paste("Linear RMSE:", sqrt(mse_linear))## [1] "Linear RMSE: 2015.98979665636"paste("Non-linear RMSE:", sqrt(mse_nonlinear))## [1] "Non-linear RMSE: 1801.15924523687"We can evaluate the model by looking at the “R Mean Squared Error” (RMSE) and comparing it to a normal linear model RMSE using the same predictors. From the results we can see that the non-linear model out performed the linear model by a difference of about 214. However if using less predictors, one could possibly obtain better results.
10.(d)
For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(fits)##
## Call: gam(formula = Outstate ~ Private + s(Accept) + s(Enroll) + s(Room.Board) +
## s(Personal) + s(Terminal) + s(perc.alumni) + s(Expend) +
## s(Grad.Rate), data = train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6562.89 -1105.46 55.91 1204.22 7257.87
##
## (Dispersion Parameter for gaussian family taken to be 3227875)
##
## Null Deviance: 8687534607 on 543 degrees of freedom
## Residual Deviance: 1646217794 on 510.0004 degrees of freedom
## AIC: 9731.806
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 2356834941 2356834941 730.151 < 2.2e-16 ***
## s(Accept) 1 612448391 612448391 189.737 < 2.2e-16 ***
## s(Enroll) 1 133668721 133668721 41.411 2.847e-10 ***
## s(Room.Board) 1 1042887706 1042887706 323.088 < 2.2e-16 ***
## s(Personal) 1 53770186 53770186 16.658 5.194e-05 ***
## s(Terminal) 1 336364664 336364664 104.206 < 2.2e-16 ***
## s(perc.alumni) 1 376312664 376312664 116.582 < 2.2e-16 ***
## s(Expend) 1 830223945 830223945 257.204 < 2.2e-16 ***
## s(Grad.Rate) 1 66909075 66909075 20.729 6.625e-06 ***
## Residuals 510 1646217794 3227875
## ---
## 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(Accept) 3 14.369 5.313e-09 ***
## s(Enroll) 3 5.854 0.0006175 ***
## s(Room.Board) 3 1.282 0.2797653
## s(Personal) 3 3.504 0.0153515 *
## s(Terminal) 3 1.961 0.1188499
## s(perc.alumni) 3 0.917 0.4323377
## s(Expend) 3 32.193 < 2.2e-16 ***
## s(Grad.Rate) 3 1.997 0.1134792
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To answer this question we can look at the “Anova for Nonparametric Effects” section. The only variables that exhibit a non-linear relationship with the response is ‘Accept’, ‘Enroll’, ‘Personal’, and ‘Expend’ with p-values of 5.313e-09, 0.0006175, 0.0153515, and <2.2e-16 respectively.