Chapter VII - Moving Beyond Linearity

All the questions are as per the ISL seventh printing.

Common

libsUsed<-c("dplyr","ggplot2","tidyverse",
            "ISLR","caret","MASS", "gridExtra",
            "pls","latex2exp","data.table")
invisible(lapply(libsUsed, library, character.only = TRUE))

Question 7.6 - Page 299

In this exercise, you will further analyze the Wage data set considered throughout this chapter.

(a) 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.

(b) 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. In this exercise, we will generate simulated data, and will then use this data to perform best subset selection.

Answer

Lets get the data.

set.seed(1984)
wageDat<-ISLR::Wage
wageDat %>% str %>% print

## 'data.frame':    3000 obs. of  11 variables:
##  $ year      : int  2006 2004 2003 2003 2005 2008 2009 2008 2006 2004 ...
##  $ age       : int  18 24 45 43 50 54 44 30 41 52 ...
##  $ maritl    : Factor w/ 5 levels "1. Never Married",..: 1 1 2 2 4 2 2 1 1 2 ...
##  $ race      : Factor w/ 4 levels "1. White","2. Black",..: 1 1 1 3 1 1 4 3 2 1 ...
##  $ education : Factor w/ 5 levels "1. < HS Grad",..: 1 4 3 4 2 4 3 3 3 2 ...
##  $ region    : Factor w/ 9 levels "1. New England",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ jobclass  : Factor w/ 2 levels "1. Industrial",..: 1 2 1 2 2 2 1 2 2 2 ...
##  $ health    : Factor w/ 2 levels "1. <=Good","2. >=Very Good": 1 2 1 2 1 2 2 1 2 2 ...
##  $ health_ins: Factor w/ 2 levels "1. Yes","2. No": 2 2 1 1 1 1 1 1 1 1 ...
##  $ logwage   : num  4.32 4.26 4.88 5.04 4.32 ...
##  $ wage      : num  75 70.5 131 154.7 75 ...
## NULL

wageDat %>% summary %>% print

##       year           age                     maritl           race
##  Min.   :2003   Min.   :18.00   1. Never Married: 648   1. White:2480
##  1st Qu.:2004   1st Qu.:33.75   2. Married      :2074   2. Black: 293
##  Median :2006   Median :42.00   3. Widowed      :  19   3. Asian: 190
##  Mean   :2006   Mean   :42.41   4. Divorced     : 204   4. Other:  37
##  3rd Qu.:2008   3rd Qu.:51.00   5. Separated    :  55
##  Max.   :2009   Max.   :80.00
##
##               education                     region               jobclass
##  1. < HS Grad      :268   2. Middle Atlantic   :3000   1. Industrial :1544
##  2. HS Grad        :971   1. New England       :   0   2. Information:1456
##  3. Some College   :650   3. East North Central:   0
##  4. College Grad   :685   4. West North Central:   0
##  5. Advanced Degree:426   5. South Atlantic    :   0
##                           6. East South Central:   0
##                           (Other)              :   0
##             health      health_ins      logwage           wage
##  1. <=Good     : 858   1. Yes:2083   Min.   :3.000   Min.   : 20.09
##  2. >=Very Good:2142   2. No : 917   1st Qu.:4.447   1st Qu.: 85.38
##                                      Median :4.653   Median :104.92
##                                      Mean   :4.654   Mean   :111.70
##                                      3rd Qu.:4.857   3rd Qu.:128.68
##                                      Max.   :5.763   Max.   :318.34
##

wageDat %>% sapply(unique) %>% sapply(length) %>% print

##       year        age     maritl       race  education     region   jobclass
##          7         61          5          4          5          1          2
##     health health_ins    logwage       wage
##          2          2        508        508

library(boot)

##
## Attaching package: 'boot'

## The following object is masked from 'package:lattice':
##
##     melanoma

a) Polynomial regression

all.deltas = rep(NA, 10)
for (i in 1:10) {
  glm.fit = glm(wage~poly(age, i), data=Wage)
  all.deltas[i] = cv.glm(Wage, glm.fit, K=10)$delta[2]
}
plot(1:10, all.deltas, xlab="Degree", ylab="CV error", type="l", pch=20, lwd=2, ylim=c(1590, 1700))
min.point = min(all.deltas)
sd.points = sd(all.deltas)
abline(h=min.point + 0.2 * sd.points, col="red", lty="dashed")
abline(h=min.point - 0.2 * sd.points, col="red", lty="dashed")
legend("topright", "0.2-standard deviation lines", lty="dashed", col="red")

