StatsPAI

简体中文 | English

StatsPAI: an Agent&Python-native Stata/R replacement for applied causal inference

    
Developed and maintained by the Stanford REAP team
Published in the Journal of Open Source Software (JOSS): Wang & Rozelle (2026), 11(125), 10604 · doi:10.21105/joss.10604

StatsPAI is for empirical researchers who would normally jump between Stata, R, and Python. Its goal is to make common Stata/R econometrics and causal-inference workflows feel native in Python: load a dataset, estimate a model, inspect diagnostics, export tables, and hand the result to an agent or notebook without leaving one API.

It is meant to be a practical replacement path for new Python-first work:

• Stata-style routines: regress, ivregress, reghdfe, csdid, rdrobust, synth, psmatch2, esttab / outreg2.
• R-style routines: lm, fixest, did, rdrobust, Synth, DoubleML, MatchIt, modelsummary, broom.
• Stata conventions where they matter: vce="robust" / vce="cluster firm" on regress, ivreg and the likelihood estimators, Stata’s small-sample factors and z / t reference distributions, test / lincom / margins, dydx() after any fit, and sp.stata("...", data=df) to run Stata command lines you already have.
• Python-native outputs: .summary(), .tidy(), .plot(), .to_latex(), .to_docx(), .to_agent_summary() where supported by the result object.
• Agent-native access: every public function is registered with a machine-readable schema (sp.list_functions(), sp.describe_function(), sp.function_schema()), and the bundled statspai-mcp server exposes the estimators to MCP clients such as Claude Code, Claude Desktop, and Cursor.
• Companion Stata tooling: our own stata-code can work with StatsPAI so agents can understand existing Stata workflows, translate them into Python, and cross-check results more smoothly.
• Companion skill repos: Auto-Empirical-Research-Skills, AER-Skills, Awesome-Journal-Skills, and Paper-WorkFlow can work alongside StatsPAI and an agent as the methods, journal, manuscript, and reproducibility skill layer.

StatsPAI is not a promise that every Stata/R command is bit-for-bit identical. The API is broad, and the numerical evidence behind it is uneven: some estimators are checked against R/Stata on identical data, others only against known-truth simulations, and many are API-stable without a numerical-parity claim yet. Every function carries a validation_status that says which case applies: validation_status distinguishes certified/validated evidence from API-stable breadth. See Validation before relying on a number for publication.

Install

pip install statspai

Python 3.10 – 3.13 on Linux and macOS. The core install covers estimation, diagnostics, the bundled datasets, .xlsx / .docx / LaTeX export, and matplotlib figures (.plot(), sp.ggdid() and the rest). seaborn / plotly styles and heavier backends are optional extras:

The interactive figure editor additionally needs ipywidgets inside Jupyter.

import statspai as sp

print(sp.datasets.list_datasets()[["name", "design", "source"]])

sp is the convention used throughout these docs, but the alias is your choice. If sp is already taken in your code (for example by import scipy as sp), import statspai as si works just as well.

StatsPAI bundles 14 datasets that load offline. Most are real published extracts (source == "bundled CSV"): Card (1995) NLSYM schooling data, LaLonde/NSW with a PSID comparison group, the U.S. Senate RD data distributed with R’s rdrobust, California Proposition 99, the castle-doctrine panel, and NHEFS, among others. A few are deterministic simulated replicas calibrated to a published design (source == "simulated"), including the Callaway–Sant’Anna mpdta panel used below; their numbers are not the numbers from the original data.

At a glance: 1,680 registered functions across 92 submodules; 756k LOC (core) + 519k LOC (tests). Run python scripts/registry_stats.py to reproduce these numbers.

If You Come From Stata Or R

