Both come from one rule that I had never seen written down anywhere I was likely to read it. In R, a missing value is a join key. NA matches NA, so every row whose key is blank matches every other row whose key is blank. SQL does the opposite, since NULL = NULL is never true there, and I suspect that is where most people’s expectation comes from.

Below I take a real pair of tables, watch both failure modes happen, measure what they cost, and end with the guards that turn a join into something that fails loudly instead.

Two tables and an obvious key

I am using two indicators from Our World in Data: life expectancy at birth and median age. Both come down as plain CSV with no key or registration, and both are shaped the same way: one row per entity per year, with an entity name and a code.

library(readr)
library(dplyr)
library(ggplot2)

owid <- function(slug) {
  read_csv(paste0("https://ourworldindata.org/grapher/", slug,
                  ".csv?v=1&csvType=full&useColumnShortNames=true"),
           show_col_types = FALSE, progress = FALSE)
}

life <- owid("life-expectancy") |>
  rename(life_exp = life_expectancy_0) |>
  filter(year == 2023) |>
  select(entity, code, life_exp)

age <- owid("median-age") |>
  rename(median_age = median_age__sex_all__age_all__variant_estimates) |>
  filter(year == 2023, !is.na(median_age)) |>
  select(entity, code, median_age)

c(life = nrow(life), age = nrow(age))
## life  age 
##  261  253

The code column is ISO 3166 alpha-3 for countries, plus a handful of OWID_ codes for things the standard does not cover. It is the column you are supposed to join on, precisely because country names are a mess across providers. Here is the part that matters:

life |> filter(is.na(code)) |> pull(entity) |> head(6)
## [1] "Americas"                               
## [2] "High-and-upper-middle-income countries" 
## [3] "Land-locked Developing Countries (LLDC)"
## [4] "Latin America and the Caribbean"        
## [5] "Least developed countries"              
## [6] "Less developed regions"

c(blank_in_life = sum(is.na(life$code)), blank_in_age = sum(is.na(age$code)))
## blank_in_life  blank_in_age 
##            15             5

Neither table is broken. Both simply carry aggregate rows next to the country rows, and an aggregate like "Least developed countries" has no ISO code to put in the column, so the field is empty. Almost every statistical source I have pulled does this somewhere: a total row, a regional subtotal, a category that predates the code list.

Related posts on DataScience+:

• Converting data from long to wide simplified: tidyverse package
• Data Manipulation with dplyr
• Proteomics Data Analysis (2/3): Data Filtering and Missing Value Imputation

The join grows

say_warnings <- function(expr) {
  withCallingHandlers(expr, warning = function(w) {
    cat("Warning:", conditionMessage(w), "\n")
    invokeRestart("muffleWarning")
  })
}

joined <- say_warnings(
  left_join(life, age, by = "code", suffix = c("_life", "_age"))
)
## Warning: Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 6 of `x` matches multiple rows in `y`.
## ℹ Row 122 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

c(rows_in = nrow(life), rows_out = nrow(joined))
##  rows_in rows_out 
##      261      321

A left join is supposed to return one row per row of life. It returned 321 from 261. Since dplyr 1.1.0 there is at least a warning, and it is a good one, but it arrives phrased as a relationship problem rather than as a missing-value problem, so it is easy to read it as "these two tables just overlap in a complicated way" and move on.

The extra rows are not complicated at all:

joined |>
  filter(is.na(code)) |>
  select(entity_life, entity_age, life_exp, median_age) |>
  head(5)
## # A tibble: 5 × 4
##   entity_life entity_age                                     life_exp median_age
##   <chr>       <chr>                                             <dbl>      <dbl>
## 1 Americas    Least developed countries                          77.3       19.2
## 2 Americas    Less developed regions                             77.3       28.4
## 3 Americas    Less developed regions, excluding China            77.3       25.6
## 4 Americas    Less developed regions, excluding least devel…     77.3       30.6
## 5 Americas    More developed regions                             77.3       41.6

15 blank-code rows on the left met 5 blank-code rows on the right and produced every combination, 75 rows in which the life expectancy of one aggregate sits beside the median age of a different one. "Americas" paired with "Least developed countries" is not a data quality issue in either source. It is a row that R invented during the join.

Why NA matches NA

The behavior is a documented default rather than an accident. Every dplyr join takes an na_matches argument, and it is set to "na":

x <- tibble(code = c("AFG", NA, NA), v = 1:3)
y <- tibble(code = c("AFG", NA),     w = c(10, 20))

left_join(x, y, by = "code")                        # default: na_matches = "na"
## # A tibble: 3 × 3
##   code      v     w
##   <chr> <int> <dbl>
## 1 AFG       1    10
## 2 <NA>      2    20
## 3 <NA>      3    20

left_join(x, y, by = "code", na_matches = "never")
## # A tibble: 3 × 3
##   code      v     w
##   <chr> <int> <dbl>
## 1 AFG       1    10
## 2 <NA>      2    NA
## 3 <NA>      3    NA

This is not a tidyverse quirk either. Base R agrees, and has its own spelling of the fix:

merge(x, y, by = "code", all.x = TRUE)
##   code v  w
## 1  AFG 1 10
## 2 <NA> 2 20
## 3 <NA> 3 20

merge(x, y, by = "code", all.x = TRUE, incomparables = NA)
##   code v  w
## 1  AFG 1 10
## 2 <NA> 2 NA
## 3 <NA> 3 NA

Which default is right depends on what the NA means to you. If it is a real category, "the unclassified group", then matching it to itself is sensible. If it means "we do not know what this is", then two unknowns are not evidence of a match, and R’s default is doing something you almost never want. In my experience with data pulled from public sources, blank keys are the second kind essentially every time.

The version that leaves no trace

The warning above only fires because both sides had several blank-code rows. Watch what happens when the lookup table has exactly one, which is what you get from any source carrying a single "World" or "Total" row.

lookup <- age |> filter(!is.na(code) | entity == "More developed regions")
sum(is.na(lookup$code))
## [1] 1

quiet <- say_warnings(
  left_join(life, lookup, by = "code", suffix = c("_life", "_age"))
)

c(rows_in = nrow(life), rows_out = nrow(quiet))
##  rows_in rows_out 
##      261      261

No warning. No change in row count. The join is now many-to-one, which is exactly the relationship a left join is meant to have, so there is nothing for dplyr to object to. And yet:

quiet |>
  filter(is.na(code)) |>
  select(entity_life, entity_age, median_age) |>
  head(5)
## # A tibble: 5 × 3
##   entity_life                             entity_age             median_age
##   <chr>                                   <chr>                       <dbl>
## 1 Americas                                More developed regions       41.6
## 2 High-and-upper-middle-income countries  More developed regions       41.6
## 3 Land-locked Developing Countries (LLDC) More developed regions       41.6
## 4 Latin America and the Caribbean         More developed regions       41.6
## 5 Least developed countries               More developed regions       41.6

Every one of the 15 blank-code rows on the left has been handed a median age of 41.6 years, the value belonging to "More developed regions". "Least developed countries" now carries the median age of the developed world. Nothing in the object tells you this happened. The row count, which is the check most of us actually run, agrees with the input and always will.

Here is what those rows look like once they are in a plot, alongside the real matches:

dsp_colors <- c("#0066CC", "#E8862D", "#159A6C", "#7D5BD6",
                "#D64580", "#2AA9B8", "#C9A227")
dsp_theme <- theme_minimal(base_size = 13) +
  theme(plot.background    = element_rect(fill = "#ECECEF", color = NA),
        panel.background   = element_rect(fill = "#ECECEF", color = NA),
        panel.grid.minor   = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color = "grey78"),
        axis.ticks         = element_blank(),
        plot.title         = element_text(face = "bold"),
        strip.text         = element_text(face = "bold"))

quiet |>
  filter(!is.na(median_age)) |>
  mutate(kind = if_else(is.na(code), "Given a value that is not theirs", "Real matches")) |>
  ggplot(aes(median_age, life_exp, color = kind)) +
  geom_point(size = 1.9, alpha = 0.85) +
  scale_color_manual(values = c("Real matches" = dsp_colors[1],
                                "Given a value that is not theirs" = dsp_colors[2])) +
  labs(title = "No warning, no extra rows, fifteen wrong values",
       subtitle = "One blank-code row in the lookup table is enough",
       x = "Median age (years)", y = "Life expectancy at birth (years)",
       color = NULL) +
  dsp_theme + theme(legend.position = "top")

The wrong rows land in the middle of the cloud, at a median age that is perfectly plausible for a rich country. Without the color there is nothing to see. This is the case I now actively check for, because the loud version at least stops you.

What the extra rows cost

Back to the first join, the one that grew. If I go on to ask how life expectancy tracks median age across the world in 2023, the invented rows come along.

is_country <- function(code) !is.na(code) & (nchar(code) == 3 | code == "OWID_KOS")

as_returned <- joined |> filter(!is.na(median_age))
countries   <- as_returned |> filter(is_country(code))

fit_summary <- function(d, label) {
  m <- lm(life_exp ~ median_age, data = d)
  tibble(table = label,
         rows  = nrow(d),
         slope = round(coef(m)[["median_age"]], 3),
         r     = round(cor(d$life_exp, d$median_age), 3),
         r2    = round(summary(m)$r.squared, 3))
}

bind_rows(fit_summary(as_returned, "join as returned"),
          fit_summary(countries,   "countries only"))
## # A tibble: 2 × 5
##   table             rows slope     r    r2
##   <chr>            <int> <dbl> <dbl> <dbl>
## 1 join as returned   317 0.504 0.727 0.528
## 2 countries only     237 0.592 0.826 0.682

80 of the 317 rows I would have modeled, 25 percent of the sample, are aggregates or invented pairs. They pull the correlation from 0.826 down to 0.727 and flatten the slope from 0.592 to 0.504 years of life expectancy per year of median age. Nothing about those numbers looks wrong on its own. A correlation of 0.727 is the sort of figure you would write into a paragraph without a second thought.

as_returned |>
  mutate(kind = if_else(is.na(code), "Invented pairs (blank code)", "Real matches")) |>
  ggplot(aes(median_age, life_exp)) +
  geom_point(aes(color = kind), size = 1.9, alpha = 0.85) +
  geom_smooth(data = countries, method = "lm", formula = y ~ x, se = FALSE,
              color = dsp_colors[1], linewidth = 0.9) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              color = dsp_colors[2], linewidth = 0.9, linetype = "22") +
  scale_color_manual(values = c("Real matches" = dsp_colors[1],
                                "Invented pairs (blank code)" = dsp_colors[2])) +
  labs(title = "Seventy-five rows that match nothing real",
       subtitle = "Dashed line: the fit on the joined table exactly as it came back",
       x = "Median age (years)", y = "Life expectancy at birth (years)",
       color = NULL) +
  dsp_theme + theme(legend.position = "top")

The orange points form vertical stripes, one per distinct blank-code value on the right-hand side, each stripe stacking every blank-code value from the left. That signature is worth memorizing. Whenever a scatter plot of joined data shows a few suspiciously straight vertical or horizontal lines of points, a key matched more rows than it should have.

Writing the join as a contract

The fix is not to remember any of this. It is to say out loud what the join is supposed to do, and let R refuse when it cannot. na_matches has been there all along, and dplyr 1.1.0 added the rest.

try_join <- function(expr) {
  out <- tryCatch(expr, error = function(e) conditionMessage(e))
  if (is.character(out)) cat("Error:", out, "\n") else cat("Returned", nrow(out), "rows\n")
}

# What I actually mean: each country on the left gets at most one match,
# blanks match nothing, and I want to hear about it if that is not true.
try_join(
  left_join(life, age, by = "code",
            na_matches   = "never",
            relationship = "many-to-one")
)
## Returned 261 rows

# The same join without the na_matches guard still cannot hold.
try_join(
  left_join(life, age, by = "code", relationship = "many-to-one")
)
## Error: Each row in `x` must match at most 1 row in `y`.
## ℹ Row 6 of `x` matches multiple rows in `y`.

Three arguments I now write by default:

• na_matches = "never" makes blank keys match nothing, which is the SQL behavior and, I would argue, the useful one for keys pulled from files.
• relationship = states the cardinality you believe in ("one-to-one", "many-to-one", "one-to-many") and errors instead of silently multiplying rows. This is the argument that catches the quiet case above, because "many-to-one" is not violated by a bad NA match, but "one-to-one" is.
• unmatched = "error" turns unmatched keys into a failure rather than a column of NA. In a left join it checks y, since rows of x are kept by definition.

The last guard is not an argument at all. Before trusting a lookup table, look at what will not match:

unmatched <- anti_join(life, age, by = "code", na_matches = "never")
unmatched |> select(entity, code) |> print(n = Inf)
## # A tibble: 19 × 2
##    entity                                                      code    
##    <chr>                                                       <chr>   
##  1 Africa                                                      OWID_AFR
##  2 Americas                                                    <NA>    
##  3 Asia                                                        OWID_ASI
##  4 Europe                                                      OWID_EUR
##  5 High-and-upper-middle-income countries                      <NA>    
##  6 Land-locked Developing Countries (LLDC)                     <NA>    
##  7 Latin America and the Caribbean                             <NA>    
##  8 Least developed countries                                   <NA>    
##  9 Less developed regions                                      <NA>    
## 10 Less developed regions, excluding China                     <NA>    
## 11 Less developed regions, excluding least developed countries <NA>    
## 12 Low-and-Lower-middle-income countries                       <NA>    
## 13 Low-and-middle-income countries                             <NA>    
## 14 Middle-income countries                                     <NA>    
## 15 More developed regions                                      <NA>    
## 16 No income group available                                   <NA>    
## 17 Northern America                                            <NA>    
## 18 Oceania                                                     OWID_OCE
## 19 Small Island Developing States (SIDS)                       <NA>

That list is the join’s own account of itself, and it is the step I would keep if I had to drop the other three. Every row on it is an aggregate, which is the answer I want: the blanks that caused all the trouble, plus four continents that one file codes and the other does not carry. No actual country is being dropped, and that is a claim I can check rather than assume.

unmatched |> filter(is_country(code)) |> nrow()
## [1] 0

When that number is not zero, the useful question is never "how many rows did I lose" but "what do the lost rows have in common". Dropped rows are hardly ever a random sample. They are the small territories, the renamed states, the entities one provider counts and the other does not.

With the guards in place the analysis is the boring one I meant to run in the first place:

panel <- inner_join(
  life |> filter(is_country(code)),
  age  |> filter(is_country(code)),
  by = "code", na_matches = "never", relationship = "one-to-one",
  suffix = c("", "_age")
)

nrow(panel)
## [1] 237

round(cor(panel$life_exp, panel$median_age), 3)
## [1] 0.826

237 countries, correlation 0.826, and the join itself is now a statement that would have failed if either file had changed shape underneath me.

What I take from this

The row count is the check everyone runs, and it is the check that the dangerous version of this bug passes. A join that returns exactly as many rows as it started with can still have filled a column with values that belong to someone else, and it will do so without a warning, because many-to-one is a perfectly legitimate relationship for a left join to have.

So the habit I would rather build is the one where every join carries its assumptions in the call. na_matches = "never" unless a blank key genuinely names a category. relationship = on every join, because writing it forces you to decide what you believe before you find out. unmatched = "error" whenever the lookup is supposed to be complete. An anti_join() first, to read the list of rows that will not match. That is four extra lines that turn a silent wrong answer into a stack trace, which is the trade I will take every time.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

The first is the test itself. R gives you two implementations of the Breusch-Pagan test, they get used interchangeably, and only one of them is really a test for non-constant variance. The other also reacts to the shape of the error distribution: on data with heavy tails and perfectly constant variance, it rejects about half the time.

