gam package
world = supervised learning
We want to model some output variable \(y\) using a set of potential predictors (\(x_1, x_2, ..., x_p\)) and the relationship could be non-linear.
task = regression
\(y\) is quantitative
a flexible algorithm
Our standard parametric models (eg: linear regression) may too rigid to represent the relationship between \(y\) and our predictors \(x\). Thus we need more flexible models.
Discuss the following questions with your group to review key concepts from recent activities.
In the previous activity, we modeled college Grad.Rate versus Expend, Enroll, and Private using data on 775 schools.
Which pre-processing steps are necessary for KNN, and why? How do we implement these pre-processing steps in R?
What is meant by the “curse of dimensionality” in KNN?
What assumptions did the KNN model make about the relationship of Grad.Rate with Expend, Enroll, and Private? Is this a pro or con?
What did the KNN model tell us about the relationship of Grad.Rate with Expend, Enroll, and Private? For example, did it give you a sense of whether grad rates are higher at private or public institutions? At institutions with higher or lower enrollments? Is this a pro or con?
Nonparametric KNN vs parametric least squares and LASSO:
LOESS: Local Regression or Locally Estimated Scatterplot Smoothing
Goal:
Build a flexible regression model of \(y\) by one quantitative predictor \(x\),
\[y = f(x) + \varepsilon\]
Idea:
Fit regression models in small localized regions, where nearby data have greater influence than far data.
Algorithm:
Define the span, aka bandwidth, tuning parameter \(h\) where \(0 \le h \le 1\). Take the following steps to estimate \(f(x)\) at each possible predictor value \(x\):
In pictures:
library(tidyverse)
# example data
plot_data <- read.csv("https://bcheggeseth.github.io/253_spring_2024/data/glucose_experiment.csv")
# plot with LOESS predictions
plot_data %>%
ggplot() +
geom_point(aes(x = time, y = glucose)) +
geom_smooth(aes(x = time, y = glucose), se=FALSE) # this step does LOESS
Open the QMD and find the LOESS discussion questions. Work on the following questions with your group.
We can plot LOESS models using geom_smooth(). Play around with the span parameter below.
library(tidyverse)
library(ISLR)
data(College)
ggplot(College, aes(x = Expend, y = Grad.Rate)) +
geom_point() +
geom_smooth(method = "loess", span = 0.5)
Note: you’ll find that you can specify span greater than 1. Use your resources to figure out what that means in terms of the algorithm.
Run the shiny app code, below. Explore the impact of the span tuning parameter h on the LOESS performance across different datasets. Continue to click the Go! button to get different datasets.
For what values of \(h\) (near 0, somewhere in the middle, near 1) do you get the following:
# Load packages & data
library(shiny)
library(tidyverse)
# Define a LOESS plotting function
plot_loess <- function(h, plot_data){
ggplot(plot_data, aes(x = x, y = y)) +
geom_point() +
geom_smooth(span = h, se = FALSE) +
labs(title = paste("h = ", h)) +
lims(y = c(-5,12))
}
# BUILD THE SERVER
# These are instructions for building the app - what plot to make, what quantities to calculate, etc
server_LOESS <- function(input, output) {
new_data <- eventReactive(input$do, {
x <- c(runif(25, -6, -2.5), runif(50, -3, 3), runif(25, 2.5, 6))
y <- 9 - x^2 + rnorm(100, sd = 0.75)
y[c(1:25, 76:100)] <- rnorm(50, sd = 0.75)
data.frame(x, y)
})
output$loesspic <- renderPlot({
plot_loess(h = input$hTune, plot_data = new_data())
})
}
# BUILD THE USER INTERFACE (UI)
# The UI controls the layout, appearance, and widgets (eg: slide bars).
ui_LOESS <- fluidPage(
sidebarLayout(
sidebarPanel(
h4("Sample 100 observations:"),
actionButton("do", "Go!"),
h4("Tune the LOESS:"),
sliderInput("hTune", "h", min = 0.05, max = 1, value = 0.05)
),
mainPanel(
h4("LOESS plot:"),
plotOutput("loesspic")
)
)
)
# RUN THE SHINY APP!
shinyApp(ui = ui_LOESS, server = server_LOESS)
Goal:
Build a flexible parametric regression model of \(y\) by one quantitative predictor \(x\),
\[y = f(x) + \varepsilon\]
Idea:
Fit polynomial models in small localized regions and constrain to make it smooth enough.
Definition:
A spline is a piecewise polynomial of degree \(k\) that is a continuous function and has continuous derivatives (of orders \(1,...,k-1\)) at its knot points (break points).
A spline function \(f\) with knot points \(t_1<...<t_m\),
More Definitions:
Spline functions can be represented using a basis and can be used within a linear regression framework. We often call this regression splines or B-splines and represent the spline function as
\[f(x) = \beta_0 + \beta_1f_1(x) + \cdots + \beta_{p}f_{p}(x)\]
where the functions \(f_1(x)\),…,\(f_p(x)\) are complex basis functions defined by knots and the degree of the polynomial.
Check out Section 7.4 in ISLR for more details!
Piecewise Polynomials in Stat 155
If you had a quadratic model (degree = 2) with 0 knots, how could you fit that model using the tools of STAT 155?
lm(y ~ x + x2, data) #if you created x2 = x^2
lm(y ~ poly(x, 2), data)If you had a piecewise linear model (degree = 1) with 1 knot at \(x=5\), how could you fit that using the tools from Stat 155?
lm(y ~ x * I(x > 5), data)These implementation methods don’t constrain them to be continuous or smooth, so we need additional tools from the splines package!
Splines in tidymodels
To build models with splines in tidymodels, we proceed with the same structure as we use for ordinary linear regression models but we’ll add some pre-processing steps to our recipe.
To work with splines, we’ll use tools from the splines package in conjunction with recipes.
step_ns() function based on ns() in that package creates the transformations needed to create a natural cubic spline function for a quantitative predictor.step_bs() function based on bs() in that package creates the transformations needed to create a basis spline (typically called B-splines) function for a quantitative predictor.# Linear Regression Model Spec
lm_spec <-
linear_reg() %>%
set_mode('regression') %>%
set_engine(engine = 'lm')
# Original Recipe
lm_recipe <- recipe(y ~ x, data = sample_data)
# Polynomial Recipe
poly_recipe <- lm_recipe %>%
step_poly(x, degree = __) # polynomial (higher degree means more flexible)
# Natural Spline Recipe
ns2_recipe <- lm_recipe %>%
step_ns(x, deg_free = __) # natural cubic spline (higher deg_free means more knots)
# Basis Spline Recipe
bs_recipe <- lm_recipe %>%
step_bs(x, options = list(knots = c(__,__))) # b-spline (cubic by default, more knots means higher deg_free)deg_free argument in step_ns() stands for degrees of freedom:
deg_free = # of knots + 1options argument in step_bs() allows you to specific specific knots:
options = list(knots = c(2, 4)) sets a knot at 2 and a knot at 4 (the choice of knots should depend on the observed range of the quantitative predictor and where you see a change in the relationship)
What is the difference between natural splines and B-splines?
See the figure below, for either natural splines or cubic splines, what values of deg_free (small, medium, large) do you get the following:
Could we use LASSO with a spline or polynomial recipe? Yes or No? Explain.
Goals
Context
Using the College data in the ISLR package, we’ll build a predictive model of Grad.Rate by 6 predictors: Private, PhD (percent of faculty with PhDs), Room.Board, Personal (estimated personal spending), Outstate (out-of-state tuition), perc.alumni (percent of alumni who donate)
Before proceeding, install the splines package by entering install.packages("splines") in the Console.
# Load packages
library(tidymodels)
library(tidyverse)
library(ISLR)
# Load data
data(College)
# Wrangle the data
college_sub <- College %>%
filter(Grad.Rate <= 100) %>%
filter((Grad.Rate > 50 | Expend < 40000)) %>%
select(Grad.Rate, Private, PhD, Room.Board, Personal, Outstate, perc.alumni)
We will model Grad.Rate as a function of the 6 predictors in college_sub.
ggplot(___, aes(___)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_classic()tidymodels to fit a linear regression model (no splines yet) with the following specifications:
Grad.Rate?set.seed(253)
# Lasso Model Spec with tune
lasso_spec <- linear_reg() %>%
set_mode('regression') %>%
set_engine(engine = 'glmnet') %>%
set_args(mixture = 1, penalty = tune())
# Recipe
lasso_recipe <- recipe(__ ~ ___, data = college_sub) %>%
step_dummy(all_nominal_predictors())
# Workflow (Recipe + Model)
lasso_wf_tune <- workflow() %>%
add_recipe(lasso_recipe) %>%
add_model(lasso_spec)
# Tune Model (trying a variety of values of Lambda penalty)
lasso_models <- tune_grid(
lasso_wf_tune, # workflow
resamples = vfold_cv(__, v = __),
metrics = metric_set(___),
grid = grid_regular(penalty(range = c(-3, 1)), levels = 30)
)
# Select parsimonious model & fit
parsimonious_penalty <- lasso_models %>%
select_by_one_std_err(metric = 'mae', desc(penalty))
# Finalize Model
ls_mod <- lasso_wf_tune %>%
finalize_workflow(parameters = parsimonious_penalty) %>%
fit(data = college_sub)
# Note which variable is the "least" important
ls_mod %>%
tidy()ls_mod and then plots of the residuals v. each of the quantitative predictors to evaluate whether the model is “right”.ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(__)) +
___ +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
We’ll compare our best LASSO model with models with spline functions for some of the quantitative predictors.
lm engine rather than glmnet) with the all predictors. Allow the quantitative predictors that had some non-linear relationships to be modeled with natural splines that have 4 degrees of freedom. Fit this model with 8 fold CV, fit_resamples (same folds as before) to compare MAE and then fit the model to the whole training data. Call this fit model ns_mod.# Model Spec
lm_spec <-
linear_reg() %>%
set_mode('regression') %>%
set_engine(engine = 'lm')
# New Recipe (remove steps needed for LASSO, add splines)
# Create Workflow (Recipe + Model)
set.seed(253)
# CV to Evaluate
ns_models <- fit_resamples(
___, # workflow
resamples = vfold_cv(__, v = 8), # cv folds
metrics = metric_set(___)
)
ns_models %>% collect_metrics()
# Fit with all data
ns_mod <- fit( ___, #workflow from above
data = college_sub
)ns_mod and then plots of the residuals v. each of the quantitative predictors to evaluate whether the model is right.# Residual plotslasso_models %>%
collect_metrics() %>%
filter(penalty == (parsimonious_penalty %>% pull(penalty)))
ns_models %>%
collect_metrics()
We can visualize the estimated nonlinear relationships between Grad.Rate and PhD & Personal, with a little bit of extra work.
Install the gam package with install.packages("gam") in the Console.
We fit the spline model with lm() for compatibility with the gam package.
lm() code below to see how we translated the spline recipe steps.We can use summary() to display the estimated coefficients.
ns(Phd...) and ns(Personal...) terms?PrivateYes and Outstate?Describing the relationship between the PhD and Personal predictors and the Grad.Rate outcome is better done with pictures than words! (Why?)
plot.Gam() creates plots of the relationship between each predictor and the outcome HOLDING THE OTHER VARIABLES CONSTANT. The thicker purple line shows this relationship, and the dotted lines show pointwise 2 standard error bands.Private, Room.Board, Outstate, and perc.alumni relate to the summary() output?PhD plot as follows: Among colleges with the same type (private vs public), room and board costs, estimated personal spending, out-of-state tuition, and percent of donating alumni, average graduation rates tend to increase with faculty PhD percentage more rapidly below 60%, and the relationship seems to flatten out above 60%.PhD example, try interpreting the Personal plot.library(gam) # For the plot.Gam() function
# Fit a spline model using lm()
ns_mod_lm <- lm(Grad.Rate ~ Private + ns(PhD, df=4) + Room.Board + ns(Personal, df=4) + Outstate + perc.alumni, data = college_sub)
# Look at the estimated coefficients
summary(ns_mod_lm)
# Plot the relationships between each predictor and
# the outcome ADJUSTING FOR THE OTHER PREDICTORS IN THE MODEL
plot.Gam(ns_mod_lm, se = TRUE, col = "orchid", lwd = 3)
Pick an algorithm
Our overall goal is to build a predictive model of Grad.Rate (\(y\)) by 6 predictors:
\[y = f(x_1,x_2,...,x_{6}) + \varepsilon\]
The table below summarizes the CV MAE for 3 possible algorithms: least squares, KNN, and splines. After examining the results, explain which model you would choose and why. NOTE: There are no typos here! All models had virtually the same CV MAE. Code is provided in the online manual if you’re curious.
| method | type | assumption about \(f(x_1,x_2,...,x_{6})\) | CV MAE |
|---|---|---|---|
| least squares | parametric | \(\beta_0 + \beta_1x_1 + \cdots + \beta_{6} x_{6}\) | 10.2 |
| LASSO w/ \(\lambda=1.2\) | parametric | \(\beta_0 + \beta_1x_1 + \cdots + \beta_{6} x_{6}\) | 10.3 |
| KNN w/ \(K = 33\) | nonparametric | average \(y\) among neighbors | 10.2 |
| splines | parametric | \(\beta_0 + f_1(x_1) + \cdots + f_{6}(x_{6})\) | 10.2 |
In these final exercises, you’ll reflect upon some algorithms we’ve learned thus far. This might feel bumpy if you haven’t reviewed the material recently, so be kind to yourself!
Some of these exercises are similar to those on Homework 3. Thus solutions are not provided. The following directions are important to deepening your learning and growth.
Convince yourself that the placement of Subset Selection (e.g. backward stepwise), LASSO, Least Squares, and Splines (which is highly related to Generalized Additive Models) make sense. (We’ll address other algorithms later this semester.)
Where would you place KNN on this graphic?
KNN, regression splines, and LOESS aren’t the only flexible regression algorithms! Chapter 7 of ISLR also discusses:
Read Chapter 7 of the book, and watch Prof Johnson’s video on these topics to learn more:
Chapter 7.7 specifically provides more information about generalized additive models (GAM). These models take the ideas of splines and loess and combine them into one framework and allow each predictor to have a non-linear relationship. They are implemented in tidymodels but can be computationally intensive. They allow the relationships to either be estimated with loess or smoothing splines, both of which are non-parameter approaches.
In general, we can think of GAM as
\[y = \beta_0 + f_1(x_1) + \cdots + f_{p}(x_{p}) + \epsilon\]
where the functions \(f_j(x)\) are estimated nonparametrically.
Chapter 7.5 provides more information about smoothing splines. In general, these techniques estimate the \(f(x)\) by penalizing wiggly behavior where wiggliness is measured by the second derivative of \(f(x)\), i.e. the degree to which the slope changes. Specifically, \(f(x)\) is estimated to minimize:
\[\sum_{i=1}^n (y_i - f(x_i))^2 + \lambda \int f''(t)^2 dt\]
That is, we want \(f(x)\) to produce small squared residuals (\(\sum_{i=1}^n (y_i - f(x_i))^2\)) and small total change in slope over the range of x (\(\int f''(t)^2 dt\)).
We need to standardize quantitative predictor variables (step_normalize(all_numeric_predictors())) so that our distance calculations aren’t dominated by variables on larger scales than others. We need to create dummy variables for categorical predictors (step_dummy(all_nominal_predictors())) so that we have a way of turning words into numbers and then quantify distances between categories. We might also use step_nzv() to remove predictors with near-zero variance prior to these steps; this removes variables that don’t contain much information and won’t be helpful for predicting our outcome. See Exercises 1–4, and their corresponding solutions, from last class for more discussion.
In a high dimensional space (when distances between cases are calculated by A LOT of predictors), nearest neighbors are much farther away and thus may not actually be “alike enough” to give a good prediction of the outcome.
The only assumption we made was that the outcome values of \(y\) should be similar if the predictor values of \(x\) are similar. No other assumptions are made.
Nothing. I’d say this is a con. There is nothing to interpret, so the model is more of a black box in terms of knowing why it gives you a particular prediction. See Exercises 7–9 from last class.
Further details are in solutions for previous activity.
Use nonparametric methods when parametric model assumptions are too rigid. Forcing a parametric method in this situation can produce misleading conclusions, bad predictions, etc.
Use parametric methods when the model assumptions hold. In such cases, parametric models provide more contextual insight (eg: meaningful coefficients) and the ability to detect which predictors are beneficial to the model.
deg_freedeg_freedeg_freeNot necessarily or at least we’d want to proceed carefully…
#a
# There might be a slight non-linear relationship
p1 <- ggplot(college_sub, aes(y = Grad.Rate, x = PhD)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_classic()
# Only deviation in the extremes of Room and Board
p2 <- ggplot(college_sub, aes(y = Grad.Rate, x = Room.Board)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_classic()
# linear relationship doesn't fit well but some extreme values in Personal
p3 <- ggplot(college_sub, aes(y = Grad.Rate, x = Personal)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_classic()
# Close to linear
p4 <- ggplot(college_sub, aes(y = Grad.Rate, x = Outstate)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_classic()
# Close to linear
p5 <- ggplot(college_sub, aes(y = Grad.Rate, x = perc.alumni)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_classic()
library(gridExtra)
grid.arrange(p1, p2, p3, p4, p5, nrow = 2)
The relationships between many of these predictors and
Grad.Ratelook reasonably linear. Perhaps two exceptions to this arePersonalandPhD.
#b
set.seed(253)
# Lasso Model Spec with tune
lasso_spec <- linear_reg() %>%
set_mode('regression') %>%
set_engine(engine = 'glmnet') %>%
set_args(mixture = 1, penalty = tune())
# Recipe
lasso_recipe <- recipe(Grad.Rate ~ ., data = college_sub) %>% # use all predictors ( ~ . )
step_dummy(all_nominal_predictors())
# Workflow (Recipe + Model)
lasso_wf_tune <- workflow() %>%
add_recipe(lasso_recipe) %>%
add_model(lasso_spec)
# Tune Model (trying a variety of values of Lambda penalty)
lasso_models <- tune_grid(
lasso_wf_tune, # workflow
resamples = vfold_cv(college_sub, v = 8), # 8-fold CV
metrics = metric_set(mae), # calculate MAE
grid = grid_regular(penalty(range = c(-5, 2)), levels = 30)
)
# Select parsimonious model & fit
parsimonious_penalty <- lasso_models %>%
select_by_one_std_err(metric = 'mae', desc(penalty))
ls_mod <- lasso_wf_tune %>%
finalize_workflow(parameters = parsimonious_penalty) %>%
fit(data = college_sub)
# Note which variable is the "least" important
ls_mod %>%
tidy()# A tibble: 7 × 3
term estimate penalty
<chr> <dbl> <dbl>
1 (Intercept) 37.1 1.17
2 PhD 0.0582 1.17
3 Room.Board 0.00102 1.17
4 Personal -0.000740 1.17
5 Outstate 0.00136 1.17
6 perc.alumni 0.285 1.17
7 Private_Yes 0 1.17
Private_Yeshas a coefficient estimate of zero, so is “least” important”.
# c
# Residuals v. Fitted/Predicted
r1 <- ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = .pred, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. PhD
r2 <- ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = PhD, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. Personal
r3 <- ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = Personal, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. Room.Board
r4 <- ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = Room.Board, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. Outstate
r5 <- ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = Outstate, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. perc.alumni
r6 <- ls_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = perc.alumni, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
grid.arrange(r1, r2, r3, r4, r5, r6, nrow = 2)
These residual plots look generally fine to me. Perhaps there is a bit of heteroskedasticity in the residual vs fitted plot. We also look to be over-estimating (residuals < 0) graduation rates for schools with low numbers of PhDs and mid-range levels (2000–5000) of personal spending.
# Model Spec
lm_spec <-
linear_reg() %>%
set_mode('regression') %>%
set_engine(engine = 'lm')
# New Recipe (remove steps needed for LASSO, add splines)
spline_recipe <- recipe(Grad.Rate ~ ., data = college_sub) %>%
step_ns(Personal, deg_free = 4) %>% # splines with 4 df for Personal
step_ns(PhD, deg_free = 4) # splines with 4 df for PhD
# Create Workflow (Recipe + Model)
spline_wf <- workflow() %>%
add_recipe(spline_recipe) %>%
add_model(lm_spec)
set.seed(253)
# CV to Evaluate
ns_models <- fit_resamples(
spline_wf, # workflow
resamples = vfold_cv(college_sub, v = 8), # cv folds
metrics = metric_set(mae) # calculate MAE
)
# check out CV metrics
ns_models %>%
collect_metrics()# A tibble: 1 × 6
.metric .estimator mean n std_err .config
<chr> <chr> <dbl> <int> <dbl> <chr>
1 mae standard 10.2 8 0.229 Preprocessor1_Model1# Fit with all data
ns_mod <- fit( spline_wf, #workflow from above
data = college_sub
)
#b. Residual plots
# Residuals v. Fitted/Predicted
sr1 <- ns_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = .pred, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. PhD
sr2 <- ns_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = PhD, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. Personal
sr3 <- ns_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = Personal, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. Room.Board
sr4 <- ns_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = Room.Board, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. Outstate
sr5 <- ns_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = Outstate, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
# Residuals v. perc.alumni
sr6 <- ns_mod %>%
augment(new_data = college_sub) %>%
mutate(.resid = Grad.Rate - .pred) %>%
ggplot(aes(x = perc.alumni, y = .resid)) +
geom_point() +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, color = "red") +
theme_classic()
grid.arrange(sr1, sr2, sr3, sr4, sr5, sr6, nrow = 2)
I still notice the same mild heteroskedasticity in the residual vs fitted plot. However, the systematic over-estimation of graduation rates for schools with low numbers of PhDs and mid-range personal spending now looks better (the “typical” residual is closer to zero over those ranges).
#c. compare CV MAE
lasso_models %>%
collect_metrics() %>%
filter(penalty == (parsimonious_penalty %>% pull(penalty)))
ns_models %>%
collect_metrics()The MAEs are very similar.
Note how we used ns(VARIABLE, df=DEG_FREE) instead of step_ns(..., deg_free).
PrivateYes interpretation: Among colleges with the same percent of faculty with a PhD, room and board costs, estimated personal spending, out-of-state tuition, and percent of donating alumni, average graduation rates at private colleges are about 9.8% higher than at public colleges.
Outstate interpretation: Among colleges that are the same in other respects, every $1000 increase in out-of-state tuition is associated with 1.3% higher average graduation rates.
Describing the relationship between the PhD and Personal predictors and the Grad.Rate outcome is better done with pictures than words because of the more complex nonlinear relationship that was estimated–we can’t sum it up in a slope for example.
The plots for Private, Room.Board, Outstate, and perc.alumni reflect the coefficients (slopes) in the summary() output.
Personal plot interpretation: among colleges that are the same in other respects, average graduation rates seem to be only weakly related to personal spending below about $1500. Average graduation rates tend to decrease more noticeably from about $1500 to $3000.
library(gam) # For the plot.Gam() function
# Fit a spline model using lm()
ns_mod_lm <- lm(Grad.Rate ~ Private + ns(PhD, df=4) + Room.Board + ns(Personal, df=4) + Outstate + perc.alumni, data = college_sub)
# Look at the estimated coefficients
summary(ns_mod_lm)
Call:
lm(formula = Grad.Rate ~ Private + ns(PhD, df = 4) + Room.Board +
ns(Personal, df = 4) + Outstate + perc.alumni, data = college_sub)
Residuals:
Min 1Q Median 3Q Max
-48.255 -8.001 -0.499 7.960 51.116
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.352e+01 6.635e+00 3.544 0.000418 ***
PrivateYes 9.809e-01 1.532e+00 0.640 0.522068
ns(PhD, df = 4)1 1.466e+01 5.096e+00 2.877 0.004130 **
ns(PhD, df = 4)2 1.302e+01 3.876e+00 3.360 0.000818 ***
ns(PhD, df = 4)3 3.177e+01 1.161e+01 2.738 0.006332 **
ns(PhD, df = 4)4 8.730e+00 4.410e+00 1.980 0.048114 *
Room.Board 1.627e-03 5.927e-04 2.746 0.006179 **
ns(Personal, df = 4)1 1.093e+00 3.097e+00 0.353 0.724235
ns(Personal, df = 4)2 -1.407e+01 4.214e+00 -3.338 0.000885 ***
ns(Personal, df = 4)3 -8.126e+00 8.593e+00 -0.946 0.344595
ns(Personal, df = 4)4 5.259e+00 1.096e+01 0.480 0.631565
Outstate 1.278e-03 2.213e-04 5.775 1.12e-08 ***
perc.alumni 3.212e-01 4.880e-02 6.582 8.61e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.24 on 762 degrees of freedom
Multiple R-squared: 0.4088, Adjusted R-squared: 0.3995
F-statistic: 43.91 on 12 and 762 DF, p-value: < 2.2e-16# Plot the relationships between each predictor and
# the outcome ADJUSTING FOR THE OTHER PREDICTORS IN THE MODEL
par(mfrow = c(2,3))
plot.Gam(ns_mod_lm, se = TRUE, col = "orchid", lwd = 3)
You can argue a couple of ways.
If you noted the non-linear relationships in the plots above: splines is more informative than KNN and less wrong than least squares.
If you didn’t think above that splines were warranted, you should pick least squares (or LASSO). It’s definitely simpler and easier to interpret.
Code
# least squares
set.seed(253)
ls_model <- lm_spec %>% # re-use lm specs from above
# don't need a recipe since we're not doing any pre-processing
fit_resamples( # get 10-fold CV MAE
Grad.Rate ~ .,
resamples = vfold_cv(college_sub, v = 10),
metrics = metric_set(mae)
)
# look at CV MAE
ls_model %>%
collect_metrics()# A tibble: 1 × 6
.metric .estimator mean n std_err .config
<chr> <chr> <dbl> <int> <dbl> <chr>
1 mae standard 10.2 10 0.365 Preprocessor1_Model1# KNN
# specifications
knn_spec <- nearest_neighbor() %>%
set_mode("regression") %>%
set_engine(engine = "kknn") %>%
set_args(neighbors = tune())
# recipe
variable_recipe <- recipe(Grad.Rate ~ ., data = college_sub) %>%
step_nzv(all_predictors()) %>%
step_dummy(all_nominal_predictors()) %>%
step_normalize(all_numeric_predictors())
# workflow
knn_workflow <- workflow() %>%
add_model(knn_spec) %>%
add_recipe(variable_recipe)
# tune K using 10-fold CV
set.seed(253)
knn_models <- knn_workflow %>%
tune_grid(
grid = grid_regular(neighbors(range = c(1, 200)), levels = 50),
resamples = vfold_cv(college_sub, v = 10),
metrics = metric_set(mae)
)
# select best K
best_K <- select_best(knn_models, metric = "mae")
best_K# A tibble: 1 × 2
neighbors .config
<int> <chr>
1 33 Preprocessor1_Model09# fit final model
knn_models %>%
collect_metrics() %>%
filter(neighbors == best_K$neighbors)# A tibble: 1 × 7
neighbors .metric .estimator mean n std_err .config
<int> <chr> <chr> <dbl> <int> <dbl> <chr>
1 33 mae standard 10.2 10 0.337 Preprocessor1_Model09Part 2 (4–7): Solutions are not provided since these questions are part of HW3.