On regress, ivreg, glm, logit, probit, poisson, nbreg, the ordered / multinomial / conditional logits, the zero-inflated and hurdle models, liml and more, vce= / robust= follow Stata’s vce() grammar (sp.feols keeps fixest’s vcov= / cluster=): True, "robust", "vce(robust)", "oim", "hc0"–"hc3", and a cluster variable written inline ("cluster firm", "vce(cluster firm)", "cl firm"). A spelling an estimator does not implement raises an error instead of quietly falling back to different standard errors. Robust and cluster SEs apply Stata’s small-sample factors, clustered OLS / IV use t(G-1), and likelihood-based fits report z. The oim / robust / cluster standard errors of 29 estimators are pinned to Stata 18 at the 1e-6 parity budget (three documented exceptions where the two optimisers stop at slightly different points). See the shared argument grammar guide.

Causal entry points also accept one shared set of argument names next to each estimator’s native spelling: id= for the panel unit, time=, first_treat= for adoption cohorts, covariates=, and running= / cutoff= for RD. sp.callaway_santanna(data=mp, y="lemp", time="year", id="countyreal", first_treat="first_treat") is the same call as the t= / i= / g= version below; a misspelt keyword names the closest one (runing= → “did you mean ‘running’?”), and sp.describe_function(name)["aliases"] lists the accepted spellings.

sp.esttab(), sp.outreg2(), and sp.modelsummary() still exist, but they are deprecated thin wrappers over sp.regtable() and emit a DeprecationWarning.

Compared With Other Python Packages

StatsPAI aims to be one broad Stata/R-style workbench. Several focused Python packages do one part of that job, often with a longer track record; if you only need that part, they are good choices.

Use StatsPAI when you want one package, one function registry, and one agent interface across the everyday Stata/R empirical workflow.

Beginner Examples With Results

The outputs below were produced with StatsPAI 1.29.0 on the bundled datasets and are pinned by tests/test_readme_examples.py, so they cannot drift from the code silently. Example 6 uses the Stata vce() grammar and post-estimation commands introduced in 1.29.0; on an older release install from source with pip install "statspai @ git+https://github.com/brycewang-stanford/StatsPAI". Long summaries are abridged (... marks omitted lines); the numbers are pinned by tests/test_readme_examples.py and tests/test_synth_placebo_pvalue.py, so they cannot silently drift from the code again.

1. OLS: the first regress / lm replacement

Question: how much higher is log wage for one more year of schooling in the Card (1995) NLSYM data?

import statspai as sp

card = sp.datasets.card_1995()
ols = sp.regress(
    "lwage ~ educ + exper + expersq + black + south + smsa",
    data=card,
    robust="hc1",
)
print(ols.summary())

Result:

Model: OLS
Method: Least Squares
Dependent Variable: lwage
...
           Coefficient  Std. Error  t-statistic  P>|t|  [0.025  0.975]
Intercept       4.7337      0.0702      67.4718 0.0000  4.5961  4.8712
educ            0.0740      0.0036      20.3208 0.0000  0.0669  0.0812
exper           0.0836      0.0067      12.4165 0.0000  0.0704  0.0968
expersq        -0.0022      0.0003      -7.0443 0.0000 -0.0029 -0.0016
black          -0.1896      0.0174     -10.8781 0.0000 -0.2238 -0.1555
south          -0.1249      0.0154      -8.1339 0.0000 -0.1550 -0.0948
smsa            0.1614      0.0152      10.6374 0.0000  0.1317  0.1912

Model Diagnostics:
--------------------
R-squared           : 0.2905
...

Read it like a Stata/R regression table: conditional on experience, race, region, and SMSA, one more year of schooling is associated with about 0.074 higher log wage (roughly 7.4%). This is a correlation, not yet a causal return: schooling is plausibly correlated with unobserved ability, which motivates the IV in example 2. The HC1 standard errors follow Stata’s vce(robust) / sandwich::vcovHC(type = "HC1") convention.

2. IV / 2SLS: replace ivregress 2sls or AER::ivreg

Question: instrument schooling with growing up near a four-year college (nearc4).

import statspai as sp