The second is what happens after. A significant test tells you the model has a problem somewhere, but it says nothing about which coefficient is damaged or by how much. On the data below, one coefficient’s standard error is off by 7% and another by 100%, from the same single p-value.

So this post is about closing that gap: what the test measures, what the fix actually changes, and which robust variance estimator to reach for. Everything runs on lmtest, sandwich and the tidyverse.

A model with a very obvious problem

I used the Ames housing data that ships with the modeldata package: every residential sale in Ames, Iowa between 2006 and 2010, with the lot and building characteristics recorded for each one. I kept the ordinary arm’s-length sales and regressed the sale price on four measurements a buyer can see.

library(tidyverse)
library(modeldata)
library(lmtest)
library(sandwich)
library(broom)

data(ames)

homes <- ames |>
  filter(Sale_Condition == "Normal") |>
  transmute(price      = Sale_Price / 1000,   # thousands of dollars
            area       = Gr_Liv_Area,         # above-grade living area, sq ft
            year_built = Year_Built,
            lot        = Lot_Area / 1000,     # thousands of sq ft
            basement   = Total_Bsmt_SF)

fit <- lm(price ~ area + year_built + lot + basement, data = homes)
nrow(homes)
## [1] 2413

That leaves 2,413 sales. The residual plot is the textbook picture, so much so that it barely needs a test:

dsp_colors <- c("#0066CC", "#E8862D", "#159A6C", "#7D5BD6",
                "#D64580", "#2AA9B8", "#C9A227")
dsp_theme <- theme_minimal(base_size = 13) +
  theme(plot.background    = element_rect(fill = "#ECECEF", color = NA),
        panel.background   = element_rect(fill = "#ECECEF", color = NA),
        panel.grid.minor   = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color = "grey78"),
        axis.ticks         = element_blank(),
        plot.title         = element_text(face = "bold"),
        strip.text         = element_text(face = "bold"))

augment(fit) |>
  ggplot(aes(.fitted, .resid)) +
  geom_hline(yintercept = 0, color = "grey40") +
  geom_point(alpha = 0.25, size = 1.1, color = dsp_colors[1]) +
  labs(title = "Residual spread grows with the predicted price",
       x = "Fitted sale price (thousands of dollars)", y = "Residual") +
  dsp_theme

(That theme snippet is reusable, so copy it once and drop it on every figure in a post.)

Cheap houses are predicted to within a few thousand dollars; expensive ones miss by a hundred thousand in either direction. The test agrees, emphatically:

bptest(fit)
## 
## 	studentized Breusch-Pagan test
## 
## data:  fit
## BP = 440.22, df = 4, p-value < 2.2e-16

Related posts on DataScience+:

• Fitting Polynomial Regression in R
• Multicollinearity in R
• How to detect heteroscedasticity and rectify it?

What the test is reacting to

Here is the part that is easy to miss. There are two Breusch-Pagan tests in common use in R, and they are not the same test.

The original 1979 statistic regresses the squared residuals on the predictors and scales the result by a variance figure that only holds if the errors are normal. bptest(fit, studentize = FALSE) computes exactly that, and car::ncvTest() computes the same normality-dependent score test (scored against the fitted values rather than the predictors, by default). The studentized version, from Koenker (1981), replaces the assumed scale with one estimated from the residuals, which makes it valid whatever the error distribution looks like. lmtest::bptest() uses the studentized version by default, which is why its printed heading says "studentized Breusch-Pagan test".

The difference matters more than the naming suggests. To see how much, I generated data where the variance really is constant, varied only the shape of the error distribution, and counted how often each version rejects at the 5% level. Every rejection in that setting is a false alarm. I then repeated it with variance that genuinely grows with x, where rejections are the desired outcome.

n <- 200
x <- rnorm(n)

errors <- function(shape) {
  e <- switch(shape,
              "normal"     = rnorm(n),
              "t(5)"       = rt(n, df = 5),
              "log-normal" = exp(rnorm(n)),
              "uniform"    = runif(n))
  (e - mean(e)) / sd(e)          # same variance, different shape
}

one_run <- function(shape, variance) {
  spread <- if (variance == "constant variance") 1 else exp(0.15 * x)
  y <- 1 + x + errors(shape) * spread
  m <- lm(y ~ x)
  tibble(shape = shape, variance = variance,
         studentized = bptest(m)$p.value,
         original    = bptest(m, studentize = FALSE)$p.value)
}

grid <- expand_grid(shape    = c("normal", "t(5)", "log-normal", "uniform"),
                    variance = c("constant variance", "variance grows with x"),
                    rep      = 1:2000)

rates <- pmap_dfr(list(grid$shape, grid$variance), one_run) |>
  pivot_longer(c(studentized, original), names_to = "test", values_to = "p") |>
  group_by(variance, shape, test) |>
  summarise(reject = mean(p < 0.05), .groups = "drop")

rates |> pivot_wider(names_from = test, values_from = reject)
## # A tibble: 8 × 4
##   variance              shape      original studentized
##   <chr>                 <chr>         <dbl>       <dbl>
## 1 constant variance     log-normal   0.525       0.05  
## 2 constant variance     normal       0.046       0.0495
## 3 constant variance     t(5)         0.187       0.046 
## 4 constant variance     uniform      0.0025      0.0495
## 5 variance grows with x log-normal   0.630       0.124 
## 6 variance grows with x normal       0.800       0.795 
## 7 variance grows with x t(5)         0.709       0.482 
## 8 variance grows with x uniform      0.91        0.990

rates |>
  mutate(shape = factor(shape, c("uniform", "normal", "t(5)", "log-normal"))) |>
  ggplot(aes(shape, reject, fill = test)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "grey30") +
  facet_wrap(~ variance) +
  scale_fill_manual(values = dsp_colors[c(2, 1)], name = NULL,
                    labels = c("original (assumes normality)",
                               "studentized (bptest default)")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "In the left panel, every rejection is a false alarm",
       subtitle = "Rejection rate at the 5% level, 2,000 simulations per bar",
       x = "Error distribution", y = "Rejections at the 5% level") +
  dsp_theme +
  theme(legend.position = "top")

Read the left panel first. With normally distributed errors both versions behave, sitting on the dashed 5% line. Swap in log-normal errors and the original test rejects 52% of the time, even though the variance never changes. Student’s t with 5 degrees of freedom, a mild amount of extra tail weight, gets it to 19%. The studentized version stays near 5% throughout.

The failure runs in the other direction too. With light-tailed uniform errors the original test rejects almost never under the null, 0.2%, which is not caution but a broken calibration that costs it detections elsewhere.

The reason is that the original statistic divides by a variance figure that is only correct when the errors are normal. Heavy tails inflate the squared residuals it is testing without inflating that denominator, so kurtosis reads as heteroscedasticity. Studentizing estimates the denominator from the residuals instead, and the sensitivity to shape disappears.

Now the right panel, where the variance really does grow with x. With normal errors the two versions have essentially identical power (80% against 80%), so the protection in the left panel is free. But look at the log-normal bars: the studentized test finds real heteroscedasticity only 12% of the time. Its higher-looking neighbor is not better, because a test that fires 52% of the time under the null cannot have its rejections interpreted.

That is the practical takeaway from this section, and it has two halves. Use the studentized version, which is what you get by default from bptest(). And do not read a non-significant result as a clean bill of health, because with skewed or heavy-tailed errors the honest test has very little power.

What robust standard errors actually change

Back to the houses, where the test result was not in doubt anyway. The usual next move is to transform the outcome, but that changes the quantity being estimated: a coefficient in a log-price model is a percentage effect, not a dollar effect, and if the dollar effect is what you wanted, you have answered a different question to fix a standard error.

The alternative keeps the model and fixes only the standard errors. sandwich::vcovHC() builds a heteroscedasticity-consistent covariance matrix, and coeftest() re-runs the coefficient table with it.

comparison <- bind_rows(
  classical = tidy(fit),
  HC3       = tidy(coeftest(fit, vcov. = vcovHC(fit))),
  .id = "vcov") |>
  filter(term != "(Intercept)") |>
  select(vcov, term, estimate, std.error, statistic) |>
  pivot_wider(names_from = vcov, values_from = c(std.error, statistic)) |>
  mutate(inflation = std.error_HC3 / std.error_classical)

print(comparison, width = Inf)
## # A tibble: 4 × 7
##   term       estimate std.error_classical std.error_HC3 statistic_classical
##   <chr>         <dbl>               <dbl>         <dbl>               <dbl>
## 1 area         0.0786             0.00159       0.00254               49.4 
## 2 year_built   0.731              0.0257        0.0275                28.5 
## 3 lot          0.685              0.0881        0.176                  7.78
## 4 basement     0.0521             0.00198       0.00303               26.4 
##   statistic_HC3 inflation
##           <dbl>     <dbl>
## 1         30.9       1.60
## 2         26.5       1.07
## 3          3.90      2.00
## 4         17.2       1.53

The point estimates never move, because ordinary least squares is still unbiased under heteroscedasticity. Only the uncertainty around them changes, and this is where a single global p-value stops being useful. The inflation factor runs from 1.07 on year_built to 2.00 on lot. The year the house was built has a standard error that was essentially fine all along. The lot size coefficient had a standard error half the size it should be, and its t statistic falls from 7.8 to 3.9.

Why that coefficient? Because heteroscedasticity only distorts a standard error where the large errors coincide with the extreme values of that predictor, and lot size is the variable with the extremes:

c(median_lot = median(homes$lot), largest_lot = max(homes$lot))
##  median_lot largest_lot 
##       9.360     215.245

# how much of the lot-size variation comes from ten properties
big <- slice_max(homes, lot, n = 10)
sum((big$lot - mean(homes$lot))^2) / sum((homes$lot - mean(homes$lot))^2)
## [1] 0.6923105

c(mean_leverage = mean(hatvalues(fit)), max_leverage = max(hatvalues(fit)))
## mean_leverage  max_leverage 
##   0.002072109   0.270229368

A handful of acreages, one of them 23 times the median lot, carry most of the information about that slope, and those same properties have the largest residuals. That is the worst possible combination, and it is invisible in the omnibus test.

One more thing worth knowing, because "robust" is often heard as "conservative": the correction is not guaranteed to make standard errors bigger. When the largest errors sit in the middle of a predictor’s range rather than at its extremes, it goes the other way.

ratio <- replicate(1000, {
  x <- rnorm(300)
  y <- 1 + x + rnorm(300, sd = exp(-0.6 * abs(x)))  # noisiest near the center
  m <- lm(y ~ x)
  sqrt(diag(vcovHC(m)))[["x"]] / sqrt(diag(vcov(m)))[["x"]]
})
c(mean_ratio = mean(ratio), share_below_one = mean(ratio < 1))
##      mean_ratio share_below_one 
##       0.6535694       1.0000000

Every simulated dataset gets a smaller robust standard error, averaging 0.65 times the classical one. The sandwich estimator is not a safety margin bolted on top of the classical one. It is a different estimate, and it can point either way.

Which sandwich

vcovHC() offers several corrections, and the choice matters far more than it appears to. HC0 is the original White estimator. HC1 applies a degrees-of-freedom factor and is what Stata’s robust option reports. HC3 divides each squared residual by (1 - h)², where h is that observation’s leverage, which approximates what you would get by leaving the point out and refitting. To see what that buys, I simulated data with real heteroscedasticity and measured how often a nominal 95% interval for the slope actually contains the true value.

covers <- function(m, V, truth = 1) {
  ci <- coefci(m, vcov. = V)["x", ]
  ci[1] <= truth & truth <= ci[2]
}

one_fit <- function(n) {
  x <- rnorm(n)
  y <- 1 + x + rnorm(n, sd = exp(0.6 * x))
  m <- lm(y ~ x)
  tibble(n = n,
         classical = covers(m, vcov(m)),
         HC0 = covers(m, vcovHC(m, type = "HC0")),
         HC1 = covers(m, vcovHC(m, type = "HC1")),
         HC3 = covers(m, vcovHC(m, type = "HC3")))
}

sizes <- c(30, 60, 120, 500)

cover <- map_dfr(rep(sizes, each = 2000), one_fit) |>
  pivot_longer(-n, names_to = "vcov", values_to = "covered") |>
  group_by(n, vcov) |>
  summarise(coverage = mean(covered), .groups = "drop")

cover |> pivot_wider(names_from = vcov, values_from = coverage)
## # A tibble: 4 × 5
##       n   HC0   HC1   HC3 classical
##   <dbl> <dbl> <dbl> <dbl>     <dbl>
## 1    30 0.898 0.912 0.953     0.828
## 2    60 0.918 0.922 0.940     0.820
## 3   120 0.934 0.936 0.942     0.815
## 4   500 0.946 0.946 0.951     0.79

cover |>
  ggplot(aes(factor(n), coverage, color = vcov, group = vcov)) +
  geom_hline(yintercept = 0.95, linetype = "dashed", color = "grey30") +
  geom_line(linewidth = 0.9) +
  geom_point(size = 2.4) +
  scale_y_continuous(labels = scales::percent, limits = c(0.75, 1)) +
  scale_color_manual(values = dsp_colors[c(5, 2, 7, 3)], name = NULL) +
  labs(title = "Coverage of a nominal 95% interval for the slope",
       subtitle = "2,000 simulations per point, variance growing with x",
       x = "Sample size", y = "Coverage") +
  dsp_theme +
  theme(legend.position = "top")

The classical interval is the flat line at the bottom, covering around 81% instead of 95%, and it gets worse as the sample grows. That is the thing to internalize about heteroscedasticity: it is not a small-sample problem that more data will wash out. The interval is converging, just to the wrong width.

Among the robust options, at n = 500 they are indistinguishable, which is the textbook claim and the reason the choice often gets waved through. At n = 30 they are not: HC0 covers 89.8% and HC1 91.2%, while HC3 reaches 95.3%. A test you believe runs at 5% is running closer to 10%. HC3 is the default in vcovHC() for exactly this reason, so plain vcovHC(fit) is already the right call, and the way to get this wrong is to type a type = argument copied from somewhere else.

The routine I would use instead

Plot the residuals, always. The picture tells you the shape of the problem, which no p-value does.

Run bptest() and leave studentize alone, and treat a non-significant result as weak evidence rather than as clearance, especially when the residuals are skewed.

Skip the test as a gate on whether to use robust standard errors. Deciding between two covariance estimators by looking at a p-value from the same data makes the reported standard error a random choice between them, and the heteroscedasticity-consistent one is valid either way, costing only a little efficiency when the variance really is constant. Report coeftest(fit, vcov. = vcovHC(fit)) and be done.

Transform the outcome when the transformed scale is the one you want to interpret, not as a repair for a standard error.

And when a coefficient’s standard error moves a lot under the correction, go look at that predictor. On the Ames data the correction was not just a technical adjustment: it pointed at ten properties that were quietly running the lot-size coefficient.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

A prototype usually focuses on one question:

Can the model produce something useful?

A production system needs to answer a harder question:

Can we trust this output enough to use it in a real workflow?

This is where human-in-the-loop validation becomes important.

Human-in-the-loop does not mean the AI system is weak. In many real-world systems, it is the difference between a risky automation and a controlled AI workflow. Human review helps teams validate uncertain outputs, correct mistakes, capture feedback, improve data quality, and create a reliable audit trail.

In this article, we will build a simple human-in-the-loop validation workflow in Python.

The goal is to create a small but practical process that:

• Loads model-generated predictions
• Assigns items for review
• Applies rule-based validation checks
• Captures human approval or correction
• Stores review decisions
• Produces a clean final dataset
• Creates a basic review summary

This pattern can be extended into a web app, Streamlit dashboard, internal review tool, or production data pipeline.

Why human-in-the-loop validation matters

