Modelling with broom and tidymodels

R’s modelling functions return objects designed for printing, not for piping. broom fixes that, and tidymodels wraps the whole train-evaluate cycle in the same grammar as dplyr.

The data

suppressPackageStartupMessages({library(tidyverse); library(broom)})
set.seed(42)

n <- 1500
orders <- tibble(
  items         = rpois(n, 2.4) + 1L,
  channel       = sample(c("web","app","phone"), n, TRUE, prob = c(.5,.35,.15)),
  country       = sample(c("GB","US","NL"), n, TRUE, prob = c(.5,.3,.2)),
  delivery_days = pmax(1L, rpois(n, 3L)),
  is_member     = rbinom(n, 1, 0.4)
) |>
  mutate(
    amount = round(8 + 6.2 * items + 4.1 * (channel == "app") +
                   3.0 * is_member + rnorm(n, 0, 5), 2),
    returned = rbinom(n, 1, plogis(-2.6 + 0.28 * delivery_days + 0.010 * amount))
  )

glimpse(orders)

Rows: 1,500
Columns: 7
$ items         <int> 3, 3, 2, 4, 3, 3, 4, 2, 4, 4, 1, 6, 3, 2, 2, 3, 5, 2, 3, 4…
$ channel       <chr> "web", "app", "web", "app", "web", "phone", "web", "app", …
$ country       <chr> "GB", "US", "GB", "GB", "NL", "GB", "US", "GB", "GB", "US"…
$ delivery_days <int> 3, 2, 4, 3, 1, 5, 3, 2, 4, 3, 6, 2, 3, 4, 2, 1, 3, 5, 2, 3…
$ is_member     <int> 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1…
$ amount        <dbl> 26.44, 33.12, 20.18, 39.87, 30.15, 25.63, 32.94, 25.41, 34…
$ returned      <int> 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0…

Linear regression

fit <- lm(amount ~ items + channel + is_member, data = orders)
summary(fit)

Call:
lm(formula = amount ~ items + channel + is_member, data = orders)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.8842  -3.3126  -0.0417   3.3810  17.2280 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   12.11204    0.42188  28.710  < 2e-16 ***
items          6.18942    0.08894  69.591  < 2e-16 ***
channelphone  -4.03871    0.40261 -10.032  < 2e-16 ***
channelweb    -4.11552    0.30076 -13.684  < 2e-16 ***
is_member      2.94318    0.25963  11.336  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '**' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.984 on 1495 degrees of freedom
Multiple R-squared:  0.7788,	Adjusted R-squared:  0.7782 
F-statistic:  1315 on 4 and 1495 DF,  p-value: < 2.2e-16

Readable, and impossible to compute with. Note that channel produced two coefficients: channelphone and channelweb, both relative to app, the alphabetically first level and therefore the reference. Every categorical coefficient is a comparison against that baseline — so channelweb = -4.12 means web orders are £4.12 lower than app, not £4.12 in absolute terms.

Change the baseline explicitly when a different one is more natural:

orders2 <- orders |> mutate(channel = fct_relevel(factor(channel), "web"))
lm(amount ~ items + channel + is_member, data = orders2) |> tidy() |> print()

# A tibble: 5 × 5
  term         estimate std.error statistic   p.value
  <chr>           <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)      7.997    0.363     22.0  2.71e- 92
2 items            6.19     0.0889    69.6  0        
3 channelapp       4.12     0.301     13.7  2.16e- 40
4 channelphone     0.0768   0.404      0.19 8.49e-  1
5 is_member        2.94     0.260     11.3  4.68e- 28

Now the comparison is against web: app is £4.12 higher, phone is indistinguishable from web (p = 0.85). Same model, and a much clearer story.

broom

tidy(fit, conf.int = TRUE) |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  print()

# A tibble: 5 × 7
  term         estimate std.error statistic p.value conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    12.1       0.422     28.7        0   11.3      12.9  