card = sp.datasets.card_1995()
iv = sp.ivreg(
    "lwage ~ (educ ~ nearc4) + exper + expersq + black + south + smsa",
    data=card,
)
print(iv.summary())

Result:

Model: IV-2SLS
Method: Two-Stage Least Squares
Dependent Variable: lwage
...
           Coefficient  Std. Error  t-statistic  P>|t|  [0.025  0.975]
...
educ            0.1323      0.0492       2.6870 0.0072  0.0358  0.2288

Model Diagnostics:
...
First-stage F (educ)        : 16.7176
...
Partial R² (educ)           : 0.0055
Hausman F-stat              : 1.5390
Hausman p-value             : 0.2149

The IV estimate (0.132) is larger than OLS but about 13 times less precise. The instrument is not strong — nearc4 explains only 0.55% of the residual variation in schooling (first-stage F ≈ 16.7) — and the Hausman test does not reject exogeneity of educ (p = 0.21). Default standard errors are the unadjusted ones with the small-sample correction used by AER::ivreg (Stata: ivregress 2sls ..., small); pass robust="hc1" for heteroskedasticity-robust errors, and the first-stage F then uses the same variance estimator, as Stata’s estat firststage does.

With a first stage this modest, report a weak-instrument-robust interval too:

ar = sp.anderson_rubin_ci(
    y="lwage", endog="educ", instruments=["nearc4"],
    exog=["exper", "expersq", "black", "south", "smsa"], data=card,
)
print(ar.summary())

Anderson-Rubin (AR) — weak-IV-robust confidence set
------------------------------------------------------------
  level                : 95%
  grid                 : 401 points on [-0.359, 0.624]
  confidence set       : [0.0384, 0.2612]

3. Staggered DiD: replace csdid or R did

Question: what is the average effect of minimum-wage increases on teen employment in the Callaway–Sant’Anna mpdta design?

import statspai as sp

mp = sp.datasets.mpdta()   # simulated replica of R did's mpdta
gt = sp.callaway_santanna(
    data=mp,
    y="lemp",
    t="year",
    i="countyreal",
    g="first_treat",
)
overall = sp.aggte(gt, type="simple", bstrap=False)
print(overall.summary())

Result:

==============================================================================
  Callaway and Sant'Anna (2021) — aggte[simple]
==============================================================================

  ATT:      -0.0330 ***
  Std. Error:  (0.0078)
  [95% CI]:    [-0.0482,  -0.0178]
  P-value:     <0.001
...
  Observations:    2,500
...

The aggregated ATT is about -0.033 log points and statistically precise. The bundled mpdta is a calibrated simulated replica, so this is not the number R reports on the original mpdta data. What is checked: the same CSV run through R did::att_gt() + aggte() and Stata csdid returns the same ATT and standard error (Track A parity module 04_csdid).

4. Regression discontinuity: replace rdrobust

Question: is there a party incumbency advantage at the zero-margin cutoff in U.S. Senate elections?

import statspai as sp

senate = sp.datasets.lee_2008_senate()  # rdrobust's Senate data: x = margin, y = vote share
rd = sp.rdrobust(data=senate, y="y", x="x", c=0)
print(rd.summary())

Result:

==============================================================================
  Sharp RD Estimation
==============================================================================

  RD Effect:       7.51 ***
  Std. Error:  (1.74)
  [95% CI]:    [4.09,  10.92]
  P-value:     <0.001

------------------------------------------------------------------------------
  Inference
------------------------------------------------------------------------------
      method  estimate     se      z  pvalue  ci_lower  ci_upper
Conventional    7.4141 1.4587 5.0826  0.0000    4.5551   10.2732
      Robust    7.5065 1.7413 4.3110  0.0000    4.0937   10.9193

------------------------------------------------------------------------------
  Observations:    1,297
...
  Bandwidth H:    17.7544
  Bandwidth B:    28.0281
...
  N Effective Left:    360
  N Effective Right:    323
...