# ANOVA
fits=list()
for (i in 1:10) {
  fits[[i]]=glm(wage~poly(age,i),data=wageDat)
}
anova(fits[[1]],fits[[2]],fits[[3]],fits[[4]],fits[[5]],
  fits[[6]],fits[[7]],fits[[8]],fits[[9]],fits[[10]])

## Analysis of Deviance Table
##
## Model  1: wage ~ poly(age, i)
## Model  2: wage ~ poly(age, i)
## Model  3: wage ~ poly(age, i)
## Model  4: wage ~ poly(age, i)
## Model  5: wage ~ poly(age, i)
## Model  6: wage ~ poly(age, i)
## Model  7: wage ~ poly(age, i)
## Model  8: wage ~ poly(age, i)
## Model  9: wage ~ poly(age, i)
## Model 10: wage ~ poly(age, i)
##    Resid. Df Resid. Dev Df Deviance
## 1       2998    5022216
## 2       2997    4793430  1   228786
## 3       2996    4777674  1    15756
## 4       2995    4771604  1     6070
## 5       2994    4770322  1     1283
## 6       2993    4766389  1     3932
## 7       2992    4763834  1     2555
## 8       2991    4763707  1      127
## 9       2990    4756703  1     7004
## 10      2989    4756701  1        3
  • The 4th degree looks the best at the moment
# 3rd or 4th degrees look best based on ANOVA test
# let's go with 4th degree fit
plot(wage~age, data=wageDat, col="darkgrey")
agelims = range(wageDat$age)
age.grid = seq(from=agelims[1], to=agelims[2])
lm.fit = lm(wage~poly(age, 4), data=wageDat)
lm.pred = predict(lm.fit, data.frame(age=age.grid))
lines(age.grid, lm.pred, col="blue", lwd=2)

b) Step function and cross-validation

# cross-validation
cv.error <- rep(0,9)
for (i in 2:10) {
  wageDat$age.cut <- cut(wageDat$age,i)
  glm.fit <- glm(wage~age.cut, data=wageDat)
  cv.error[i-1] <- cv.glm(wageDat, glm.fit, K=10)$delta[1]  # [1]:std, [2]:bias-corrected
}
cv.error

## [1] 1732.337 1682.978 1636.736 1635.600 1624.174 1610.688 1604.081 1612.005
## [9] 1607.022

cv.error

## [1] 1732.337 1682.978 1636.736 1635.600 1624.174 1610.688 1604.081 1612.005
## [9] 1607.022

plot(2:10, cv.error, type="b")

cut.fit <- glm(wage~cut(age,8), data=wageDat)
preds <- predict(cut.fit, newdata=list(age=age.grid), se=TRUE)
se.bands <- preds$fit + cbind(2*preds$se.fit, -2*preds$se.fit)
plot(wageDat$age, wageDat$wage, xlim=agelims, cex=0.5, col="darkgrey")
title("Fit with 8 Age Bands")
lines(age.grid, preds$fit, lwd=2, col="blue")
matlines(age.grid, se.bands, lwd=1, col="blue", lty=3)

Question 7.8 - Page 299

Fit some of the non-linear models investigated in this chapter to the Auto data set. Is there evidence for non-linear relationships in this data set? Create some informative plots to justify your answer.

Answer

autoDat<-ISLR::Auto

autoDat %>% pivot_longer(-c(mpg,name),names_to="Params",values_to="Value") %>% ggplot(aes(x=mpg,y=Value)) +
  geom_point() +
  facet_wrap(~ Params, scales = "free_y")

Very clearly there is a lot of non-linearity in the mpg data, especially for acceleration, weight, displacement, horsepower.

rss = rep(NA, 10)
fits = list()
for (d in 1:10) {
    fits[[d]] = lm(mpg ~ poly(displacement, d), data = autoDat)
    rss[d] = deviance(fits[[d]])
}
rss %>% print

##  [1] 8378.822 7412.263 7392.322 7391.722 7380.838 7270.746 7089.716 6917.401
##  [9] 6737.801 6610.190

