Validated Analyses

What golden validation means

The engine is checked against a committed fixture set generated once by R 4.5.3 (survival 3.8.6) by a test suite that runs as part of the engine's test run. The fixture file at python/tests/golden_r/fixtures.json contains the reference numbers produced by R on deterministic, seed-fixed datasets; the parametrized test suite at python/tests/test_golden_r.py loads those fixtures and compares 62 engine-output cases key-by-key at explicit tolerances. Coverage spans descriptive tests, regression and GLM families (including Gamma GLM and Poisson with exposure/offset), panel estimators, time-series and forecasting (seasonal decomposition, DFM backtest), PCA, and distribution fitting (MLE + Kolmogorov-Smirnov, vs MASS::fitdistr). Because the fixtures are committed — not regenerated on each run — neither a working R installation nor internet access is required to run the suite. The regeneration command (run from python/tests/golden_r/) is:

Rscript generate_fixtures.R

A further golden suite cross-validates the specialist latent-variable methods against their R reference package, using the same committed-fixture approach (generate_sem_lavaan_reference.R produces the numbers; the test bakes them in with the R command cited, so CI still needs no R):

• SEM family vs lavaan 0.6.21 (test_golden_r_sem_lavaan.py): confirmatory factor analysis, the structural model, and saturated path analysis. semopy and lavaan both fit by maximum likelihood with marker-variable scaling, so estimates agree to optimizer slack (raw estimates ~3×10⁻⁴, standardized ~1×10⁻⁴, fit indices ~1×10⁻⁴; model df exact).

SHAP is validated against the reference shap library — Shapley attributions are model-specific and the browser path uses a sampled KernelExplainer, so R is not the appropriate cross-validation anchor.

The suite demonstrably fails when formulas are perturbed: mutating the pooled-SD denominator in parametric.py (from n1 + n2 - 2 to n1 + n2) causes cohens_d_independent to fail; multiplying a fitted coefficient in regression.py's output path by 1.001 causes both linreg_simple and linreg_multiple to fail.

Tolerance disclosure

Survival and Poisson outputs are rounded to 4 decimal places by the engine; agreement with R is therefore verified to within 5×10⁻⁵ (half a unit in the last reported digit). All other analyses agree with R to a relative tolerance of 10⁻⁶ unless noted.

Three cases carry explicit scope limitations:

1. Tukey HSD adjusted p-values below approximately 1×10⁻⁸ are verified only as "both systems report p < 1×10⁻⁸". Both scipy and R approximate the studentized-range CDF (ptukey) and diverge in relative terms deep in the tail while remaining in full agreement that the result is far below any standard significance threshold. Mean differences are verified at rtol = 10⁻⁶.

