Chapter 7: Multivariate Adaptive Regression Splines

Note: Some results may differ from the hard copy book due to the changing of sampling procedures introduced in R 3.6.0. See http://bit.ly/35D1SW7 for more details. Access and run the source code for this notebook here.

Hidden chapter requirements used in the book to set the plotting theme and load packages used in hidden code chunks:

# Set the graphical theme
ggplot2::theme_set(ggplot2::theme_light())

# Set global knitr chunk options
knitr::opts_chunk$set(
  fig.align = "center",
  fig.height = 3.5,
  warning = FALSE, 
  message = FALSE
)

# Load required packages
library(dplyr)

Prerequisites

For this chapter we will use the following packages:

# Helper packages
library(dplyr)     # for data wrangling
library(ggplot2)   # for awesome plotting

# Modeling packages
library(earth)     # for fitting MARS models
library(caret)     # for automating the tuning process

# Model interpretability packages
library(vip)       # for variable importance
library(pdp)       # for variable relationships

To illustrate various concepts we’ll continue with the ames_train and ames_test data sets:

# Stratified sampling with the rsample package
set.seed(123)
ames <- AmesHousing::make_ames()
split  <- rsample::initial_split(ames, prop = 0.7, strata = "Sale_Price")
ames_train  <- rsample::training(split)
ames_test   <- rsample::testing(split)

The basic idea

Figure 7.1:

set.seed(123)  # for reproducibility
x <- seq(from = 0, to = 2 * pi, length = 500)
y <- sin(x) + rnorm(length(x), sd = 0.3)
df <- data.frame(x, y) %>%
  filter(x < 6)

