EDI v1.0.0.post3 Python Baseline Benchmarks

Generated: 2026-08-17 14:51:03 UTC

Python analog of benchmark_model_fits_R.html — same three tables (point-estimate, Wald/full-inference, utility math kernels), same table shape, same three-color row coding. Produced by R/benchmark/benchmark_model_fits_python.py.

Status: EDI Python bindings (built fresh on every run)

EDI's C++ model-fitting kernels are bound to Python in the edi_kernels package under python/ (released on PyPI — see R/package_metadata/python_bindings_package_spec.md). This script does not own or edit that package; at the top of every run it configures and builds its _core pybind11 module (via CMake, into an isolated /tmp directory) and imports whatever compiles, so the EDI Time column reflects the real state of the bindings at run time rather than a stale snapshot — different runs of this script can show different EDI coverage as the bindings evolve, and a build failure there simply leaves every EDI column NA (not a crash of this script). Point-estimate rows with a working, wired binding call it directly; rows where the Wald table needs standard errors the bindings don't yet expose stay NA. Same-input, same-process discipline: for every wired row, both the EDI and canonical closures are built from the exact same generated arrays (never redrawn separately). Each side is timed at its own fastest available point-estimate-only entry point — the same convention the R report already uses (R's canonical glm.fit() point-estimate row likewise never forms the explicit variance-covariance matrix; that only happens in the separate Wald table). EDI rows therefore use estimate_only=True where that argument exists. Where statsmodels' public API has no equivalent cheaper mode — .fit() is the only entry point it exposes, full stop — the canonical column still reflects .fit()'s real, honest cost: that is genuinely the fastest a Python user can get this model fit today, not a handicap imposed on canonical Python. (Checked directly: statsmodels' .fit() doesn't secretly race ahead by skipping something EDI bothers to compute — accessing .bse after .fit() costs no measurable additional time, since IRLS's per-iteration weighted-least-squares solve already forms the same (X'WX)⁻¹ a covariance estimate would need; the two APIs just draw the "point estimate vs. full inference" line in different places.)

Benchmark Dataset Specification

• Sample size (N): 1,000 subjects for most families; 500 for survival families and the Wilcoxon Hodges-Lehmann row (O(n²) pairwise-difference computation, matching the R harness's scale = 0.5); 5,000 elements for the utility-function vectors.
• Predictors (p): intercept + a balanced binary treatment + 4 continuous covariates ~ Normal(0, 1); covariate coefficients ~ Normal(0, 0.5); treatment coefficient fixed at 0.5.
• Same family-generation formulas as benchmark_model_fits.R (logistic/log-binomial/identity-binomial links, Poisson mean model, Beta(mu*phi, (1-mu)*phi) proportions, exponential survival times with ~20% censoring, 3-level ordinal construction).
• Stratified Cox exception: a low-cardinality (2x3) strata grid is injected before outcome generation, matching the R harness, so the row exercises a genuinely stratified fit.

Methodology

• Bare-metal canonical timing: each canonical row constructs the model object and calls its lowest-level .fit()/equivalent directly on pre-built NumPy arrays (or, where the package requires it — lifelines, statsmodels.OrderedModel internals — a pre-built pandas.DataFrame), inside the timed region only; data generation happens once, outside the timed closure.
• Point-estimate vs. Wald tables: the point-estimate table times a bare fit (matching each row's fastest available canonical entry point — lstsq instead of OLS().fit() where that's a real, distinct fast path); the Wald table times a full fit that also produces the treatment coefficient's standard error and two-sided p-value (.bse/.pvalues off the same fitted result for most statsmodels families; a package switch to lifelines.CoxPHFitter for the unstratified-Cox row specifically, since scikit-survival's bare-metal CoxPHSurvivalAnalysis doesn't expose a variance; a finite-difference delta-method SE for the four G-computation rows, off .cov_params()). The Wald table's row set follows R's wald_specs rather than reusing the point-estimate table's — it adds several nonparametric-test-only classes with no point-estimate-table row at all (Lin's estimator, Miettinen-Nurminen/Newcombe risk-difference CIs, Jonckheere-Terpstra, Ridit, Gehan-Wilcoxon, KM median difference) and omits several point-estimate-table rows that R's own Wald table never covers (identity-binomial, modified Poisson, fractional logit, ordered probit, cauchit, cloglog). Classes with no EDI C++ kernel at all — the EDI R class delegates to base R or a CRAN package (quantile regression via quantreg::rq, the pooled-variance t-test, the Wilcoxon rank-sum test itself, Fisher exact) — are omitted from both tables entirely: with no EDI side to time, a permanently-NA row carries no information. One row, restricted-mean-survival-time difference, is excluded from the Wald table for the same reason R's own Wald table excludes it — see below.
• Utility-function timing: each row calls the scipy/numpy vectorized function once over a fixed-length input vector, mirroring benchmark/fast_math_utils_bench.cpp's apples-to-apples vectorized-vs-vectorized discipline (the EDI side likewise calls the vectorized elementwise wrappers exposed by python/cpp/bindings_fast_math.cpp, never a Python-level scalar loop).
• Averaging: medians over 30 cold timing samples via an adaptive-batch time.perf_counter harness (mirrors the R harness's adaptive system.time split, target 200ms/batch).
• Significance: Welch's two-sample t-test (scipy.stats.ttest_ind(..., equal_var=False)) between the EDI and canonical timing replicate distributions — real for any row with a wired, working EDI binding; NA for rows where no binding is wired or available.
• Row highlighting: light green = Speedup > 1 and Timing Pval < 0.05; light grey = NA timing comparison (EDI not bound yet, or a fit failed); light blue = no canonical Python implementation exists at all for this model family/function.
• Package versions this report was run against: numpy 2.5.1, scipy 1.18.0, pandas 3.0.5, statsmodels 0.14.6, sksurv 0.28.0, lifelines 0.30.0.

Point-Estimate Results (32 of 46 rows timed, 14 Baseline Gaps)

EDI faster, significant   NA timing comparison   no canonical Python implementation

Families with no clean, actively-maintained Python canonical equivalent — an absent/NA comparison here is more honest than a mismatched substitute baseline (same discipline python_bindings_package_spec.md and the R report apply):

• InferenceContinKKGLMM (future kernel fast_gaussian_lmm): Gaussian LMM with a random intercept, fit via MLE: statsmodels' MixedLM uses REML/a different estimation path by default and isn't a direct like-for-like timing comparison -- see python_bindings_package_spec.md Baseline Gaps.
• InferenceCountKKCondPoissonOneLik (future kernel fast_cpoisson_combined): KK combined (matched-pair + reservoir) joint-likelihood estimator: no canonical analog in either R or Python.
• InferenceCountKKGLMM (future kernel fast_poisson_glmm): Poisson GLMM with a random intercept, fit via adaptive Gauss-Hermite quadrature + MLE: no pure-Python package offers a comparable ML (not variational/Bayesian) GLMM fitter -- see python_bindings_package_spec.md Baseline Gaps.
• InferenceCountKKHurdlePoissonOneLik (future kernel fast_hurdle_poisson_glmm): Hurdle-Poisson GLMM with a random intercept: no pure-Python package combines a hurdle count model with adaptive-quadrature GLMM fitting.
• InferenceIncidKKCondLogitGLMMOneLik (future kernel fast_clogit_plus_glmm): KK combined (matched-pair + reservoir) joint-likelihood estimator: no canonical analog in either R or Python (see python_bindings_package_spec.md Baseline Gaps).
• InferenceOrdinalAdjCatLogitRegr (future kernel fast_adjacent_category_logit): No identified Python package implements the adjacent-category logit link (R uses VGAM::vglm(acat())).
• InferenceOrdinalCauchitRegr (future kernel fast_ordinal_cauchit_regression): Same as cloglog: OrderedModel's distr= only documents 'probit'/'logit'.
• InferenceOrdinalCloglogRegr (future kernel fast_ordinal_cloglog_regression): statsmodels.miscmodels.ordinal_model.OrderedModel's distr= argument only officially supports 'probit'/'logit' strings; cloglog is not a documented/tested option, so this is treated as a gap rather than an unverified custom-distribution hack.
• InferenceOrdinalContRatioRegr (future kernel fast_continuation_ratio_regression): No identified Python package implements the continuation-ratio link (R uses VGAM::vglm(cratio())).
• InferenceOrdinalKKCLMM (future kernel fast_ordinal_clmm): Same as InferenceOrdinalKKGLMM, generalized to logit/probit/cauchit/cloglog links.
• InferenceOrdinalKKGLMM (future kernel fast_ordinal_glmm): Proportional-odds ordinal GLMM with a random intercept: no pure-Python package offers ML ordinal-GLMM fitting with adaptive quadrature.
• InferenceOrdinalStereotypeLogitRegr (future kernel fast_stereotype_logit): Stereotype logit ordinal regression: R uses VGAM::vglm(multinomial(...)) style fitting; no identified Python package implements the stereotype-logit link specifically.
• InferencePropZeroOneInflatedBetaRegr (future kernel fast_zero_one_inflated_beta): Zero-one-inflated beta regression: no canonical package in either R or Python.
• InferenceSurvivalKKWeibullFrailtyOneLik (future kernel fast_weibull_frailty): Weibull AFT with shared log-normal frailty: no clean Python package (R side notes only a partial/PH-parameterized frailtypack analog).

32 of 46 rows have a working canonical Python timing; 0 canonical fit(s) failed on this run (see console log); 14 are documented Baseline Gaps.

Wald Test Performance / Full Inference (30 of 43 rows timed, 13 Baseline Gaps)

EDI faster, significant   NA timing comparison   no canonical Python implementation

Same Baseline Gap families as the point-estimate table above, plus 1 row(s) excluded outright (no variance computation, not a Baseline Gap in the usual sense):

• InferenceSurvivalRestrictedMeanDiff: Not in R's own Wald table either: lifelines' restricted_mean_survival_time (like R's survival:::survmean) returns only a point value; a variance requires a separate bootstrap-based computation, a different computational profile than every other Wald row here.

30 of 43 rows have a working canonical Python timing; 0 canonical fit(s) failed on this run (see console log); 13 are documented Baseline Gaps; 1 rows from the point-estimate table are excluded entirely (see above).

Utility / Math Kernel Performance (14 of 14 functions timed)

EDI faster, significant   NA timing comparison   no canonical Python implementation

EDI's internal fast_* scalar math kernels — every one that exists in R/EDI/src (fast_digamma, fast_trigamma, fast_lgamma, fast_lbeta, fast_qnorm, fast_log_pnorm, fast_log_dnorm, fast_dnbinom_mu, fast_pchisq_upper, fast_erfc, pnorm_fast, dnorm_fast, fast_atan, fast_log1pexp) — vs. their scipy/numpy vectorized equivalents, over a length-5000 vector. All 14 are bound in python/cpp/bindings_fast_math.cpp as vectorized elementwise wrappers around the same portable inline scalar functions EDI's model-fitting kernels call internally, so every row times a real EDI call against its canonical vectorized counterpart.

14 of 14 functions have a working canonical timing; 0 failed on this run.

See R/package_metadata/python_bindings_package_spec.md for the full kernel-binding plan.