Regression conventions

Part of the numerical-validation suite. The recorded conventions behind the regression and GLM numbers.

GLM quantile residuals — the per-family definitions

Decided 2026-08-09, before any code, as step 1a of docs/releases/release_0.16.0.md item 1 (kept private). The diagnostic panel for linear regression shipped in 0.12.0; extending it to the GLM families is new statistics, not plumbing, because "residual of what?" has a different answer per family. These are those answers.

The residual: randomized quantile residuals, every family

For an observation y with fitted distribution F(· ; θ̂), the quantile residual (Dunn & Smyth 1996) is

r = Φ⁻¹(u),    u ~ Uniform( F(y⁻ ; θ̂), F(y ; θ̂) ]

where F(y⁻) is the left limit — F(y−1) for an integer-valued outcome, and F(y) itself for a continuous one, which collapses the interval to a point and makes the residual exact and non-random.

Why not deviance or Pearson. For a discrete outcome they are misleading by construction: the raw, deviance and Pearson residuals of a Bernoulli fit each take one value per outcome level at a given fitted mean, so the scatter forms arcs — one per level — and a user who recognises the Residuals-vs-Fitted layout reads structure that is an artefact of the outcome being binary, not evidence about the model. Poisson fits at small μ band the same way. Quantile residuals have no such artefact and are exactly N(0, 1) under a correctly specified model (treating θ̂ as known), for every family — so one reading of the panel applies across the whole catalog.

Gamma is the deliberate uniformity choice, and it departs from R. The outcome is continuous, so deviance residuals do not band and R's plot.glm would draw them. We use the exact (non-randomized) quantile residual anyway, so that one residual concept, one reference line and one sentence in the User Guide cover all six families. The panel already departs from R's default for the discrete families — that departure is the entire reason the panel exists — and being internally consistent is worth more than matching a default we have already, correctly, left behind. The User Guide states which residual is drawn so a user cross-checking against R is never surprised.

The CDF bracket per family

