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).

Scale invariance of the singularity verdict — the 0.22.0 convention

Recorded 2026-08-24 for release 0.22.0 item 1a. It governs analyze::correlation::invert, which sits under 76 call sites in 31 analyze modules — every GLM family, CFA, IRT, mediation, panel, survival, factor and VAR path that forms an information or cross-product matrix.

The rule. Whether a matrix is called singular must not depend on the units its variables are measured in. invert therefore two-sided-equilibrates first — every row and every column rescaled to max-abs 1, as m = diag(row) · scaled · diag(col) — runs Gauss-Jordan with partial pivoting on scaled, and undoes the scaling on the result. The pivot floor is 1e-13 on that equilibrated scale, matching the |R_kk| dependence floor the QR solver uses on unit-norm columns (above), and for the same reason: exact dependence leaves a pivot at rounding level, 1e-16 … 1e-14.

Two-sided rather than row-only equilibration is deliberate. Row scaling alone makes the verdict invariant under a uniform rescaling of the whole matrix; only two-sided scaling makes it invariant under per-variable rescaling, which is the case that actually arises — a covariance mixing variables in micrograms and tonnes.

What this replaced, and why it was wrong. Until 0.22.0 the floor was an absolute 1e-12 tested against the raw matrix, so the verdict moved with the units in both directions:

• an information matrix with entries near 1e6 was tested at an effective 1e-18 relative — accepting matrices with no reliable digits left and inverting them into printed statistics;
• one with entries near 1e-7 was tested at an effective 1e-5 relative — refusing perfectly well-conditioned matrices as "singular".

A statistic this was silently corrupting. The change was caught by multiple_imputation_test.dart, whose FMI fixture regressed a noiseless y = 2x + 3. Under the old inverter the twenty imputations' slopes differed in the 16th significant digit (1.9999999999999964 … 2.0000000000000018); Rubin's rules take the ratio of between- to within-imputation variance, both ~1e-28 here, and turned that pure rounding noise into a printed FMI of 0.638 — "64% of the information is missing", on data where the predictor determines the outcome exactly. The equilibrated inverter returns the slope exactly, every imputation is identical, and the FMI is 0, which is the correct answer: when the observed data determine the missing values, no information is missing. Both cases are now pinned — a fixture with residual scatter for FMI > 0, and the noiseless one for FMI = 0.

The conditioning bound. invert_with_condition additionally returns the ratio of the largest to the smallest pivot magnitude met during elimination — a lower bound on κ₂ of the equilibrated matrix, on the same "never overstates" footing as the QR bound above.