Many AI systems produce outputs that look reasonable most of the time.

For example, an AI system might:

• Classify support tickets
• Label mobile apps by category
• Detect risky documents
• Score leads
• Summarise customer conversations
• Extract fields from invoices
• Recommend product tags
• Identify policy violations

The challenge is that "mostly correct" is not always good enough.

A wrong classification may affect reporting.
A wrong risk score may affect compliance.
A wrong label may affect user experience.
A wrong extracted field may affect business operations.

In production, we need to know when an AI output should be accepted automatically, when it should be reviewed, and how corrections should be stored.

A good human-in-the-loop system should answer:

• Which predictions need review?
• Who reviewed them?
• What decision did they make?
• What was corrected?
• Why was it corrected?
• Was the final value written back to the dataset?
• Can we analyse review patterns later?

If these answers are stored properly, human review becomes more than manual checking. It becomes a feedback loop.

Related posts on DataScience+:

• Machine Learning for Diabetes with Python
• Topic Modeling in Python with NLTK and Gensim
• Top Python Libraries for Machine Learning

Sample dataset

For this tutorial, we will use a small CSV file called ai_predictions.csv.

Each row contains an item that has been classified by an AI model.

item_id,item_name,predicted_category,predicted_audience,confidence,source
1,Math Games for Kids,Education,Children,0.91,android
2,Fast Racing Challenge,Racing,Teens,0.74,android
3,Invoice Scanner Pro,Productivity,Adults,0.86,ios
4,Monster Battle Arena,Action,Children,0.58,android
5,ABC Learning App,Education,Children,0.96,ios
6,Crypto Trading Signals,Finance,Teens,0.62,android
7,Daily Meditation Guide,Health,Adults,0.82,ios
8,Princess Coloring Book,Education,Children,0.88,android
9,Real Car Parking 3D,Racing,Children,0.67,android
10,Legal Contract Reader,Productivity,Adults,0.79,ios

In a real system, this data could come from:

* A machine learning model
* An LLM classification workflow
* An internal AI service
* A batch prediction pipeline
* A database table
* An API response
* A spreadsheet export

For this example, we will keep the dataset simple.

## Step 1: Load the predictions

First, load the CSV file with pandas.

import pandas as pd

df = pd.read_csv("ai_predictions.csv")

print(df.head())

Now check the structure of the data.

print(df.info())
print(df["predicted_category"].value_counts())

At this stage, we have raw model predictions. These are not yet approved. They are simply outputs that need either automatic acceptance or human review.

Step 2: Define review rules

Not every prediction needs human review.

A simple production workflow might auto-approve high-confidence predictions and send lower-confidence or higher-risk predictions to a review queue.

For this tutorial, we will review an item if:

• confidence is below 0.80
• the predicted audience is Children and confidence is below 0.90
• the category is Finance
• the category is Action and audience is Children

These rules are only examples. In a real system, review rules should be based on business risk, compliance needs, and model performance.

def needs_review(row):
    if row["confidence"] < 0.80:
        return True

    if row["predicted_audience"] == "Children" and row["confidence"] < 0.90:
        return True

    if row["predicted_category"] == "Finance":
        return True

    if row["predicted_category"] == "Action" and row["predicted_audience"] == "Children":
        return True

    return False

df["needs_review"] = df.apply(needs_review, axis=1)

print(df[["item_id", "item_name", "predicted_category", "predicted_audience", "confidence", "needs_review"]])

This creates a clear separation between items that can be auto-approved and items that require review.

Step 3: Create a review queue

Now we can create a review queue containing only items that require human validation.

review_queue = df[df["needs_review"] == True].copy()

review_queue = review_queue.sort_values(
    by=["confidence"],
    ascending=True
)

print(review_queue[[
    "item_id",
    "item_name",
    "predicted_category",
    "predicted_audience",
    "confidence"
]])

Sorting by confidence helps reviewers focus on the most uncertain items first.

In a larger system, the review queue could also be prioritised by:

• business value
• customer impact
• compliance risk
• model confidence
• number of previous errors
• source system
• deadline or SLA

Step 4: Simulate reviewer decisions

In a real application, reviewer decisions would be captured through a form, dashboard, or web application.

For this tutorial, we will simulate reviewer input using a small dictionary.

Each reviewer decision includes:

• item_id
• review_status
• final_category
• final_audience
• reviewer_notes

review_decisions = [
    {
        "item_id": 2,
        "review_status": "approved",
        "final_category": "Racing",
        "final_audience": "Teens",
        "reviewer_notes": "Prediction looks correct."
    },
    {
        "item_id": 4,
        "review_status": "corrected",
        "final_category": "Action",
        "final_audience": "Teens",
        "reviewer_notes": "Game is not suitable for children based on content."
    },
    {
        "item_id": 6,
        "review_status": "corrected",
        "final_category": "Finance",
        "final_audience": "Adults",
        "reviewer_notes": "Finance-related app should not be labelled for teens."
    },
    {
        "item_id": 9,
        "review_status": "approved",
        "final_category": "Racing",
        "final_audience": "Children",
        "reviewer_notes": "Child-friendly racing game."
    },
    {
        "item_id": 10,
        "review_status": "approved",
        "final_category": "Productivity",
        "final_audience": "Adults",
        "reviewer_notes": "Prediction looks correct."
    }
]

reviews = pd.DataFrame(review_decisions)

print(reviews)

This review table is important because it becomes part of the audit trail.

It tells us not only what the final answer was, but also what the human reviewer changed.

Step 5: Merge predictions with review decisions

Now merge the original predictions with the review decisions.

validated = df.merge(
    reviews,
    on="item_id",
    how="left"
)

print(validated.head())

For items that were not reviewed, we can use the original model prediction as the final approved value.

validated["final_category"] = validated["final_category"].fillna(
    validated["predicted_category"]
)

validated["final_audience"] = validated["final_audience"].fillna(
    validated["predicted_audience"]
)

validated["review_status"] = validated["review_status"].fillna(
    "auto_approved"
)

validated["reviewer_notes"] = validated["reviewer_notes"].fillna(
    "Auto-approved based on validation rules."
)

print(validated[[
    "item_id",
    "item_name",
    "predicted_category",
    "final_category",
    "predicted_audience",
    "final_audience",
    "review_status"
]])

Now every row has a final value.

Some were auto-approved.
Some were approved by a human.
Some were corrected by a human.

This is much better than keeping raw model outputs with no validation context.

Step 6: Identify changed predictions

It is useful to know where the human reviewer disagreed with the model.

validated["category_changed"] = (
    validated["predicted_category"] != validated["final_category"]
)

validated["audience_changed"] = (
    validated["predicted_audience"] != validated["final_audience"]
)

validated["any_change"] = (
    validated["category_changed"] | validated["audience_changed"]
)

print(validated[[
    "item_id",
    "item_name",
    "predicted_category",
    "final_category",
    "predicted_audience",
    "final_audience",
    "any_change"
]])

This gives us a simple correction dataset.

In a real system, this can help answer:

• Which categories are most often corrected?
• Which audience labels are most unreliable?
• Which sources produce more review failures?
• Which rules are sending too many or too few items for review?
• Where should the model or prompt be improved?

Step 7: Create a review summary

Now we can summarise the review process.

summary = {
    "total_items": len(validated),
    "auto_approved": (validated["review_status"] == "auto_approved").sum(),
    "human_reviewed": (validated["review_status"] != "auto_approved").sum(),
    "approved_by_reviewer": (validated["review_status"] == "approved").sum(),
    "corrected_by_reviewer": (validated["review_status"] == "corrected").sum(),
    "changed_predictions": validated["any_change"].sum()
}

for key, value in summary.items():
    print(f"{key}: {value}")

We can also calculate review rates.

total_items = len(validated)
human_reviewed = summary["human_reviewed"]
changed_predictions = summary["changed_predictions"]

review_rate = human_reviewed / total_items
correction_rate = changed_predictions / total_items

print(f"Review rate: {review_rate:.1%}")
print(f"Correction rate: {correction_rate:.1%}")

These metrics are useful because they show how much human effort is required and how often the model is corrected.

A high correction rate may indicate:

• the model is weak for certain categories
• the prompt needs improvement
• the training data is outdated
• the review rules are catching genuinely risky cases
• the category definitions are unclear

Step 8: Analyse corrections by category

Now let us see which predicted categories had the most corrections.

corrections_by_category = (
    validated[validated["any_change"] == True]
    .groupby("predicted_category")
    .agg(
        corrections=("item_id", "count")
    )
    .reset_index()
    .sort_values("corrections", ascending=False)
)

print(corrections_by_category)

If a category is corrected often, it may need closer attention.

We can also look at corrections by source.

corrections_by_source = (
    validated.groupby("source")
    .agg(
        total_items=("item_id", "count"),
        corrections=("any_change", "sum")
    )
    .reset_index()
)

corrections_by_source["correction_rate"] = (
    corrections_by_source["corrections"] / corrections_by_source["total_items"]
)

print(corrections_by_source)

This can help identify whether errors are linked to a specific source system, platform, or data pipeline.

Step 9: Export the final dataset

Once review is complete, export the validated dataset.

final_output = validated[[
    "item_id",
    "item_name",
    "final_category",
    "final_audience",
    "review_status",
    "reviewer_notes"
]]

final_output.to_csv("validated_ai_outputs.csv", index=False)

print(final_output)

This final output is the dataset that downstream systems should use.

The raw model output is still preserved, but the business process uses the validated output.

Step 10: Save the review audit trail

It is also important to save the review history.

audit_trail = validated[[
    "item_id",
    "item_name",
    "predicted_category",
    "predicted_audience",
    "confidence",
    "needs_review",
    "review_status",
    "final_category",
    "final_audience",
    "reviewer_notes",
    "category_changed",
    "audience_changed",
    "any_change"
]]

audit_trail.to_csv("ai_review_audit_trail.csv", index=False)

print(audit_trail)

The audit trail is useful for:

• debugging
• compliance
• reviewer training
• model evaluation
• future retraining
• stakeholder trust

Without an audit trail, human review becomes a manual activity that disappears after the correction is made.

With an audit trail, review becomes structured feedback data.

Step 11: Visualise review outcomes

A simple chart can show how many items were auto-approved, approved, or corrected.

import matplotlib.pyplot as plt

status_counts = validated["review_status"].value_counts()

plt.figure(figsize=(8, 5))
plt.bar(status_counts.index, status_counts.values)
plt.xlabel("Review status")
plt.ylabel("Number of items")
plt.title("AI validation outcomes")
plt.tight_layout()
plt.show()

We can also visualise correction rate by source.

plt.figure(figsize=(8, 5))
plt.bar(
    corrections_by_source["source"],
    corrections_by_source["correction_rate"]
)
plt.xlabel("Source")
plt.ylabel("Correction rate")
plt.title("Correction rate by source")
plt.tight_layout()
plt.show()

These simple visuals can help teams understand how the validation process is performing.

What this workflow gives us

This small workflow creates a basic but useful human-in-the-loop validation layer.

It helps answer:

• Which AI outputs were auto-approved?
• Which outputs required human review?
• Which predictions were corrected?
• What final values should downstream systems use?
• Which categories or sources produce more corrections?
• What evidence exists for audit and improvement?

This is the foundation of a more reliable AI system.

How to extend this into production

This tutorial uses CSV files and simulated reviewer decisions, but the same logic can be extended into a real production workflow.

Possible improvements include:

Use a database

Store predictions, reviews, and final outputs in PostgreSQL, MySQL, BigQuery, Snowflake, or another database.

Add reviewer identities

Track who reviewed each item and when.

Add timestamps

Store prediction time, review time, approval time, and update time.

Add role-based access

Different users may need different permissions:

• reviewer
• validator
• admin
• analyst
• auditor

Add review queues

A real system should assign work based on priority, risk, category, source, or reviewer capacity.

Add model feedback

Corrected outputs can be reused for future evaluation, prompt tuning, or model retraining.

Add alerts

If correction rates increase, the system should alert the team.

Add data quality monitoring

Track whether errors are coming from model weakness, source data problems, or unclear business definitions.

Common mistakes to avoid

Treating human review as temporary

Some teams assume human review is only needed until the model gets better. In reality, many production systems need review permanently for high-risk or uncertain cases.

Not storing corrections

If corrections are made but not stored, the system cannot learn from them.

Reviewing too much

If every item needs review, the AI system may not be reducing work. Review rules should focus attention where it matters most.

Reviewing too little

If the system auto-approves risky or low-confidence outputs, mistakes can quietly enter production.

Not defining final ownership

Someone must own the final approved value. Is it the model, the reviewer, the validator, or the business team?

Final thoughts

Human-in-the-loop validation is not a sign that an AI system has failed.

It is often the process that makes the AI system safe enough to use.

The important difference is structure.

Unstructured review creates bottlenecks, confusion, and lost feedback. Structured review creates validated outputs, audit trails, correction datasets, and continuous improvement.

For production AI systems, the goal is not to remove humans from every decision. The goal is to use human judgment where it adds the most value, and then capture that judgment as data.

That is how human review becomes part of the AI system, not a workaround for it.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

Can we make this work?

A developer builds a prototype. A product manager tests a few prompts. The results look useful. Soon the team wants to connect the model to documents, internal tools, customer data, support tickets, reports, or dashboards.

That is usually when the harder questions appear:

• How much will this cost when more users start using it?
• Which workflows are generating the most token usage?
• Are we sending sensitive data into the model?
• Which prompts are risky?
• Are we tracking failures, retries, and expensive requests?
• Can we explain where the data went?

Many AI prototypes fail not because the model is weak, but because the team has no visibility around cost, risk, and data movement.

In this article, we will build a simple Python-based LLM cost and risk visibility tool. The goal is not to create a full governance platform. The goal is to show how a small team can start monitoring AI usage in a practical way.

We will create a lightweight workflow that:

• loads sample LLM request logs;
• estimates token-based cost;
• flags risky prompts;
• groups usage by user and workflow;
• identifies expensive or sensitive requests;
• produces a simple summary report.

This pattern can be extended later into a dashboard, database table, Streamlit app, or internal monitoring tool.

Why cost and risk visibility matters

LLM usage can grow quietly.

One user testing a chatbot may not cost much. But once an AI feature is used across a team, costs can increase through:

• long prompts;
• large uploaded documents;
• repeated retries;
• agent loops;
• unnecessary model calls;
• verbose responses;
• background automations;
• multiple tools calling the model without tracking.

Risk grows in a similar way.

At the prototype stage, users may paste test data. In production, they may paste customer details, financial information, confidential documents, employee data, or business-sensitive material.

If the application does not log what is happening, the team cannot manage it.

A simple monitoring layer can help answer:

• Which workflows use the most tokens?
• Which users or teams generate the highest cost?
• Which prompts contain sensitive terms?
• Which model calls failed or retried?
• Which requests need review?

Related posts on DataScience+:

• Machine Learning for Diabetes with Python
• Topic Modeling in Python with NLTK and Gensim
• Top Python Libraries for Machine Learning

Sample dataset

For this example, we will use a small CSV file called llm_requests.csv.

Each row represents one LLM request.

request_id,user_id,workflow,model,prompt,response,status
1,u001,customer_support,gpt-4,"Summarise this customer complaint: My card number is 4111 1111 1111 1111 and I was charged twice.","The customer reports a duplicate charge.",success
2,u002,sales_email,gpt-4,"Write a follow-up email for a SaaS lead interested in pricing.","Here is a professional follow-up email...",success
3,u003,hr_policy,gpt-4,"Can you review this employee performance note and suggest improvements?","The note can be rewritten as follows...",success
4,u001,customer_support,gpt-3.5,"Classify this ticket as billing, technical, or account access.","Billing",success
5,u004,legal_review,gpt-4,"Review this contract clause and identify risks related to liability.","The clause may create liability exposure...",success
6,u005,analytics,gpt-3.5,"Generate SQL to calculate monthly active users from the events table.","SELECT DATE_TRUNC('month', event_time)...",success
7,u006,finance,gpt-4,"Analyse this invoice and bank account number 12345678 for payment validation.","The invoice appears valid...",success