p1 <- ggplot(df, aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle("(A) Assumed linear relationship")

p2 <- ggplot(df, aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  stat_smooth( method = "lm", se = FALSE, formula = y ~ poly(x, 2, raw = TRUE)) +
  ggtitle("(B) Degree-2 polynomial regression")

p3 <- ggplot(df, aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  stat_smooth( method = "lm", se = FALSE, formula = y ~ poly(x, 3, raw = TRUE)) +
  ggtitle("(C) Degree-3 polynomial regression")

# fit step function model (6 steps)
step_fit <- lm(y ~ cut(x, 5), data = df)
step_pred <- predict(step_fit, df)

p4 <- ggplot(cbind(df, step_pred), aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  geom_line(aes(y = step_pred), size = 1, color = "blue") +
  ggtitle("(D) Step function regression")

gridExtra::grid.arrange(p1, p2, p3, p4, nrow = 2)

Multivariate adaptive regression splines

Figure 7.2:

# one knot
mars1 <- mda::mars(df$x, df$y, nk = 3, prune = FALSE)
p1 <- df %>%
  mutate(predicted = as.vector(mars1$fitted.values)) %>%
  ggplot(aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  geom_line(aes(y = predicted), size = 1, color = "blue") +
  ggtitle("(A) One knot")

# two knots
mars2 <- mda::mars(df$x, df$y, nk = 5, prune = FALSE)
p2 <- df %>%
  mutate(predicted = as.vector(mars2$fitted.values)) %>%
  ggplot(aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  geom_line(aes(y = predicted), size = 1, color = "blue") +
  ggtitle("(B) Two knots")

mars3 <- mda::mars(df$x, df$y, nk = 7, prune = FALSE)
p3 <- df %>%
  mutate(predicted = as.vector(mars3$fitted.values)) %>%
  ggplot(aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  geom_line(aes(y = predicted), size = 1, color = "blue") +
  ggtitle("(C) Three knots")

mars4 <- mda::mars(df$x, df$y, nk = 9, prune = FALSE)
p4 <- df %>%
  mutate(predicted = as.vector(mars4$fitted.values)) %>%
  ggplot(aes(x, y)) +
  geom_point(size = 1, alpha = .2) +
  geom_line(aes(y = predicted), size = 1, color = "blue") +
  ggtitle("(D) Four knots")


gridExtra::grid.arrange(p1, p2, p3, p4, nrow = 2)

Fitting a basic MARS model

# Fit a basic MARS model
mars1 <- earth(
  Sale_Price ~ .,  
  data = ames_train   
)

# Print model summary
print(mars1)

Selected 36 of 40 terms, and 28 of 307 predictors
Termination condition: RSq changed by less than 0.001 at 40 terms
Importance: Gr_Liv_Area, Year_Built, Total_Bsmt_SF, Overall_QualExcellent, ...
Number of terms at each degree of interaction: 1 35 (additive model)
GCV 547654257    RSS 1.047912e+12    GRSq 0.9150216    RSq 0.9207205

summary(mars1) %>% .$coefficients %>% head(10)

                             Sale_Price
(Intercept)                234803.84945
h(2787-Gr_Liv_Area)           -50.67882
h(Year_Built-2004)           3595.39940
h(2004-Year_Built)           -373.40103
h(Total_Bsmt_SF-1298)          56.30424
h(1298-Total_Bsmt_SF)         -29.95470
h(Bsmt_Unf_SF-536)            -24.47153
h(536-Bsmt_Unf_SF)             16.28784
Overall_QualExcellent       79769.47075
Overall_QualVery_Excellent 117138.64127

plot(mars1, which = 1)

# Fit a basic MARS model
mars2 <- earth(
  Sale_Price ~ .,  
  data = ames_train,
  degree = 2
)

# check out the first 10 coefficient terms
summary(mars2) %>% .$coefficients %>% head(10)

                                           Sale_Price
(Intercept)                              2.568344e+05
h(Gr_Liv_Area-2787)                     -1.449341e+02
h(2787-Gr_Liv_Area)                     -4.966604e+01
h(Year_Built-2004)                       4.548491e+03
h(2004-Year_Built)                      -7.249201e+02
h(Total_Bsmt_SF-1298)                    8.060668e+01
h(1298-Total_Bsmt_SF)                   -4.196757e+01
h(Bsmt_Unf_SF-1017)*h(2787-Gr_Liv_Area) -2.352446e-02
h(1017-Bsmt_Unf_SF)*h(2787-Gr_Liv_Area)  8.535269e-03
Condition_1Norm*h(Gr_Liv_Area-2787)      2.785045e+02

Tuning

# create a tuning grid
hyper_grid <- expand.grid(
  degree = 1:3, 
  nprune = seq(2, 100, length.out = 10) %>% floor()
)

head(hyper_grid)

# Cross-validated model
set.seed(123)  # for reproducibility
cv_mars <- train(
  x = subset(ames_train, select = -Sale_Price),
  y = ames_train$Sale_Price,
  method = "earth",
  metric = "RMSE",
  trControl = trainControl(method = "cv", number = 10),
  tuneGrid = hyper_grid
)

# View results
cv_mars$bestTune

ggplot(cv_mars)

Table 7.1:

Cross-validated RMSE results for tuned MARS and regression models.

	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.	NA's
LM	20844.33	22581.04	24947.45	26098.00	27695.65	39521.49	0
PCR	28983.85	31243.24	31723.22	34042.15	37152.29	42528.87	0
PLS	20918.60	22324.61	24548.09	25459.51	25448.77	39084.91	0
ENET	20204.50	21264.34	23761.18	24721.83	25594.41	38721.64	0
MARS	20698.98	21730.78	22878.15	27899.28	24879.49	70420.90	0

Feature interpretation

# variable importance plots
p1 <- vip(cv_mars, num_features = 40, bar = FALSE, value = "gcv") + ggtitle("GCV")
p2 <- vip(cv_mars, num_features = 40, bar = FALSE, value = "rss") + ggtitle("RSS")

gridExtra::grid.arrange(p1, p2, ncol = 2)

# extract coefficients, convert to tidy data frame, and
# filter for interaction terms
cv_mars$finalModel %>%
  coef() %>%  
  broom::tidy() %>%  
  filter(stringr::str_detect(names, "\\*")) 

# Construct partial dependence plots
p1 <- partial(cv_mars, pred.var = "Gr_Liv_Area", grid.resolution = 10) %>% 
  autoplot()

p2 <- partial(cv_mars, pred.var = "Year_Built", grid.resolution = 10) %>% 
  autoplot()

p3 <- partial(cv_mars, pred.var = c("Gr_Liv_Area", "Year_Built"), 
              grid.resolution = 10) %>% 
  plotPartial(levelplot = FALSE, zlab = "yhat", drape = TRUE, colorkey = TRUE, 
              screen = list(z = -20, x = -60))

# Display plots side by side
gridExtra::grid.arrange(p1, p2, p3, ncol = 3)

Attrition data

# get attrition data
df <- rsample::attrition %>% mutate_if(is.ordered, factor, ordered = FALSE)

# Create training (70%) and test (30%) sets for the rsample::attrition data.
# Use set.seed for reproducibility
set.seed(123)
churn_split <- rsample::initial_split(df, prop = .7, strata = "Attrition")
churn_train <- rsample::training(churn_split)
churn_test  <- rsample::testing(churn_split)


# for reproducibiity
set.seed(123)

# cross validated model
tuned_mars <- train(
  x = subset(churn_train, select = -Attrition),
  y = churn_train$Attrition,
  method = "earth",
  trControl = trainControl(method = "cv", number = 10),
  tuneGrid = hyper_grid
)

# best model
tuned_mars$bestTune

# plot results
ggplot(tuned_mars)

Table 7.2:

# train logistic regression model
set.seed(123)
glm_mod <- train(
  Attrition ~ ., 
  data = churn_train, 
  method = "glm",
  family = "binomial",
  preProc = c("zv", "center", "scale"),
  trControl = trainControl(method = "cv", number = 10)
  )

# train regularized logistic regression model
set.seed(123)
penalized_mod <- train(
  Attrition ~ ., 
  data = churn_train, 
  method = "glmnet",
  family = "binomial",
  preProc = c("zv", "center", "scale"),
  trControl = trainControl(method = "cv", number = 10),
  tuneLength = 10
  )

# extract out of sample performance measures
summary(resamples(list(
  Logistic_model = glm_mod, 
  Elastic_net = penalized_mod,
  MARS_model = tuned_mars
  )))$statistics$Accuracy %>%
  kableExtra::kable(caption = "Cross-validated accuracy results for tuned MARS and regression models.") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))

Cross-validated accuracy results for tuned MARS and regression models.

	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.	NA's
Logistic_model	0.8365385	0.8495146	0.8792476	0.8757893	0.8907767	0.9313725	0
Elastic_net	0.8446602	0.8759280	0.8834951	0.8835759	0.8915469	0.9411765	0
MARS_model	0.8155340	0.8578463	0.8780697	0.8708500	0.8907767	0.9029126	0