anova(fits[[1]],fits[[2]],fits[[3]],fits[[4]],fits[[5]],
  fits[[6]],fits[[7]],fits[[8]],fits[[9]],fits[[10]])

## Analysis of Variance Table
##
## Model  1: mpg ~ poly(displacement, d)
## Model  2: mpg ~ poly(displacement, d)
## Model  3: mpg ~ poly(displacement, d)
## Model  4: mpg ~ poly(displacement, d)
## Model  5: mpg ~ poly(displacement, d)
## Model  6: mpg ~ poly(displacement, d)
## Model  7: mpg ~ poly(displacement, d)
## Model  8: mpg ~ poly(displacement, d)
## Model  9: mpg ~ poly(displacement, d)
## Model 10: mpg ~ poly(displacement, d)
##    Res.Df    RSS Df Sum of Sq       F    Pr(>F)
## 1     390 8378.8
## 2     389 7412.3  1    966.56 55.7108 5.756e-13 ***
## 3     388 7392.3  1     19.94  1.1494  0.284364
## 4     387 7391.7  1      0.60  0.0346  0.852549
## 5     386 7380.8  1     10.88  0.6273  0.428823
## 6     385 7270.7  1    110.09  6.3455  0.012177 *
## 7     384 7089.7  1    181.03 10.4343  0.001344 **
## 8     383 6917.4  1    172.31  9.9319  0.001753 **
## 9     382 6737.8  1    179.60 10.3518  0.001404 **
## 10    381 6610.2  1    127.61  7.3553  0.006990 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Confirming our visual indications, we see that the second degree models work well.

library(glmnet)

## Loading required package: Matrix

##
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
##
##     expand, pack, unpack

## Loaded glmnet 3.0-2

library(boot)

cv.errs = rep(NA, 15)
for (d in 1:15) {
    fit = glm(mpg ~ poly(displacement, d), data = Auto)
    cv.errs[d] = cv.glm(Auto, fit, K = 15)$delta[2]
}
which.min(cv.errs)

## [1] 10

Strangely, we seem to have ended up with a ten variable model here.

# Step functions
cv.errs = rep(NA, 10)
for (c in 2:10) {
    Auto$dis.cut = cut(Auto$displacement, c)
    fit = glm(mpg ~ dis.cut, data = Auto)
    cv.errs[c] = cv.glm(Auto, fit, K = 10)$delta[2]
}
which.min(cv.errs) %>% print

## [1] 9

library(splines)
cv.errs = rep(NA, 10)
for (df in 3:10) {
    fit = glm(mpg ~ ns(displacement, df = df), data = Auto)
    cv.errs[df] = cv.glm(Auto, fit, K = 10)$delta[2]
}
which.min(cv.errs) %>% print

## [1] 10

library(gam)

## Loading required package: foreach

##
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
##
##     accumulate, when

## Loaded gam 1.16.1

# GAMs
fit = gam(mpg ~ s(displacement, 4) + s(horsepower, 4), data = Auto)

## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored

summary(fit)

##
## Call: gam(formula = mpg ~ s(displacement, 4) + s(horsepower, 4), data = Auto)
## Deviance Residuals:
##      Min       1Q   Median       3Q      Max
## -11.2982  -2.1592  -0.4394   2.1247  17.0946
##
## (Dispersion Parameter for gaussian family taken to be 15.3543)
##
##     Null Deviance: 23818.99 on 391 degrees of freedom
## Residual Deviance: 5880.697 on 382.9999 degrees of freedom
## AIC: 2194.05
##
## Number of Local Scoring Iterations: 2
##
## Anova for Parametric Effects
##                     Df  Sum Sq Mean Sq F value  Pr(>F)
## s(displacement, 4)   1 15254.9 15254.9 993.524 < 2e-16 ***
## s(horsepower, 4)     1  1038.4  1038.4  67.632 3.1e-15 ***
## Residuals          383  5880.7    15.4
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
##                    Npar Df Npar F     Pr(F)
## (Intercept)
## s(displacement, 4)       3 13.613 1.863e-08 ***
## s(horsepower, 4)         3 15.606 1.349e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Question 7.9 - Pages 299-300