In a real system, these logs could come from:

• application logs;
• API gateway logs;
• database tables;
• prompt logging middleware;
• SaaS AI tool exports;
• internal workflow automation systems.

For this tutorial, the CSV keeps things simple.

Step 1: Load the data

import pandas as pd

df = pd.read_csv("llm_requests.csv")

print(df.head())

We should also check the basic structure.

print(df.info())
print(df["workflow"].value_counts())

This gives us a starting point: who is using the system, which workflows are active, and what kind of model calls are being made.

Step 2: Estimate token usage

In production, the best option is to use actual token counts returned by the model provider.

For this example, we will use a rough approximation:

• 1 token is around 4 characters in English text.
• Total tokens = prompt tokens + response tokens.

This is not exact, but it is good enough for a basic monitoring prototype.

def estimate_tokens(text):
    if pd.isna(text):
        return 0
    return max(1, round(len(str(text)) / 4))

df["prompt_tokens"] = df["prompt"].apply(estimate_tokens)
df["response_tokens"] = df["response"].apply(estimate_tokens)
df["total_tokens"] = df["prompt_tokens"] + df["response_tokens"]

print(df[["request_id", "workflow", "model", "total_tokens"]])

In a real application, replace this approximation with the token usage returned by the API.

Step 3: Add model pricing

Different models have different costs. We can create a simple pricing table.

The numbers below are only example prices. You should replace them with current pricing from your model provider.

pricing = {
    "gpt-4": {
        "input_per_1k": 0.03,
        "output_per_1k": 0.06
    },
    "gpt-3.5": {
        "input_per_1k": 0.0015,
        "output_per_1k": 0.002
    }
}

Now we can estimate cost per request.

def estimate_cost(row):
    model = row["model"]

    if model not in pricing:
        return 0

    input_cost = (row["prompt_tokens"] / 1000) * pricing[model]["input_per_1k"]
    output_cost = (row["response_tokens"] / 1000) * pricing[model]["output_per_1k"]

    return input_cost + output_cost

df["estimated_cost_usd"] = df.apply(estimate_cost, axis=1)

print(df[["request_id", "workflow", "model", "total_tokens", "estimated_cost_usd"]])

This gives us an estimated cost for each LLM request.

Step 4: Summarise cost by workflow

Cost is more useful when grouped by business workflow.

workflow_cost = (
    df.groupby("workflow")
    .agg(
        requests=("request_id", "count"),
        total_tokens=("total_tokens", "sum"),
        estimated_cost_usd=("estimated_cost_usd", "sum")
    )
    .reset_index()
    .sort_values("estimated_cost_usd", ascending=False)
)

print(workflow_cost)

This helps answer:

• Which workflow is most expensive?
• Which workflows may need prompt optimisation?
• Which use cases are growing fastest?

A workflow with high cost is not automatically bad. It may be valuable. But without this visibility, the team cannot make informed decisions.

Step 5: Summarise cost by user

User-level monitoring can help detect unusual behaviour.

user_cost = (
    df.groupby("user_id")
    .agg(
        requests=("request_id", "count"),
        total_tokens=("total_tokens", "sum"),
        estimated_cost_usd=("estimated_cost_usd", "sum")
    )
    .reset_index()
    .sort_values("estimated_cost_usd", ascending=False)
)

print(user_cost)

This can be useful for internal AI tools where different teams use the same system.

For example, if one user or team generates unusually high usage, it may indicate:

• a genuine high-value use case;
• a prompt that is too long;
• repeated retries;
• misuse;
• an automation loop;
• lack of user training.

Step 6: Flag risky prompts

Now we can add a simple risk scanner.

This is not a full data-loss-prevention system. It is a basic first layer that flags prompts containing sensitive terms or patterns.

import re

risk_patterns = {
    "payment_card": r"b(?:d[ -]*?){13,16}b",
    "bank_account": r"bd{8}b",
    "confidential_terms": r"b(confidential|private|secret|internal only)b",
    "employee_data": r"b(employee|performance note|salary|disciplinary)b",
    "legal_content": r"b(contract|liability|clause|legal)b"
}

Now apply the patterns to each prompt.

def detect_risks(text):
    if pd.isna(text):
        return []

    detected = []
    text = str(text).lower()

    for risk_name, pattern in risk_patterns.items():
        if re.search(pattern, text, flags=re.IGNORECASE):
            detected.append(risk_name)

    return detected

df["risk_flags"] = df["prompt"].apply(detect_risks)
df["risk_count"] = df["risk_flags"].apply(len)
df["has_risk"] = df["risk_count"] > 0

print(df[["request_id", "workflow", "prompt", "risk_flags"]])

This gives us a basic view of which prompts may need review.

Step 7: Create a risk summary

We can now summarise risk by workflow.

risk_summary = (
    df.groupby("workflow")
    .agg(
        total_requests=("request_id", "count"),
        risky_requests=("has_risk", "sum"),
        estimated_cost_usd=("estimated_cost_usd", "sum")
    )
    .reset_index()
)

risk_summary["risk_rate"] = (
    risk_summary["risky_requests"] / risk_summary["total_requests"]
)

risk_summary = risk_summary.sort_values("risk_rate", ascending=False)

print(risk_summary)

This helps identify which workflows are most likely to involve sensitive or high-risk content.

For example:

• HR workflows may contain employee data.
• Finance workflows may contain bank details.
• Legal workflows may contain contracts.
• Customer support workflows may contain personal information.

The goal is not to block every request. The goal is to understand which workflows need stronger controls.

Step 8: Add a simple priority score

A useful monitoring tool should help teams decide what to review first.

We can create a simple priority score based on:

• risk count;
• estimated cost;
• model used;
• workflow type.

high_risk_workflows = ["finance", "legal_review", "hr_policy"]

def calculate_priority(row):
    score = 0

    # Risk flags
    score += row["risk_count"] * 3

    # Expensive model
    if row["model"] == "gpt-4":
        score += 2

    # High-risk workflow
    if row["workflow"] in high_risk_workflows:
        score += 3

    # Higher token usage
    if row["total_tokens"] > 100:
        score += 1

    return score

df["priority_score"] = df.apply(calculate_priority, axis=1)

review_queue = df.sort_values("priority_score", ascending=False)

print(review_queue[[
    "request_id",
    "workflow",
    "model",
    "total_tokens",
    "estimated_cost_usd",
    "risk_flags",
    "priority_score"
]])

This creates a lightweight review queue.

Requests with higher scores may need:

• manual review;
• prompt rewriting;
• workflow restrictions;
• user training;
• model downgrade;
• stronger data controls.

Step 9: Generate a simple report

Now we can generate a short summary.

total_requests = len(df)
total_cost = df["estimated_cost_usd"].sum()
risky_requests = df["has_risk"].sum()
risk_rate = risky_requests / total_requests

print("LLM Usage Summary")
print("-----------------")
print(f"Total requests: {total_requests}")
print(f"Estimated cost: ${total_cost:.4f}")
print(f"Risky requests: {risky_requests}")
print(f"Risk rate: {risk_rate:.1%}")

print("nTop workflows by cost:")
print(workflow_cost.head())

print("nHighest priority requests:")
print(review_queue[[
    "request_id",
    "workflow",
    "model",
    "risk_flags",
    "priority_score"
]].head())

This kind of report can be run daily or weekly.

It can also be exported to CSV.

workflow_cost.to_csv("workflow_cost_summary.csv", index=False)
risk_summary.to_csv("workflow_risk_summary.csv", index=False)
review_queue.to_csv("llm_review_queue.csv", index=False)

Step 10: Visualise cost by workflow

A simple bar chart can make the result easier to understand.

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.bar(workflow_cost["workflow"], workflow_cost["estimated_cost_usd"])
plt.xlabel("Workflow")
plt.ylabel("Estimated cost in USD")
plt.title("Estimated LLM cost by workflow")
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
plt.show()

You can also visualise risky requests.

plt.figure(figsize=(10, 6))
plt.bar(risk_summary["workflow"], risk_summary["risky_requests"])
plt.xlabel("Workflow")
plt.ylabel("Risky requests")
plt.title("Risky LLM requests by workflow")
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
plt.show()

These charts are simple, but they are enough to start a conversation with product, data, security, or finance teams.

What this prototype shows

This small project gives us a basic visibility layer around LLM usage.

It helps answer:

• Which workflows use the most tokens?
• Which workflows cost the most?
• Which users generate the most usage?
• Which prompts contain risky content?
• Which requests should be reviewed first?

This is not a complete AI governance system, but it is a practical starting point.

Many organisations do not need a complex platform on day one. They need a simple way to see what is happening.

How to improve this further

This prototype can be extended in many ways.

Use real token counts

Instead of estimating tokens by character length, collect actual token usage from the model provider.

Add user and team metadata

Join request logs with team, department, or cost-centre data.

Track success and failure rates

Add fields such as:

• error type;
• retry count;
• latency;
• timeout;
• fallback model;
• user feedback.

Add data classification

Use a more advanced scanner to detect:

• personally identifiable information;
• financial data;
• health data;
• legal data;
• customer records;
• source code;
• credentials.

Build a Streamlit dashboard

A simple Streamlit interface could show:

• total cost;
• cost by workflow;
• cost by user;
• risky prompts;
• review queue;
• model usage;
• trend over time.

Store logs in a database

Instead of using CSV files, store logs in PostgreSQL, BigQuery, Snowflake, or another analytics database.

Add governance actions

For high-risk requests, the system could:

• flag for review;
• block the request;
• redact sensitive values;
• require approval;
• route to a safer model;
• warn the user before submission.

Final thoughts

AI systems do not become production-ready just because the model works.

They need visibility.

Teams need to understand how data moves, which workflows create risk, where costs are coming from, and which requests require human review.

A simple Python monitoring layer can provide that first level of visibility. It does not need to be perfect. It just needs to make invisible problems visible.

For teams adopting LLMs, that is often the difference between an exciting prototype and a system that can actually be trusted in production.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

I kept wondering how much of that dark interior anyone actually lives in, so I built the map myself, for Nevada, and then measured the thing the shading hides. Nevada really does look like one of the worst hospital deserts in the country, even though almost everyone in it lives in Las Vegas or Reno, a few minutes from an emergency room. The fix is to weight the same result by population instead of by area, and the difference turns out to be the whole story. And because a gap measured in one state could be waved off as just how the measurement behaves, I run the identical pipeline on Massachusetts, a state with the opposite geography, and let the two results face each other.

Along the way you get a reusable toolkit: pulling a facility list from an API, geocoding each facility from its street address, measuring distance to the nearest one with sf, and weighting a spatial result by population instead of by area. Everything is open data and every package is on CRAN, and the whole pipeline needs no API key at all.

library(httr2)
library(dplyr)
library(tidyr)
library(zipcodeR)
library(tidygeocoder)
library(sf)
library(tigris)
library(ggplot2)

Nevada hospitals data

The Centers for Medicare & Medicaid Services publishes a Hospital General Information file listing every hospital in the country. It is queryable as JSON, and its datastore API takes a conditions filter, so I ask for Nevada only and get the whole state’s list in one page, no paging loop.

url <- "https://data.cms.gov/provider-data/api/1/datastore/query/xubh-q36u/0"

hosp_raw <- request(url) |>
  req_url_query(
    limit = 500,
    `conditions[0][property]` = "state",
    `conditions[0][operator]` = "=",
    `conditions[0][value]` = "NV"
  ) |>
  req_perform() |>
  resp_body_json() |>
  (\(x) bind_rows(lapply(x$results, as_tibble)))()

nrow(hosp_raw)
## [1] 46

That is every hospital in Nevada. But "hospital" is a broad label, and I want emergency care. The file carries an emergency_services flag and a hospital_type, so I keep the emergency-capable set: general acute-care, critical-access, and rural-emergency hospitals that report an ER.

er_types <- c(
  "Acute Care Hospitals",
  "Critical Access Hospitals",
  "Rural Emergency Hospital"
)

er_hosp <- hosp_raw |>
  filter(hospital_type %in% er_types, emergency_services == "Yes") |>
  mutate(zip5 = substr(zip_code, 1, 5))

nrow(er_hosp)
## [1] 34

Related posts on DataScience+:

• How Happy is Your Country? — Happy Planet Index Visualized
• Analysis and Visualization of Blue Bikes Sharing in Boston
• Machine Learning Results in R: one plot to rule them all! (Part 1 – Classification Models)

Geocoding hospital addresses

The file gives street addresses but no coordinates, so I have to geocode them. The US Census runs a free, keyless batch geocoder, and tidygeocoder talks to it straight from R. With only a few dozen Nevada hospitals this is a single quick call, no batching.

geo <- er_hosp |>
  geocode(
    street = address,
    city = citytown,
    state = state,
    postalcode = zip5,
    method = "census",
    quiet = TRUE
  )

The Census matcher resolves most of the addresses; the few it misses usually have a PO-box or otherwise unmatchable street line. For those I fall back to the ZIP centroid, which is fine as a minority backstop, with a coalesce().

zip_ll <- zip_code_db |> transmute(zip5 = zipcode, zlat = lat, zlng = lng)

er_geo <- geo |>
  left_join(zip_ll, by = "zip5") |>
  mutate(
    long = coalesce(long, zlng), # Census miss -> ZIP centroid
    lat = coalesce(lat, zlat)
  ) |>
  filter(!is.na(lat))

nrow(er_geo)
## [1] 34

Now I turn the coordinates into a spatial object and project it. Modern sf, through its s2 backend, measures true distances on the sphere even from raw longitude and latitude, so the old rule that you must never compute distance on lon/lat no longer bites. I still project to EPSG:5070, an equal-area continental projection in meters, because the area and centroid calculations later behave in a flat equal-area space and the map draws in the right shape.

er <- er_geo |>
  st_as_sf(coords = c("long", "lat"), crs = 4326) |>
  st_transform(5070)

The desert map

I pull Nevada’s census tracts with tigris (geometry only, no survey data). The tract geometry carries ALAND, its land area, which I will need later.

tr <- tracts("Nevada", cb = TRUE, year = 2022, progress_bar = FALSE) |>
  st_transform(5070)

The measurement is two sf calls. For each tract centroid, st_nearest_feature() finds the closest emergency-capable hospital, and st_distance(..., by_element = TRUE) returns that one distance. I divide by 1609.34 for miles.

ctr <- st_centroid(tr)
i <- st_nearest_feature(ctr, er)
tr$er_mi <- as.numeric(st_distance(ctr, er[i, ], by_element = TRUE)) / 1609.34

Now the desert map. I shade each tract by its distance on a square-root color scale, so a few very remote tracts do not stretch it and flatten everything else, and I drop a dot on every emergency-capable hospital so you can see where the care actually sits. I use one editorial theme, tinted light gray with the axes stripped away; copy it into your own maps.

rng <- range(tr$er_mi, na.rm = TRUE)