The data are the extract of Cattaneo, Frandsen & Titiunik (2015, doi:10.1515/jci-2013-0010) that ships with R’s rdrobust (the loader keeps its historical name): x is the party’s vote-share margin in the election at time t and y its vote share (0–100) in the election at t+2, following rdrobust’s own illustration. Barely winning at t raises the vote share at t+2 by about 7.4 percentage points (conventional), 7.5 with robust bias correction; the headline line reports the robust bias-corrected estimate and CI. On this data the default MSE-optimal bandwidths, estimates, and standard errors match R rdrobust::rdrobust() and Stata rdrobust (Track A parity module 06_rd). Observations are 1,297 because 93 rows have a missing outcome.

5. Synthetic control: replace Stata/R synth

Question: how did California’s Proposition 99 affect cigarette sales?

import statspai as sp

prop99 = sp.datasets.california_prop99()
sc = sp.synth(
    data=prop99,
    outcome="cigsale",
    unit="state",
    time="year",
    treated_unit="California",
    treatment_time=1989,
)
print(sc.summary())

Result:

==============================================================================
  Synthetic Control Method
==============================================================================

  ATT:      -19.5 *
  Std. Error:  (10.8)
  [95% CI]:    [-58.6,  19.6]
  P-value:     0.077

------------------------------------------------------------------------------
  Detailed Estimates
------------------------------------------------------------------------------
         unit  weight
         Utah  0.3939
      Montana  0.2318
       Nevada  0.2049
  Connecticut  0.1091
New Hampshire  0.0454
     Colorado  0.0148
...

The estimate says California consumed about 20 fewer packs per capita per year after the intervention. The p-value is the in-space placebo rank: California’s post/pre RMSPE ratio ranks 3rd of the 39 states, so p = 3/39 ≈ 0.077. This default matches on pre-treatment outcomes only; pass covariates= (e.g. ["lnincome", "retprice", "age15to24", "beer"]) for an ADH-style predictor specification. That path re-solves the nested V-W problem for every placebo state and is much slower: pass n_jobs=-1 to fit the placebos in parallel (bit-identical results), or placebo=False while iterating on the specification.

Read this number with its caveat, which the full summary also prints: classical SCM weights are often not uniquely identified on empirical data, and different correct solvers can land on different donor weights. StatsPAI’s native solver is certified on uniquely identified designs and labelled identification-dependent elsewhere. On this specification R Synth reaches an ATT of about -19.59 rather than -19.51; pass backend="synth" (needs a local R with the Synth package; outcome-lag specification only) when you need R’s exact numbers.

6. Logit, clustered SEs and post-estimation: Stata’s logit + margins

Question: in Thornton’s Malawi experiment, how much did a randomly offered cash incentive raise the probability that people collected their HIV test results? Villages are the clusters.

import statspai as sp

hiv = sp.datasets.thornton_hiv(complete_case=True)
fit = sp.logit("got ~ C(any) + distvct + male + age", data=hiv,
               vce="cluster villnum")
print(fit.summary())
print(sp.margins(fit, data=hiv, variables=["any"]).round(4))   # margins, dydx(any)

Result:

Model: Logit
Method: Maximum Likelihood (Newton-Raphson)
Dependent Variable: got
...
               Coefficient  Std. Error  z-statistic  P>|z|  [0.025  0.975]
Intercept          -0.6370      0.2010      -3.1691 0.0015 -1.0310 -0.2431
C(any)[T.1.0]       2.0178      0.0994      20.3029 0.0000  1.8230  2.2126
distvct            -0.1696      0.0408      -4.1610 0.0000 -0.2496 -0.0897
male               -0.0530      0.1051      -0.5046 0.6139 -0.2590  0.1529
age                 0.0100      0.0034       2.9046 0.0037  0.0032  0.0167
...

  variable  dy/dx      se        z  pvalue  ci_lower  ci_upper
0    1.any   0.45  0.0217  20.7367     0.0    0.4075    0.4926