μ is the fitted mean, α the NB2 dispersion, φ the gamma dispersion (the engine's Pearson statistic per residual df), π the structural-zero probability and θ the hurdle-crossing probability.

with the two mixture CDFs written in terms of the count CDF F_c and its mass at zero f_c(0):

F_ZI(k) = π + (1−π)·F_c(k)                            k ≥ 0,  F_ZI(−1) = 0
F_H(0)  = 1−θ
F_H(k)  = (1−θ) + θ·(F_c(k) − f_c(0)) / (1 − f_c(0))  k ≥ 1,  F_H(−1)  = 0

The three CDFs are stats.rs's poisson_cdf, negbin_cdf and gamma_cdf, each built on the regularized incomplete gamma / beta functions the engine already carries, and each pinned against SciPy in stats.rs's own tests.

The seed policy

Randomization collides with the reproducibility contract — the same journaled command must produce the same output — so the draw is seeded, from the engine's existing convention: the integer seed option, default 1 (DEFAULT_SEED in resample.rs), through the same SplitMix64 generator the bootstrap uses, drawn once per case in listwise order. Same command + same data ⇒ the same plot, always. The seed is named in the table note exactly as the bootstrap names its own.

Two rejected alternatives, recorded so they are not re-litigated. A clock or entropy seed breaks the contract outright. A data-derived seed (hashing the response column, say) looks attractive but is worse than a constant: it changes when any case changes, so two runs on nearly-identical data get unrelated realizations and cannot be compared — the opposite of what a diagnostic is for.

The seed is deliberately user-settable, because Dunn & Smyth's own advice is to inspect several realizations: any single one carries randomization noise, and a feature that appears under one seed and not another is noise, not misfit. Varying seed is how that check is made, and because it is journaled a report can say which realization it shows.

The reference lines — where the linear panel's convention does not transfer

• Quantile Residuals vs Linear Predictor — a horizontal line at zero. Quantile residuals have mean zero under a correct model, so this transfers unchanged from the OLS panel.
• Q-Q Plot of Quantile Residual — the identity line (slope 1, intercept 0), not the first/third-quartile fit the linear panel draws (the R qqline convention qq_plot uses). Quantile residuals are standard normal by construction under a correct model, so the theoretical line is known a priori; fitting a line through the sample quartiles would absorb exactly the location and scale misfit the plot exists to reveal. This is the one place in the panel where the OLS convention is deliberately not carried over, and the reference-line table records slope = 1, intercept = 0 as constants rather than as a fit.

The x axis: the linear predictor, not the fitted mean

Both scatters are drawn against the linear predictor η (the link scale). This is R's plot.glm convention — stats:::plot.lm takes predict(x), which for a glm object is the link-scale prediction, and carries the source comment ## != fitted() for glm. On the response scale the fitted values of a binomial or count fit compress against their boundary and the plot loses its resolution. For the zero-inflated and hurdle models it is the count component's linear predictor, since that is the component the count CDF describes.

The panel names

The OLS titles are never reused: "Residuals vs Fitted", "Scale-Location" and "Q-Q Plot of Std. Residual" name different quantities, and a returning user must not be told they are looking at the three plots they already know.

Each note names the family, says whether the residual is randomized, and gives the seed when it is.

Undefined residuals are skipped, not nulled

The 0.12.0 item-2d convention. A fitted probability at the boundary — perfect separation in a logistic fit — puts u at 0 or 1, so Φ⁻¹(u) is ±∞; that case is dropped from the table rather than emitted as a null. With no finite residual anywhere the table is empty and the chart layer's empty-check suppresses the plot, which is the honest outcome for a degenerate fit.

Out of scope for v1

No influence diagnostics for the GLM families (Cook's distance, leverage, DFBETA), no half-normal plots with simulated envelopes, and no binned-residual plot. Each is a defensible follow-up and each is its own validation problem.

Mediation follow-ups — the 0.21.0 conventions

Recorded 2026-08-21, before any code, for docs/releases/release_0.21.0.md item 5 (kept private): the three follow-ups the 0.15.0 moderated-mediation scope line left behind. Each sub-item adds its subsection here when its row a lands; the pins live in test/mediation_followups_validation_test.dart against mediation_followups_reference_values.json from generate_mediation_followups_references.R.

The fixture is a deterministic synthetic dataset (mediation_followups.csv, n = 200, fixed seed, create-if-absent) built so every follow-up has something to find: a real X·W interaction on the a-path and a real M₁·W interaction on the b-path (so the Johnson–Neyman quadratics have finite roots inside the observed range — asserted at generation), a second moderator Z with a three-way X·W·Z term, and a serial chain X → M₁ → M₂ → Y. The references are R lm for every OLS quantity (the same algebra statsmodels reproduced for 0.15.0) and the interactions package for the Johnson–Neyman boundaries.

Johnson–Neyman regions (5·JN)

• What is solved. For each moderated path, the conditional effect θ(w) = β₁ + β₃·w with Var θ(w) = V₁₁ + w²V₃₃ + 2wV₁₃ (the 0.15.0 conditional-slope convention); the boundaries are the roots of |θ(w)| = t_crit·SE(w), i.e. of (β₃² − t²V₃₃)·w² + 2(β₁β₃ − t²V₁₃)·w + (β₁² − t²V₁₁) = 0, with t_crit Student's t on the model's residual df at the interval's α (the confidence argument; two-sided). A vanishing quadratic coefficient falls back to the linear root; a negative discriminant means no crossing. PROCESS and interactions::johnson_neyman define the boundaries identically (same t, same df); the latter is the pin.
• Only boundaries inside the observed range of W are reported (PROCESS's convention; interactions solves the same roots but plots over a widened range). The effect at each boundary equals ±t_crit·SE by construction — the table prints both so the reader can see the identity — with the percent of analyzed cases below and above. Which side of a boundary is significant is read from the quadratic's sign at the interval midpoint, and the note spells the region out in words ("significant for W < a and W > b", "for W > a", "over the whole observed range", "nowhere in the observed range").
• The plot is the conditional effect and its confidence% band over a 41-point grid of the observed range (the boundaries inserted), with the zero line and the boundaries as vertical reference lines — a band-line chart (the forecast chart's geometry with a moderator axis). The grid table is emitted beside it so the chart replays from the journal.
• Not for the indirect effect. Johnson–Neyman is analytic for the OLS conditional effect only; the conditional indirect effect's inference is bootstrap-only, so no region is drawn for it (PROCESS makes the same choice). Recorded as a non-goal.
• Pins. The a-path (m1 ~ x*w) boundaries (both inside the range) with the effect, SE and percent-below at the lower one; the b-path (y ~ x + m1*w) boundaries and their count inside the range — against interactions::johnson_neyman(…, alpha = .05)$bounds and lm/vcov.

Three-way interactions (5·3W)

• The model shapes. A second moderator Z (moderator2, numeric, distinct from every other role) joins W in the moderated model(s) only with every product: a-path M ~ X + W + X·W + Z + X·Z + W·Z + X·W·Z (+ C) (PROCESS model 11's first stage), b-path Y ~ X + M + W + M·W + Z + M·Z + W·Z + M·W·Z (+ C) per mediator (model 18's second stage), both both. The coefficients are labelled a1…a7 / b1…b7 in that order (a5 = X·Z, a7 = X·W·Z; b5 = M·Z, b7 = M·W·Z). Without Z the output is byte-identical to 0.15.0's.
• Probes and the grid. Each moderator is probed by the 0.15.0 rule on its own (≤ 5 distinct values → each value; else mean ∓/± 1 SD, SD with n − 1); the grid is Z-outer, W-inner, so a table reads "at Z = z: W = …" down the page, and the interaction plot is drawn once per Z probe (titles "… at Z = label"). Probe values are held fixed across bootstrap resamples.
• Conditional effects at (w, z): the contrast a1 + a3·w + a5·z + a7·wz (b analogue) with Var = cᵀVc from the coefficient covariance, t on the residual df — pinned against lm/vcov.
• Indices (Hayes 2018). With exactly one path moderated: the index of moderated moderated mediation a7·b / a·b7 (how W's moderation of the indirect effect changes with Z) and the conditional index of moderated mediation by W at each Z probe (a3 + a7·z)·b / a·(b3 + b7·z); both bootstrapped like the conditional indirect effects. With both paths moderated neither exists (the indirect effect is bilinear in W), and the output says so.
• Johnson–Neyman with Z solves the regions of W at each Z probe on the combined coefficients θ(w | z) = (β₁ + β₅z) + (β₃ + β₇z)·w with the covariance of those two combinations — tables "… at Z = label".
• Pins. lm(m1 ~ x*w*z): a7 and its SE; the conditional slope and SE at (w_hi, z_lo); the conditional indirect effect there (b from lm(y ~ x + m1)); a7·b; (a3 + a7·z_lo)·b. lm(y ~ x + m1*w*z): b7 and its SE; the conditional b-slope and SE at (w_lo, z_hi).

Serial mediation (5·SER)

• The model. serial with exactly two mediators in causal order (the listed order: M₁ affects M₂): M₁ ~ X (+C), M₂ ~ X + M₁ (+C) — the d₂₁ path — Y ~ X + M₁ + M₂ (+C) and the total model Y ~ X (+C), all OLS (PROCESS model 6). Specific indirect effects a₁b₁, a₂b₂ and the serial a₁d₂₁b₂, their sum the total indirect; c = c′ + Σ holds exactly (asserted). The path rows are labelled a1, a2, d21, b1, b2, c′, c; the indirect rows X → M₁ → Y, X → M₂ → Y, X → M₁ → M₂ → Y, total. Pairwise contrasts between the specific indirects are not reported (a recorded cut).
• With one moderator. W enters the a-paths (X·W in both mediator models — M₂'s keeps M₁ after the moderated terms), the b-paths (M₁·W, M₂·W in the outcome model) or both; d₂₁ is not moderated (recorded; PROCESS model 83 / 87 shapes less the direct-path moderation). Each specific indirect effect is conditional on W — the serial one a₁(w)·d₂₁·b₂(w) — and, with one path moderated, linear in W, so each carries an index of moderated mediation: the serial path's is a₁₃·d₂₁·b₂ under a-moderation, a₁·d₂₁·b₂₃ under b-moderation. The serial effect is the third block of the conditional-indirect table and the third row of the index table. A second moderator is not combinable with serial (errors as data). Bootstrap CIs follow the v1 contract (not pinned).
• Pins. lm paths on the fixture: d21 and its SE, a2, the three specific indirects and their total, the total effect c; the a-moderated serial model's d21, the conditional serial indirect at the W mean and its index a3·d21·b2.

Regression through the origin — the 0.21.0 conventions

Recorded 2026-08-21, before any code, for docs/releases/release_0.21.0.md item 12: the intercept argument on linear_regression (default true; false fits the model without a constant — regression through the origin, y = b₁x₁ + … + bₖxₖ + e). The summary statistics change meaning without a constant and packages have differed historically, so the definitions are fixed here first; the pins are NIST's NoInt1 / NoInt2 in nist_strd_test.dart and a statsmodels set (reg_linear/noint_*) in numerical_validation_test.dart.

Why the statistics change. With a constant in the model the residuals sum to zero and the total sum of squares about the mean splits exactly into regression and residual parts (SST = SSR + SSE); R² = SSR/SST is then the share of the variance explained. Without a constant neither holds — the residuals need not sum to zero, the centred decomposition fails, and the centred R² can come out negative. Every reference that certifies or prints a through-origin fit (NIST StRD NoInt1/NoInt2; R's summary(lm(y ~ x - 1)); statsmodels OLS without add_constant, which switches to uncentered_tss when it detects no constant) therefore uses the uncentred decomposition about zero, and so does the engine:

n is the sum of the case weights throughout, as in the intercept model.

The sibling paths, decided per path:

• Robust (HC0/HC1/HC3) standard errors — allowed, unchanged formulas over the no-constant design (HC1's N/(N − p) uses p = k; HC3's leverage is the no-constant hat diagonal).
• Stepwise entry (forward/backward) — allowed; the selection oracle fits each candidate model through the origin too, so the p-values it ranks by are the ones the final table prints.
• Casewise diagnostics, diagnostic plots, saved columns — allowed; leverage and Cook's distance come from the no-constant hat matrix with p = k.
• Survey design — refused "(v1)" as errors-as-data, like casewise, diagnostics and save under a design (the design-based fit keeps its constant).
• Journal — command_for_request emits intercept: false only when the option is off, so every existing .cks script and every default run stays byte-identical; the dispatcher passes the option to the engine only when false.

References. NIST StRD NoInt1 (n = 11) and NoInt2 (n = 3), both average difficulty, certify B1, its standard deviation, the residual standard deviation, R² and the ANOVA table with F — all in exactly the uncentred form above (NoInt1's certified R² = 0.999365… is SSR/(SSR + SSE) = 200457.73/(200457.73 + 127.27); a centred R² on the same data would be 1.0). statsmodels OLS(bp, [age, chol, bmi]).fit() on clinic.csv pins R², adjusted R², F and its p, and the age coefficient and SE against the cross-package suite; R's lm(y ~ . - 1) prints the same definitions (summary.lm computes r.squared with df.int = 0).

Linear least squares by QR — the 0.21.0 conventions

Recorded 2026-08-21 for docs/releases/release_0.21.0.md item 14, raised from the StRD coverage work: the regression_linear solver moved from the normal equations to Householder QR, and the output gained a conditioning note. No statistic changed definition — the pins are the existing linear StRD datasets (Filip joining them) and the reg_linear/* cross-package set, all unchanged except that Filip now passes.

The solver. wls_qr forms A = √W·X (the design row with or without the constant), equilibrates its columns to unit norm, reduces it by Householder reflections to R, and back-substitutes R β = Qᵀ√W·y; the covariance bread is (XᵀWX)⁻¹ = D⁻¹ R⁻¹R⁻ᵀ D⁻¹ (D the column norms). Equilibration leaves the solution unchanged and makes the two checks below unit-free. Every downstream statistic reads β and the bread exactly as it read the inverted XᵀWX — model-based SEs, HC0/1/3 sandwiches, leverage, Cook's distance, the stepwise oracle's p-values, the survey-design Taylor sandwich — so there is one linear solver, not two.

Dependence. A column whose |R_kk| (on unit-norm columns) falls below 1e-13 is declared linearly dependent on those before it and the fit is refused with the existing "predictors are collinear" error-as-data — exact dependence leaves |R_kk| at rounding level (1e-16 … 1e-14), while NIST's Filip, the worst-conditioned design the engine is certified on, stays well above (its κ bound is 1.9 × 10⁷). n < p and an identically zero column are refused the same way.

The conditioning note (the gate-1 surface the Filip record asked for). The solver reports a lower bound on κ₂ of the equilibrated design — the ratio of the largest to the smallest |R_kk|, which never overstates the conditioning and is exact for an orthogonal design. At or beyond 1e8 the Model Summary note adds: "The design is ill-conditioned (condition number at least κ; from the QR of the unit-scaled predictors): the coefficients carry correspondingly fewer reliable digits — consider centring or rescaling the predictors." The threshold is the point where double precision leaves at most ~8 of its 16 digits in the coefficients. It is an accuracy note, not a collinearity diagnostic — Tolerance/VIF remain the statistical ones — and it is deliberately phrased as a bound ("at least"). Pinned by an ffi.rs test (present on degree-9 raw powers over [1, 1.2], absent on a plain predictor, exact dependence still refused).