ggplot(tr) +
  geom_sf(aes(fill = er_mi), color = NA) +
  # a dot on every ER hospital: they cluster in the pale near-ER areas
  # (Las Vegas, Reno) and vanish from the dark empty center
  geom_sf(
    data = er,
    inherit.aes = FALSE,
    shape = 21,
    size = 1.1,
    stroke = 0.25,
    fill = "#111111",
    color = "white"
  ) +
  scale_fill_viridis_c(
    option = "magma",
    direction = -1,
    transform = "sqrt",
    limits = rng,
    breaks = c(5, 10, 20, 40),
    name = "Miles to nearest ER  "
  ) +
  guides(
    fill = guide_colorbar(
      barwidth = 14,
      barheight = 0.5,
      title.position = "top",
      title.hjust = 0.5
    )
  ) +
  labs(
    title = "Nevada: how far is the nearest ER?",
    subtitle = "Distance from each census tract to the nearest hospital with a general ER",
    caption = "Data: CMS Hospital General Information"
  ) +
  theme_void(base_size = 12) +
  theme(
    plot.background = element_rect(fill = "#ECECEF", color = NA),
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 17),
    plot.subtitle = element_text(color = "grey30")
  )

The pattern looks damning: a vast dark center hundreds of miles from an ER, with the hospital dots huddled into a couple of small clusters. Read as "where hospitals are out of reach," the map is misleading. The next section shows why.

The problem: the map colors land, not people

Look at the dark tracts. They are real, they genuinely are far from an ER, but they are enormous and nearly empty. The people live in the small pale specks: Las Vegas and Reno. A choropleth gives every tract visual weight in proportion to its area, so a thousand square miles of empty desert dominates the image while a dense city takes up a few pixels. The map answers "how much land is far from an ER," when the question we care about is "how many people are."

I can answer both with the data already in hand. Census tracts are drawn to hold roughly equal population, so counting tracts is a good stand-in for counting people, and ALAND gives me area. I compute the share of Nevada within 20 miles of an ER weighted two ways: by people (each tract equal) and by land (each tract weighted by area).

tr_tab <- st_drop_geometry(tr)

by_people <- mean(tr_tab$er_mi <= 20) * 100
by_land <- sum(tr_tab$ALAND[tr_tab$er_mi <= 20]) / sum(tr_tab$ALAND) * 100

c(by_people = round(by_people), by_land = round(by_land))
## by_people   by_land 
##        94        17

There is the whole story in two numbers. 94% of Nevadans live within 20 miles of an ER, but only 17% of the state’s land does, a gap of 78 points. The desert is real as geography and almost empty as a matter of people.

A control state with the opposite geography

Two numbers from one state are suggestive, not conclusive. Maybe a gap like that is just what this measurement produces everywhere, and every state’s map overstates its desert by about the same amount. The way to rule that out is a control with the opposite geography, so I picked Massachusetts: small, dense, and with emergency care spread across the whole state. Every step above is the same handful of calls with a different state code, so I fold the pipeline into one function and run it once more.

desert_tracts <- function(abbr, name) {
  er <- request(url) |>
    req_url_query(
      limit = 500,
      `conditions[0][property]` = "state",
      `conditions[0][operator]` = "=",
      `conditions[0][value]` = abbr
    ) |>
    req_perform() |>
    resp_body_json() |>
    (\(x) bind_rows(lapply(x$results, as_tibble)))() |>
    filter(hospital_type %in% er_types, emergency_services == "Yes") |>
    mutate(zip5 = substr(zip_code, 1, 5)) |>
    geocode(
      street = address,
      city = citytown,
      state = state,
      postalcode = zip5,
      method = "census",
      quiet = TRUE
    ) |>
    left_join(zip_ll, by = "zip5") |>
    mutate(long = coalesce(long, zlng), lat = coalesce(lat, zlat)) |>
    filter(!is.na(lat)) |>
    st_as_sf(coords = c("long", "lat"), crs = 4326) |>
    st_transform(5070)

  tr <- tracts(name, cb = TRUE, year = 2022, progress_bar = FALSE) |>
    st_transform(5070)
  ctr <- st_centroid(tr)
  i <- st_nearest_feature(ctr, er)
  tr$er_mi <- as.numeric(st_distance(ctr, er[i, ], by_element = TRUE)) / 1609.34
  tr$state <- name

  list(er = er, tracts = tr)
}

ma <- desert_tracts("MA", "Massachusetts")

First, the same desert map for Massachusetts, drawn on the same square-root color scale with Nevada’s limits, so a shade of color means the same distance in both maps.

ggplot(ma$tracts) +
  geom_sf(aes(fill = er_mi), color = NA) +
  geom_sf(
    data = ma$er,
    inherit.aes = FALSE,
    shape = 21,
    size = 1.1,
    stroke = 0.25,
    fill = "#111111",
    color = "white"
  ) +
  scale_fill_viridis_c(
    option = "magma",
    direction = -1,
    transform = "sqrt",
    limits = rng,
    breaks = c(5, 10, 20, 40),
    name = "Miles to nearest ER  "
  ) +
  guides(
    fill = guide_colorbar(
      barwidth = 14,
      barheight = 0.5,
      title.position = "top",
      title.hjust = 0.5
    )
  ) +
  labs(
    title = "Massachusetts: how far is the nearest ER?",
    subtitle = "Same measurement and color scale as the Nevada map",
    caption = "Data: CMS Hospital General Information"
  ) +
  theme_void(base_size = 12) +
  theme(
    plot.background = element_rect(fill = "#ECECEF", color = NA),
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 17),
    plot.subtitle = element_text(color = "grey30")
  )

On Nevada’s scale Massachusetts never gets near the dark end: its deepest shade is the outer tip of Cape Cod, about 27 miles from an ER, while Nevada’s remotest tract sits 3.9 times farther out. The 57 emergency-capable hospitals dot the state end to end instead of huddling in two corners. Nothing here tempts a "desert" headline, which is what a control should look like.

Now the same 20-mile shares, for both states side by side.

both <- bind_rows(
  tibble(state = "Nevada", er_mi = tr_tab$er_mi, ALAND = tr_tab$ALAND),
  ma$tracts |> st_drop_geometry() |> select(state, er_mi, ALAND)
)

shares <- both |>
  group_by(state) |>
  summarise(
    by_people = mean(er_mi <= 20) * 100,
    by_land = sum(ALAND[er_mi <= 20]) / sum(ALAND) * 100
  ) |>
  mutate(gap = by_people - by_land)

shares
## # A tibble: 2 × 4
##   state         by_people by_land    gap
##   <chr>             <dbl>   <dbl>  <dbl>
## 1 Massachusetts      99.8    99.3  0.486
## 2 Nevada             94.4    16.5 77.8

Massachusetts barely registers a gap: 99.8% of its people and 99.3% of its land sit within 20 miles of an ER, 0.5 points apart, against Nevada’s 78. So the gap is not something the measurement manufactures on its own. It is a property of Nevada’s geography, and it took the control to earn that sentence.

The cumulative curves put the contrast in one picture. For each state I plot the share within a given distance weighted by people next to the same share weighted by land.

curve <- both |>
  group_by(state) |>
  reframe(
    mi = seq(0, 60, 1),
    People = sapply(mi, \(m) mean(er_mi <= m) * 100),
    Land = sapply(mi, \(m) sum(ALAND[er_mi <= m]) / sum(ALAND) * 100)
  ) |>
  pivot_longer(c(People, Land), names_to = "weight", values_to = "pct")

ggplot(curve, aes(mi, pct, color = weight)) +
  geom_line(linewidth = 1.1) +
  facet_wrap(~ factor(state, c("Nevada", "Massachusetts"))) +
  scale_color_manual(values = c(People = "#0066CC", Land = "#E8862D")) +
  labs(
    title = "Same measurement, two states, counted two ways",
    subtitle = "Share of each state within a distance of an ER, weighted by population vs by land area",
    x = "Miles to nearest ER",
    y = "% within",
    color = NULL,
    caption = "Data: CMS Hospital General Information; tract land area from US Census TIGER"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.background = element_rect(fill = "#ECECEF", color = NA),
    panel.background = element_rect(fill = "#ECECEF", color = NA),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(color = "grey78"),
    axis.ticks = element_blank(),
    strip.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(color = "grey30"),
    legend.position = "bottom"
  )

In the Nevada panel the blue people line snaps up to near 100% within a few miles while the orange land line crawls: the signature of a place that looks empty on a map but lives in a couple of cities. In Massachusetts the two lines climb together and are finished by about 20 miles, because when people and hospitals spread over the land together there is nothing for the two weightings to disagree about. The divergence, not the dark shading, is what a desert of land rather than a desert of people looks like, and it is exactly what the choropleth hides.

(Two caveats worth naming. First, equal-population tracts are an approximation, and rural tracts run a little smaller, so if anything they understate how concentrated the population really is: weighting by an actual population count only widens the gap. Second, each state sees only its own hospitals, so a few border towns that actually rely on an out-of-state ER, Mesquite near St. George, Utah, or Laughlin near Bullhead City, Arizona, read as farther from care than they are. That, too, only overstates the desert.)

Make it your own

The pipeline is general: desert_tracts() already takes any state, or point the same steps at any table of facilities, pharmacies, clinics, grocery stores, and st_nearest_feature() measures access to whatever you give it. The last figure is the piece worth carrying past this dataset. Any time you shade a map by a rate or a distance over regions that vary wildly in population, from rural tracts to whole countries, the map is weighting by area, and your eye reads it as importance. Before you call something a desert, weight it by the people who actually live there, and check a place where you expect no desert at all, so you know what agreement looks like. Sometimes the desert stays. In Nevada, I found, almost nobody lives in it.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

That bothered me, because k-means will happily partition anything. Hand it a featureless cloud of noise and it returns tidy, well-separated groups. So I wanted to answer a narrow question: can the standard validation toolkit tell the difference between data that really contains groups and data that does not? I ran the experiment, and the answer surprised me enough to write it up.

Short version: bootstrap stability, the check most often held up as the rigorous one, could not tell them apart at all.

The data

I used the Coffee Quality Institute ratings distributed through TidyTuesday: professional cuppers scoring arabica lots on aroma, flavor, aftertaste, acidity, body and balance. It is a good test case because "are there distinct coffee flavor profiles?" is a real question someone would reach for k-means to answer.

library(tidyverse)
library(cluster)
library(fpc)
library(diptest)

url <- paste0("https://raw.githubusercontent.com/rfordatascience/tidytuesday/",
              "master/data/2020/2020-07-07/coffee_ratings.csv")

coffee <- read_csv(url, show_col_types = FALSE) |>
  filter(species == "Arabica") |>
  select(aroma, flavor, aftertaste, acidity, body, balance) |>
  drop_na() |>
  filter(if_all(everything(), ~ .x > 0))   # a few lots scored all zeros

X <- scale(coffee)
dim(X)
## [1] 1310    6

That leaves 1310 coffees on 6 sensory scores. One thing to note before clustering, because it matters later:

round(cor(coffee), 2)
##            aroma flavor aftertaste acidity body balance
## aroma       1.00   0.74       0.69    0.61 0.55    0.61
## flavor      0.74   1.00       0.86    0.74 0.66    0.73
## aftertaste  0.69   0.86       1.00    0.71 0.67    0.76
## acidity     0.61   0.74       0.71    1.00 0.61    0.64
## body        0.55   0.66       0.67    0.61 1.00    0.67
## balance     0.61   0.73       0.76    0.64 0.67    1.00

pca <- prcomp(X)
round(summary(pca)$importance[2, 1:3], 3)
##   PC1   PC2   PC3 
## 0.739 0.079 0.063

Every score correlates with every other score between 0.55 and 0.86, and the first principal component absorbs 73.9% of the variance with near-identical loadings on all six variables. Coffees that score well on aroma score well on everything. This looks like one general quality dimension rather than distinct profiles, but let us proceed exactly as a tutorial would and see what the diagnostics say.

Related posts on DataScience+:

• K Means Clustering in R
• Predict Customer Churn – Logistic Regression, Decision Tree and Random Forest
• A Gentle Introduction on Market Basket Analysis — Association Rules

The standard pipeline

dsp_colors <- c("#0066CC", "#E8862D", "#159A6C", "#7D5BD6",
                "#D64580", "#2AA9B8", "#C9A227")
dsp_theme <- theme_minimal(base_size = 13) +
  theme(plot.background    = element_rect(fill = "#ECECEF", color = NA),
        panel.background   = element_rect(fill = "#ECECEF", color = NA),
        panel.grid.minor   = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color = "grey78"),
        axis.ticks         = element_blank(),
        plot.title         = element_text(face = "bold"),
        strip.text         = element_text(face = "bold"))

elbow <- map_dfr(1:10, ~ tibble(
  k   = .x,
  wss = kmeans(X, .x, nstart = 25, iter.max = 50)$tot.withinss
))

ggplot(elbow, aes(k, wss)) +
  geom_line(color = dsp_colors[1], linewidth = 0.9) +
  geom_point(color = dsp_colors[1], size = 2.4) +
  scale_x_continuous(breaks = 1:10) +
  labs(title = "The elbow says three",
       x = "Number of clusters (k)",
       y = "Total within-cluster sum of squares") +
  dsp_theme

A textbook elbow at k = 3. The drop in within-cluster sum of squares is 3,354 going from one cluster to two, 1,107 going to three, then only 490 going to four. So k = 3 it is.

km3 <- kmeans(X, 3, nstart = 25, iter.max = 50)

coffee |>
  mutate(cluster = km3$cluster) |>
  group_by(cluster) |>
  summarise(n = n(), across(everything(), ~ round(mean(.x), 2))) |>
  arrange(flavor)
## # A tibble: 3 × 8
##   cluster     n aroma flavor aftertaste acidity  body balance
##     <int> <dbl> <dbl>  <dbl>      <dbl>   <dbl> <dbl>   <dbl>
## 1       2   298  7.23   7.09       6.97    7.19  7.22    7.1 
## 2       1   706  7.58   7.54       7.42    7.53  7.52    7.54
## 3       3   306  7.88   7.89       7.79    7.9   7.83    7.89

Three clean groups. But look at the summary rows: the clusters are ordered on every single variable at once. Cluster means rise together across aroma, flavor, aftertaste, acidity, body and balance. These are not flavor profiles, they are quality tiers, which is what you get when you cut a one-dimensional gradient into three pieces. The cluster means on PC1 confirm it:

round(tapply(pca$x[, 1], km3$cluster, mean), 2)
##     1     2     3 
## -0.06  2.82 -2.61

The three checks, and how they do

Now the validation. I ran the three diagnostics that show up most often.

Silhouette width measures how much closer each point sits to its own cluster than to the next nearest one, averaged over all points. Kaufman and Rousseeuw’s rule of thumb is that below 0.25 there is no substantial structure, and 0.25 to 0.50 is weak.

sil_k <- function(M, k) {
  km <- kmeans(M, k, nstart = 25, iter.max = 50)
  mean(silhouette(km$cluster, dist(M))[, 3])
}
coffee_sil <- sil_k(X, 3)
round(coffee_sil, 3)
## [1] 0.303

The gap statistic compares the within-cluster dispersion to what you would get on a structureless reference distribution, and cluster::clusGap picks a k for you.

gap_coffee <- clusGap(X, FUN = kmeans, nstart = 25, K.max = 8, B = 50)
maxSE(gap_coffee$Tab[, "gap"], gap_coffee$Tab[, "SE.sim"], method = "firstSEmax")
## [1] 3

round(gap_coffee$Tab[, "gap"], 3)
## [1] 1.138 1.200 1.248 1.250 1.251 1.245 1.235 1.223

Bootstrap stability is the one usually presented as the serious check. fpc::clusterboot resamples the data, re-clusters, matches the new clusters to the original ones and reports a mean Jaccard similarity per cluster. Above 0.85 is conventionally called highly stable.

jac_k <- function(M, k, seed) {
  cb <- clusterboot(M, B = 100, clustermethod = kmeansCBI, krange = k,
                    runs = 25, count = FALSE, seed = seed)
  mean(cb$bootmean)
}
coffee_jac <- jac_k(X, 3, seed = 7)
round(coffee_jac, 3)
## [1] 0.963