As in Stata, the 9 rows with a missing age or village id leave the estimation sample (N = 2,825, 119 villages); the 4 dropped for a missing cluster variable are reported in a StatsPAIWarning and in fit.model_info["n_missing_cluster_dropped"].

The logit coefficient is on the log-odds scale; sp.margins reports the effect on the probability, with a delta-method standard error. any is a 0/1 treatment, so it enters as C(any) (Stata’s i.any) and the margin is the discrete change from 0 to 1. Being offered any incentive raised the probability of collecting results by about 45 percentage points, close to the raw difference in means of 45.1 points. Without C() the column counts as continuous, as in Stata without i., and margins returns the average slope (0.3558). That slope is not the effect of the offer. The coefficients, the village-clustered standard errors, and the marginal effect match Stata 18 logit got i.any distvct male age, vce(cluster villnum) followed by margins, dydx(any) (0.4500223, SE 0.0217017).

Tests and linear combinations use the full covariance matrix and the fit’s own reference distribution (χ² / z here, F / t after OLS):

fit.test("distvct = 0")    # chi2(1) = 17.31, p < 0.001  (Stata: test distvct)
fit.lincom("1.any + male") # 1.9648, SE 0.1412           (Stata: lincom 1.any + male)

If you already have the Stata lines, run them as they are:

ame = sp.stata("""
logit got i.any distvct male age, vce(cluster villnum)
margins, dydx(any)
""", data=hiv)

sp.stata returns the output of the last line — here the same marginal-effects table. A line it cannot translate faithfully (for example xtreg, fe without the panel id that xtset would have supplied) raises instead of running a different model; sp.from_stata(line) shows the Python call without running it.

Export Results

sp.regtable() is the single table builder behind every format. Build the table once, then write it wherever your co-authors need it:

import statspai as sp

card = sp.datasets.card_1995()
m1 = sp.regress("lwage ~ educ", data=card, robust="hc1")
m2 = sp.regress("lwage ~ educ + exper + expersq", data=card, robust="hc1")
m3 = sp.regress("lwage ~ educ + exper + expersq + black + south + smsa",
                data=card, robust="hc1")
m4 = sp.ivreg("lwage ~ (educ ~ nearc4) + exper + expersq + black + south + smsa",
              data=card, robust="hc1")

tbl = sp.regtable(
    m1, m2, m3, m4,
    model_labels=["OLS (1)", "OLS (2)", "OLS (3)", "2SLS (4)"],
    coef_labels={"educ": "Years of schooling", "exper": "Experience",
                 "expersq": "Experience squared", "black": "Black",
                 "south": "South", "smsa": "SMSA"},
    drop=["Intercept"],
    title="Returns to Schooling (Card 1995)",
    notes=["HC1 robust SE. Column (4) instruments schooling with nearc4."],
)
print(tbl)                    # terminal
tbl.to_excel("table1.xlsx")   # Excel
tbl.to_word("table1.docx")    # Word
tbl.to_latex()                # LaTeX source; also .to_markdown(), .to_html()

The images are the .xlsx files written by tbl.to_excel(), rendered with LibreOffice: the Card (1995) table above, and a LaLonde/NSW table regressing 1978 earnings on NSW treatment with a PSID comparison group. The LaLonde table is also a warning about observational comparisons: the treatment coefficient moves from -635 to +1,551 once pre-treatment earnings and demographics are controlled for. See the export guide for journal templates, standard-error formats, and single-model exports.

Interactive Plot Editing

If you miss Stata’s Graph Editor, use sp.interactive(fig) on any matplotlib figure returned by StatsPAI. In Jupyter it opens an editing panel next to a live preview, so beginners can adjust a figure without learning every matplotlib option first. Requires pip install ipywidgets.

What it is for:

• change titles, labels, fonts, colors, markers, line widths, grids, legends (including entry names), axis limits, figure size, and export DPI;
• recolor a series together with its error bars and confidence band, or apply a palette such as okabe_ito (color-blind safe) or grayscale;
• switch among StatsPAI’s publication themes (academic, aea, minimal, cn_journal) and the built-in matplotlib and seaborn styles, for this figure only;
• keep the data layer locked (protect_data=True by default): axis limits never hide a point, error bar or band, and data cannot be made invisible;
• export reproducible Python code for the edits. Tests check, pixel for pixel, that the code reproduces the edited figure and that undo restores the original;
• save journal-ready PDFs with embedded TrueType fonts, and run editor.check() for text that is too small or series that look alike in black and white.

Outside a notebook (scripts, servers, agents) pass the edits as a list: sp.interactive(fig, actions=[{"action": "apply_theme", "theme": "aea"}]). See the figure editor guide.

import statspai as sp

mp = sp.datasets.mpdta()
gt = sp.callaway_santanna(data=mp, y="lemp", t="year",
                          i="countyreal", g="first_treat")
agg = sp.aggte(gt, type="dynamic", bstrap=False)
fig, ax = sp.ggdid(agg)

editor = sp.interactive(fig)   # edit the plot in Jupyter
print(editor.generate_code())  # copy reproducible matplotlib edits

The screenshot above shows the intended workflow: preview on one side, editing controls on the other, and code export for reproducibility.

Everyday Workflow

import statspai as sp

card = sp.datasets.card_1995()
r1 = sp.regress(
    "lwage ~ educ + exper + expersq + black + south + smsa",
    data=card,
    robust="hc1",
)
r2 = sp.ivreg("lwage ~ (educ ~ nearc4) + exper + expersq + black + south + smsa", data=card)

print(r1.summary())                          # human-readable table
print(r1.tidy().head())                      # broom-style dataframe
print(r1.test("black = south"))              # Wald test, like Stata's `test`
print(r1.lincom("black - south"))            # like Stata's `lincom`
tbl = sp.regtable(r1, r2, model_labels=["OLS", "2SLS"])
tbl.to_word("table.docx")                    # Word table
tbl.to_excel("results.xlsx")                 # Excel table

Useful docs:

• Getting started and Cookbook
• Choosing an estimator: DiD, IV, RD, matching, synthetic control
• Migrating: from R, Stata/R command translators, shared argument grammar
• Exporting regression tables
• Agents: agent API, MCP workflow for economists
• Evidence: stability and validation tiers, parity matrix, tier census

Using StatsPAI From An Agent

The same registry that powers sp.help() is exposed three ways.

In Python — discover functions and their schemas without reading source:

import statspai as sp

sp.list_functions(core=True)                   # the ~30 everyday verbs, in order
sp.list_functions(category="causal")[:5]      # names
sp.describe_function("rdrobust")               # parameters, aliases, validation
sp.function_schema("rdrobust")                 # JSON schema for tool calling
sp.from_stata("reghdfe y x, absorb(id year) vce(cluster id)")
# {'tool': 'feols', 'python_code': "sp.feols('y ~ x | id + year', data=df, cluster='id')", ...}
sp.stata("regress y x, vce(cluster id)", data=df)   # translate and run

From the shell — statspai list, statspai describe rdrobust, statspai search "synthetic control"; and estimation with the same payload an MCP client gets: statspai run callaway_santanna --data panel.csv --arg y=lemp --arg g=first_treat --arg t=year --arg i=countyreal, or the schema-driven shortcuts statspai did --data panel.csv --y lemp --treat treated --time year --id id. statspai route did --answer design=staggered routes a question to estimators without data; statspai mcp starts the server.