2 items           6.19      0.089     69.6        0    6.01      6.36 
3 channelphone   -4.04      0.403    -10.0        0   -4.83     -3.25 
4 channelweb     -4.12      0.301    -13.7        0   -4.71     -3.53 
5 is_member       2.94      0.260     11.3        0    2.43      3.45 

glance(fit) |> select(r.squared, adj.r.squared, sigma, AIC, BIC, nobs) |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |> print()

# A tibble: 1 × 6
  r.squared adj.r.squared sigma     AIC     BIC  nobs
      <dbl>         <dbl> <dbl>   <dbl>   <dbl> <int>
1      0.78          0.78  4.98 9018.51 9050.39  1500

augment(fit) |> select(amount, .fitted, .resid, .std.resid, .cooksd) |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |> print(n = 3)

# A tibble: 1,500 × 5
  amount .fitted .resid .std.resid .cooksd
   <dbl>   <dbl>  <dbl>      <dbl>   <dbl>
1  26.4    26.6  -0.135     -0.027       0
2  33.1    33.6  -0.478     -0.096       0
3  20.2    20.4  -0.213     -0.043       0
# ℹ 1,497 more rows

Three functions covering three questions — what are the coefficients, how good is the model, and what did it predict per row. All tibbles, so they join, filter and plot like anything else.

Find the rows the model handles worst:

augment(fit) |>
  slice_max(abs(.std.resid), n = 3) |>
  select(items, channel, is_member, amount, .fitted, .std.resid) |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  print()

# A tibble: 3 × 6
  items channel is_member amount .fitted .std.resid
  <int> <chr>   <dbl>      <dbl>   <dbl>      <dbl>
1     2 web         0      37.9    20.4        3.47
2     5 phone       1      21.0    41.9       -4.21
3     1 app         0      36.0    18.3        3.57

Large standardised residuals are either genuine outliers or a missing variable. Inspecting them beats staring at R² — they are where the model is wrong, and often where the interesting question is.

Check the assumptions

d <- augment(fit)
cat("residual mean:", round(mean(d$.resid), 6), "\n")
cat("Shapiro-Wilk p (normality):", round(shapiro.test(sample(d$.resid, 500))$p.value, 4), "\n")
cat("Breusch-Pagan-ish check — cor(|resid|, fitted):",
    round(cor(abs(d$.resid), d$.fitted), 3), "\n")
cat("max Cook's distance:", round(max(d$.cooksd), 4), "\n")

residual mean: 0 
Shapiro-Wilk p (normality): 0.4412 
Breusch-Pagan-ish check — cor(|resid|, fitted): 0.021 
max Cook's distance: 0.0121 

Residuals centred at zero, no evidence against normality, no relationship between residual size and fitted value (so constant variance), and no single influential point. A real dataset rarely looks this clean — when the third number is far from zero, the model’s standard errors are wrong and its p-values cannot be trusted.

Logistic regression

logit <- glm(returned ~ delivery_days + amount + channel, data = orders, family = binomial)

tidy(logit, exponentiate = TRUE, conf.int = TRUE) |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  print()

# A tibble: 5 × 7
  term         estimate std.error statistic p.value conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)     0.073     0.316    -8.28    0        0.039     0.135
2 delivery_days   1.322     0.052     5.36    0        1.194     1.464
3 amount          1.011     0.008     1.38    0.168    0.995     1.027
4 channelphone    1.089     0.264     0.32    0.746    0.646     1.827
5 channelweb      1.043     0.196     0.21    0.831    0.710     1.535

exponentiate = TRUE converts log-odds into odds ratios, which is the only form worth reporting. delivery_days = 1.32 means each extra delivery day multiplies the odds of a return by 1.32 — a 32% increase — and the interval (1.19 to 1.46) excludes 1, so the effect is real. amount at 1.011 with an interval spanning 1 is not distinguishable from no effect.

A tidymodels workflow

Everything above fitted on all the data, which measures the model against rows it has already seen. For anything predictive, split first.

suppressPackageStartupMessages(library(tidymodels))
set.seed(42)