This question uses the variables dis (the weighted mean of distances to five Boston employment centers) and nox (nitrogen oxides concentration in parts per 10 million) from the Boston data. We will treat dis as the predictor and nox as the response.

(a) Use the poly() function to fit a cubic polynomial regression to predict nox using dis. Report the regression output, and plot the resulting data and polynomial fits.

(b) Plot the polynomial fits for a range of different polynomial degrees (say, from 1 to 10), and report the associated residual sum of squares.

(c) Perform cross-validation or another approach to select the optimal degree for the polynomial, and explain your results.

(d) Use the bs() function to fit a regression spline to predict nox using dis. Report the output for the fit using four degrees of freedom. How did you choose the knots? Plot the resulting fit.

(e) Now fit a regression spline for a range of degrees of freedom, and plot the resulting fits and report the resulting RSS. Describe the results obtained.

(f) Perform cross-validation or another approach in order to select the best degrees of freedom for a regression spline on this data. Describe your results.

Answer

boston<-MASS::Boston
boston %>% str %>% print

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
## NULL

boston %>% summary %>% print

##       crim                zn             indus            chas
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000
##       nox               rm             age              dis
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127
##       rad              tax           ptratio          black
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90
##      lstat            medv
##  Min.   : 1.73   Min.   : 5.00
##  1st Qu.: 6.95   1st Qu.:17.02
##  Median :11.36   Median :21.20
##  Mean   :12.65   Mean   :22.53
##  3rd Qu.:16.95   3rd Qu.:25.00
##  Max.   :37.97   Max.   :50.00

boston %>% sapply(unique) %>% sapply(length) %>% print

##    crim      zn   indus    chas     nox      rm     age     dis     rad     tax
##     504      26      76       2      81     446     356     412       9      66
## ptratio   black   lstat    medv
##      46     357     455     229

a) Polynomial

fit.03 <- lm(nox~poly(dis,3), data=boston)
dislims <- range(boston$dis)
dis.grid <- seq(dislims[1], dislims[2], 0.1)
preds <- predict(fit.03, newdata=list(dis=dis.grid), se=TRUE)
se.bands <- preds$fit + cbind(2*preds$se.fit, -2*preds$se.fit)
par(mfrow=c(1,1), mar=c(4.5,4.5,1,1), oma=c(0,0,4,0))
plot(boston$dis, boston$nox, xlim=dislims, cex=0.5, col="darkgrey")
title("Degree 3 Polynomial Fit")
lines(dis.grid, preds$fit, lwd=2, col="blue")
matlines(dis.grid, se.bands, lwd=1, col="blue", lty=3)

summary(fit.03)

##
## Call:
## lm(formula = nox ~ poly(dis, 3), data = boston)
##
## Residuals:
##       Min        1Q    Median        3Q       Max
## -0.121130 -0.040619 -0.009738  0.023385  0.194904
##
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)    0.554695   0.002759 201.021  < 2e-16 ***
## poly(dis, 3)1 -2.003096   0.062071 -32.271  < 2e-16 ***
## poly(dis, 3)2  0.856330   0.062071  13.796  < 2e-16 ***
## poly(dis, 3)3 -0.318049   0.062071  -5.124 4.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06207 on 502 degrees of freedom
## Multiple R-squared:  0.7148, Adjusted R-squared:  0.7131
## F-statistic: 419.3 on 3 and 502 DF,  p-value: < 2.2e-16

b) Multiple Polynomials

rss.error <- rep(0,10)
for (i in 1:10) {
  lm.fit <- lm(nox~poly(dis,i), data=boston)
  rss.error[i] <- sum(lm.fit$residuals^2)
}
rss.error

##  [1] 2.768563 2.035262 1.934107 1.932981 1.915290 1.878257 1.849484 1.835630
##  [9] 1.833331 1.832171

plot(rss.error, type="b")

c) Cross validation and polynomial selection

require(boot)
set.seed(1)
cv.error <- rep(0,10)
for (i in 1:10) {
  glm.fit <- glm(nox~poly(dis,i), data=boston)
  cv.error[i] <- cv.glm(boston, glm.fit, K=10)$delta[1]  # [1]:std, [2]:bias-corrected
}
cv.error

##  [1] 0.005558263 0.004085706 0.003876521 0.003863342 0.004237452 0.005686862
##  [7] 0.010278897 0.006810868 0.033308607 0.004075599