Over MCP — the package installs a statspai-mcp stdio server (pure Python, no extra dependencies). By default it lists the hand-curated estimators, workflow and pipeline tools plus three discovery meta-tools (search_functions → describe_function → call_function) that reach every registered function; statspai-mcp --profile full advertises every auto-generated tool as its own entry (several hundred, about 2 MB, more than a client context window holds; for schema dumps and debugging). Clients that connect to a URL use statspai-mcp --transport http (binds 127.0.0.1:8765/mcp by default; any other address needs STATSPAI_MCP_HTTP_TOKEN). It also ships workflow prompts (for example audit_did_result, stata_command_workflow) and resources such as statspai://catalog. Tools take a data_path (CSV, Stata .dta, and other formats pandas can read) and return structured JSON with data provenance, captured runtime warnings and any unsupported arguments. For Claude Code:

claude mcp add statspai -- statspai-mcp

For Claude Desktop, Cursor, and other clients:

{
  "mcpServers": {
    "statspai": { "command": "statspai-mcp", "args": [] }
  }
}

See the MCP workflow guide for data handoff, result handles, and the recommended detect → estimate → audit loop.

As a Claude Code skill — statspai skill install copies the packaged statspai-analysis skill (applied-econ / epidemiology / ML-causal pipelines with paper-ready exports) into ~/.claude/skills; statspai skill validate re-checks every API claim in it against the installed package, and with --run executes its pipelines on synthetic data. The installed skill is a copy and does not follow package upgrades: statspai skill status reports whether it is current, statspai skill install --force refreshes it. Agents using the library start at AGENTS.md and llms.txt.

Validation: What Has Been Checked, And What Has Not

StatsPAI has a large API surface, so validation status matters.

import statspai as sp

print(sp.describe_function("ivreg")["validation_status"])   # 'certified'
print(sp.list_functions(validation_status="certified")[:5])

Every registered function carries one of these tiers (counts at the last release; sp.list_functions(validation_status=...) gives the live numbers):

In other words, 1,025 of the 1,680 registered functions carry numerical evidence. Breadth is not the same as validation; check the tier of the functions you depend on.

The evidence scorecard puts the four checks a specialist makes on one page, each regenerated from committed artifacts: agreement with R / Stata on the same data, confidence interval coverage in Monte Carlo, speed against the reference on the same task, and the estimators that still have no evidence.

Cross-language parity, made queryable

The tiers above are derived from an auditable parity index: every verified function records what it was aligned against, to what tolerance, on which test, and how closely it matched. Each row traces to a committed test artifact (the pinned StatsPAI ↔ R ↔ Stata harness, version-locked via renv.lock + per-run provenance) — nothing is asserted from memory.

import statspai as sp

s = sp.parity_status("feols")
print(s)
# feols: bit-exact vs fixest::feols [py/R/Stata] (headline rel_est 5.2e-15 within rel_est<=1e-06, rel_se<=1e-06)
s["reference_versions"]          # {'R': 'R version 4.5.2 (2025-10-31)', 'fixest': '0.14.0'}

sp.parity_summary()              # coverage counts, including the unverified gap
sp.parity_matrix(status="bit-exact")

Grades: bit-exact (headline relative error ≤ 1e-6 against a named R/Stata reference), aligned (a documented, pre-registered looser tolerance), analytical-only (recovers a known DGP truth or closed-form identity), external-replication (reproduces published-paper numbers), and unverified (registered but no parity evidence attached yet — the honest gap). The full, auto-generated matrix is published at docs/parity.md.

For your own data, sp.cross_validate re-runs one estimand through every independent engine installed locally and reports whether they agree:

card = sp.datasets.card_1995()
cv = sp.cross_validate(card, "iv", y="lwage", endog=["educ"], instruments=["nearc4"],
                       covariates=["exper", "expersq", "black", "south", "smsa"])
print(cv.summary())

Engine              Estimate     Std.Err                95% CI    status
------------------------------------------------------------------------
statspai             0.13229     0.04923      [0.0358, 0.2288]        ok
pyfixest             0.13229     0.04923      [0.0358, 0.2288]        ok
linearmodels         0.13229     0.04923      [0.0358, 0.2288]        ok
R::fixest            0.13229     0.04923      [0.0358, 0.2288]        ok
------------------------------------------------------------------------
VERDICT: ✓ AGREE   (4/4 engines ran)

