StatsPAI Grammar — learn the keywords once¶

Stata's real genius is not that it has many commands — it's that its grammar is consistent: , vce(cluster id), , robust, weights, margins, estat mean the same thing across two hundred commands. R is more fragmented, but each package is internally self-consistent.

StatsPAI exposes 1,100+ flat functions under a single import statspai as sp. That breadth is a strength only if the grammar stays consistent — otherwise the same idea spelled five different ways (robust / vcov / vce / se_type) makes a thousand functions feel like a thousand dialects. This page is the contract that keeps the grammar honest, and the two CI gates that enforce it.

The pitch vs. R: in R you hand-assemble sandwich + clubSandwich + fwildclusterboot + lmtest to get a full SE menu, and each has its own calling convention. StatsPAI's goal is that the entire SE menu (cluster / wild bootstrap / CR2–CR3 / multiway / Conley) is the same, consistently spelled parameter on every estimator. We are not all the way there yet (see the matrix below) — and we track exactly how far, in code.

1. Canonical keyword vocabulary¶

The single source of truth is statspai/_house_style.py. The canonical spellings, ratified 2026-06-30:

Concept	Canonical	Accepted aliases (forwarded)	Notes
Standard-error / variance type	`vce`	`robust`, `vcov`, `se_type`	Matches Stata `vce()`. `vcov` stays a permanent alias on the `fixest` family (pyfixest/R parity).
Outcome / dependent variable	`y`	`outcome`, `depvar`, `dependent`	EconML-mirrored `Y` is kept verbatim on the CATE/forest family.
Treatment variable	`treat`	`treatment`, `treatvar`	EconML-mirrored `T`/`W` kept verbatim on the CATE/forest family.
Clustering variable	`cluster`	`cluster_var`	`clusters` (plural) is reserved for multiway helpers that take a list.
Observation weights	`weights`	`weight`, `sample_weight`	`w`/`W` are not weights (spatial weight matrix / GMM weighting matrix).
Input dataframe	`data`	`df`, `frame`	Already near-universal (486 uses).
Panel identifier	`id`	`unit`, `i`, `group` (where it means the unit)	Ratified 2026-09-15. Matches `sp.did(id=)` and Stata `xtset id`.
Time / period	`time`	`t`
Adjustment covariates	`covariates`	`controls`, `covs`, `X`, `x` (where it means covariates)	Ratified 2026-09-15. In the RD family `x` stays the running variable.
Adoption cohort	`first_treat`	`g`, `cohort`	The period a unit is first treated — not a 0/1 indicator.

On the causal entry points the native spelling (which mirrors the reference package: i/t/g/x in R did, group/controls in Stata did_imputation, x/c/covs in rdrobust) stays in the signature, and the house-style spelling is accepted as well. Both spellings in one call raise TypeError; a misspelt keyword names the closest real parameter. The accepted aliases of any function are listed by sp.describe_function(name)["aliases"].

sp.callaway_santanna(df, y="y", id="county", time="year", first_treat="g0",
                     covariates=["lpop"])       # == i=, t=, g=, x=
sp.bjs(df, y="y", id="county", time="year", first_treat="g0")  # == group=
sp.rdrobust(df, y="vote", running="margin", cutoff=0)          # == x=, c=
sp.psm(df, y="earn", treat="train", covariates=["age", "educ"])  # == d=, X=

Where to start, and running Stata lines as they are¶

sp.list_functions(core=True) lists the ~30 everyday verbs (regression, test / lincom / margins, regtable, and the DiD / RD / synth / matching / DML headliners) — the rest of the 1,100+ names are specialised.

sp.stata runs the Stata lines you already know and returns the same object as the equivalent sp.* call:

fit = sp.stata("regress y x1 x2, vce(cluster firm)", data=df)
sp.stata("""
    logit yb x1 x2, vce(robust)
    margins, dydx(*)
""", data=df)                          # margins of the logit above
sp.stata("lincom x1 - x2", result=fit)