plot(cv.error, type="b")
  • I feel like the second degree fit would be the most reasonable, though the fourth degree seems to be doing well.

d) Regression spline

fit.sp <- lm(nox~bs(dis, df=4), data=boston)
pred <- predict(fit.sp, newdata=list(dis=dis.grid), se=T)
plot(boston$dis, boston$nox, col="gray")
lines(dis.grid, pred$fit, lwd=2)
lines(dis.grid, pred$fit+2*pred$se, lty="dashed")
lines(dis.grid, pred$fit-2*pred$se, lty="dashed")

# set df to select knots at uniform quantiles of `dis`
attr(bs(boston$dis,df=4),"knots")  # only 1 knot at 50th percentile

##     50%
## 3.20745

e) Range of regression splines

rss.error <- rep(0,7)
for (i in 4:10) {
  fit.sp <- lm(nox~bs(dis, df=i), data=boston)
  rss.error[i-3] <- sum(fit.sp$residuals^2)
}
rss.error

## [1] 1.922775 1.840173 1.833966 1.829884 1.816995 1.825653 1.792535

plot(4:10, rss.error, type="b")
  • As the model gains more degrees of freedom, it tends to over fit to the training data better

f) Cross validation for best spline

cv.error <- rep(0,7)
for (i in 4:10) {
  glm.fit <- glm(nox~bs(dis, df=i), data=boston)
  cv.error[i-3] <- cv.glm(boston, glm.fit, K=10)$delta[1]
}

## Warning in bs(dis, degree = 3L, knots = c(`50%` = 3.1523), Boundary.knots =
## c(1.1296, : some 'x' values beyond boundary knots may cause ill-conditioned
## bases

## Warning in bs(dis, degree = 3L, knots = c(`50%` = 3.1523), Boundary.knots =
## c(1.1296, : some 'x' values beyond boundary knots may cause ill-conditioned
## bases

## Warning in bs(dis, degree = 3L, knots = c(`50%` = 3.2157), Boundary.knots =
## c(1.137, : some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`50%` = 3.2157), Boundary.knots =
## c(1.137, : some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`33.33333%` = 2.35953333333333, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`33.33333%` = 2.35953333333333, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`33.33333%` = 2.38403333333333, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`33.33333%` = 2.38403333333333, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`25%` = 2.07945, `50%` = 3.1323, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`25%` = 2.07945, `50%` = 3.1323, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`25%` = 2.1103, `50%` = 3.2797, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`25%` = 2.1103, `50%` = 3.2797, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`20%` = 1.9682, `40%` = 2.7147, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`20%` = 1.9682, `40%` = 2.7147, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`20%` = 1.95434, `40%` = 2.59666, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`20%` = 1.95434, `40%` = 2.59666, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`16.66667%` = 1.82203333333333, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`16.66667%` = 1.82203333333333, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`16.66667%` = 1.8226, `33.33333%` =
## 2.3817, : some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`16.66667%` = 1.8226, `33.33333%` =
## 2.3817, : some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`14.28571%` = 1.7936, `28.57143%`
## = 2.16972857142857, : some 'x' values beyond boundary knots may cause ill-
## conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`14.28571%` = 1.7936, `28.57143%`
## = 2.16972857142857, : some 'x' values beyond boundary knots may cause ill-
## conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`12.5%` = 1.754625, `25%` = 2.10215, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`12.5%` = 1.754625, `25%` = 2.10215, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`12.5%` = 1.751575, `25%` = 2.08755, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

## Warning in bs(dis, degree = 3L, knots = c(`12.5%` = 1.751575, `25%` = 2.08755, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

cv.error

## [1] 0.003898810 0.003694675 0.003732665 0.003766202 0.003716389 0.003723126
## [7] 0.003727358

plot(4:10, cv.error, type="b")
  • A fifth degree polynomial is clearly indicated

Question 10 - Page 300

This question relates to the College data set.