So: silhouette 0.303 (weak, but tutorials publish worse), gap statistic endorsing k = 3, and a bootstrap Jaccard of 0.963, comfortably in "highly stable" territory. Two of the three checks pass, and the one that is lukewarm would not stop most people.

The control

Here is the part I actually wanted to run. I generated a dataset with the same number of rows, the same number of variables and the same covariance matrix as the coffee scores, drawn from a single multivariate normal. By construction it contains zero clusters. Then I put it through the identical pipeline.

null_like <- function(M) {
  N <- MASS::mvrnorm(nrow(M), mu = rep(0, ncol(M)), Sigma = cov(M))
  colnames(N) <- colnames(M)
  N
}
X_null <- null_like(X)

project <- function(M, label) {
  pc <- prcomp(M)$x
  tibble(PC1 = pc[, 1], PC2 = pc[, 2], data = label,
         cluster = factor(kmeans(M, 3, nstart = 25, iter.max = 50)$cluster))
}

bind_rows(project(X, "Coffee scores (real data)"),
          project(X_null, "Random data, no clusters at all")) |>
  ggplot(aes(PC1, PC2, color = cluster)) +
  geom_point(alpha = 0.55, size = 1.2) +
  facet_wrap(~ data) +
  scale_color_manual(values = dsp_colors) +
  labs(title = "Same pipeline, same picture", x = "PC1", y = "PC2") +
  dsp_theme +
  theme(legend.position = "none")

Both panels show three crisp, well-separated, roughly equal bands. The right panel is pure noise. Nothing in that picture distinguishes the real data from the fake, because in both cases k-means is doing the same thing: slicing a single elongated blob into three along its longest axis.

The diagnostics on the noise:

c(silhouette = round(sil_k(X_null, 3), 3),
  jaccard    = round(jac_k(X_null, 3, seed = 7), 3))
## silhouette    jaccard 
##      0.273      0.925

A structureless cloud scores about as well as the coffee data on silhouette and just as well on bootstrap stability.

Calibrating against the null

One draw could be luck, so I repeated it. Twenty null datasets, each matched to the real data’s size and covariance, each pushed through the same two checks. I did the same for a positive control: the palmerpenguins body measurements, where three real species are known to be present.

penguins_X <- palmerpenguins::penguins |>
  drop_na(bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) |>
  select(bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) |>
  scale()

calibrate <- function(M, k, R = 20) {
  map_dfr(1:R, function(i) {
    N <- null_like(M)
    tibble(rep = i, silhouette = sil_k(N, k), jaccard = jac_k(N, k, seed = i))
  })
}

null_coffee   <- calibrate(X, 3)
null_penguins <- calibrate(penguins_X, 3)

observed <- tibble(
  data       = c("Coffee", "Penguins"),
  silhouette = c(coffee_sil, sil_k(penguins_X, 3)),
  jaccard    = c(coffee_jac, jac_k(penguins_X, 3, seed = 7))
)
observed |> mutate(across(where(is.numeric), ~ round(.x, 3)))
## # A tibble: 2 × 3
##   data     silhouette jaccard
##   <chr>         <dbl>   <dbl>
## 1 Coffee        0.303   0.963
## 2 Penguins      0.447   0.965

nulls <- bind_rows(mutate(null_coffee, data = "Coffee"),
                   mutate(null_penguins, data = "Penguins")) |>
  pivot_longer(c(silhouette, jaccard), names_to = "index", values_to = "null")

obs_long <- observed |>
  pivot_longer(c(silhouette, jaccard), names_to = "index", values_to = "obs")

labels <- c(silhouette = "Silhouette width", jaccard = "Bootstrap Jaccard")

ggplot(nulls, aes(y = data)) +
  geom_point(aes(x = null), color = "grey55", alpha = 0.6, size = 1.8) +
  geom_point(data = obs_long, aes(x = obs), color = dsp_colors[2], size = 4) +
  facet_wrap(~ index, scales = "free_x", labeller = labeller(index = labels)) +
  labs(title = "Observed value (orange) against 20 null datasets (grey)",
       x = NULL, y = NULL) +
  dsp_theme +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(color = "grey78"))

This is the result worth taking away. On bootstrap Jaccard, the orange dots are buried in the grey. The coffee continuum scores 0.963, squarely inside its own null range of 0.911 to 0.977. The penguins, which genuinely contain three species, score 0.965 against a null range of 0.82 to 0.956, clearing the noise by 0.009. Real groups, an artificial cut through a gradient, and pure noise all land in the same narrow band above 0.9, and the ordering between them is meaningless.

That is not a bug in clusterboot. It is what the statistic measures. Jaccard stability asks whether the algorithm reproduces the same partition on resampled data, and a smooth elongated cloud has extremely reproducible cut points, precisely because there is nothing there to make the boundary wobble. High stability is evidence that your k-means run is deterministic, not that your clusters exist.

Silhouette does better, but only once it is calibrated. Raw, the coffee value of 0.303 and the penguin value of 0.447 are both "weak" by the usual thresholds. Against their own nulls the picture separates: coffee clears its null range (0.263 to 0.295) by only 0.008, while the penguins beat their null range (0.254 to 0.286) by 0.161, a margin roughly 19 times larger. The absolute number was uninformative. The comparison was informative.

A note on the gap statistic, which I left out of that plot because it behaves differently:

gap_k <- function(M) {
  g <- clusGap(M, FUN = kmeans, nstart = 25, K.max = 8, B = 50)
  maxSE(g$Tab[, "gap"], g$Tab[, "SE.sim"], method = "firstSEmax")
}
c(noise = gap_k(X_null), coffee = gap_k(X), penguins = gap_k(penguins_X))
##    noise   coffee penguins 
##        1        3        5

Run on the pure noise it correctly returns k = 1, so it is not useless. But it endorsed k = 3 on the coffee continuum, and on the penguins it misses the three species. Its reference distribution is a uniform box, so it responds to any departure from uniformity, skewness included, and not to clustering specifically.

What did work

If the question is "are there groups", then a better thing to test is whether the data is multimodal along the direction that carries the structure. Hartigan and Hartigan’s dip test does exactly that, and diptest::dip.test is one line.

dip_coffee   <- dip.test(prcomp(X)$x[, 1])
dip_penguins <- dip.test(prcomp(penguins_X)$x[, 1])
c(coffee = dip_coffee$p.value, penguins = dip_penguins$p.value)
##       coffee     penguins 
## 8.544241e-01 2.465773e-05

bind_rows(
  tibble(pc1 = scale(prcomp(X)$x[, 1])[, 1],
         data = "Coffee: one bump"),
  tibble(pc1 = scale(prcomp(penguins_X)$x[, 1])[, 1],
         data = "Penguins: two bumps")
) |>
  ggplot(aes(pc1)) +
  geom_density(fill = dsp_colors[1], color = NA, alpha = 0.75) +
  facet_wrap(~ data, scales = "free_y") +
  labs(title = "The check that separated them",
       x = "First principal component (scaled)", y = "Density") +
  dsp_theme

The coffee scores give p = 0.85, entirely consistent with a single bump. The penguins give p = 2.5e-05. That is the separation the other two indices failed to make, from a test that took a fraction of a second.

The honest caveat: the dip test here looks only at PC1. Groups could separate along some other direction and this would miss them, so on a real problem run it on the first few components, or on the pairwise distance distribution.

What I would do instead

Three things came out of this that I will apply from now on.

Always run a matched null. null_like() above is four lines. Generate data with your data’s covariance and no clusters, run your entire pipeline on it, and report your index next to that reference. Any index without a null is a number without a scale, and as this showed, the thresholds people quote from memory can be met by noise.

Do not treat stability as evidence of existence. Bootstrap Jaccard is a genuinely useful statistic for a different question: given that groups exist, which of my clusters are solid and which are an artifact of a few points. Reported as proof that the groups are real, it will agree with you no matter what.

Set runs when you call clusterboot. kmeansCBI defaults to runs = 1, a single random start per bootstrap replicate. On the penguins that default made the Jaccard swing between roughly 0.63 and 0.86 depending only on the seed, because k-means kept landing in local optima. With runs = 25 it is stable across seeds. If you have ever seen a stability score change when you changed the seed, this is likely why.

And when the data really is a continuum, as the coffee scores are, cutting it into three is still a perfectly reasonable thing to do. Quality tiers are useful. Just describe them as what they are, thresholds on a gradient that you chose, rather than as types you discovered. The difference matters to whoever reads the result next and assumes the groups were out there waiting.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

The R code needs three packages:

library(haven) # read SAS .xpt files
library(dplyr) # data wrangling
library(survey) # complex-survey models

Import the data

Each NHANES component is a SAS transport file (.xpt), which both languages read directly from the CDC for this post. I pull four files from the 2017 to 2018 cycle: demographics (which also holds the survey-design columns), body measures, blood pressure, and glycohemoglobin (HbA1c).

Python uses pandas.read_sas:

import pandas as pd
import numpy as np
import urllib.request
import io

base = "https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/"

def nhanes(f):
    req = urllib.request.Request(base + f + ".xpt",
                                 headers={"User-Agent": "Mozilla/5.0"})
    raw = urllib.request.urlopen(req).read()
    return pd.read_sas(io.BytesIO(raw), format="xport")

demo = nhanes("DEMO_J")
bmx  = nhanes("BMX_J")
bpx  = nhanes("BPX_J")
ghb  = nhanes("GHB_J")
demo.shape
## (9254, 46)

R uses haven::read_xpt:

base <- "https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/"
nhanes <- function(f) read_xpt(paste0(base, f, ".xpt"))

demo <- nhanes("DEMO_J")
bmx <- nhanes("BMX_J")
bpx <- nhanes("BPX_J")
ghb <- nhanes("GHB_J")
dim(demo)
## [1] 9254   46

Related posts on DataScience+:

• How to Incorporate ML.Net With Algorithmic Trading
• Building A Logistic Regression in Python, Step by Step
• Logistic Regression with Python using Titanic data

Merge and prepare the variables

Every file shares an id SEQN. I average the four blood-pressure readings, left-join the components onto demographics, rename the predictors, code sex as a label, and derive a diabetes outcome (HbA1c of 6.5% or higher). A missing HbA1c stays missing rather than counting as "no diabetes".

Python:

bpx = bpx.assign(SBP=bpx[["BPXSY1", "BPXSY2", "BPXSY3", "BPXSY4"]].mean(axis=1))

df = (demo[["SEQN", "RIAGENDR", "RIDAGEYR",
            "WTMEC2YR", "SDMVPSU", "SDMVSTRA"]]
      .merge(bmx[["SEQN", "BMXBMI"]], on="SEQN", how="left")
      .merge(bpx[["SEQN", "SBP"]],    on="SEQN", how="left")
      .merge(ghb[["SEQN", "LBXGH"]],  on="SEQN", how="left"))

df = df.rename(columns={"RIDAGEYR": "age", "BMXBMI": "bmi", "SBP": "sbp"})
df = df.assign(
    sex      = np.where(df["RIAGENDR"] == 1, "Male", "Female"),
    diabetes = np.where(df["LBXGH"] >= 6.5, 1.0,
               np.where(df["LBXGH"].notna(), 0.0, np.nan)),
)
adult = df[df["age"] >= 20]
adult[["age", "sex", "bmi", "sbp", "diabetes"]].head()
##      age     sex   bmi         sbp  diabetes
## 2   66.0  Female  31.7  200.000000       0.0
## 5   66.0  Female  23.7  142.000000       0.0
## 6   75.0  Female  38.9  118.666667       0.0
## 8   56.0    Male  21.3  101.333333       0.0
## 10  67.0    Male  23.5  104.666667       0.0

R:

bpx <- bpx |>
  mutate(
    SBP = rowMeans(across(c(BPXSY1, BPXSY2, BPXSY3, BPXSY4)), na.rm = TRUE)
  )

df <- demo |>
  select(SEQN, RIAGENDR, RIDAGEYR, WTMEC2YR, SDMVPSU, SDMVSTRA) |>
  left_join(select(bmx, SEQN, BMXBMI), by = "SEQN") |>
  left_join(select(bpx, SEQN, SBP), by = "SEQN") |>
  left_join(select(ghb, SEQN, LBXGH), by = "SEQN") |>
  rename(age = RIDAGEYR, bmi = BMXBMI, sbp = SBP) |>
  mutate(
    sex = if_else(RIAGENDR == 1, "Male", "Female"),
    diabetes = case_when(LBXGH >= 6.5 ~ 1, !is.na(LBXGH) ~ 0, TRUE ~ NA_real_)
  )
adult <- filter(df, age >= 20)
head(adult[c("age", "sex", "bmi", "sbp", "diabetes")])
## # A tibble: 6 × 5
##     age sex      bmi   sbp diabetes
##   <dbl> <chr>  <dbl> <dbl>    <dbl>
## 1    66 Female  31.7  200         0
## 2    66 Female  23.7  142         0
## 3    75 Female  38.9  119.        0
## 4    56 Male    21.3  101.        0
## 5    67 Male    23.5  105.        0
## 6    54 Female  39.9  162         1

Linear regression

I model systolic blood pressure on age, BMI, and sex. Both languages use a formula with a categorical sex (Female is the reference).

Python (statsmodels, missing="drop" drops rows with missing data):

import statsmodels.formula.api as smf

lin = smf.ols("sbp ~ age + bmi + C(sex)", data=adult, missing="drop").fit()
lin.params.round(3)
## Intercept         87.734
## C(sex)[T.Male]     1.691
## age                0.528
## bmi                0.378
## dtype: float64

R (lm drops incomplete rows by default):

lin <- lm(sbp ~ age + bmi + sex, data = adult)
round(coef(lin), 3)
## (Intercept)         age         bmi     sexMale 
##      87.734       0.528       0.378       1.691

lin.summary() in Python and summary(lin) in R print the full table with standard errors and p-values.

Both languages report about 0.53 mmHg per year of age, 0.38 per BMI point, and 1.7 higher for men.

Logistic regression

Same formula for the binary diabetes outcome. The coefficients are on the log-odds scale, so I exponentiate them to read odds ratios.

Python (smf.logit):

logit = smf.logit("diabetes ~ age + bmi + C(sex)",
                  data=adult, missing="drop").fit(disp=0)
np.exp(logit.params).round(3)
## Intercept         0.001
## C(sex)[T.Male]    1.283
## age               1.052
## bmi               1.068
## dtype: float64

R (glm with family = binomial):

logit <- glm(diabetes ~ age + bmi + sex, data = adult, family = binomial)
round(exp(coef(logit)), 3)
## (Intercept)         age         bmi     sexMale 
##       0.001       1.052       1.068       1.283

Both languages report an odds ratio of about 1.05 per year of age, 1.07 per BMI point, and 1.28 for men versus women.

Add the survey weights

The fits above treat every row equally. NHANES is a complex survey: it oversamples some groups and samples in clusters, so any estimate meant for the US population needs the sampling weight (WTMEC2YR), the stratum (SDMVSTRA), and the cluster (SDMVPSU). In both languages I describe the design once, then reuse it, and I build it on the full data and restrict to adults afterward so the design stays intact for correct standard errors.

Python uses the svy package (built on polars). where= restricts to the adult domain, and family is "gaussian" for the linear model and "binomial" for the logistic one:

import polars as pl
import svy

design = svy.Design(stratum="SDMVSTRA", psu="SDMVPSU", wgt="WTMEC2YR")
sample = svy.Sample(pl.from_pandas(df), design)

lin_w = sample.glm.fit(y="sbp", x=["age", "bmi", svy.Cat("sex", ref="Female")],
                       family="gaussian", where=svy.col("age") >= 20)
{c.term: round(float(c.est), 3) for c in lin_w.coefs}
## {'_intercept_': 87.602, 'age': 0.458, 'bmi': 0.418, 'sex_Male': 3.1}