Non-matrices are data, not panics. A non-square argument, a non-finite entry, or an identically zero row or column all return None (the caller's existing "singular matrix" error-as-data) rather than propagating NaN into a printed statistic, which the old code did.

The SPD route: Cholesky for information matrices — the 0.22.0 convention

Recorded 2026-08-24 for release 0.22.0 item 1b.

Most of invert's call sites pass a symmetric positive-definite matrix — an information matrix, a cross-product XᵀX, a model-implied covariance. For those, analyze::correlation::invert_spd is the better route: half the arithmetic of Gauss-Jordan, no pivot search (an SPD matrix needs none), and a conditioning bound that falls out of the factorization for free.

The scaling is symmetric, d_i = 1/√m_ii, so the factored matrix is the correlation matrix of whatever m is the covariance of. That is the SPD-specific counterpart of the two-sided equilibration above and buys the same thing: the positive-definiteness verdict does not depend on the variables' units.

The two routes agree by construction. The Cholesky pivot at step k is L_kk², which is what Gauss-Jordan's k-th pivot reduces to on an SPD matrix, so both test the same 1e-13 floor on the same scale: a matrix one route calls singular, the other does too. This is pinned rather than asserted — invert_spd_agrees_with_gauss_jordan_on_every_spd_fixture runs both routes over an identity, a correlation matrix, a cubic design's XᵀX, Hilbert-5 and a badly scaled covariance, and requires agreement to 1e-9 relative to each row's largest entry (the small entries of an ill-conditioned inverse carry no digits of their own; demanding they agree absolutely would be testing noise).

The conditioning bound is the ratio of the largest to the smallest Cholesky pivot. For triangular T, σ_max ≥ max|T_kk| and σ_min ≤ min|T_kk|, so κ₂(L) ≥ max L_kk / min L_kk; squaring gives a rigorous lower bound on κ₂(LLᵀ). Like the QR bound it never overstates — pinned by requiring the identity to report exactly 1 and Hilbert-5 to report at or below its true κ₂ of ~4.8 × 10⁵.

Where it is surfaced. Cox regression's covariance and standard errors now come through this route, and the "Variables in the Equation" note carries the accuracy warning at or beyond the same 1e8 threshold the linear model uses, worded for its provenance: "The information matrix is ill-conditioned (condition number at least κ; from the Cholesky of the unit-scaled information): …". The note text is now produced by one shared regression::condition_note(condition, subject, provenance) for both the linear and the Cox surfaces, so the two cannot drift apart.

An absent bound is silent. A caller whose SPD route declines (a final information matrix that is not quite positive definite — a flat likelihood) records the bound as NaN, and the note treats NaN as "no bound was computed" rather than "infinitely ill-conditioned". This needs saying because NaN < threshold is false in IEEE arithmetic, so the guard has to be explicit or the note would print the word NaN at the reader.

Not yet adopted elsewhere. The other SPD call sites — CFA, IRT, factor, panel, mediation, the multinomial and ordinal families — still go through invert. That is a deliberate stopping point, not an oversight: each adoption moves the digits of a printed standard error and so needs its own pin check, and invert is correct for them (the two routes agree). The GLM families are the exception and are handled by item 1c, which takes them through the QR path rather than either inverse.

The normal quantile, and the far tail of the p-value — the 0.22.0 conventions

Recorded 2026-08-24 for release 0.22.0 item 1d. The second half of this section was not on the plan: it is what looking hard at the first half turned up.

The normal quantile Φ⁻¹

stats::inverse_normal_cdf is Acklam's rational approximation (absolute error ~1.2e-9) followed by one Halley step against normal_cdf. The step costs one CDF and one density evaluation and lifts the result to ~3e-16 relative, pinned against scipy.stats.norm.ppf from the pinned reference environment.

The step is taken in the smaller tail. Refining p = 0.999… directly would form the residual Φ(z) − p by subtracting two numbers agreeing in their leading digits, throwing away exactly the precision the step exists to gain; the reflected problem Φ(−z) = 1 − p is solved instead. The far tail, where the density has underflowed, keeps the unrefined value.

The pin was 5e-15 at first, not the 1e-15 the plan asked for, and the reason was recorded rather than papered over: the Halley step cannot be more accurate than the CDF it reads, and the incomplete-gamma routines then under normal_cdf floor at ~1–2e-15 (measured branch by branch against a 50-digit mpmath oracle). ln_gamma was tested as the suspect and is not — substituting an exact ln Γ(½) moved nothing. Getting below that needed a dedicated erfc for the normal case rather than a route through the incomplete gamma, which was deferred because it would move every normal p-value in the engine and wanted its own pass.

The floor moved (2026-08-26, item 12). With normal_cdf on the dedicated erfc described in the next section, the observed worst is 2.7e-16 relative (at p = 0.05, where the old route was weakest) and the pin is now 1e-15 — the figure originally asked for, carrying ~4× headroom because exp and ln may differ by an ULP between the three platforms in the tier-2 matrix. It is twenty million times tighter than the 1e-6 the old unit test asserted, which was itself six orders looser than the unrefined function achieved and so could not have caught a regression at all.

The far tail of the p-value

upper_gamma(a, x) was computing 1 − (1 − Q). gamma_continued_fraction already returns Q directly; lower_gamma wrapped it as 1 − Q for x ≥ a + 1; upper_gamma was 1 − lower_gamma. The accurate tail was computed and then destroyed by the round trip through 1. For Q below ~1e-16 the inner subtraction rounds to exactly 1.0 and the function returned exactly zero; above it, Q came back quantized to multiples of machine epsilon.

Since normal_two_tailed_p(z) is upper_gamma(0.5, z²/2) and chi_square_p(χ², df) is upper_gamma(df/2, χ²/2), that put a floor under every normal and chi-square p-value the engine printed: zero beyond |z| ≈ 8.3, and only ~6 correct digits by |z| ≈ 7. Nothing in the suite noticed, because a p-value of 0 and a p-value of 1e-20 both print as "< .001" — the defect was invisible at the surface and total underneath. upper_gamma now reads the continued fraction directly.

The crossover is a, not a + 1. The branch is chosen by which of P and Q is the small one, since that is what decides whether a subtraction from 1 costs digits — and Q becomes the small one once x passes the distribution's mean a. lower_gamma keeps a + 1 for the mirror-image reason. Sharing a + 1 would leave the band a ≤ x < a + 1 computing a smallish Q as 1 − P: at a = ½ that band holds z ≈ 1.0 … 1.7, so the 5% two-tailed critical value sat inside it, and Φ⁻¹(0.05) came back with ~9× the rounding error of Φ⁻¹(0.025) beside it. That asymmetry is what led to the finding.

Both are pinned against scipy at 1e-12 relative from z = 3 out to z = 37 (p ≈ 1e-299) and for chi-square to df = 20, plus a blunt companion assertion that a far-tail p-value is not zero — the shape the regression would take.

The normal tail: a dedicated erfc

Recorded 2026-08-26 for release 0.22.0 item 12 — the remainder item 1d measured and deferred.

normal_cdf(z) is ½·erfc(−z/√2) and normal_two_tailed_p(z) is erfc(|z|/√2), over stats::erfc — W. J. Cody's rational Chebyshev approximations (Math. Comp. 23, 1969; coefficients as in his SPECFUN CALERF), three ranges joined at |x| = 0.46875 and |x| = 4, each rational good to better than 1e-18 so that only the arithmetic is left. The exponential is formed as exp(−h²)·exp(−(x−h)(x+h)) with h = ⌊16x⌋/16, so the large exponent is the square of a number with four fractional bits and is exact. Before this the normal tail was upper_gamma(½, z²/2) — the same quantity through the general incomplete gamma, whose series and continued fraction both floor at ~1–2e-15. The chi-square, gamma and Poisson tails keep the incomplete gamma; they need it, and its floor is recorded above.

Pinned against a 50-digit mpmath oracle at 1e-15 relative, over 413 points in [−40, 40]: dense on [−6, 6], both tails, the three joins bracketed on both sides, and the tiny-argument region. The table (rust/src/erfc_mpmath_reference.txt, generated by tool/validation/generate_erfc_references.py, mpmath pinned in tool/validation/requirements.txt) prints each reference to 20 significant digits so the literal's own half-ULP rounding is the only oracle-side error. Observed worst on Linux/glibc: 5.2e-16. Past x ≈ 26.54 the true value is below the smallest normal double; those rows stay in the table and assert that the underflow is graceful — never negative, never NaN, in the same subnormal neighbourhood — rather than skipping the region.

Nothing moved. The engine-wide sweep this change was deferred to carry — 549 Rust tests and the full Dart suite, including every gate-4 reference pin under the ~100 call sites of the two functions — passed with no pin moved at any tolerance. That is the expected outcome, stated so it is not mistaken for an untested one: the old route was already correct to ~1e-15, and the pins on printed statistics sit at 1e-6 to 1e-12, so a change in the 16th digit of Φ cannot reach them. What it reaches is the quantile's Halley step (above) and the far tail below.

The far-tail pin tightened 1e-12 → 1e-13, and why not further. The observed worst against scipy is now 2.6e-14 at z = 37. Against the mpmath truth it is 8.8e-14 — and scipy itself is 1.1e-13 from that truth. The residual is not in erfc (a few ULPs) but in its argument: z/√2 rounds by half an ULP, and erfc's relative sensitivity to its argument is 2x², which at x ≈ 26 is ~1.5e-13. Every double-precision implementation shares that; scipy multiplies by 1/√2 rather than dividing by √2, as the engine now does, which is why the two agree more closely with each other than either does with the truth (the division form was measured at 2.1e-13 against scipy). A pin below the argument's own rounding would be pinning noise.

GLM: IRLS through the QR — the 0.22.0 convention

Recorded 2026-08-24 for release 0.22.0 item 1c.

fit_glm takes each Fisher-scoring step through wls_qr on the working problem, instead of forming XᵀWX and inverting it. The two are the same step in exact arithmetic; they are not in double precision, because forming the information squares the condition number of the working design — and the inverse it hands back is what every standard error, Wald statistic and confidence interval in the GLM output is built from. The covariance now comes off the QR factor as R⁻¹R⁻ᵀ.

The regressand is the step, not the textbook working response. Regressing (y − μ)/(dμ/dη) gives (XᵀWX)⁻¹·score — the Fisher step — while the textbook η + (y − μ)/(dμ/dη) gives the whole of β. Both are correct; the step form keeps the accumulated β in the accumulator rather than requiring the solve to reproduce its magnitude, which matters most when η has grown large, i.e. exactly when the fit is fragile.

What it bought, measured rather than asserted. On the well-conditioned validation fixtures the two routes are indistinguishable — the LRE against statsmodels is identical to two decimals for Poisson, negative binomial, probit and gamma, and within noise for logistic (14.58 vs 14.94). This is recorded because it is the honest result: squaring a small condition number costs nothing. The difference appears where the linear case's did — degree-8 raw powers over [1, 1.1], κ ≈ 6.9 × 10¹², where the formed information is singular to the inverter and the legacy route cannot take a step at all, while the QR takes a finite one (glm_qr_step_survives_an_ill_conditioned_design_the_normal_equations_lose).

The per-family LRE row, and the oracle that does not exist. NIST StRD certifies univariate summaries, linear least squares and nonlinear least squares — all 58 datasets are committed under test/validation/data/nist/, and none is a GLM. No standards body publishes certified logistic or Poisson coefficients. The row is therefore a cross-implementation agreement against statsmodels, not a certified accuracy, and is labelled as such:

Negative binomial sits lowest because its dispersion is estimated alongside the coefficients and the two implementations do not use the same inner criterion — a modelling difference, not a precision one. Floors sit under the achieved values with room for cross-platform libm variation.

A behaviour change this exposed. Under quasi-complete separation the logistic MLE does not exist and β diverges. The old path had an accidental brake: as μ saturated, XᵀWX shrank until it tripped the inverter's absolute pivot floor, the loop broke, and an arbitrary finite β was returned as though it were a fit. The QR route has no such brake, so β now runs to the iteration cap and the fitted probabilities saturate — which is what makes the separation visible to the guards that look for it. The propensity-score fixture turned out to be exactly such a case (2z − x is exactly 1 for all five treated and at most 1 for every control); it was replaced with a non-separable one, and the separating frame is kept as its own fixture behind quasi_separation_is_reported_not_silently_fitted. Reporting separation is the correct outcome; returning a truncated β was not.

Tolerances: what was tightened, and what stays loose — the 0.22.0 sweep

Recorded 2026-08-24 for release 0.22.0 item 1e. The standing rule for that item is that no tolerance is ever loosened, and every one that stays loose says why.

Measured first, then tightened. Every one of the 267 pins was run and its achieved error compared with its tolerance, rather than tightening by guesswork. Nine were tightened, each to at least 100× above what it actually achieves so cross-platform libm variation cannot trip them:

The tolerances are authored in generate_references.py and mirrored into the committed JSON; both were changed together and no reference value was touched, which is why the JSON diff is exactly nine lines.

The glm_resid family is not loose — it is not a tolerance. The plan listed "20 pins at tol ≥ 1e-2, glm_resid/logit_case0 at 0.46" as loose tolerances to tighten. Twelve of those twenty are containment brackets. For a discrete family the randomized quantile residual only satisfies Φ(r) ∈ [F(y⁻), F(y)], so the generator pins that interval as midpoint ± half-width and closeTo tests containment. logit_case0's 0.46 is half the width of [1 − μ, 1] at μ = 0.925 — not slack, the check itself. Tightening them would not sharpen anything; it would break them. (The gamma rows are different: a continuous family gives Φ(r) = F(y) exactly, and they are pinned at 1e-9.)

What stays loose, and why. nls/b0 and nls/b1 remain at 1e-2, and it is not a precision limit: in b0 + b1·exp(b2·(age−50)) the two trade off along a valley, so a second optimizer reaches the same minimum at a slightly different point — an identifiability property of the model, which no amount of arithmetic care removes. The certified nonlinear oracle is not this row at all: it is the NIST StRD suite, which runs all 27 nonlinear datasets against certified values (Misra1a LRE 9.3, Thurber 5.7).

The rule is now enforced, not promised. every tolerance at or above 1e-2 has a recorded reason fails if any pin sits at tol ≥ 1e-2 without an entry, and equally if an entry lingers after its pin is tightened. Containment brackets are recognised from their own source string rather than by an id list — whether a bracket's half-width lands above or below 1e-2 depends on the data, so the rule is what is stable — and their count is pinned so a regeneration cannot quietly turn them back into what look like sloppy tolerances.

The Gamma log-likelihood's dispersion — the 0.23.0 convention

Found by the second oracle (0.23.0 item 8): reg_gamma/m2ll is the one GLM pin on which R and statsmodels disagree, by 3.6 digits, while the coefficients, the deviance and the Pearson dispersion of the same fit agree to 10+.

The Gamma family's log-likelihood depends on a dispersion φ that the IRLS fit does not estimate, so every implementation has to pick one to evaluate the likelihood at, and two are in print:

• The Pearson estimate, φ̂ = χ²_P / (n − p) — what statsmodels' GLM reports as scale and evaluates llf at.
• The deviance estimate, φ̂ = D / n — what R's logLik.glm uses (through the family's aic function, which is documented to do so).

