EDI 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 are in-progress (built fresh on every run)

EDI's C++ model-fitting kernels are being bound to Python in an active, uncommitted, in-progress scaffold under python/ (a separate effort from this benchmark script — see R/package_metadata/python_bindings_package_spec.md). This script does not own or edit that scaffold; at the top of every run it configures and builds whatever is currently there (via CMake, into an isolated /tmp directory) and imports whatever compiles, so the EDI Time column reflects the real state of that work at run time rather than a stale snapshot — different runs of this script can show different EDI coverage as that scaffold evolves, 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 without one, or where the Wald table needs standard errors these 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 matches R's wald_specs exactly (41 classes) rather than reusing the point-estimate table's 40 — it adds several nonparametric-test-only classes with no point-estimate-table row at all (pooled-variance t-test, Wilcoxon rank-sum, Lin's estimator, Fisher exact, 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). 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 (a scalar-loop binding would unfairly penalize a C++ side that hasn't been written yet).
• 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 (33 of 47 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).

33 of 47 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 (34 of 47 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 three rows 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.

34 of 47 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. Only fast_pchisq_upper is declared in the python/ fast_math stub (python/cpp/bindings_fast_math.cpp), and even that build currently fails to import (an unrelated undefined symbol from other in-progress bindings in the same compiled module — python/build_verify/), so no row has a usable EDI binding yet and every row is grey today.

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.