the-stats-duck v0.7 — a regression in every GROUP BY

the-stats-duck — our open-source DuckDB extension for statistics — just shipped v0.7.0, and it's now a one-line install from the DuckDB community extensions (see below).

v0.6 taught DuckDB to fit a regression: lm took an R-style formula and gave you back a coefficient table for the whole query. Good for one model. But clinical trial work is rarely just one model — it's usually one model per site, per arm, per cohort — and the standard errors you report had better survive a reviewer who knows the residuals aren't textbook-tidy.

So v0.7's headline is lm_fit: ordinary least squares as a DuckDB aggregate. One regression per GROUP BY group, with robust and cluster-robust standard errors on request.

One regression per group

Because lm_fit is an aggregate, GROUP BY does the fanning-out for free. Here it fits body_mass_g ~ flipper_length_mm + bill_length_mm separately for each penguin species, in a single query — live, in your browser:

Three species, three fitted models, one pass over the data. Change GROUP BY species to GROUP BY island and you get a model per island instead — the aggregate doesn't care.

The shape of it

lm_fit takes the response and the predictor row, and returns a STRUCT:

lm_fit(y, x [, vcov [, cluster [, add_intercept]]])

y — the response value (one per row).
x — the predictor row as a LIST(DOUBLE), e.g. [flipper_length_mm, bill_length_mm]. An intercept is prepended automatically unless add_intercept is set to false.
The result STRUCT carries a coefficients list — one {term, estimate, std_error, t_statistic, p_value} per predictor — alongside n, k, df_residual, r_squared, adj_r_squared, sigma, f_statistic, f_p_value, and the vcov_type actually used.

Predictors come back by position — (Intercept), then x1, x2, … in list order — so you unnest the coefficients to get a tidy row per term. On a tiny dataset that looks like this:

SELECT (u).term, (u).estimate, (u).std_error
FROM (
  SELECT unnest((lm_fit(y, [x1])).coefficients) AS u
  FROM (VALUES
      (2.1, 1.0), (3.9, 2.0), (6.2, 3.0), (7.8, 4.0),
      (10.1, 5.0), (12.2, 6.0), (13.8, 7.0), (16.1, 8.0)
  ) AS t(y, x1)
);
--   term         estimate  std_error
--   (Intercept)  0.0357    0.1404
--   x1           1.9976    0.0278

lm_fit is the aggregate, positional sibling of v0.6's formula-driven lm/lm_summary. Reach for lm when you want a formula over one table; reach for lm_fit when you want a model per group, robust standard errors, or both.

Standard errors that survive review

The point estimates are just OLS. What changes with the vcov argument is the standard errors — and with them the t-statistics and p-values. That's usually the number a methodologist actually argues about.

-- Heteroskedasticity-consistent (HC3) — the conservative small-sample default
SELECT (u).term, (u).std_error
FROM (SELECT unnest((lm_fit(y, [x1, x2], 'HC3')).coefficients) AS u FROM mytable);
--   term         std_error
--   (Intercept)  3.2802
--   x1           0.4749
--   x2           0.3839

The full menu:

'const' — classical OLS (the default).
'HC0'–'HC3' — heteroskedasticity-consistent. 'HC1' is Stata's , robust; 'HC3' is the conservative choice when n is small.
'CR0' / 'CR1' — cluster-robust, for data with repeated measurement inside a subject, site, or center. Pass a per-row cluster key ('cluster' is accepted as an alias for 'CR1', the Stata / statsmodels default), and inference uses a t(G − 1) reference with G clusters:

SELECT (u).term, (u).estimate, (u).std_error
FROM (
  SELECT unnest(
    (lm_fit(response, [dose], 'CR1', subject_id::VARCHAR)).coefficients
  ) AS u
  FROM observations
);

Only the covariance estimator moves — the coefficients are identical across all of them. Every estimator is validated against statsmodels' cov_type='HC*' and 'cluster'.

One nicety for the per-group case: a degenerate group — fewer rows than predictors, a singular design, or fewer than two clusters — returns NULL for that group instead of aborting the whole query. You get every model that can be fit and NULLs where one can't, which is exactly what you want when you've just fanned out across 200 sites.

bootstrap(), now identical on every machine

A quieter but important fix: a seeded bootstrap() now returns the exact same resamples on every operating system. Before v0.7, the same seed could draw a different set of resamples on Windows, macOS, and Linux — so a bootstrap result you shared with a colleague wouldn't reproduce number-for-number on their machine. The estimates were always statistically sound; they just weren't bit-for-bit identical across platforms.

Now they are: give bootstrap() a seed and anyone who reruns it, on any OS, gets the same resamples and the same interval. If you like this kind of bug, we wrote up two related RNG traps from the sampler work, and the full reasoning behind this one is in the engineering note.

Under the hood: a shared numerical core

lm_fit is built on a new, well-tested numerical core — the matrix machinery that regression and robust standard errors actually run on (the least-squares solve, the matrix inverses, and the "sandwich" that turns residuals into robust standard errors). It's deliberately self-contained and reusable, so the same carefully-checked code can back more of the extension over time: lm_fit is the first to use it, and the existing lm / lm_summary move onto it next. More of the statistics now stands on one foundation, validated the same way.

Go play

the-stats-duck is MIT-licensed and now lives in the DuckDB community extensions — so installing it is one line in any DuckDB session, no building from source:

INSTALL stats_duck FROM community;
LOAD stats_duck;

It's the statistics engine under Bedevere and KoliLang; the code and full v0.7.0 release notes live on GitHub. Drop it into any DuckDB session and your GROUP BY grows a regression department.

And if there's an estimator or test you need that isn't there yet, get in touch — commissioned functions ship in the open-source extension, so everyone gets them.

Not dead. What? Doing your regressions.