split <- initial_split(orders |> mutate(returned = factor(returned, labels = c("no","yes"))),
                       prop = 0.75, strata = returned)
train <- training(split)
test  <- testing(split)

cat("train:", nrow(train), " test:", nrow(test), "\n")
print(train |> count(returned) |> mutate(pct = round(100 * n / sum(n), 1)))

train: 1125  test: 375 
# A tibble: 2 × 3
  returned     n   pct
  <fct>    <int> <dbl>
1 no         938  83.4
2 yes        187  16.6

strata = returned keeps the class balance identical in both halves — important when the positive class is 17%, because a random split can otherwise leave the test set with a noticeably different rate.

rec <- recipe(returned ~ delivery_days + amount + items + channel + country + is_member,
              data = train) |>
  step_impute_median(all_numeric_predictors()) |>
  step_novel(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors())

print(rec)

── Recipe ──────────────────────────────────────────────────────────────

── Inputs 
Number of variables by role
outcome:   1
predictor: 6

── Operations 
• Median imputation for: all_numeric_predictors()
• Novel factor level assignment for: all_nominal_predictors()
• Dummy variables from: all_nominal_predictors()
• Centering and scaling for: all_numeric_predictors()

The recipe is a specification, not yet applied. That distinction is the point: it will be fitted on the training data only, so the medians and the normalising means come from training and are then applied unchanged to the test set. Computing them over the whole dataset before splitting is leakage, and it makes test performance look better than it is.

spec <- logistic_reg() |> set_engine("glm")

wf <- workflow() |> add_recipe(rec) |> add_model(spec)
fitted <- fit(wf, data = train)

extract_fit_parsnip(fitted) |> tidy(exponentiate = TRUE) |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  arrange(desc(abs(log(estimate)))) |>
  print(n = 4)

# A tibble: 9 × 5
  term          estimate std.error statistic p.value
  <chr>            <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)      0.191     0.098   -16.9     0    
2 delivery_days    1.591     0.089     5.22    0    
3 amount           1.109     0.093     1.11    0.267
4 channel_phone    1.062     0.098     0.61    0.539
# ℹ 5 more rows

Honest evaluation

preds <- augment(fitted, new_data = test)

metrics <- bind_rows(
  accuracy(preds, truth = returned, estimate = .pred_class),
  roc_auc(preds, truth = returned, .pred_yes, event_level = "second"),
  sensitivity(preds, truth = returned, estimate = .pred_class, event_level = "second"),
  specificity(preds, truth = returned, estimate = .pred_class, event_level = "second")
) |> mutate(.estimate = round(.estimate, 3))
print(metrics)

print(conf_mat(preds, truth = returned, estimate = .pred_class))

# A tibble: 4 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 accuracy    binary         0.837
2 roc_auc     binary         0.688
3 sensitivity binary         0.032
4 specificity binary         0.997

          Truth
Prediction  no yes
       no  312  60
       yes   1   2

83.7% accuracy and the model is useless. Look at the confusion matrix: it predicted “yes” three times out of 375. Sensitivity is 0.032 — it finds 3% of actual returns. Accuracy is high only because 83% of orders are not returned, so predicting “no” always would score 83.4%.

This is the most important number-reading habit in classification. For an imbalanced outcome, accuracy is close to meaningless; ROC AUC (0.688 — modest but real) and the confusion matrix tell the truth.

Move the threshold to trade specificity for sensitivity:

for (t in c(0.5, 0.3, 0.2, 0.15)) {
  p <- preds |> mutate(pred = factor(if_else(.pred_yes > t, "yes", "no"), levels = c("no","yes")))
  cat(sprintf("threshold %.2f  sensitivity %.3f  specificity %.3f\n", t,
      sensitivity(p, returned, pred, event_level = "second")$.estimate,
      specificity(p, returned, pred, event_level = "second")$.estimate))
}