(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.

(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.

(c) Evaluate the model obtained on the test set, and explain the results obtained.

(d) For which variables, if any, is there evidence of a non-linear relationship with the response?

Answer

colDat<-ISLR::College
colDat %>% str %>% print

## 'data.frame':    777 obs. of  18 variables:
##  $ Private    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Apps       : num  1660 2186 1428 417 193 ...
##  $ Accept     : num  1232 1924 1097 349 146 ...
##  $ Enroll     : num  721 512 336 137 55 158 103 489 227 172 ...
##  $ Top10perc  : num  23 16 22 60 16 38 17 37 30 21 ...
##  $ Top25perc  : num  52 29 50 89 44 62 45 68 63 44 ...
##  $ F.Undergrad: num  2885 2683 1036 510 249 ...
##  $ P.Undergrad: num  537 1227 99 63 869 ...
##  $ Outstate   : num  7440 12280 11250 12960 7560 ...
##  $ Room.Board : num  3300 6450 3750 5450 4120 ...
##  $ Books      : num  450 750 400 450 800 500 500 450 300 660 ...
##  $ Personal   : num  2200 1500 1165 875 1500 ...
##  $ PhD        : num  70 29 53 92 76 67 90 89 79 40 ...
##  $ Terminal   : num  78 30 66 97 72 73 93 100 84 41 ...
##  $ S.F.Ratio  : num  18.1 12.2 12.9 7.7 11.9 9.4 11.5 13.7 11.3 11.5 ...
##  $ perc.alumni: num  12 16 30 37 2 11 26 37 23 15 ...
##  $ Expend     : num  7041 10527 8735 19016 10922 ...
##  $ Grad.Rate  : num  60 56 54 59 15 55 63 73 80 52 ...
## NULL

colDat %>% summary %>% print

##  Private        Apps           Accept          Enroll       Top10perc
##  No :212   Min.   :   81   Min.   :   72   Min.   :  35   Min.   : 1.00
##  Yes:565   1st Qu.:  776   1st Qu.:  604   1st Qu.: 242   1st Qu.:15.00
##            Median : 1558   Median : 1110   Median : 434   Median :23.00
##            Mean   : 3002   Mean   : 2019   Mean   : 780   Mean   :27.56
##            3rd Qu.: 3624   3rd Qu.: 2424   3rd Qu.: 902   3rd Qu.:35.00
##            Max.   :48094   Max.   :26330   Max.   :6392   Max.   :96.00
##    Top25perc      F.Undergrad     P.Undergrad         Outstate
##  Min.   :  9.0   Min.   :  139   Min.   :    1.0   Min.   : 2340
##  1st Qu.: 41.0   1st Qu.:  992   1st Qu.:   95.0   1st Qu.: 7320
##  Median : 54.0   Median : 1707   Median :  353.0   Median : 9990
##  Mean   : 55.8   Mean   : 3700   Mean   :  855.3   Mean   :10441
##  3rd Qu.: 69.0   3rd Qu.: 4005   3rd Qu.:  967.0   3rd Qu.:12925
##  Max.   :100.0   Max.   :31643   Max.   :21836.0   Max.   :21700
##    Room.Board       Books           Personal         PhD
##  Min.   :1780   Min.   :  96.0   Min.   : 250   Min.   :  8.00
##  1st Qu.:3597   1st Qu.: 470.0   1st Qu.: 850   1st Qu.: 62.00
##  Median :4200   Median : 500.0   Median :1200   Median : 75.00
##  Mean   :4358   Mean   : 549.4   Mean   :1341   Mean   : 72.66
##  3rd Qu.:5050   3rd Qu.: 600.0   3rd Qu.:1700   3rd Qu.: 85.00
##  Max.   :8124   Max.   :2340.0   Max.   :6800   Max.   :103.00
##     Terminal       S.F.Ratio      perc.alumni        Expend
##  Min.   : 24.0   Min.   : 2.50   Min.   : 0.00   Min.   : 3186
##  1st Qu.: 71.0   1st Qu.:11.50   1st Qu.:13.00   1st Qu.: 6751
##  Median : 82.0   Median :13.60   Median :21.00   Median : 8377
##  Mean   : 79.7   Mean   :14.09   Mean   :22.74   Mean   : 9660
##  3rd Qu.: 92.0   3rd Qu.:16.50   3rd Qu.:31.00   3rd Qu.:10830
##  Max.   :100.0   Max.   :39.80   Max.   :64.00   Max.   :56233
##    Grad.Rate
##  Min.   : 10.00
##  1st Qu.: 53.00
##  Median : 65.00
##  Mean   : 65.46
##  3rd Qu.: 78.00
##  Max.   :118.00

colDat %>% sapply(unique) %>% sapply(length) %>% print

##     Private        Apps      Accept      Enroll   Top10perc   Top25perc
##           2         711         693         581          82          89
## F.Undergrad P.Undergrad    Outstate  Room.Board       Books    Personal
##         714         566         640         553         122         294
##         PhD    Terminal   S.F.Ratio perc.alumni      Expend   Grad.Rate
##          78          65         173          61         744          81

plotLEAP=function(leapObj){
  par(mfrow = c(2,2))
  bar2=which.max(leapObj$adjr2)
  bbic=which.min(leapObj$bic)
  bcp=which.min(leapObj$cp)
  plot(leapObj$rss,xlab="Number of variables",ylab="RSS",type="b")
  plot(leapObj$adjr2,xlab="Number of variables",ylab=TeX("Adjusted R^2"),type="b")
  points(bar2,leapObj$adjr2[bar2],col="green",cex=2,pch=20)
  plot(leapObj$bic,xlab="Number of variables",ylab=TeX("BIC"),type="b")
  points(bbic,leapObj$bic[bbic],col="blue",cex=2,pch=20)
  plot(leapObj$cp,xlab="Number of variables",ylab=TeX("C_p"),type="b")
  points(bcp,leapObj$cp[bcp],col="red",cex=2,pch=20)
}

a) Train test

train_ind = sample(colDat %>% nrow,100)
test_ind = setdiff(seq_len(colDat %>% nrow), train_ind)

Best subset selection

train_set<-colDat[train_ind,]
test_set<-colDat[-train_ind,]

library(leaps)

modelFit<-regsubsets(Outstate~.,data=colDat,nvmax=20)
modelFit %>% summary %>% print

## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = colDat, nvmax = 20)
## 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 17
## Selection Algorithm: exhaustive
##           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 ) "*"        "*"  "*"    "*"    "*"       "*"       "*"
## 17  ( 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 ) " "         "*"        "*"   "*"      "*" "*"      "*"
## 17  ( 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 ) "*"         "*"    "*"
## 17  ( 1 ) "*"         "*"    "*"

We might want to take a look at these.

par(mfrow=c(2,2))
plot(modelFit)
plot(modelFit,scale='Cp')
plot(modelFit,scale='r2')
plot(modelFit,scale='adjr2')

plotLEAP(modelFit %>% summary)
  • So we like 14 variables, namely
coefficients(modelFit,id=14) %>% print

##   (Intercept)    PrivateYes          Apps        Accept        Enroll
## -1.817040e+03  2.256946e+03 -2.999022e-01  8.023519e-01 -5.372545e-01
##     Top10perc   F.Undergrad    Room.Board      Personal           PhD
##  2.365529e+01 -9.569936e-02  8.741819e-01 -2.478418e-01  1.269506e+01
##      Terminal     S.F.Ratio   perc.alumni        Expend     Grad.Rate
##  2.297296e+01 -4.700560e+01  4.195006e+01  2.003912e-01  2.383197e+01
  • But five seems like a better bet.
coefficients(modelFit,id=5)

##   (Intercept)    PrivateYes    Room.Board           PhD   perc.alumni
## -2864.6325619  2936.7416766     1.0677573    40.5334088    61.3147684
##        Expend
##     0.2253945

b) GAM