Engines that are not installed (pyfixest, or R with fixest) are skipped. Every engine is asked for the same variance estimator (including cluster= and vcov=) and the same small-sample convention, so the standard errors are compared as well as the point estimates.

Beyond point-parity, a Track-B coverage study runs B=1000 Monte Carlo replications per estimator and checks that 95% confidence intervals hit their nominal rate on known-truth DGPs, against a 99% Wilson acceptance band of [0.935, 0.967]. The thirteen materialized nominal rows (twelve known-truth DGPs) — OLS on an RCT (0.952), a 2×2 DiD (0.955), strong-instrument IV (0.962), Callaway–Sant’Anna staggered ATT (0.947), Sun–Abraham overall ATT (0.950), a two-way FE panel through sp.panel (0.948) and through sp.fast.feols (0.955), entropy balancing (0.945), a causal-forest AIPW ATE at 2,000 trees (0.959), DML IRM ATE (0.968), SDID with placebo SEs (0.928), sharp RD with the robust CI (0.934), and DML PLR with the default learners (0.883) — each also record bias, Monte Carlo SD and SE calibration. The SDID and RD shortfalls come with calibrated SEs (RD’s intervals equal R rdrobust’s draw by draw); the PLR shortfall is regularisation bias of the default gradient-boosting nuisances (0.95 with the true nuisances); see tests/coverage_monte_carlo/FINDINGS.md. The committed artifacts live under tests/coverage_monte_carlo/results_b1000/.

Life Cycle

StatsPAI is in active development (PyPI status “Beta”) and uses semantic versioning. Releases add functions and optional parameters, but a stable function is never broken without notice. The old behaviour first warns (DeprecationWarning, or FutureWarning for a changing default) for at least one minor release, and every such change is listed in MIGRATION.md. A fix that changes returned numbers is flagged as a correctness fix there and in the changelog. The 32 core entry points (sp.list_functions(core=True)) have snapshot-tested signatures. A function marked experimental may change in any minor release and makes no numerical-validation claim; the tiers above say how far each function’s numbers have been checked. The glossary defines these terms.

Changelog

Release notes live outside the README:

• CHANGELOG.md for the full version history.
• MIGRATION.md for deprecations and correctness fixes that change numbers.
• Docs changelog page for the rendered documentation site.

The README is intentionally focused on first-time users.

Paper

StatsPAI is described in a peer-reviewed paper in the Journal of Open Source Software (2026, 11(125), 10604): https://doi.org/10.21105/joss.10604. The reviewer-facing material prepared for that review remains the quickest way to audit the package:

• JOSS reviewer guide
• JOSS validation dossier
• Design rationale and FAQ
• Examples
• Contributing
• Support

Citation

If you use StatsPAI in research, cite the JOSS paper (preferred) and the underlying method papers for each estimator. sp.citation() returns the paper citation, sp.citation(which="software") the versioned software entry, and many result objects expose estimator-level citation helpers.

_{The JOSS article page, published 3 September 2026}

@article{wang2026statspaijoss,
  author  = {Wang, Biaoyue and Rozelle, Scott},
  title   = {StatsPAI: A Unified, Agent-Native Python Toolkit for
             Causal Inference and Applied Econometrics},
  journal = {Journal of Open Source Software},
  year    = {2026},
  volume  = {11},
  number  = {125},
  pages   = {10604},
  doi     = {10.21105/joss.10604},
  url     = {https://doi.org/10.21105/joss.10604}
}

@software{wang2026statspai,
  author  = {Wang, Biaoyue and Rozelle, Scott},
  title   = {StatsPAI: A Unified, Agent-Native Python Toolkit for
             Causal Inference and Applied Econometrics},
  year    = {2026},
  version = {1.29.0},
  doi     = {10.5281/zenodo.19933900},
  url     = {https://doi.org/10.5281/zenodo.19933900},
  license = {MIT}
}

License

MIT. See LICENSE.