A line whose translation is incomplete raises rather than running a different model — xtreg y x, fe needs the panel id from an earlier xtset, so it raises with that note instead of fitting pooled OLS; write it as sp.feols("y ~ x | id", data=df). sp.from_stata(line) shows the translation without running it.

False friends (deliberately not unified)¶

Some parameters look like an inference keyword but mean something orthogonal — unifying them would be a correctness bug, not a consistency win:

w / W — spatial weight matrix (spatial econometrics), GMM weighting matrix, do-calculus node set, network adjacency matrix.
cov_type — random-effects covariance structure in mixed models ('unstructured'), not an SE type.
group — a comparison/cohort group in decompositions and DiD, not clustering.

The lint excludes these via an explicit allowlist so its signal stays trustworthy.

Accepting the canonical spelling today (additive, reversible)¶

Estimators opt in via the accepts_aliases decorator, which forwards the canonical spelling to the existing parameter without renaming it — old call sites keep working unchanged:

import statspai as sp

sp.regress("y ~ x", df, vce="hc1")      # canonical — accepted today
sp.regress("y ~ x", df, robust="hc1")   # existing spelling — unchanged
# coefficients and SEs are bit-identical between the two

This is strictly additive: no defaults change and no warnings are raised.

Signatures are not renamed (decision, 2026-09-17). Where a parameter mirrors the reference implementation — i / t / g / x in R did, group / controls in Stata did_imputation, x / c / covs in rdrobust, robust= on Stata-style estimators — the signature keeps that name, because it is what a user porting code from R or Stata already writes, and a rename would break every existing call. The house-style spelling is accepted alongside it through @accepts_aliases. Revisit only if a parameter has no reference-package counterpart, and then through a DeprecationWarning and a MIGRATION.md entry with at least one minor release of overlap.

The reviewer's sharpest question: is wild cluster bootstrap usable on any estimator, or only on feols? The honest answer, tracked in scripts/se_menu_matrix.py:

estimator	classical	hc_robust	cluster	twoway	cr2_cr3	wild_cluster_boot	conley	jackknife
`feols`	✓	✓	✓	✓	✓	✓	✓	✓
`hdfe_ols`	✓	✓	✓	✓	✓	✓	✓	✓
`fepois`	✓	✓	✓	✓	✓	✓	✓	✓
`feglm`	✓	✓	✓	✓	✓	✓	✓	✓
`regress`	✓	✓	✓	✓	✓	✓	✓	✓
`ivreg`	✓	✓	✓	✓	✓	✓	✓	✓
`ppmlhdfe`	✓	✓	✓	✓	✓	✓	✓	✓
`panel`	✓	✓	✓	✓	✓	✓	✓	✓
`callaway_santanna`	·	·	✓	·	·	✓	·	·
`did`	·	·	✓	·	·	✓	·	·
`dml`	✓	·	·	·	·	·	·	·
`rdrobust`	·	✓	✓	·	·	·	·	·
`aipw`	✓	✓	✓	·	·	·	·	·
`ipw`	·	✓	✓	·	·	·	·	·
`tobit`	✓	✓	✓	·	·	·	·	·
`heckman`	✓	✓	✓	·	·	·	·	·
`synth`	·	·	·	·	·	·	·	✓

Legend: ✓ native (estimator parameter, correct for its FE/IV/GLM structure) · ○ standalone (correct, via an sp.* SE function reading stored design/residuals) · ⚠ standalone-unsafe (a standalone function that re-parses the formula and refits plain OLS — wrong for FE/IV) · · n/a.

What the matrix makes explicit today:

regress has a full, native SE menu (8/8) — all the standalone SE options are now vce= parameters on sp.regress itself:

sp.regress("y ~ x", df, vce="CR2", cluster="firm")       # bias-reduced CR2
sp.regress("y ~ x", df, vce="CR3", cluster="firm")       # bias-reduced CR3 (== jackknife)
sp.regress("y ~ x", df, vce="wild", cluster="firm")      # WCR cluster bootstrap
sp.regress("y ~ x", df, cluster=["firm","year"])          # CGM-2011 two-way
sp.regress("y ~ x", df, vce="conley",                    # Conley spatial HAC
               conley_lat="lat", conley_lon="lon", conley_cutoff=200.0)

