Applied Machine Learning for Marketing
By the end of this notebook you will be able to:
step_poly() to capture nonlinear relationships between predictors and the outcomeyardsticktune(), tune_grid(), and select_best() to choose the optimal penaltyTextbook references: TMWR Chapters 6 (Fitting Models with parsnip), 7 (A Model Workflow), 9 (Judging Model Effectiveness); ISLR Chapters 3 (Linear Regression), 6 (Linear Model Selection and Regularization)
#install.packages(
# "https://cran.r-project.org/src/contrib/Archive/vip/vip_0.4.6.tar.gz",
# repos = NULL,
# type = "source"
#)
library(tidymodels)
library(tidyverse)
library(glmnet) # Ridge, Lasso, Elastic Net engine
library(patchwork) # combining plots
library(vip) # variable importance plots
tidymodels_prefer()
set.seed(2025)We simulate 2,500 customers from a retail subscription service. The outcome variable is annual_spend — total customer spending in the past year — a continuous, right-skewed variable typical of real marketing data.
Predicting annual spend is one of the most common tasks in Customer Lifetime Value (CLV) modeling: knowing which customers will spend the most allows marketers to allocate acquisition budgets, design loyalty programs, and prioritize retention efforts.
set.seed(2025)
n <- 2500
clv_data <- tibble(
customer_id = paste0("CLV", str_pad(1:n, 5, pad = "0")),
age = round(pmin(pmax(rnorm(n, mean = 42, sd = 12), 18), 85)),
income = round(rlnorm(n, meanlog = 10.8, sdlog = 0.55)),
tenure_months = round(pmax(1, rnorm(n, mean = 30, sd = 18))),
num_products = sample(
1:6,
n,
replace = TRUE,
prob = c(0.25, 0.30, 0.20, 0.12, 0.08, 0.05)
),
freq_purchases = rpois(n, lambda = 8),
recency_days = round(pmax(0, rexp(n, rate = 1 / 30))), # mean = 1/rate
avg_basket = round(rlnorm(n, meanlog = 4.0, sdlog = 0.45), 2), # $
channel = sample(
c("online", "in_store", "both"),
n,
replace = TRUE,
prob = c(0.45, 0.25, 0.30)
),
loyalty_tier = sample(
c("Bronze", "Silver", "Gold", "Platinum"),
n,
replace = TRUE,
prob = c(0.35, 0.30, 0.22, 0.13)
),
region = sample(
c("West", "South", "Midwest", "Northeast"),
n,
replace = TRUE
),
promo_response = sample(
c("yes", "no"),
n,
replace = TRUE,
prob = c(0.40, 0.60)
)
) |>
mutate(
# annual_spend — continuous outcome
# driven by income, frequency, basket size, products owned, loyalty
annual_spend = round(
pmax(
10,
200 +
0.004 * income +
18 * freq_purchases +
0.85 * avg_basket +
45 * num_products -
0.8 * recency_days +
2.0 * tenure_months +
case_when(
loyalty_tier == "Platinum" ~ 350,
loyalty_tier == "Gold" ~ 200,
loyalty_tier == "Silver" ~ 80,
TRUE ~ 0
) +
if_else(channel == "both", 120, 0) +
if_else(promo_response == "yes", 75, 0) +
rnorm(n, 0, 180) # noise
),
2
),
# Introduce missingness
income = if_else(runif(n) < 0.07, NA_real_, income),
avg_basket = if_else(runif(n) < 0.05, NA_real_, avg_basket),
# Encode factors
channel = factor(channel),
loyalty_tier = factor(
loyalty_tier,
levels = c("Bronze", "Silver", "Gold", "Platinum")
),
region = factor(region),
promo_response = factor(promo_response)
)
glimpse(clv_data)Rows: 2,500
Columns: 13
$ customer_id <chr> "CLV00001", "CLV00002", "CLV00003", "CLV00004", "CLV000…
$ age <dbl> 49, 42, 51, 57, 46, 40, 47, 41, 38, 50, 37, 21, 37, 51,…
$ income <dbl> 47226, 76957, 28248, 75347, 97159, 68629, 55455, 65584,…
$ tenure_months <dbl> 70, 72, 22, 1, 2, 48, 34, 41, 27, 36, 33, 38, 34, 10, 2…
$ num_products <int> 1, 2, 5, 1, 4, 4, 1, 1, 2, 4, 5, 3, 2, 1, 2, 5, 1, 2, 4…
$ freq_purchases <int> 5, 16, 9, 7, 12, 8, 5, 10, 9, 5, 11, 10, 15, 7, 12, 6, …
$ recency_days <dbl> 51, 35, 7, 45, 3, 34, 18, 9, 25, 36, 12, 2, 34, 21, 11,…
$ avg_basket <dbl> 25.72, 33.65, 59.96, 42.05, 74.31, 126.22, 77.90, 23.18…
$ channel <fct> both, both, online, both, in_store, online, both, both,…
$ loyalty_tier <fct> Bronze, Gold, Silver, Silver, Gold, Silver, Platinum, B…
$ region <fct> Midwest, Midwest, Midwest, Northeast, Midwest, West, We…
$ promo_response <fct> no, yes, no, yes, no, yes, yes, no, no, no, yes, no, no…
$ annual_spend <dbl> 532.60, 1517.24, 662.31, 1006.08, 1129.05, 1106.61, 117…library(corrr)
clv_data |>
select(
age,
income,
tenure_months,
num_products,
freq_purchases,
recency_days,
avg_basket,
promo_response,
annual_spend
) |>
correlate() |>
shave() |>
fashion(2)# Distribution of annual_spend
p1 <- ggplot(clv_data, aes(x = annual_spend)) +
geom_histogram(bins = 40, fill = "#154734", color = "white", alpha = 0.85) +
labs(
title = "Distribution of Annual Spend",
x = "Annual spend ($)",
y = "Count"
) +
theme_minimal(base_size = 11)
# Log-transformed spend
p2 <- ggplot(clv_data, aes(x = log(annual_spend))) +
geom_histogram(bins = 40, fill = "#B08D57", color = "white", alpha = 0.85) +
labs(
title = "Log-Transformed Annual Spend",
x = "log(Annual spend)",
y = "Count"
) +
theme_minimal(base_size = 11)
p1 + p2
# Spend by loyalty tier
clv_data |>
group_by(loyalty_tier) |>
summarise(
n = n(),
mean_spend = round(mean(annual_spend)),
median_spend = round(median(annual_spend)),
sd_spend = round(sd(annual_spend))
)| loyalty_tier | n | mean_spend | median_spend | sd_spend |
|---|---|---|---|---|
| Bronze | 920 | 855 | 847 | 265 |
| Silver | 726 | 913 | 903 | 246 |
| Gold | 529 | 1035 | 1012 | 252 |
| Platinum | 325 | 1198 | 1180 | 279 |
# Spend vs. key continuous predictors
p3 <- ggplot(clv_data, aes(x = freq_purchases, y = annual_spend)) +
geom_point(alpha = 0.2, color = "#154734") +
geom_smooth(method = "lm", color = "#B08D57", se = FALSE) +
labs(
title = "Spend vs. Purchase Frequency",
x = "Number of purchases",
y = "Annual spend ($)"
) +
theme_minimal(base_size = 11)
p4 <- ggplot(clv_data, aes(x = avg_basket, y = annual_spend)) +
geom_point(alpha = 0.2, color = "#154734") +
geom_smooth(method = "lm", color = "#B08D57", se = FALSE) +
labs(
title = "Spend vs. Average Basket Size",
x = "Average basket ($)",
y = "Annual spend ($)"
) +
theme_minimal(base_size = 11)
p3 + p4
annual_spend is right-skewed — a small group of high spenders pulls the mean above the median. This is characteristic of customer spending distributions in retail: most customers spend modestly, while a small segment of “whales” drives disproportionate revenue.
Throughout this module we predict annual_spend directly (on the dollar scale). Note that in production, some practitioners log-transform the outcome before modeling and back-transform predictions — a design choice with implications for how we interpret RMSE.
Before adding regularization, we establish a baseline using ordinary least squares (OLS) linear regression. This gives us a benchmark to compare against Ridge, Lasso, and Elastic Net.
OLS minimizes the residual sum of squares (RSS):
\[RSS = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2\]
It fits coefficients that produce the smallest total squared prediction error on the training data — with no constraint on how large the coefficients can be.
lm_spec <- linear_reg() |>
set_engine("lm")
lm_wf <- workflow() |>
add_recipe(clv_rec) |>
add_model(lm_spec)
lm_fit <- fit(lm_wf, data = clv_train)
# Coefficient table
lm_fit |>
extract_fit_parsnip() |>
tidy() |>
filter(term != "(Intercept)") |>
arrange(desc(abs(estimate)))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| loyalty_tier_Platinum | 355.8851970 | 13.858359 | 25.6801839 | 0.0000000 |
| loyalty_tier_Gold | 191.0371125 | 11.717718 | 16.3032696 | 0.0000000 |
| channel_online | -131.7541609 | 10.097418 | -13.0483025 | 0.0000000 |
| income | 129.5213415 | 4.326477 | 29.9369063 | 0.0000000 |
| channel_in_store | -126.5855237 | 11.624506 | -10.8895399 | 0.0000000 |
| promo_response_yes | 82.3956197 | 8.707265 | 9.4628584 | 0.0000000 |
| num_products | 70.4320278 | 4.306286 | 16.3556305 | 0.0000000 |
| loyalty_tier_Silver | 68.7282590 | 10.699641 | 6.4234175 | 0.0000000 |
| freq_purchases | 57.2852024 | 4.315193 | 13.2752342 | 0.0000000 |
| tenure_months | 34.0349477 | 4.305643 | 7.9047307 | 0.0000000 |
| recency_days | -23.0654099 | 4.311924 | -5.3492158 | 0.0000001 |
| avg_basket | 13.8197731 | 4.307926 | 3.2079877 | 0.0013580 |
| region_South | -10.5914643 | 12.114347 | -0.8742910 | 0.3820657 |
| age | -2.8038609 | 4.306990 | -0.6510024 | 0.5151204 |
| region_West | 0.3545659 | 12.220944 | 0.0290130 | 0.9768572 |
| region_Northeast | 0.0277501 | 12.221718 | 0.0022706 | 0.9981886 |
lm_preds <- augment(lm_fit, new_data = clv_test)
lm_metrics <- lm_preds |>
metric_set(rmse, rsq, mae)(truth = annual_spend, estimate = .pred)
lm_metrics| .metric | .estimator | .estimate |
|---|---|---|
| rmse | standard | 194.6322969 |
| rsq | standard | 0.5135786 |
| mae | standard | 153.8205246 |
ggplot(lm_preds, aes(x = .pred, y = annual_spend - .pred)) +
geom_point(alpha = 0.3, color = "#154734") +
geom_hline(yintercept = 0, linetype = "dashed", color = "#B08D57") +
labs(
title = "Linear Regression — Residuals vs. Fitted",
subtitle = "Patterns here suggest the model is missing something",
x = "Predicted annual spend ($)",
y = "Residual ($)"
) +
theme_minimal(base_size = 12)
OLS produces unbiased coefficient estimates — but with many predictors, it can overfit the training data by assigning large coefficients to noise variables. The solution is regularization: adding a penalty to the loss function that discourages overly large coefficients. The three methods that follow (Ridge, Lasso, Elastic Net) all use the glmnet engine and differ only in how they penalize coefficients.
Ridge regression adds an L2 penalty to the OLS loss function:
\[\text{Ridge loss} = RSS + \lambda \sum_{j=1}^{p} \beta_j^2\]
where \(\lambda\) (the penalty parameter) controls the strength of regularization. As \(\lambda\) increases, all coefficients shrink toward — but never exactly reach — zero.
ridge_spec <- linear_reg(penalty = 0.1, mixture = 0) |>
set_engine("glmnet")
ridge_wf <- workflow() |>
add_recipe(clv_rec) |>
add_model(ridge_spec)
ridge_fit <- fit(ridge_wf, data = clv_train)
# Coefficients
ridge_fit |>
extract_fit_parsnip() |>
tidy() |>
filter(term != "(Intercept)") |>
arrange(desc(abs(estimate)))| term | estimate | penalty |
|---|---|---|
| loyalty_tier_Platinum | 333.8281897 | 0.1 |
| loyalty_tier_Gold | 174.9886601 | 0.1 |
| income | 123.6629653 | 0.1 |
| channel_online | -121.9869035 | 0.1 |
| channel_in_store | -116.5302104 | 0.1 |
| promo_response_yes | 79.0380181 | 0.1 |
| num_products | 67.3757718 | 0.1 |
| loyalty_tier_Silver | 56.7448909 | 0.1 |
| freq_purchases | 54.3468457 | 0.1 |
| tenure_months | 32.3753840 | 0.1 |
| recency_days | -21.6341098 | 0.1 |
| avg_basket | 13.3707191 | 0.1 |
| region_South | -10.7833167 | 0.1 |
| age | -2.9827022 | 0.1 |
| region_Northeast | 0.9131551 | 0.1 |
| region_West | 0.3367024 | 0.1 |
ridge_preds <- augment(ridge_fit, new_data = clv_test)
ridge_metrics <- ridge_preds |>
metric_set(rmse, rsq, mae)(truth = annual_spend, estimate = .pred)
ridge_metrics| .metric | .estimator | .estimate |
|---|---|---|
| rmse | standard | 194.8722078 |
| rsq | standard | 0.5113965 |
| mae | standard | 154.1712337 |
In linear_reg() with the glmnet engine, the mixture parameter controls the blend of L1 (Lasso) and L2 (Ridge) penalties:
mixture = 0 → pure Ridge (L2 only)mixture = 1 → pure Lasso (L1 only)0 < mixture < 1 → Elastic Net (blend of both)The penalty parameter is \(\lambda\) — the strength of regularization. We set it manually here (penalty = 0.1) as a starting point; in Section 7 we will search for the optimal value using cross-validation.
Lasso (Least Absolute Shrinkage and Selection Operator) adds an L1 penalty:
\[\text{Lasso loss} = RSS + \lambda \sum_{j=1}^{p} |\beta_j|\]
The key difference from Ridge: because the L1 penalty is proportional to the absolute value of the coefficient (not its square), Lasso can set coefficients to exactly zero — performing automatic variable selection.
lasso_spec <- linear_reg(penalty = 0.1, mixture = 1) |>
set_engine("glmnet")
lasso_wf <- workflow() |>
add_recipe(clv_rec) |>
add_model(lasso_spec)
lasso_fit <- fit(lasso_wf, data = clv_train)
# Coefficients — note which are exactly zero
lasso_fit |>
extract_fit_parsnip() |>
tidy() |>
filter(term != "(Intercept)") |>
arrange(desc(abs(estimate))) |>
mutate(selected = if_else(estimate != 0, "✓ kept", "✗ zeroed"))| term | estimate | penalty | selected |
|---|---|---|---|
| loyalty_tier_Platinum | 353.987476 | 0.1 | ✓ kept |
| loyalty_tier_Gold | 189.254114 | 0.1 | ✓ kept |
| channel_online | -130.470300 | 0.1 | ✓ kept |
| income | 129.157344 | 0.1 | ✓ kept |
| channel_in_store | -125.142151 | 0.1 | ✓ kept |
| promo_response_yes | 81.775982 | 0.1 | ✓ kept |
| num_products | 70.108144 | 0.1 | ✓ kept |
| loyalty_tier_Silver | 67.049408 | 0.1 | ✓ kept |
| freq_purchases | 56.880115 | 0.1 | ✓ kept |
| tenure_months | 33.708878 | 0.1 | ✓ kept |
| recency_days | -22.711488 | 0.1 | ✓ kept |
| avg_basket | 13.539281 | 0.1 | ✓ kept |
| region_South | -10.072894 | 0.1 | ✓ kept |
| age | -2.515086 | 0.1 | ✓ kept |
| region_Northeast | 0.000000 | 0.1 | ✗ zeroed |
| region_West | 0.000000 | 0.1 | ✗ zeroed |
lasso_preds <- augment(lasso_fit, new_data = clv_test)
lasso_metrics <- lasso_preds |>
metric_set(rmse, rsq, mae)(truth = annual_spend, estimate = .pred)
lasso_metrics| .metric | .estimator | .estimate |
|---|---|---|
| rmse | standard | 194.618812 |
| rsq | standard | 0.513421 |
| mae | standard | 153.808602 |
# Coefficient path — how coefficients change as penalty increases
lasso_fit |>
extract_fit_engine() |>
autoplot() +
labs(
title = "Lasso Coefficient Path",
subtitle = "As penalty (λ) increases, more coefficients are driven to zero",
x = "Log(λ)",
y = "Coefficient value"
) +
theme_minimal(base_size = 12)
Each line in the path plot traces one predictor’s coefficient as \(\lambda\) increases from left (small penalty, near-OLS) to right (large penalty, all coefficients near zero). Predictors whose lines hit zero earliest are the least important — Lasso is performing variable selection automatically as the penalty grows.
Predictors whose coefficients are driven to zero earlier tend to have weaker or more redundant contributions to the model, after accounting for the other predictors. However, this should not be interpreted as a perfect ranking of real-world importance, because correlated predictors, scaling, and sampling variation can affect which variables Lasso keeps or removes.
Elastic Net is a hybrid regularization method that combines the penalties of both Ridge (\(L_2\)) and Lasso (\(L_1\)) regression. It is particularly useful when dealing with highly correlated features (multicollinearity) or when the number of predictors (\(p\)) is much larger than the number of observations (\(n\)).
The objective of Elastic Net is to find the coefficient vector \(\beta\) that minimizes the following penalized loss function:
\[\min_{\beta} \left\{ \frac{1}{2n} \sum_{i=1}^{n} \left( y_i - \beta_0 - \sum_{j=1}^{p} x_{ij} \beta_j \right)^2 + \lambda \left[ \frac{1 - \alpha}{2} \sum_{j=1}^{p} \beta_j^2 + \alpha \sum_{j=1}^{p} |\beta_j| \right] \right\}\]
Elastic Net blends the L2 (Ridge) and L1 (Lasso) penalties:
\[\text{Elastic Net loss} = RSS + \lambda \left[ \frac{1-\alpha}{2} \sum_j \beta_j^2 + \alpha \sum_j |\beta_j| \right]\]
where \(\alpha\) is the mixture parameter in tidymodels. Setting mixture = 0.5 gives equal weight to Ridge and Lasso penalties.
penalty (\(\lambda\)) and mixture (\(\alpha\))enet_spec <- linear_reg(penalty = 0.1, mixture = 0.5) |>
set_engine("glmnet")
enet_wf <- workflow() |>
add_recipe(clv_rec) |>
add_model(enet_spec)
enet_fit <- fit(enet_wf, data = clv_train)
# Coefficients
enet_fit |>
extract_fit_parsnip() |>
tidy() |>
filter(term != "(Intercept)") |>
arrange(desc(abs(estimate))) |>
mutate(selected = if_else(estimate != 0, "✓ kept", "✗ zeroed"))| term | estimate | penalty | selected |
|---|---|---|---|
| loyalty_tier_Platinum | 353.79686 | 0.1 | ✓ kept |
| loyalty_tier_Gold | 189.17576 | 0.1 | ✓ kept |
| channel_online | -130.45930 | 0.1 | ✓ kept |
| income | 129.08666 | 0.1 | ✓ kept |
| channel_in_store | -125.14741 | 0.1 | ✓ kept |
| promo_response_yes | 81.80668 | 0.1 | ✓ kept |
| num_products | 70.09379 | 0.1 | ✓ kept |
| loyalty_tier_Silver | 67.04226 | 0.1 | ✓ kept |
| freq_purchases | 56.88143 | 0.1 | ✓ kept |
| tenure_months | 33.72669 | 0.1 | ✓ kept |
| recency_days | -22.73849 | 0.1 | ✓ kept |
| avg_basket | 13.57720 | 0.1 | ✓ kept |
| region_South | -10.19479 | 0.1 | ✓ kept |
| age | -2.56791 | 0.1 | ✓ kept |
| region_Northeast | 0.00000 | 0.1 | ✗ zeroed |
| region_West | 0.00000 | 0.1 | ✗ zeroed |
enet_preds <- augment(enet_fit, new_data = clv_test)
enet_metrics <- enet_preds |>
metric_set(rmse, rsq, mae)(truth = annual_spend, estimate = .pred)
enet_metrics| .metric | .estimator | .estimate |
|---|---|---|
| rmse | standard | 194.6192854 |
| rsq | standard | 0.5134054 |
| mae | standard | 153.8075839 |
Linear regression assumes the relationship between each predictor and the outcome is a straight line. Real marketing data rarely behaves that way — the effect of freq_purchases on annual_spend may accelerate at higher frequencies, or recency_days may have diminishing negative returns.
Polynomial regression extends linear regression by adding powered versions of a predictor:
\[\hat{y} = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \ldots\]
In tidymodels, step_poly() handles this transformation automatically inside the recipe.
Polynomial regression is still a linear model — linear in the coefficients \(\beta\), even though the predictor enters nonlinearly. This means we can still use linear_reg(), interpret coefficients, and apply regularization. The polynomial terms are just new features created by step_poly().
# Does the relationship look linear or curved?
p5 <- ggplot(clv_train, aes(x = freq_purchases, y = annual_spend)) +
geom_point(alpha = 0.2, color = "#154734") +
geom_smooth(
method = "lm",
color = "#B08D57",
se = FALSE,
linetype = "dashed",
linewidth = 0.9
) +
geom_smooth(
method = "loess",
color = "#4F8A8B",
se = FALSE,
linewidth = 0.9
) +
labs(
title = "Spend vs. Purchase Frequency",
subtitle = "Dashed = linear fit | Solid = LOESS (flexible)",
x = "Number of purchases",
y = "Annual spend ($)"
) +
theme_minimal(base_size = 11)
p6 <- ggplot(clv_train, aes(x = tenure_months, y = annual_spend)) +
geom_point(alpha = 0.2, color = "#154734") +
geom_smooth(
method = "lm",
color = "#B08D57",
se = FALSE,
linetype = "dashed",
linewidth = 0.9
) +
geom_smooth(
method = "loess",
color = "#4F8A8B",
se = FALSE,
linewidth = 0.9
) +
labs(
title = "Spend vs. Tenure",
subtitle = "Does the relationship flatten out over time?",
x = "Tenure (months)",
y = "Annual spend ($)"
) +
theme_minimal(base_size = 11)
p5 + p6
poly_rec <- recipe(annual_spend ~ ., data = clv_train) |>
update_role(customer_id, new_role = "ID") |>
step_impute_median(all_numeric_predictors()) |>
step_log(income, avg_basket, base = 10, offset = 1) |>
step_poly(freq_purchases, tenure_months, degree = 2) |>
step_normalize(all_numeric_predictors()) |>
step_dummy(all_nominal_predictors()) |>
step_zv(all_predictors())
poly_wf <- workflow() |>
add_recipe(poly_rec) |>
add_model(lasso_spec) # use Lasso to handle the extra polynomial terms
poly_fit <- fit(poly_wf, data = clv_train)
poly_fit |>
extract_fit_parsnip() |>
tidy() |>
filter(term != "(Intercept)", estimate != 0) |>
arrange(desc(abs(estimate)))| term | estimate | penalty |
|---|---|---|
| loyalty_tier_Platinum | 353.8744300 | 0.1 |
| loyalty_tier_Gold | 189.2216828 | 0.1 |
| channel_online | -130.3069998 | 0.1 |
| income | 129.0249911 | 0.1 |
| channel_in_store | -125.0171441 | 0.1 |
| promo_response_yes | 81.7204315 | 0.1 |
| num_products | 70.0676683 | 0.1 |
| loyalty_tier_Silver | 66.8774991 | 0.1 |
| freq_purchases_poly_1 | 56.8906581 | 0.1 |
| tenure_months_poly_1 | 33.5495223 | 0.1 |
| recency_days | -22.7454483 | 0.1 |
| avg_basket | 13.5596879 | 0.1 |
| region_South | -10.0160372 | 0.1 |
| freq_purchases_poly_2 | -2.6920383 | 0.1 |
| age | -2.4520127 | 0.1 |
| tenure_months_poly_2 | -0.4284127 | 0.1 |
poly_preds <- augment(poly_fit, new_data = clv_test)
poly_metrics <- poly_preds |>
metric_set(rmse, rsq, mae)(truth = annual_spend, estimate = .pred)
poly_metrics| .metric | .estimator | .estimate |
|---|---|---|
| rmse | standard | 194.5807686 |
| rsq | standard | 0.5135506 |
| mae | standard | 153.9285791 |
Adding polynomial terms increases the number of predictors. With degree = 2 on two variables, we add 2 new columns (freq_purchases^2, tenure_months^2). With degree = 3 on many variables, the explosion of terms can lead to overfitting. Combining step_poly() with Lasso regularization lets the model automatically zero out polynomial terms that don’t contribute — giving us flexibility without the overfitting risk.
So far we have set penalty = 0.1 for all regularized models — a reasonable starting value, but essentially a guess. The optimal penalty depends on the data, the number of predictors, and the signal-to-noise ratio. The solution is to search over a range of penalty values and use cross-validation to find the one that minimizes prediction error on held-out data.
This is the same fit_resamples() machinery you used in Module 3 — the difference is that now we run it across many penalty values simultaneously using tune_grid(). The penalty becomes a tune() placeholder rather than a fixed value, and tidymodels handles the rest.
# Use tune() as a placeholder for the penalty
lasso_tune_spec <- linear_reg(penalty = tune(), mixture = 1) |>
set_engine("glmnet")
lasso_tune_wf <- workflow() |>
add_recipe(clv_rec) |>
add_model(lasso_tune_spec)
# Cross-validation folds
set.seed(2025)
clv_folds <- vfold_cv(clv_train, v = 10)
# Penalty grid — 30 values on a log scale from 10^-4 to 10^1
penalty_grid <- grid_regular(
penalty(range = c(-4, 1)),
levels = 30
)
penalty_grid| penalty |
|---|
| 0.0001000 |
| 0.0001487 |
| 0.0002212 |
| 0.0003290 |
| 0.0004894 |
| 0.0007279 |
| 0.0010826 |
| 0.0016103 |
| 0.0023950 |
| 0.0035622 |
| 0.0052983 |
| 0.0078805 |
| 0.0117210 |
| 0.0174333 |
| 0.0259294 |
| 0.0385662 |
| 0.0573615 |
| 0.0853168 |
| 0.1268961 |
| 0.1887392 |
| 0.2807216 |
| 0.4175319 |
| 0.6210169 |
| 0.9236709 |
| 1.3738238 |
| 2.0433597 |
| 3.0391954 |
| 4.5203537 |
| 6.7233575 |
| 10.0000000 |
set.seed(2025)
lasso_tune_results <- tune_grid(
lasso_tune_wf,
resamples = clv_folds,
grid = penalty_grid,
metrics = metric_set(rmse, rsq, mae)
)
# Plot tuning results
autoplot(lasso_tune_results) +
labs(
title = "Lasso Tuning — RMSE and R² Across Penalty Values",
subtitle = "Each point = 10-fold CV mean; shaded band = ± 1 SE"
) +
theme_minimal(base_size = 12)
# Best penalty by lowest RMSE
best_penalty <- select_best(lasso_tune_results, metric = "rmse")
best_penalty| penalty | .config |
|---|---|
| 0.9236709 | pre0_mod24_post0 |
# One-standard-error rule — simpler model within 1 SE of the best
best_penalty_1se <- select_by_one_std_err(
lasso_tune_results,
metric = "rmse",
desc(penalty)
)
best_penalty_1se| penalty | .config |
|---|---|
| 6.723357 | pre0_mod29_post0 |
select_best() picks the penalty that minimizes mean RMSE — but this penalty may be on the edge of the performance plateau, where a slightly larger (simpler) penalty gives nearly identical RMSE but a more parsimonious model.
select_by_one_std_err() picks the largest penalty whose mean RMSE is within one standard error of the minimum. This gives a simpler, more regularized model with no meaningful sacrifice in performance — often preferred in practice.
# Plug the best penalty back into the workflow
final_lasso_wf <- finalize_workflow(lasso_tune_wf, best_penalty)
# Fit on the full training set with the tuned penalty
final_lasso_fit <- fit(final_lasso_wf, data = clv_train)
# Evaluate on the test set — this is the honest final estimate
final_preds <- augment(final_lasso_fit, new_data = clv_test)
final_metrics <- final_preds |>
metric_set(rmse, rsq, mae)(truth = annual_spend, estimate = .pred)
final_metrics| .metric | .estimator | .estimate |
|---|---|---|
| rmse | standard | 194.625871 |
| rsq | standard | 0.513041 |
| mae | standard | 153.857371 |
# Which predictors did the tuned Lasso keep?
final_lasso_fit |>
extract_fit_parsnip() |>
tidy() |>
filter(term != "(Intercept)") |>
arrange(desc(abs(estimate))) |>
mutate(selected = if_else(estimate != 0, "✓ kept", "✗ zeroed"))| term | estimate | penalty | selected |
|---|---|---|---|
| loyalty_tier_Platinum | 350.369019 | 0.9236709 | ✓ kept |
| loyalty_tier_Gold | 185.898280 | 0.9236709 | ✓ kept |
| income | 128.462717 | 0.9236709 | ✓ kept |
| channel_online | -128.039251 | 0.9236709 | ✓ kept |
| channel_in_store | -122.446318 | 0.9236709 | ✓ kept |
| promo_response_yes | 80.576446 | 0.9236709 | ✓ kept |
| num_products | 69.491345 | 0.9236709 | ✓ kept |
| loyalty_tier_Silver | 63.926006 | 0.9236709 | ✓ kept |
| freq_purchases | 56.110525 | 0.9236709 | ✓ kept |
| tenure_months | 33.081464 | 0.9236709 | ✓ kept |
| recency_days | -22.032671 | 0.9236709 | ✓ kept |
| avg_basket | 12.996948 | 0.9236709 | ✓ kept |
| region_South | -8.819330 | 0.9236709 | ✓ kept |
| age | -1.965201 | 0.9236709 | ✓ kept |
| region_Northeast | 0.000000 | 0.9236709 | ✗ zeroed |
| region_West | 0.000000 | 0.9236709 | ✗ zeroed |
Now that we have a tuned Lasso, let us compare all five models on equal footing using 10-fold cross-validation.
# Helper: fit a workflow and return CV metrics with a model label
cv_metrics <- function(wf, label) {
fit_resamples(
wf,
resamples = clv_folds,
metrics = metric_set(rmse, rsq, mae)
) |>
collect_metrics() |>
mutate(model = label)
}
lm_cv <- cv_metrics(lm_wf, "OLS Linear")
ridge_cv <- cv_metrics(ridge_wf, "Ridge (λ=0.1)")
lasso_cv <- cv_metrics(lasso_wf, "Lasso (λ=0.1)")
enet_cv <- cv_metrics(enet_wf, "Elastic Net (λ=0.1)")
poly_cv <- cv_metrics(poly_wf, "Polynomial + Lasso")
all_cv <- bind_rows(lm_cv, ridge_cv, lasso_cv, enet_cv, poly_cv)
# Summary table
all_cv |>
filter(.metric == "rmse") |>
select(model, mean, std_err) |>
arrange(mean)| model | mean | std_err |
|---|---|---|
| Lasso (λ=0.1) | 192.6698 | 3.359255 |
| Elastic Net (λ=0.1) | 192.6776 | 3.359883 |
| OLS Linear | 192.7114 | 3.356660 |
| Polynomial + Lasso | 192.9355 | 3.385298 |
| Ridge (λ=0.1) | 192.9722 | 3.385736 |
all_cv |>
filter(.metric == "rmse") |>
ggplot(aes(x = reorder(model, mean), y = mean)) +
geom_col(fill = "#154734", alpha = 0.85) +
geom_errorbar(
aes(ymin = mean - std_err, ymax = mean + std_err),
width = 0.3,
color = "#B08D57",
linewidth = 0.8
) +
coord_flip() +
labs(
title = "10-Fold CV RMSE — All Regression Models",
subtitle = "Error bars show ± 1 standard error | Lower RMSE = better",
x = NULL,
y = "Mean RMSE ($)"
) +
theme_minimal(base_size = 12)
bind_rows(
lm_metrics |> mutate(model = "OLS Linear"),
ridge_metrics |> mutate(model = "Ridge (λ=0.1)"),
lasso_metrics |> mutate(model = "Lasso (λ=0.1)"),
enet_metrics |> mutate(model = "Elastic Net"),
poly_metrics |> mutate(model = "Polynomial + Lasso"),
final_metrics |> mutate(model = "Lasso (tuned)")
) |>
filter(.metric %in% c("rmse", "rsq")) |>
select(model, .metric, .estimate) |>
pivot_wider(names_from = .metric, values_from = .estimate) |>
arrange(rmse)| model | rmse | rsq |
|---|---|---|
| Polynomial + Lasso | 194.5808 | 0.5135506 |
| Lasso (λ=0.1) | 194.6188 | 0.5134210 |
| Elastic Net | 194.6193 | 0.5134054 |
| Lasso (tuned) | 194.6259 | 0.5130410 |
| OLS Linear | 194.6323 | 0.5135786 |
| Ridge (λ=0.1) | 194.8722 | 0.5113965 |
A few things to watch for:
After fitting the best model, it is natural to ask: which predictors matter most? The vip package computes variable importance from the magnitude of standardized coefficients in regularized models.
final_lasso_fit |>
extract_fit_parsnip() |>
vip(
num_features = 12,
aesthetics = list(fill = "#154734", color = "white", alpha = 0.85)
) +
labs(
title = "Variable Importance — Tuned Lasso",
subtitle = "Importance based on absolute standardized coefficient magnitude",
x = "Importance"
) +
theme_minimal(base_size = 12)
Variable importance from a Lasso coefficient plot is a simple form of model explainability — you can see exactly which customer characteristics drive spend predictions. In Module 10 we will explore more sophisticated explainability tools (SHAP values, partial dependence plots) from TMWR Chapter 18, which work for any model type including the non-linear ones in Module 6.
| Method | Penalty | Zeroes Coefficients? | Best When |
|---|---|---|---|
| OLS | None | No | Predictors >> n is not a problem; interpretability priority |
| Ridge | L2: \(\lambda \sum \beta_j^2\) | No | Multicollinearity; keep all predictors |
| Lasso | L1: \(\lambda \sum |\beta_j|\) | Yes | Many predictors; want automatic selection |
| Elastic Net | L1 + L2 blend | Yes (partial) | Correlated predictors; want selection + stability |
| Polynomial | Any of the above | Depends on model | Nonlinear relationships visible in scatter plots |
# 1. Replace fixed penalty with tune()
spec <- linear_reg(penalty = tune(), mixture = 1) |> set_engine("glmnet")
# 2. Build a grid of candidate penalties
grid <- grid_regular(penalty(range = c(-4, 1)), levels = 30)
# 3. Search with cross-validation
results <- tune_grid(wf, resamples = folds, grid = grid,
metrics = metric_set(rmse, rsq))
# 4. Select best penalty
best <- select_best(results, metric = "rmse")
# 5. Finalize and fit on full training data
final_fit <- fit(finalize_workflow(wf, best), data = train)Fit Ridge (mixture = 0) and Lasso (mixture = 1) with an identical, large penalty of penalty = 10. Extract the coefficients from both models.
# Your code hereFit three Elastic Net models with penalty = 0.1 and mixture values of 0.25, 0.50, and 0.75. Compare the number of nonzero coefficients and test RMSE across all three.
mixture increases from 0.25 toward 0.75, what happens to the number of nonzero coefficients? Why?mixture value gives the best test RMSE? Is the difference meaningful given the standard errors you saw in the CV comparison?# Your code hereRefit the polynomial recipe with degree = 3 for both freq_purchases and tenure_months. Compare cross-validated RMSE for degree 1 (linear), degree 2, and degree 3.
mixture = 0 and penalty = 0.001).# Your code hereThe tuning example in Section 7 only tuned the penalty for Lasso (mixture = 1 fixed). Now tune both penalty and mixture for an Elastic Net model using a regular grid with 5 levels each.
penalty and mixture combination by RMSE? Is the optimal mixture closer to Ridge (0) or Lasso (1)?# Hint: use grid_regular(penalty(), mixture(), levels = 5)
# Your code hereUsing the variable importance plot and coefficient table from the tuned Lasso:
# Your code here