log_w = sample.glm.fit(y="diabetes", x=["age", "bmi", svy.Cat("sex", ref="Female")],
                       family="binomial", where=svy.col("age") >= 20)
{c.term: round(float(np.exp(c.est)), 3) for c in log_w.coefs}
## {'_intercept_': 0.0, 'age': 1.057, 'bmi': 1.08, 'sex_Male': 1.552}

R uses the survey package. svydesign sets up the design, subset restricts to adults, and svyglm fits the models (quasibinomial for logistic to avoid a harmless weights warning):

options(survey.lonely.psu = "adjust")
des <- svydesign(
  ids = ~SDMVPSU,
  strata = ~SDMVSTRA,
  weights = ~WTMEC2YR,
  nest = TRUE,
  data = df
)
des_ad <- subset(des, age >= 20)

lin_w <- svyglm(sbp ~ age + bmi + sex, design = des_ad)
round(coef(lin_w), 3)
## (Intercept)         age         bmi     sexMale 
##      87.602       0.458       0.418       3.100

log_w <- svyglm(
  diabetes ~ age + bmi + sex,
  design = des_ad,
  family = quasibinomial()
)
round(exp(coef(log_w)), 3)
## (Intercept)         age         bmi     sexMale 
##       0.000       1.057       1.080       1.552

The weighted estimates are population estimates with design-based standard errors, and the two toolkits agree: the male-versus-female blood-pressure gap is about 3.1 mmHg and the diabetes odds ratio for men is about 1.55, in both Python and R.

That’s it. I hope you find this tutorial useful, and if you have questions, leave a comment below.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

Let’s answer all three on real data. First I replace glm(), predict(), and tidy() with their tidymodels equivalents one for one and see what is different (mostly the output, not the math). Then I show what tidymodels adds around the model, the part base R has no direct function for. Finally I weigh whether it is worth using for everyday epidemiological research.

The data: predicting diabetes in NHANES

I use NHANES , the US National Health and Nutrition Examination Survey, which ships as an R package. NHANES is a complex survey with sampling weights, and any valid population estimate has to carry those weights (using the survey  package and survey::svyglm() in place of glm()). For this tutorial, I keep the example simple and unweighted.

library(tidymodels)
library(NHANES)
tidymodels_prefer() # solves naming conflicts in R by making tidymodels functions take priority over other packages

data(NHANES)
nh <- NHANES |>
  filter(Age >= 18) |>
  select(
    Diabetes,
    Age,
    BMI,
    Gender,
    PhysActive,
    BPSysAve,
    TotChol,
    DirectChol,
    Poverty,
    Smoke100
  ) |>
  filter(!is.na(Diabetes)) |>
  mutate(Diabetes = factor(Diabetes, levels = c("No", "Yes"))) |>
  distinct() # NHANES resampled some individuals; drop exact duplicate rows

head(nh)
## # A tibble: 6 × 10
##   Diabetes   Age   BMI Gender PhysActive BPSysAve TotChol DirectChol Poverty Smoke100
##   <fct>    <int> <dbl> <fct>  <fct>         <int>   <dbl>      <dbl>   <dbl> <fct>   
## 1 No          34  32.2 male   No              113    3.49       1.29    1.36 Yes     
## 2 No          49  30.6 female No              112    6.7        1.16    1.91 Yes     
## 3 No          45  27.2 female Yes             118    5.82       2.12    5    No      
## 4 No          66  23.7 male   Yes             111    4.99       0.67    2.2  Yes     
## 5 No          58  23.7 male   Yes             104    4.24       0.96    5    No      
## 6 No          54  26.0 male   Yes             134    6.41       1.16    2.2  No

That leaves 4,834 distinct adults, 11.3% with diagnosed diabetes. The distinct() matters more than it looks: the NHANES package is a teaching sample that resampled some individuals, so the raw extract carries many exact duplicate rows. Leave them in and copies of one person would land on both sides of the train/test split I make later, quietly flattering any out-of-sample estimate, so I keep one row per person.

Note the factor(..., levels = c("No", "Yes")): it is not cosmetic. glm() models the probability of the last level, so with No first I am modeling P(diabetes), which is what I want. Reverse it and every odds ratio silently inverts.

Related posts on DataScience+:

• Exploring Vitamin D deficiency in the United States: NHANES 2001-2010
• Weight loss in the U.S. – An analysis of NHANES data with tidyverse
• How to Perform a Logistic Regression in R

Replacing the base functions, one for one

glm() and broom, and their parsnip equivalent

Here is the base-R fit, tidied into odds ratios.

fit_glm <- glm(
  Diabetes ~ Age +
    BMI +
    Gender +
    PhysActive +
    BPSysAve +
    TotChol +
    DirectChol +
    Poverty +
    Smoke100,
  data = nh,
  family = binomial
)

tidy(fit_glm, exponentiate = TRUE, conf.int = TRUE) |>
  select(term, OR = estimate, conf.low, conf.high, p.value)
## # A tibble: 10 × 5
##    term               OR conf.low conf.high  p.value
##    <chr>           <dbl>    <dbl>     <dbl>    <dbl>
##  1 (Intercept)   0.00218 0.000670   0.00693 7.77e-25
##  2 Age           1.06    1.05       1.06    1.85e-46
##  3 BMI           1.08    1.06       1.09    7.43e-22
##  4 Gendermale    1.14    0.901      1.45    2.68e- 1
##  5 PhysActiveYes 0.965   0.766      1.22    7.61e- 1
##  6 BPSysAve      1.01    1.00       1.02    2.52e- 3
##  7 TotChol       0.780   0.700      0.868   6.34e- 6
##  8 DirectChol    0.531   0.376      0.743   2.70e- 4
##  9 Poverty       0.880   0.820      0.945   4.26e- 4
## 10 Smoke100Yes   1.21    0.975      1.51    8.32e- 2

Now the same model through parsnip, the tidymodels interface: pick a model type, set an engine, fit().

fit_psn <- logistic_reg() |>
  set_engine("glm") |>
  fit(
    Diabetes ~ Age +
      BMI +
      Gender +
      PhysActive +
      BPSysAve +
      TotChol +
      DirectChol +
      Poverty +
      Smoke100,
    data = nh
  )

tidy(fit_psn, exponentiate = TRUE, conf.int = TRUE) |>
  select(term, OR = estimate, conf.low, conf.high, p.value)
## # A tibble: 10 × 5
##    term               OR conf.low conf.high  p.value
##    <chr>           <dbl>    <dbl>     <dbl>    <dbl>
##  1 (Intercept)   0.00218 0.000670   0.00693 7.77e-25
##  2 Age           1.06    1.05       1.06    1.85e-46
##  3 BMI           1.08    1.06       1.09    7.43e-22
##  4 Gendermale    1.14    0.901      1.45    2.68e- 1
##  5 PhysActiveYes 0.965   0.766      1.22    7.61e- 1
##  6 BPSysAve      1.01    1.00       1.02    2.52e- 3
##  7 TotChol       0.780   0.700      0.868   6.34e- 6
##  8 DirectChol    0.531   0.376      0.743   2.70e- 4
##  9 Poverty       0.880   0.820      0.945   4.26e- 4
## 10 Smoke100Yes   1.21    0.975      1.51    8.32e- 2

Identical coefficients, because parsnip called glm() for you and tidy() is the same broom. So far tidymodels has changed nothing.

predict(): the one replacement that actually behaves differently

Prediction is where the swap is not cosmetic. Let’s score three people with each.

newpeople <- nh[c(10, 200, 900), ]
newpeople
## # A tibble: 3 × 10
##   Diabetes   Age   BMI Gender PhysActive BPSysAve TotChol DirectChol Poverty Smoke100
##   <fct>    <int> <dbl> <fct>  <fct>         <int>   <dbl>      <dbl>   <dbl> <fct>   
## 1 No          60  25.8 male   No              152    6.39       1.34    1.03 Yes     
## 2 No          22  52.1 male   No              128    3.67       0.88    0.98 No      
## 3 No          45  22.0 female No              114    4.78       1.55    2.03 No

# you must remember type = "response" to get a probability
predict(fit_glm, newpeople, type = "response")
##          1          2          3 
## 0.14843669 0.21032214 0.03234337

With type = "response", predict() returns the three people as probabilities between 0 and 1, the modeled diabetes risk for each. Omit the argument and you get log-odds instead, which is the default trap the comment warns about.

# TIDYMODELS: always a tibble, one row per input row, standardized names
predict(fit_psn, newpeople, type = "prob") # .pred_No / .pred_Yes
## # A tibble: 3 × 2
##   .pred_No .pred_Yes
##      <dbl>     <dbl>
## 1    0.852    0.148 
## 2    0.790    0.210 
## 3    0.968    0.0323

tidymodels returns the same risks, but as a tibble with one column per outcome level: .pred_No and .pred_Yes, which sum to 1 in each row. There is no scale to remember, type = "prob" always means probabilities, and the column names say which level each one belongs to. .pred_Yes is the diabetes risk, the same number the type = "response" vector gave above.

augment() goes one step further and binds the predictions onto the data:

augment(fit_psn, newpeople) |>
  select(Diabetes, Age, BMI, .pred_class, .pred_Yes)
## # A tibble: 3 × 5
##   Diabetes   Age   BMI .pred_class .pred_Yes
##   <fct>    <int> <dbl> <fct>           <dbl>
## 1 No          60  25.8 No             0.148 
## 2 No          22  52.1 No             0.210 
## 3 No          45  22.0 No             0.0323

So the whole base-R workflow has a direct translation:

Conclusion: for a single fit, tidymodels computes exactly what glm(), predict(), and broom already give you. The output is tidier and the prediction defaults are safer, but the math is the same, so the swap on its own is not a reason to switch.

What tidymodels adds

The framework is not really about replacing glm(); it is about the workflow around the model.

Preprocessing that travels with predict()

In base R, if you scale or impute your predictors, you have to redo that transformation, identically, on every new dataset you predict on, using the training statistics. Forget to, or use the new data’s own mean instead, and you get wrong or leaky predictions. A recipe attaches the preprocessing to the model so predict() applies it automatically.

set.seed(2026)
split <- initial_split(nh, prop = 0.75, strata = Diabetes)
train <- training(split)
test <- testing(split)

rec <- recipe(Diabetes ~ ., data = train) |>
  step_impute_median(all_numeric_predictors()) |> # learn medians on train
  step_impute_mode(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors()) # center/scale on train stats

wf_fit <- workflow() |>
  add_recipe(rec) |>
  add_model(logistic_reg() |> set_engine("glm")) |>
  fit(train)

The payoff shows up when a new row is incomplete. Here is a test patient whose BMI is missing:

patient <- test[1, ]
patient$BMI <- NA_real_

# base glm cannot score a row with a gap
predict(fit_glm, patient, type = "response")
##  1 
## NA

# the workflow imputes BMI (with the training median) before predicting
predict(wf_fit, patient, type = "prob")
## # A tibble: 1 × 2
##   .pred_No .pred_Yes
##      <dbl>     <dbl>
## 1    0.945    0.0551

The base model returns NA; the workflow silently imputes the missing BMI with the median it learned from the training data and returns a real risk.

Model performance

The question a risk model has to answer is how it does on data it has never seen, measured the way risk models are judged: discrimination (ROC AUC) and calibration. Fitting was the easy part; the work is estimating those two honestly, on data the fit did not touch. That is the piece base R leaves you to assemble, and where tidymodels earns its keep.

vfold_cv() splits the training data into 10 folds; the model is fit on 9 of them and scored on the tenth, ten times over, so every performance number comes from data the fit never saw. Doing this by hand means writing the fold loop and re-learning the recipe inside each fold, so the imputation and scaling use only that fold’s training portion. Get that last part wrong (impute once on the whole set) and the held-out fold has quietly seen the training data through a shared median.

set.seed(2026)
folds <- vfold_cv(train, v = 10, strata = Diabetes)

log_spec <- logistic_reg() |> set_engine("glm")
mset <- metric_set(roc_auc, brier_class) # discrimination + calibration

res_log <- fit_resamples(
  workflow() |> add_recipe(rec) |> add_model(log_spec),
  folds,
  metrics = mset
)

fit_resamples() fits the workflow on every fold, re-learns the recipe inside each one, and collects the two metrics without ever touching the test set.

collect_metrics(res_log) |>
  select(.metric, cv_estimate = mean, std_err)
## # A tibble: 2 × 3
##   .metric     cv_estimate std_err
##   <chr>             <dbl>   <dbl>
## 1 brier_class      0.0860 0.00172
## 2 roc_auc          0.816  0.0105

Cross-validated, the logistic model reaches an AUC of 0.816, and the number is trustworthy precisely because no fold was scored on data it trained on. You could get the same value in base R, but you would own the whole fold loop and the per-fold preprocessing described above. tidymodels does that bookkeeping for you: the recipe travels into every resample automatically, so the shortcut that leaks is not even on the table.

Discrimination and calibration plots

lf_log <- last_fit(
  workflow() |> add_recipe(rec) |> add_model(log_spec),
  split,
  metrics = mset
)
preds <- collect_predictions(lf_log)

dsp_colors <- c(
  "#0066CC",
  "#E8862D",
  "#159A6C",
  "#7D5BD6",
  "#D64580",
  "#2AA9B8",
  "#C9A227"
)
dsp_theme <- theme_minimal(base_size = 13) +
  theme(
    plot.background = element_rect(fill = "#ECECEF", color = NA),
    panel.background = element_rect(fill = "#ECECEF", color = NA),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "grey78"),
    axis.ticks = element_blank(),
    plot.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold")
  )

preds |>
  roc_curve(truth = Diabetes, .pred_Yes, event_level = "second") |>
  ggplot(aes(1 - specificity, sensitivity)) +
  geom_abline(linetype = "dashed", color = "grey60") +
  geom_path(linewidth = 1.1, color = dsp_colors[1]) +
  coord_equal() +
  labs(
    title = "Discrimination: who is at risk of diabetes?",
    subtitle = "ROC curve on held-out test data",
    x = "False positive rate",
    y = "True positive rate"
  ) +
  dsp_theme

On the held-out test set the model’s AUC is 0.816, close to the cross-validated estimate, which is the reassurance you want: it did not fall apart on data it had never seen. But discrimination is only half the story. A risk model also has to be calibrated: when it predicts a 30% risk, close to 30% of the people carrying that prediction should actually have diabetes. The next plot checks that by binning people into risk deciles and comparing predicted against observed rates.

preds |>
  mutate(bin = ntile(.pred_Yes, 10)) |>
  group_by(bin) |>
  summarise(
    pred = mean(.pred_Yes),
    obs = mean(Diabetes == "Yes"),
    .groups = "drop"
  ) |>
  ggplot(aes(pred, obs)) +
  geom_abline(linetype = "dashed", color = "grey60") +
  geom_line(linewidth = 0.9, color = dsp_colors[1]) +
  geom_point(size = 2.4, color = dsp_colors[1]) +
  labs(
    title = "Calibration: are the predicted risks honest?",
    subtitle = "Predicted vs observed rate by risk decile (dashed line = perfect)",
    x = "Mean predicted risk",
    y = "Observed diabetes rate"
  ) +
  dsp_theme

The points track the diagonal closely through the low and middle deciles and wobble only at the high-risk end, where the bins are smallest and noisiest and predicted risk climbs toward 40%. The Brier score, 0.087 (lower is better), is a single number that folds discrimination and calibration together.

Is any of this useful in epidemiological and clinical research?

It depends on the aim:

Association and etiologic papers report adjusted odds or hazard ratios with confidence intervals, usually from a survey-weighted or matched design. For that, base glm() or survey::svyglm() with broom::tidy() is exactly right and tidymodels adds almost nothing.

Prediction and prognostic papers are the other case: developing or validating a clinical risk score, the kind of work that reporting guidelines like TRIPOD govern. There the deliverable is discrimination and calibration on unseen data, with internal validation by resampling. That is what tidymodels streamlines.

So it is not glm() versus tidymodels, and broom lives inside both. Replace your base functions with tidymodels when you want the workflow it enables, not for the sake of the swap: use glm() or svyglm() and broom to explain associations, and reach for the full framework when your focus is a risk model you have to validate.

That’s it. I hope you find this tutorial useful, and if you have questions, leave a comment below.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.

I use a handful of packages, all on CRAN. tidycensus fetches ACS estimates (and, handily, their map geometry in one call); the rest are for wrangling and plotting.

library(tidycensus)
library(dplyr)
library(stringr)
library(forcats)
library(ggplot2)
library(sf)
library(ggrepel)
library(patchwork)

tidycensus needs a free Census API key. Request one at api.census.gov/data/key_signup.html, then run census_api_key("YOUR_KEY", install = TRUE) once and restart R.

Finding the right variable

The ACS reports commute times across dozens of tables. There is no single "mean commute" variable, but there is an easy way to build one. Table B08013 gives the aggregate travel time to work (every commuter’s minutes added together) for an area, and B08012 gives the number of workers who commute. Divide one by the other and you get the mean. Every ACS table has an _001 row for the total, which is all I need here.

I pull both for every county in the country at once. Leaving out the state argument returns all ~3,200 counties. I use the wide output so each variable becomes its own column.

cty <- get_acs(
  geography = "county",
  variables = c(agg = "B08013_001", workers = "B08012_001"),
  year = 2022, survey = "acs5", output = "wide"
) |>
  mutate(mean_commute = aggE / workersE)

nat_mean <- sum(cty$aggE, na.rm = TRUE) / sum(cty$workersE, na.rm = TRUE)
round(nat_mean, 1)
## [1] 26.7

Averaged across every commuter in the country, the mean trip is 26.7 minutes each way. Before mapping anything, it is worth looking at the extremes to check the numbers make sense.

cty |>
  filter(workersE > 50000) |>
  arrange(desc(mean_commute)) |>
  transmute(NAME, mean_commute = round(mean_commute, 1)) |>
  head(6)
## # A tibble: 6 × 2
##   NAME                        mean_commute
##   <chr>                              <dbl>
## 1 Bronx County, New York              44.6
## 2 Charles County, Maryland            44.2
## 3 Richmond County, New York           43.8
## 4 Queens County, New York             43.6
## 5 Kings County, New York              42.7
## 6 Monroe County, Pennsylvania         39.8

The longest commutes among populous counties are the outer boroughs of New York and the exurban counties feeding Washington and New York. That is exactly what you would expect, so the pipeline is working.

Related posts on DataScience+:

• Building Heatmaps in R with ggplot2 package
• Visualising Thefts using Heatmaps in ggplot2
• How Happy is Your Country? — Happy Planet Index Visualized

The national map

To map this I ask tidycensus for the same data with geometry by adding geometry = TRUE. That returns an sf object I can hand straight to ggplot2. The one extra touch is tigris::shift_geometry(), which relocates Alaska, Hawaii, and Puerto Rico under the lower 48 so the whole country fits one frame.

I reuse a small editorial theme across every figure in this post. It tints the entire canvas a light gray and strips the chart down to the data. For a map I also drop the axes.

dsp_theme <- theme_minimal(base_size = 13) +
  theme(plot.background    = element_rect(fill = "#ECECEF", color = NA),
        panel.background   = element_rect(fill = "#ECECEF", color = NA),
        panel.grid.minor   = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color = "grey78"),
        axis.ticks         = element_blank(),
        plot.title         = element_text(face = "bold"),
        plot.subtitle      = element_text(color = "grey30"))

map_theme <- dsp_theme +
  theme(panel.grid.major = element_blank(),
        axis.text = element_blank(), axis.title = element_blank(),
        legend.position = "bottom")

Now the map itself. The data comes down with geometry, I shift it, and geom_sf() fills each county by its mean commute. County outlines would only clutter 3,200 shapes, so I draw them borderless and overlay a single thin white state layer to give the eye something to orient by. I pull state boundaries from tigris and shift them the same way, so they line up with the counties. I also drop labels on a few large metros so the geography is easy to read, using ggrepel to nudge the text clear of the coastline.

cty_geo <- get_acs(
  geography = "county",
  variables = c(agg = "B08013_001", workers = "B08012_001"),
  year = 2022, survey = "acs5", output = "wide", geometry = TRUE
) |>
  mutate(mean_commute = aggE / workersE) |>
  tigris::shift_geometry()

states <- tigris::states(cb = TRUE, year = 2022) |>
  filter(!STUSPS %in% c("VI", "GU", "MP", "AS")) |>
  tigris::shift_geometry()

# a few metros to anchor the eye; transformed into the shifted map's CRS
cities <- tibble::tribble(
  ~lab,            ~lon,     ~lat,
  "New York",      -74.00,   40.71,
  "San Francisco", -122.42,  37.77,
  "Washington",    -77.04,   38.90,
  "Boston",        -71.06,   42.36,
  "Chicago",       -87.63,   41.88
) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326) |>
  st_transform(st_crs(cty_geo))

city_labels <- geom_text_repel(
  data = cities, aes(label = lab, geometry = geometry),
  stat = "sf_coordinates", size = 3.2, fontface = "bold",
  color = "grey10", bg.color = "#ECECEF", bg.r = 0.15,
  min.segment.length = 0, segment.color = "grey30",
  segment.size = 0.25, box.padding = 0.5, seed = 1
)

ggplot(cty_geo) +
  geom_sf(aes(fill = mean_commute), color = NA) +
  geom_sf(data = states, fill = NA, color = "white", linewidth = 0.2) +
  city_labels +
  scale_fill_viridis_c(option = "magma", direction = -1,
                       name = "Mean commute (min)  ",
                       breaks = c(15, 25, 35, 45)) +
  labs(title = "How long is the drive to work?",
       subtitle = "Mean travel time to work by county, ACS 2018-2022",
       caption = "Data: US Census Bureau, ACS 5-year") +
  guides(fill = guide_colorbar(barwidth = 12, barheight = 0.5,
                               title.position = "top", title.hjust = 0.5)) +
  map_theme

The geography tells a clear story. Commutes are longest in the dark bands of the East, the Appalachians and the Deep South, and around every major metro. They are shortest across the light expanse of the Great Plains and rural Alaska, where towns are small and the drive across one is a matter of minutes. Long commutes are a feature of density and distance-to-a-big-city, not of open country.

Which metros have it worst

The county map hints at metros but does not rank them. For that I switch the geography to "cbsa", the Census term for metropolitan and micropolitan areas, and pull population too so I can keep only the large metros where the number is stable.

metro <- get_acs(
  geography = "cbsa",
  variables = c(agg = "B08013_001", workers = "B08012_001", pop = "B01003_001"),
  year = 2022, survey = "acs5", output = "wide"
) |>
  filter(str_detect(NAME, "Metro Area"), popE > 500000) |>
  mutate(mean_commute = aggE / workersE,
         label = paste0(str_extract(NAME, "^[^-,]+"), ", ",
                        str_remove(str_extract(NAME, ", [A-Z]{2}"), ", ")))

A horizontal bar chart is the natural way to compare the top of the list. I sort with fct_reorder() so the bars line up longest-first and add the value at the end of each.

metro |>
  slice_max(mean_commute, n = 15) |>
  ggplot(aes(mean_commute, fct_reorder(label, mean_commute))) +
  geom_col(fill = "#0066CC", width = 0.72) +
  geom_text(aes(label = sprintf("%.1f", mean_commute)), hjust = -0.15, size = 3.6) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.10))) +
  labs(title = "The metros with the longest commutes",
       subtitle = "Mean travel time to work, metros over 500k people, ACS 2018-2022",
       x = NULL, y = NULL, caption = "Data: US Census Bureau, ACS 5-year") +
  dsp_theme +
  theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank())

New York tops the list, no surprise. What is more interesting is the company it keeps: Stockton, Riverside, and Modesto are not big job centers themselves, they are the affordable edges of the Bay Area and Los Angeles, where people move for cheaper housing and pay for it in driving time. The longest commutes are half about big cities and half about the sprawl around them.

The result I did not expect: by mode

Here is where the data surprised me. I assumed transit riders, gliding past traffic, would have the shortest commutes. The opposite is true. Table B08136 breaks aggregate travel time down by how people get to work, and B08301 gives the matching worker counts, so the same divide-and-average trick gives a mean commute per mode. The mode categories are split across several rows, so I pull the whole tables and pick the pieces I need.

agg <- get_acs("us", table = "B08136", year = 2022, survey = "acs5")
cnt <- get_acs("us", table = "B08301", year = 2022, survey = "acs5")

A <- function(code) agg$estimate[agg$variable == code]
N <- function(...) sum(sapply(c(...), function(x) cnt$estimate[cnt$variable == x]))

modes <- data.frame(
  mode = c("Drove alone", "Carpooled", "Bus", "Subway / rail",
           "Commuter rail / ferry", "Walked"),
  minutes = c(A("B08136_003") / N("B08301_003"),
              A("B08136_004") / N("B08301_004"),
              A("B08136_008") / N("B08301_011"),
              A("B08136_009") / N("B08301_012", "B08301_014"),
              A("B08136_010") / N("B08301_013", "B08301_015"),
              A("B08136_011") / N("B08301_019"))
)

ggplot(modes, aes(minutes, fct_reorder(mode, minutes))) +
  geom_col(fill = "#0066CC", width = 0.68) +
  geom_text(aes(label = sprintf("%.0f min", minutes)), hjust = -0.15, size = 4) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.13))) +
  labs(title = "Transit takes far longer than driving",
       subtitle = "Mean travel time to work by commute mode, US, ACS 2018-2022",
       x = NULL, y = NULL, caption = "Data: US Census Bureau, ACS 5-year") +
  dsp_theme +
  theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank())

Someone who drives alone averages about 26 minutes. Someone on commuter rail or a ferry averages 71, nearly three times as long. This is not because trains are slow. It is a selection effect: transit is worth taking mostly for long trips into dense downtowns, so the people who use it are the ones with the longest journeys to begin with, and their door-to-door time includes walking, waiting, and transfers. The mode does not cause the long commute; the long commute is why they chose the mode. It is a good reminder that a difference in a group average is rarely a clean cause-and-effect.

A map makes the selection effect obvious. If I map the share of workers who commute by public transportation, table B08301 again, I can see where those long transit trips actually happen. Transit share is tiny in almost every county, so I put the color scale on a square-root transform to keep the busy metros from washing everything else out.

ts <- get_acs(
  geography = "county",
  variables = c(transit = "B08301_010", total = "B08301_001"),
  year = 2022, survey = "acs5", output = "wide", geometry = TRUE
) |>
  mutate(transit_share = 100 * transitE / totalE) |>
  tigris::shift_geometry()

ggplot(ts) +
  geom_sf(aes(fill = transit_share), color = NA) +
  geom_sf(data = states, fill = NA, color = "grey35", linewidth = 0.15) +
  city_labels +
  scale_fill_viridis_c(option = "mako", direction = -1, transform = "sqrt",
                       breaks = c(1, 5, 15, 30),
                       name = "% commuting by transit  ") +
  labs(title = "Transit is a big-metro phenomenon",
       subtitle = "Share of workers commuting by public transportation, ACS 2018-2022",
       caption = "Data: US Census Bureau, ACS 5-year") +
  guides(fill = guide_colorbar(barwidth = 12, barheight = 0.5,
                               title.position = "top", title.hjust = 0.5)) +
  map_theme

The country is almost empty. Transit collapses to a handful of dark islands: New York above all, then the Bay Area, Washington, Boston, and Chicago. These are the same metros that topped the commute-time chart, and that is the whole point. The long transit average is not a national fact about trains; it is a few very large, very dense metros where a long ride beats an even longer drive, showing up in the one national number.

Inside a single metro

A metro-wide average hides as much as it shows. To see the texture I drop to the census-tract level, the finest geography the ACS maps well, for two metros with very different shapes. tidycensus takes a vector of counties, and I keep only tracts with enough commuters to be reliable.

get_metro <- function(state, counties, name) {
  get_acs("tract", state = state, county = counties,
          variables = c(agg = "B08013_001", workers = "B08012_001"),
          year = 2022, survey = "acs5", output = "wide", geometry = TRUE) |>
    mutate(mean_commute = aggE / workersE, metro = name) |>
    filter(workersE > 200)
}

nyc <- get_metro("NY", c("New York", "Kings", "Queens", "Bronx", "Richmond"),
                 "New York City")
atl <- get_metro("GA", c("Fulton", "DeKalb", "Cobb", "Gwinnett", "Clayton"),
                 "Atlanta")

I draw both on a shared color scale with patchwork so the two panels are directly comparable, then collect the single legend.

rng <- range(c(nyc$mean_commute, atl$mean_commute), na.rm = TRUE)

one_map <- function(d, ttl) {
  ggplot(d) +
    geom_sf(aes(fill = mean_commute), color = NA) +
    scale_fill_viridis_c(option = "magma", direction = -1, limits = rng, name = "min") +
    labs(title = ttl) +
    theme_void(base_size = 12) +
    theme(plot.background = element_rect(fill = "#ECECEF", color = NA),
          plot.title = element_text(face = "bold", hjust = 0.02))
}

(one_map(nyc, "New York City") | one_map(atl, "Atlanta")) +
  plot_layout(guides = "collect") +
  plot_annotation(
    title = "Same average, two different shapes",
    subtitle = "Mean travel time to work by census tract, ACS 2018-2022",
    caption = "Data: US Census Bureau, ACS 5-year",
    theme = theme(plot.background = element_rect(fill = "#ECECEF", color = NA),
                  plot.title = element_text(face = "bold", size = 16),
                  plot.subtitle = element_text(color = "grey30")))

The two cities could not be more different. New York is dark almost everywhere, a median tract commute of 42.7 minutes, and the longest trips are in the outer boroughs, not the core: people commute into Manhattan, not out of it. Atlanta shows the classic ring: a bright, short-commute core in Midtown and Buckhead fading to dark exurban edges, with a tract median of just 30 minutes. Same country, same survey, two completely different geographies of getting to work.

One more number worth knowing

The averages are one thing, but the tail is where the human cost sits. Table B08303 bins commuters by how long they travel, so I can count the "super-commuters" who spend 60 minutes or more each way.

b <- get_acs("us", table = "B08303", year = 2022, survey = "acs5")
total <- b$estimate[b$variable == "B08303_001"]
over60 <- sum(b$estimate[b$variable %in% c("B08303_012", "B08303_013")])
round(100 * over60 / total, 1)
## [1] 8.9

About 8.9 percent of American commuters, roughly 12.3 million people, spend at least an hour getting to work, one way. That is two-plus hours of every working day spent in transit.

Make it your own

Everything here comes from three or four ACS variables and the divide-an-aggregate-by-a-count trick. Swap the table codes and you can map median rent, income, insurance coverage, or the share of people who work from home, at any geography from the whole nation down to a single tract. Change the state and county arguments and the tract maps redraw for your own city. The ACS has thousands of variables; load_variables(2022, "acs5") lets you search them all.

That’s it. I hope you find this tutorial useful, and if you have questions, leave a comment below.

This article was first published on DataScience+, a community of R and Python tutorial authors. Have a data-science technique worth sharing? Write for us — no pitch required.