library(gam)

fit = gam(Outstate ~ Private+s(Apps,3)+Accept+Enroll+
            Top10perc+F.Undergrad+Room.Board+
            Personal+PhD+Terminal+S.F.Ratio+
            perc.alumni+Expend+Grad.Rate
        , data = colDat)

## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored

fit2 = gam(Outstate ~ Private+s(Room.Board,2)+s(PhD,3)+s(perc.alumni)+Expend
        , data = colDat)

## Warning in model.matrix.default(mt, mf, contrasts): non-list contrasts argument
## ignored

par(mfrow=c(2,2))
plot(fit,se=TRUE)


par(mfrow=c(2,2))
plot(fit2,se=TRUE)

c) Evaluate

pred <- predict(fit, test_set)
mse.error <- mean((test_set$Outstate - pred)^2)
mse.error %>% print

## [1] 3691891

gam.tss = mean((test_set$Outstate - mean(test_set$Outstate))^2)
test.rss = 1 - mse.error/gam.tss
test.rss %>% print

## [1] 0.7731239

pred2 <- predict(fit2, test_set)
mse.error2 <- mean((test_set$Outstate - pred2)^2)
mse.error2 %>% print

## [1] 4121902

gam.tss2 = mean((test_set$Outstate - mean(test_set$Outstate))^2)
test.rss2 = 1 - mse.error2/gam.tss2
test.rss2 %>% print