2. Wilcoxon signed-rank: both the rank statistic (W) and the p-value are validated. The engine exposes a method="auto" | "exact" | "approx" option that mirrors R's wilcox.test semantics — auto uses the exact distribution when there are fewer than 50 non-zero differences with no ties or zeros, and the continuity-corrected normal approximation otherwise; an impossible exact request falls back to the approximation with a note, exactly as R warns and falls back. The golden suite pins both quantities against R: case_wilcoxon_signed_rank forces method="approx" to match the fixture's exact = FALSE, correct = TRUE, and W and the p-value both agree at rtol = 10⁻⁶. Additional cases pin the auto/exact/approx branches against R fixtures including tied and zero differences. One caveat: in the exact branch the engine reports no z-statistic (scipy's exact result carries none, and a two-sided exact p has no signed z to invert), so the z-based effect size r = |Z| / √n_eff is reported only when the approximation branch is used.

3. AIC and BIC are excluded from all comparisons. R and statsmodels/lifelines use different additive constants in their information-criterion formulas; a matching AIC would not validate the model and a mismatching one would not indicate a model error. Log-likelihood (llf) is compared instead for linear and logistic regression, where both systems agree to rtol = 10⁻⁶ (verified exactly on linreg_simple).

Validated analyses

Analysis	Engine function	Quantities compared	Tolerance
One-sample t-test	parametric.ttest_one_sample	t, df, p, mean	rtol 10⁻⁶
Independent t-test (pooled)	parametric.ttest_independent	t, df, p	rtol 10⁻⁶
Independent t-test (Welch)	parametric.ttest_independent	t, df, p	rtol 10⁻⁶
Paired t-test	parametric.ttest_paired	t, df, p	rtol 10⁻⁶
Cohen's d (independent)	parametric.ttest_independent	d (pooled-SD formula)	rtol 10⁻⁶
Mann-Whitney U	nonparametric.mann_whitney	U, p	rtol 10⁻⁶
Wilcoxon signed-rank	nonparametric.wilcoxon_signed_rank	W, p (approx branch; see note 2)	rtol 10⁻⁶
Shapiro-Wilk normality	diagnose.test_normality	W, p	rtol 10⁻⁶
Kolmogorov-Smirnov normality	diagnose.test_normality	D, p	rtol 10⁻⁶
Bartlett homogeneity	diagnose.check_homogeneity	statistic, p	rtol 10⁻⁶
Levene homogeneity (median)	diagnose.check_homogeneity	F, p	rtol 10⁻⁶
One-way ANOVA	parametric.anova_oneway	F, p, df, SS between, SS within	rtol 10⁻⁶
Tukey HSD post-hoc	parametric.anova_oneway	mean differences, adjusted p (see note 1)	rtol 10⁻³; atol 10⁻⁸ in tail
Kruskal-Wallis	nonparametric.kruskal_wallis	H, df, p	rtol 10⁻⁶
Friedman	nonparametric.friedman_test	χ², df, p	rtol 10⁻⁶
Chi-square independence	nonparametric.chi_square_independence	χ², df, p	rtol 10⁻⁶
Chi-square goodness-of-fit	nonparametric.chi_square_goodness_of_fit	χ², df, p	rtol 10⁻⁶
Pearson correlation	correlation.correlation_pair	r, p, 95% CI	rtol 10⁻⁶; CI at rtol 10⁻⁴ (1.96 approximation)
Spearman correlation	correlation.correlation_pair	ρ, p	rtol 10⁻⁶
Kendall correlation	correlation.correlation_pair	τ, p	rtol 10⁻⁶
Linear regression (simple)	regression.linear_regression	coef, SE, t, p, R², adj-R², F, log-likelihood	rtol 10⁻⁶
Linear regression (multiple)	regression.linear_regression	coef, SE, t, p, R², adj-R², F, log-likelihood	rtol 10⁻⁶
PCA	factor_cluster.pca_analysis	variance ratios, corrected sdev	rtol 10⁻⁶
Logistic regression	regression.logistic_regression	coef, SE, z, p, log-likelihood	rtol 10⁻⁴ (iterative MLE)
Poisson regression	regression_ext.poisson_regression	coef, SE, p, deviance	atol 5×10⁻⁵ (engine rounds to 4 dp)
Kaplan-Meier survival	survival.kaplan_meier	survival probabilities (first 10 times), median	atol 5×10⁻⁵ (engine rounds to 4 dp); median at rtol 10⁻⁶
Cox proportional hazards	survival.cox_regression	coef, SE, z, p, HR	atol 5×10⁻⁵ + rtol 10⁻⁴ (iterative, Efron ties)
Log-rank test	survival.kaplan_meier	χ², p	atol 5×10⁻⁵ (engine rounds to 4 dp)
Cronbach's alpha	factor_cluster.cronbach_alpha	α (classic k/(k-1) formula)	rtol 10⁻⁶
Confirmatory factor analysis	latent.cfa	free loadings (raw/SE/std.all), factor correlation, χ²/df/CFI/TLI/RMSEA	vs lavaan; estimates atol 5×10⁻⁴, df exact
Structural equation model	latent.sem	structural coefficient + measurement loadings (raw/SE/std.all), fit indices	vs lavaan; estimates atol 5×10⁻⁴, df exact
Path analysis (saturated)	latent.path_analysis	path coefficients (raw/SE/std.all), R² per endogenous	vs lavaan; estimates atol 5×10⁻⁴, R² atol 10⁻⁶
R environment: R version 4.5.3 (2026-03-11 ucrt), survival 3.8.6, jsonlite 2.0.0; latent-variable suite uses lavaan 0.6.21.

NIST StRD numerical-accuracy benchmark

Beyond agreement with R, the engine is benchmarked against the NIST Statistical Reference Datasets (StRD) — the recognised certified-value standard for the numerical accuracy of statistical software. Committed tests (test_nist_strd_univariate.py, test_nist_strd_anova.py, test_nist_strd_linreg.py) report the correct-digit count (LRE = −log₁₀(|engine − certified| / |certified|)) against NIST's multiple-precision certified values, with honest floors:

• Univariate summary statistics (PiDigits, Lew, Mavro, NumAcc1–4): mean and sample SD reach 13–15+ correct digits on well-conditioned data. On the NumAcc stress sets (SD/|mean| ≈ 10⁻⁷–10⁻⁸) the SD degrades gracefully to ~9.5 (NumAcc3) and ~8.3 (NumAcc4) digits — the documented cost of float64 cancellation; the engine's stable two-pass variance avoids collapse.
• One-way ANOVA (SiRstv, AtmWtAg, SmLs01/04/07): F, SS-between and SS-within reach 12–15 digits at ordinary data scale; df always exact. Known limitation: the textbook deviation-formula sum-of-squares loses precision on extreme-scale data — SmLs07 (values ≈ 10¹²) drops to ~2.7–4.4 correct digits (catastrophic cancellation). Documented, not hidden.
• Linear least squares (Norris, Pontius, Longley, Wampler1, Filip): Norris and the classically-multicollinear Longley reach ≥10 digits (the SVD solver shrugs off collinearity); polynomial fits 6–14 digits. Known limitation: on Filip (degree-10, condition number ≈ 10¹⁵) the engine's default pinv/SVD OLS solver returns 0 correct digits — a total failure under extreme ill-conditioning. statsmodels' QR solver recovers ~7.9 digits on the same data; moving the engine to QR (with a pinv fallback for rank-deficient designs) is the recommended fix.

A full per-method validation map — what is NIST-certified vs R-cross-validated vs library-faithful vs internally rigor-audited vs unverified — is maintained at validation-matrix.md. It is honest by construction (adversarially audited; a 2026-06-17 audit caught and corrected three over-tier claims and a procedure-counting error). Of 415 distinct analyses: 3 NIST-certified, 214 R-cross-validated, 18 library-faithful (parity-tested), 50 internally rigor-audited, 130 unverified beyond a smoke contract — most of the last being thin calls into peer-used libraries (statsmodels/scipy/lifelines/etc.).

The academic-validation campaign (2026-06)

The R-cross-validated tier grew from 42 to 214 after the committed academic-validation campaign — the suite python/tests/test_acadval_wave1..21_*.py, which drives ~224 engine methods against R across 17 families (GLM, ordered/limited-dependent, penalized/system regression, mixed models, panel estimators, diagnostics, post-estimation margins, time-series core + advanced, nonparametric/categorical tests, ANOVA/MANOVA, factor/scaling, IRT/SEM/latent-class, survival, meta-analysis/resampling, complex survey, missing data, and modern causal designs). Each wave's reference numbers were produced once by a deterministic, committed R script (R 4.5.3 with the relevant CRAN packages — MASS, pscl, betareg, plm, pgmm, lmtest, sandwich, car, forecast, urca, vars, tsDyn, rugarch, ARDL, mirt, eRm, difR, CTT, poLCA, mclust, Gifi, MASS::lda, cancor, survival, flexsurv, cmprsk, survRM2, metafor, boot, survey, mice, naniar, norm, Matching, AER, …) and are baked into the test files as documented constants — CI runs the comparison with no R installation required.

172 methods were newly R-validated. Of the 214 total, 49 carry an honest "R (partial: …)" marker: the engine agrees with R on its primary outputs (point estimates / closed-form statistics) but a named secondary quantity legitimately differs — a documented backend or algorithm disagreement (e.g. ZIP/ZINB standard errors vs pscl, FGLS intercept vs nlme::gls, System-GMM Windmeijer SE vs pgmm, Hotelling-Lawley F-approx vs stats::manova, HAC bandwidth/kernel SE vs sandwich::vcovHAC, bootstrap Monte-Carlo SE/CI/p with a different RNG, or AIC/BIC additive-constant conventions). These disagreements are pinned in the tests at the magnitude at which they actually occur, never hidden behind a loosened tolerance, and are inventoried per-method in the validation matrix. Methods that showed no clean R agreement on any primary quantity (Ridge/ElasticNet vs glmnet objective normalization, panel RE/between θ, DCC-GARCH grid-vs-MLE, IRT CML-vs-ML, MGWR, k-means/MDS/MCA algorithm mismatches, etc.) were deliberately not upgraded — they remain library-faithful or internally rigor-audited with the disagreement documented.