Validated vs Stata reghdfe vce(cluster firm year), acreg ... spatial dist(), and R sandwich::vcovCL(HC2/3) to machine precision on the same 400-obs panel (see tests/reference_parity/test_ols_se_external_parity.py). - Wild cluster bootstrap is native on feols via vce="wild" (and on the panel hdfe_ols path). It runs the WCR bootstrap on the FE-absorbed within design:

sp.feols("y ~ x | firm", data=df, vce="wild", cluster="firm")
# point estimates = feols; p-values & CIs = wild cluster bootstrap

Externally validated against Stata boottest (David Roodman's canonical implementation): on an identical 600-obs / 15-cluster panel the point estimate and CRV1 cluster SE match reghdfe to ~1e-9, and the wild p-value matches boottest's exact 2¹⁵ Rademacher enumeration to Monte-Carlo error (0.0265 vs 0.0263). See tests/reference_parity/test_feols_wild_boottest_parity.py. - ivreg has native wild via vce="wild" — the WRE bootstrap (Davidson-MacKinnon 2010), the IV-correct procedure that resamples both the structural and reduced-form residuals and refits 2SLS:

sp.ivreg("y ~ w + (d ~ z1 + z2)", data=df, vce="wild", cluster="firm")
# the endogenous coefficient gets a wild-bootstrap p-value & CI

Validated against Stata boottest after ivreg2 across strong-IV (p 0.2016 vs 0.20155) and weak-IV (p 0.3415 vs 0.3412) regimes — the weak-IV case rules out the naive reduced form (which gives 0.426). Also validated for two endogenous regressors (p 0.2101/0.0151 vs boottest 0.2108/0.0141). See tests/reference_parity/test_iv_wild_boottest_parity.py. - ivreg has native two-way clustering via cluster=["a", "b"] — the CGM (2011) inclusion-exclusion sandwich on the projected regressors, matching Stata ivreg2, cluster(a b) small. - ivreg has native bias-reduced cluster SEs via vce="CR2" (Bell- McCaffrey) and vce="CR3" (== vce="jackknife") — the Pustejovsky-Tipton (2018) adjustment on the projected 2SLS regressors, matching R clubSandwich::vcovCR(ivreg, type=...) to machine precision. - ivreg has native Conley spatial HAC via vce="conley", conley_lat=, conley_lon=, conley_cutoff= — the spatial kernel on the projected 2SLS scores with Stata acreg's planar distance (111 km/deg, cos(lat) longitude), matching acreg ... spatial exactly. - ivreg's SE menu is now complete (8/8 native): every cell — classical / HC / cluster / two-way / CR2-CR3 / wild bootstrap / Conley / jackknife — is a native, externally-validated option. No ⚠ cells remain in the whole matrix. - feols has native bias-reduced / spatial SEs via vce="CR2", vce="CR3" (== vce="jackknife") and vce="conley" on the FE-absorbed within design. The within-transform's leverage adjustment reproduces R clubSandwich::vcovCR(plm, model="within", type=...) to machine precision, so the entire feols row is 8/8 native (tests/reference_parity/test_feols_bias_reduced_parity.py). - panel(method="fe") mirrors that menu. The one-way entity fixed-effects estimator accepts vce="CR2"/"CR3"/"jackknife", vce="conley" (with conley_lat=/conley_lon=/conley_cutoff=) and two-way clustering via cluster=["a", "b"]. Since OLS on the entity-demeaned design reproduces the linearmodels FE coefficients, the CR2/CR3 SEs equal the identical clubSandwich::vcovCR(plm) anchor as feols, and Conley / two-way match sp.regress on the hand-demeaned data (Stata acreg / reghdfe conventions) — see tests/reference_parity/test_panel_bias_reduced_parity.py.

sp.panel(df, "y ~ x", entity="firm", time="t", method="fe",
         vce="CR2", cluster="firm")                 # bias-reduced CR2
sp.panel(df, "y ~ x", entity="firm", time="t", method="fe",
         cluster=["firm", "region"])                # two-way cluster

- fepois / feglm have native bias-reduced cluster SEs via vce="CR2", vce="CR3" (== vce="jackknife"). The GLM adjustment is the IRLS-weighted generalisation of CR2 (d = dμ/dη, V = Var(μ), working weight w = d²/V) on the fixed-effects-as-dummies design, matching R clubSandwich::vcovCR(glm, type=...) for the Poisson / logit / probit / gaussian families. (Unlike OLS, the weighted projection does not carry the CR2 leverage through FE absorption, so the dummy design is required — it is guarded against high-dimensional FE, which should use cluster= or the wild bootstrap instead.) See tests/reference_parity/test_feglm_bias_reduced_parity.py.

sp.fepois("y ~ x1 + x2 | firm", data=df, vce="CR2", cluster="clu")
sp.feglm("y ~ x1 | firm", data=df, family="logit", vce="CR3", cluster="clu")

- fepois / feglm also have a native wild cluster bootstrap via vce="wild" — the restricted score wild cluster bootstrap (Kline-Santos 2012), the method Stata boottest runs after poisson / logit, implemented with boottest's exact studentization (cluster-share-centered CRVE denominator + strict exceedance counting, reverse-engineered from its enumerated bootstrap distribution and corroborated against boottest.mata). In the enumerated regime (2^G <= wild_reps) the p-values are bit-exact matches to Stata boottest, verified for Poisson and logit. It returns wild-bootstrap p-values per coefficient (SE/CI stay cluster-robust). So the answer to "is wild cluster bootstrap usable on any estimator or only feols?" is now: regress, feols, ivreg (WRE), fepois, feglm all have a native, externally-validated vce="wild".

sp.fepois("y ~ x1 + x2 | firm", data=df, vce="wild", cluster="clu")

- ppmlhdfe has native two-way clustering via cluster=["a", "b"] — the CGM (2011) inclusion-exclusion sandwich on the FE-residualised PPML design with the single G_min/(G_min-1) small-sample factor, byte-identical to Stata ppmlhdfe ..., cluster(a b) (its one-way path already matches Stata exactly). See tests/reference_parity/test_ppmlhdfe_twoway_parity.py.

sp.ppmlhdfe("y ~ x1 + x2 | o + d", data=df, cluster=["origin", "dest"])

- ppmlhdfe gains the boottest-convention wild bootstrap at scale, plus CR2/CR3. Stata's boottest cannot run after ppmlhdfe at all (no constraints() support — verified empirically), so ppmlhdfe(vce="wild") is a beyond-Stata capability: the score bootstrap touches only per-cluster one-step contributions (never per-observation leverage), and their weighted Frisch-Waugh-Lovell reduction onto the FE-absorbed design is exact (verified to 1e-17) — so it runs on 500-level FE in under a second and is byte-identical to fepois(vce="wild") (itself bit-exact vs boottest) on low-dimensional FE. vce="CR2"/"CR3" use the reference-matching dummy design with the same high-dim guard as fepois. - fepois / feglm gain vce="conley" — GLM spatial HAC referenced to R conleyreg (which natively supports Poisson/logit/probit; Stata acreg is OLS/2SLS-only). The implementation reproduces conleyreg's spherical uniform kernel (haversine, earth radius 6371.01 km, read from its C++ source) to ~1e-7. Note the documented convention split: the OLS menu follows acreg's planar distance, the GLM menu follows conleyreg's spherical distance — each the convention of the reference that exists. - hdfe_ols and panel complete their rows (8/8). hdfe_ols gains the canonical vce= menu on its absorber's within design: vce="robust"/"hc1" (HC1 with reghdfe's N/(N-k-df_a) factor — matches Stata reghdfe, vce(robust) exactly), vce="CR2"/"CR3"/"jackknife" (same frozen clubSandwich plm anchor as feols/panel), vce="conley", and vce="wild" as shorthand for its native WCR bootstrap. panel(method="fe") gains vce="wild" — the same WCR engine as sp.regress(vce="wild") on the entity-within design, byte-identical p-values to regress on the hand-demeaned data with the same seed. Every regression-family estimator (regress, feols, ivreg, panel, hdfe_ols) is now 8/8 native.