## [1] 0.7466987

This is pretty good model, all told.

d) Summary

summary(fit) %>% print

##
## Call: gam(formula = Outstate ~ Private + s(Apps, 3) + Accept + Enroll +
##     Top10perc + F.Undergrad + Room.Board + Personal + PhD + Terminal +
##     S.F.Ratio + perc.alumni + Expend + Grad.Rate, data = colDat)
## Deviance Residuals:
##       Min        1Q    Median        3Q       Max
## -6641.083 -1262.806    -5.698  1270.911  9965.901
##
## (Dispersion Parameter for gaussian family taken to be 3749048)
##
##     Null Deviance: 12559297426 on 776 degrees of freedom
## Residual Deviance: 2849276343 on 760 degrees of freedom
## AIC: 13985.3
##
## Number of Local Scoring Iterations: 2
##
## Anova for Parametric Effects
##              Df     Sum Sq    Mean Sq  F value    Pr(>F)
## Private       1 4034912907 4034912907 1076.250 < 2.2e-16 ***
## s(Apps, 3)    1 1344548030 1344548030  358.637 < 2.2e-16 ***
## Accept        1   90544274   90544274   24.151 1.091e-06 ***
## Enroll        1  144471570  144471570   38.535 8.838e-10 ***
## Top10perc     1 1802244831 1802244831  480.721 < 2.2e-16 ***
## F.Undergrad   1   45230645   45230645   12.065 0.0005430 ***
## Room.Board    1 1110285773 1110285773  296.151 < 2.2e-16 ***
## Personal      1   47886988   47886988   12.773 0.0003738 ***
## PhD           1  220249039  220249039   58.748 5.476e-14 ***
## Terminal      1   66366007   66366007   17.702 2.892e-05 ***
## S.F.Ratio     1  190811028  190811028   50.896 2.274e-12 ***
## perc.alumni   1  225293653  225293653   60.094 2.904e-14 ***
## Expend        1  258162295  258162295   68.861 4.805e-16 ***
## Grad.Rate     1   57947219   57947219   15.457 9.214e-05 ***
## Residuals   760 2849276343    3749048
## ---
## 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(Apps, 3)        2  8.571 0.0002085 ***
## Accept
## Enroll
## Top10perc
## F.Undergrad
## Room.Board
## Personal
## PhD
## Terminal
## S.F.Ratio
## perc.alumni
## Expend
## Grad.Rate
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit2) %>% print

##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, 2) + s(PhD,
##     3) + s(perc.alumni) + Expend, data = colDat)
## Deviance Residuals:
##       Min        1Q    Median        3Q       Max
## -8676.030 -1345.678    -8.409  1265.524  9590.459
##
## (Dispersion Parameter for gaussian family taken to be 4175193)
##
##     Null Deviance: 12559297426 on 776 degrees of freedom
## Residual Deviance: 3194023899 on 765.0002 degrees of freedom
## AIC: 14064.05
##
## Number of Local Scoring Iterations: 2
##
## Anova for Parametric Effects
##                   Df     Sum Sq    Mean Sq F value    Pr(>F)
## Private            1 3751107814 3751107814  898.43 < 2.2e-16 ***
## s(Room.Board, 2)   1 2913770756 2913770756  697.88 < 2.2e-16 ***
## s(PhD, 3)          1 1149711330 1149711330  275.37 < 2.2e-16 ***
## s(perc.alumni)     1  556759894  556759894  133.35 < 2.2e-16 ***
## Expend             1  554812125  554812125  132.88 < 2.2e-16 ***
## Residuals        765 3194023899    4175193
## ---
## 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, 2)       1 4.9853 0.0258517 *
## s(PhD, 3)              2 9.1614 0.0001171 ***
## s(perc.alumni)         3 0.8726 0.4548496
## Expend
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1