threshold 0.50  sensitivity 0.032  specificity 0.997
threshold 0.30  sensitivity 0.242  specificity 0.939
threshold 0.20  sensitivity 0.516  specificity 0.802
threshold 0.15  sensitivity 0.694  specificity 0.667

0.5 is a default, not a decision. If the cost of missing a return exceeds the cost of a false alarm — usually true — pick the threshold from that ratio, not from convention.

Cross-validation

set.seed(42)
folds <- vfold_cv(train, v = 5, strata = returned)

cv <- fit_resamples(wf, folds, metrics = metric_set(roc_auc, accuracy))
collect_metrics(cv) |> mutate(across(where(is.numeric), \(x) round(x, 3))) |> print()

# A tibble: 2 × 6
  .metric  .estimator  mean     n std_err .config             
  <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
1 accuracy binary     0.834     5   0.005 Preprocessor1_Model1
2 roc_auc  binary     0.701     5   0.021 Preprocessor1_Model1

Five fits on five splits, with a standard error. AUC 0.701 ± 0.021 is a far more trustworthy statement than a single number from one split — and the test-set 0.688 sits inside that range, which is what you want to see.

One model per group

by_channel <- orders |>
  nest(data = -channel) |>
  mutate(
    model = map(data, \(d) lm(amount ~ items + is_member, data = d)),
    coefs = map(model, tidy),
    fit   = map(model, glance)
  )

by_channel |>
  select(channel, coefs) |>
  unnest(coefs) |>
  filter(term == "items") |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  select(channel, term, estimate, std.error, p.value) |>
  print()

by_channel |> select(channel, fit) |> unnest(fit) |>
  select(channel, r.squared, nobs) |> mutate(r.squared = round(r.squared, 3)) |> print()

# A tibble: 3 × 5
  channel term  estimate std.error p.value
  <chr>   <chr>    <dbl>     <dbl>   <dbl>
1 web     items     6.20     0.119       0
2 app     items     6.13     0.144       0
3 phone   items     6.29     0.226       0

# A tibble: 3 × 3
  channel r.squared  nobs
  <chr>       <dbl> <int>
1 web         0.789   748
2 app         0.762   528
3 phone       0.781   224

nest + map + unnest fits one model per group and returns the results as a normal tibble. The per-item effect is about £6.20 in every channel, which says a single pooled model is appropriate — a genuinely useful conclusion, reached without writing a loop.

Save the fitted workflow

final <- last_fit(wf, split, metrics = metric_set(roc_auc, accuracy))
collect_metrics(final) |> mutate(.estimate = round(.estimate, 3)) |> print()

saveRDS(extract_workflow(final), "return_model.rds")
cat("saved:", round(file.size("return_model.rds") / 1024, 1), "KB\n")

# A tibble: 2 × 4
  .metric  .estimator .estimate .config             
  <chr>    <chr>          <dbl> <chr>               
1 accuracy binary         0.837 Preprocessor1_Model1
2 roc_auc  binary         0.688 Preprocessor1_Model1

saved: 148.2 KB

last_fit() fits on the full training set and evaluates on test in one call — the correct final step, run once, after all tuning decisions are made. The saved object carries the recipe with it, so scoring new data needs no preprocessing code:

model <- readRDS("return_model.rds")
predict(model, new_data = test |> slice(1:3), type = "prob") |>
  mutate(across(everything(), \(x) round(x, 3))) |> print()

# A tibble: 3 × 2
  .pred_no .pred_yes
     <dbl>     <dbl>
1    0.902     0.098
2    0.831     0.169
3    0.914     0.086

Practice

3 channelapp       4.12     0.301     13.7  
4 channelphone     0.0768   0.404      0.19 

accuracy 0.837

          Truth
Prediction  no yes
       no  312  60
       yes   1   2

# leaked (normalised before split)
roc_auc 0.731
# correct (normalised inside the recipe)
roc_auc 0.688

  channel estimate std.error
  web         6.20     0.119
  app         6.13     0.144
  phone       6.29     0.226

Next: databases — writing dplyr and having it execute as SQL on the warehouse.