sp.hdfe_ols("y ~ x | firm + year", data=df, vce="robust")        # reghdfe HC1
sp.panel(df, "y ~ x", entity="firm", time="t", method="fe",
         vce="wild", cluster="firm", seed=42)                    # WCR bootstrap

DiD gets the canonical few-clusters inference. did(method="2x2", vce="wild", cluster=...) runs the WCR wild cluster bootstrap on the interaction regression — the setting the wild-bootstrap literature was built for (MacKinnon-Webb 2017) — validated against Stata boottest after reg y d t dt, vce(cluster ...) and byte-identical to sp.regress(vce="wild") on the same design. classical/hc_robust for DiD stay deliberately absent (unclustered DiD SEs are the Bertrand-Duflo-Mullainathan mistake). For Callaway-Sant'Anna, the Mammen multiplier bootstrap on unit-level influence functions (sp.aggte) is the wild bootstrap for that estimator (CS 2021 §4.1; R did::mboot), so its cell is native by construction. ppmlhdfe also completes vce="conley" via the conleyreg-referenced GLM spatial HAC. Remaining blanks (synth, rdrobust, dml rows) are documented in scripts/se_menu_matrix.py as methodologically legitimate n/a.
Treatment effects and limited dependent variables take Stata's weights and clusters (1.32). sp.aipw / sp.ipw accept weights= (Stata [pw=]) and cluster=; sp.tobit and sp.heckman(method="ml") accept vce=, cluster= and weights=. Each is pinned against Stata 18 (teffects aipw / teffects ipw / tobit / heckman) to 1e-10 or better in tests/reference_parity/test_teffects_design_stata_parity.py and test_ldv_design_stata_parity.py:

sp.aipw(df, y="y", treat="d", covariates=["x1", "x2"], weights="w", cluster="g",
        cross_fit=False, se_method="sandwich")      # teffects aipw [iw=w], vce(cluster g)
sp.ipw(df, "y", "d", ["x1", "x2"], estimand="ATT", weights="w",
       se_method="sandwich")                        # teffects ipw ... [pw=w], atet
sp.tobit(df, "y", ["x1", "x2"], ll=0, vce="cluster g")
sp.heckman(df, y="y", x=["x1"], select="s", z=["x1", "z"], method="ml", weights="w")

One convention to know: teffects aipw refuses [pw=], and its [iw=], vce(robust) treats the weights as frequencies (effective N = sum of the weights). sp.aipw(weights=) is the sampling-weight SE, which Stata reports under vce(cluster).

This is the gap the SE-menu wiring work closes, estimator by estimator. The matrix is the scoreboard: the CI gate ratchets the native count up and the unsafe count down, so coverage can only improve.

3. The CI gates¶

Two --check gates keep the grammar from drifting:

python scripts/signature_house_style.py --check   # keyword-spelling ratchet
python scripts/se_menu_matrix.py --check          # SE-menu coverage ratchet

signature_house_style introspects every public callable and counts legacy-spelling sites per theme against a frozen baseline. A new function using robust= where vce= is canonical raises the count and fails the gate. False friends are excluded by allowlist, and a legacy-named parameter whose canonical spelling the function accepts through @accepts_aliases counts as converged — the gate measures what a caller can type, so removing an alias fails it too.
se_menu_matrix validates that every estimator and SE function it references still resolves under sp.* (drift guard), and ratchets the native-cell count up / unsafe-cell count down.

Both are ratchets, not hard zeros: the surface has historical drift, so the gates block regressions while letting convergence land incrementally — the counts only move in the improving direction.