ChakataStat prints the Pearson form. The reported dispersion is the Pearson estimate (the Dispersion row and the standard errors use it), so the −2 log-likelihood, AIC and any likelihood-ratio comparison are evaluated at the same φ̂ the rest of the table is built on — one dispersion per fit, not two. It is also the number the primary pin holds. R's value is recorded beside it as a divergence, not a discrepancy: same fit, a different φ̂ in the likelihood.

The unbalanced panel's variance components — the 0.23.0 convention

Found by the second oracle (0.23.0 item 8, lane 5): on the balanced Grunfeld panel every Swamy–Arora quantity agrees between plm and linearmodels to 15 digits; on the unbalanced grunfeld_ub.csv five pins — θ's range, the RE slope and its SE, the Hausman statistic — disagree by 3–4 digits, all downstream of one number.

The random-effects estimator needs σ²ᵤ, the between-entity variance, and Swamy & Arora (1972) estimate it from the between regression's residual variance minus σ²ₑ divided by the common T. An unbalanced panel has no common T, and the literature offers more than one replacement:

• The harmonic-mean T̄ — Stata's xtreg, re and linearmodels: the balanced formula with T replaced by the harmonic mean of the entity lengths.
• Baltagi & Chang (1994) — plm's random.method = "swar": a trace-based degrees-of-freedom correction on the between regression (Baltagi 2013, §9.2) that reduces to the same formula when the panel is balanced.

ChakataStat uses the harmonic-mean T̄ (recorded in rust/src/analyze/panel.rs and pinned against linearmodels). Both are legitimate Swamy–Arora estimators; the choice is the one a Stata user will reproduce. plm's values are recorded beside the five affected pins as a divergence, not a discrepancy.