statspai.fixest

fixest

High-dimensional fixed effects estimation via pyfixest.

This module provides thin wrappers around pyfixest's estimation functions, converting results into StatsPAI's EconometricResults for seamless integration with outreg2 and the rest of the ecosystem.

Requires: pip install pyfixest

Examples:

>>> import os
>>> import tempfile
>>> import numpy as np
>>> import pandas as pd
>>> from statspai.fixest import feols, fepois
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> df = pd.DataFrame({'firm': rng.integers(0, 20, n),
...                    'year': rng.integers(2000, 2006, n),
...                    'experience': rng.normal(10, 3, n),
...                    'rd_spending': rng.normal(5, 1, n)})
>>> df['wage'] = 1 + 0.5 * df['experience'] + rng.normal(size=n)
>>> df['patents'] = rng.poisson(np.exp(0.2 * df['rd_spending'] - 0.5))
>>>
>>> # Two-way fixed effects with clustered SEs
>>> # (pyfixest's first call JIT-compiles for ~20 s, hence the skips)
>>> result = feols("wage ~ experience | firm + year",
...               data=df, vcov={"CRV1": "firm"})
>>> report = result.summary()
>>>
>>> # Poisson regression
>>> result = fepois("patents ~ rd_spending | firm", data=df)
>>>
>>> # Works with outreg2
>>> from statspai import outreg2
>>> path = os.path.join(tempfile.mkdtemp(), "table.xlsx")
>>> outreg2(result, filename=path)
Regression results exported to: ...table.xlsx

etable

etable(*results: EconometricResults, **kwargs: Any) -> Any

Display a pyfixest-style regression table for StatsPAI results.

If every result carries a _pyfixest_fit reference, uses pyfixest's native etable. Otherwise returns one row per term with "coef*** (se)" cells, formatted by :func:sp.regtable so both tables report identical numbers.

Parameters:

Name	Type	Description	Default
*results	EconometricResults	One or more fitted results.	()
**kwargs	Any	Passed to pyfixest.etable().	{}

Returns:

Type	Description
str or DataFrame

Examples:

Side-by-side comparison of nested specifications:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n = 200
>>> df = pd.DataFrame({"x1": rng.normal(size=n), "x2": rng.normal(size=n),
...                    "firm": rng.integers(0, 10, n)})
>>> df["wage"] = 1.0 + 0.5 * df["x1"] - 0.3 * df["x2"] + rng.normal(0, 0.5, n)
>>> r1 = sp.regress("wage ~ x1", data=df)
>>> r2 = sp.regress("wage ~ x1 + x2", data=df)
>>> tab = sp.etable(r1, r2)  # DataFrame: one coefficient column per model
>>> list(tab.columns)
['(1)', '(2)']

With sp.feols fits (requires the fixest extra), pyfixest's native styled table is returned instead:

>>> f1 = sp.feols("wage ~ x1 | firm", data=df)
>>> f2 = sp.feols("wage ~ x1 + x2 | firm", data=df)
>>> tab = sp.etable(f1, f2)

A list works too: sp.etable([r1, r2]).

feglm

feglm(fml: str, data: DataFrame, family: str = 'gaussian', vcov: Optional[Union[str, Dict[str, str]]] = None, *, cluster: Optional[str] = None, wild_reps: int = 9999, wild_weight_type: str = 'rademacher', seed: Optional[int] = None, conley_lat: Optional[str] = None, conley_lon: Optional[str] = None, conley_cutoff: Optional[float] = None, backend: str = 'auto', bias_correction: bool = False, **kwargs: Any) -> Union[EconometricResults, List[EconometricResults]]

Estimate GLM (logit, probit, negative binomial, Gaussian) with high-dimensional fixed effects.

Parameters:

Name	Type	Description	Default
fml	str	pyfixest formula.	required
data	DataFrame	Input dataset.	required
family	str	GLM family: "gaussian", "logit", "probit", or "negbin": the negative binomial (NB2, variance mu + mu^2 / theta) with the fixed effects estimated as parameters, as fixest::fenegbin fits it (not Stata's conditional xtnbreg, fe). theta is estimated jointly; it and its standard error are in model_info['theta'] / model_info['theta_se'] (fixest lists it as the coefficient .theta). Rows its effects fit perfectly (groups with only zero outcomes, singletons) are dropped, as fixest does; weights= are likelihood weights. Native engine only: the vce= options that go through pyfixest are refused for it.	"gaussian"
vcov	str or dict	Variance-covariance estimator (vce= is the canonical alias). Also accepts vce="CR2"/"CR3"/"jackknife" (with cluster=) for the clubSandwich bias-reduced cluster-robust SEs, and vce="wild" (with cluster=) for the restricted score wild cluster bootstrap (Kline-Santos 2012; bit-exact vs Stata boottest in the enumerated regime).	None
cluster	str	Cluster id column for the extended vce= menu (also a shorthand for one-way {"CRV1": cluster}).	None
wild_reps	int	Replications for vce="wild". When 2**G <= wild_reps the full Rademacher grid is enumerated (deterministic).	9999
wild_weight_type	str	Wild weight distribution ("rademacher" or "webb").	"rademacher"
seed	int	RNG seed for sampled (non-enumerated) vce="wild" draws.	None
backend	(auto, pyfixest, native)	'auto' fits family='logit' / 'probit' on StatsPAI's own IRLS engine whenever the call is y ~ x1 + ... | fe1 + ... with numeric regressors, vcov in iid / hetero / {'CRV1': 'a'} / {'CRV1': 'a + b'} and at most weights= among the extra keywords, and through pyfixest otherwise (family='gaussian' always). The native fit follows fixest::feglm: observations in a fixed-effect group whose outcomes are all 0 or all 1 are dropped, p-values are normal, and IRLS runs to a relative deviance change of 1e-10. 'pyfixest' always delegates; 'native' refuses calls the native engine cannot fit.	'auto'
bias_correction	bool	With family='logit' / 'probit' and unit (and period) effects, correct the slopes for the incidental-parameter bias (order 1/T from unit effects, 1/N from period effects) by the analytical correction for strictly exogenous regressors, as R alpaca::biasCorr(L = 0) does, and report the variance at the corrected slopes without small-sample factors (alpaca's convention): vcov 'iid', 'hetero' or one-way {'CRV1': c}. The uncorrected slopes and the shift are in model_info['coef_uncorrected'] / ['bias_shift']. Unweighted fits with one or two absorbed effects; native engine. The correction assumes regressors uncorrelated with past and future errors (no lagged outcomes). Bib key: fernandezval2016individual.	False
**kwargs	Any	Additional arguments passed to pyfixest.feglm().	{}

Returns:

Type	Description
EconometricResults or list of EconometricResults

Examples:

Logit GLM with firm fixed effects:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 500
>>> firm = rng.integers(0, 8, n)
>>> x1 = rng.normal(size=n)
>>> p = 1.0 / (1.0 + np.exp(-(0.2 + 0.8 * x1)))
>>> y = (rng.random(n) < p).astype(int)
>>> df = pd.DataFrame({"y": y, "x1": x1, "firm": firm})
>>> res = sp.feglm("y ~ x1 | firm", data=df, family="logit")
>>> "x1" in res.params.index
True

feols

feols(fml: str, data: DataFrame, vcov: Optional[Union[str, Dict[str, str]]] = None, *, weights: Optional[str] = None, ssc: Optional[Any] = None, fixef_rm: Optional[str] = 'singleton', collin_tol: float = 1e-06, lean: bool = False, cluster: Optional[str] = None, wild_reps: int = 999, wild_weight_type: str = 'rademacher', seed: Optional[int] = None, conley_lat: Optional[str] = None, conley_lon: Optional[str] = None, conley_cutoff: Optional[float] = None, backend: str = 'auto', weights_as: Optional[str] = None, **kwargs: Any) -> Union[EconometricResults, List[EconometricResults]]

Estimate OLS / IV with high-dimensional fixed effects.

Uses the Frisch-Waugh-Lovell theorem for fast absorption of high-dimensional fixed effects.

Parameters:

Name	Type	Description	Default
fml	str	A pyfixest-style formula. Examples: "Y ~ X1 + X2" — plain OLS "Y ~ X1 | firm + year" — two-way fixed effects "Y ~ 1 | firm | X1 ~ Z1" — IV with fixed effects "Y ~ X1 | csw0(firm, year)" — multiple estimations "Y ~ X1 | firm + state[year]" — varying slopes: state fixed effects plus state-specific year trends (state[[year]] for the slopes alone). pyfixest does not absorb these, so such formulas run on StatsPAI's own HDFE kernel (model_info['backend'] == 'statspai-native'), matching R fixest for vcov iid / hetero / {'CRV1': col} and weights=.	required
data	DataFrame	Input dataset.	required
vcov	str or dict	Variance-covariance estimator (vce= is the canonical alias). "iid" — classical "HC1", "HC2", "HC3" — heteroskedasticity-robust {"CRV1": "firm"} — cluster-robust {"CRV1": "firm + year"} — two-way clustering other spellings of the same request give identical results: dict keys CL1 / CR1 / cluster (CL3 / CR3 for CRV3), list values, and the strings "cluster firm", "vce(cluster firm)", "cluster:firm", "cluster(firm)", "~firm". An unknown spelling raises MethodIncompatibility. vce="CR2" / "CR3" / "jackknife" (with cluster=) — Pustejovsky-Tipton bias-reduced cluster-robust on the FE-absorbed within design; matches R clubSandwich::vcovCR(plm). vce="wild" (with cluster=) — WCR wild cluster bootstrap (Cameron-Gelbach-Miller 2008); validated against Stata boottest. vce="conley" (with conley_lat=/conley_lon=/conley_cutoff=) — Conley spatial HAC (Stata acreg planar-distance convention). "DK" / "NW" with vcov_kwargs={'time_id': 't', 'panel_id': 'id', 'lag': L} — Driscoll-Kraay / Newey-West panel HAC through pyfixest (R fixest DK(L) / NW(L) with panel.id = ~id + t); without the ids they raise.	None
weights	str	Column name for regression weights. Alone they are analytic weights (Stata [aw=], R fixest weights=).	None
weights_as	(analytic, sampling)	What the weights are. 'analytic' is what weights= alone already means. 'sampling' (Stata [pw=]) gives the same coefficients with vcov='hetero', unless vcov=, vce= or cluster= asks for another. Saying which also silences the warning about dispersed weights under the classical variance. 'frequency' is refused.	'analytic'
ssc	optional	Small-sample correction via pyfixest.ssc().	None
fixef_rm	(singleton, none)	Singleton fixed-effect groups: "singleton" drops the observations in a fixed-effect group of one, repeatedly until none is left, as R fixest, Stata reghdfe and pyfixest do by default; "none" keeps them. Coefficients are the same either way; the number of observations and the degrees of freedom behind the standard errors are not. .. versionchanged:: 1.40.0 The default was "none" (singletons kept) before 1.40.	'singleton'
collin_tol	float	Collinearity tolerance.	1e-6
lean	bool	If True, drop large intermediate arrays to save memory.	False
cluster	str	Cluster id column for the extended vce= menu; also a shorthand for one-way {"CRV1": cluster}.	None
wild_reps	int	Bootstrap replications for vce="wild".	999
wild_weight_type	str	Wild weight distribution ("rademacher", "webb", "mammen").	"rademacher"
seed	int	RNG seed for vce="wild".	None
conley_lat	str	Coordinate columns (decimal degrees) for vce="conley".	None
conley_lon	str	Coordinate columns (decimal degrees) for vce="conley".	None
conley_cutoff	float	Conley distance cutoff in km for vce="conley".	None
backend	(auto, pyfixest, native)	'auto' fits on StatsPAI's own engine (sp.fast.feols, on the Rust kernels when statspai-hdfe is installed) whenever the call is y ~ x1 + ... | fe1 + ... with numeric regressors and i(factor[, var], ref=value) terms, or the IV form y ~ x | fe | d1 + ... ~ z1 + ... with absorbed effects; vcov in iid / hetero / {'CRV1': 'a'} / {'CRV1': 'a + b'}; no ssc= or extra pyfixest keyword. Multiple estimation (y1 + y2 ~, sw(), csw(), ...) is native when every model is. Other calls go through pyfixest (also when regressors are absorbed or collinear, which fixest drops). model_info['backend'] says which ran. The native IV first-stage F is fixest's ivwald1. The native fit follows fixest's conventions: same coefficients, SEs within the demeaning tolerance, rows with missing values dropped first. 'pyfixest' always delegates. 'native' uses StatsPAI's engine for everything it supports, including varying slopes (sp.hdfe_ols), and refuses the rest.	'auto'
**kwargs	Any	Additional arguments passed to pyfixest.feols().	{}

Returns:

Type	Description
EconometricResults or list of EconometricResults	Single result for simple formulas, list for multiple estimations (e.g. csw0/sw/sw0 syntax).

Examples:

Two-way fixed effects with clustered SEs:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> firm = rng.integers(0, 8, n)
>>> year = rng.integers(0, 5, n)
>>> x1 = rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 + 0.1 * firm + 0.2 * year + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({"y": y, "x1": x1, "firm": firm, "year": year})
>>> res = sp.feols("y ~ x1 | firm + year", data=df, vcov={"CRV1": "firm"})
>>> "x1" in res.params.index
True

Multiple estimation (csw0 returns a list, one fit per FE set):

>>> results = sp.feols("y ~ x1 | csw0(firm, year)", data=df)
>>> summaries = [r.summary() for r in results]

Use with outreg2 to build a regression table:

>>> r1 = sp.feols("y ~ x1", data=df)
>>> r2 = sp.feols("y ~ x1 | firm", data=df)
>>> sp.outreg2(r1, r2, filename="table.xlsx")

fepois

fepois(fml: str, data: DataFrame, vcov: Optional[Union[str, Dict[str, str]]] = None, *, weights: Optional[str] = None, ssc: Optional[Any] = None, fixef_rm: Optional[str] = 'singleton', collin_tol: float = 1e-06, iwls_tol: float = 1e-08, iwls_maxiter: int = 25, cluster: Optional[str] = None, wild_reps: int = 9999, wild_weight_type: str = 'rademacher', seed: Optional[int] = None, conley_lat: Optional[str] = None, conley_lon: Optional[str] = None, conley_cutoff: Optional[float] = None, backend: str = 'auto', **kwargs: Any) -> Union[EconometricResults, List[EconometricResults]]

Estimate Poisson regression with high-dimensional fixed effects.

Parameters:

Name	Type	Description	Default
fml	str	pyfixest formula. E.g. "Y ~ X1 | firm".	required
data	DataFrame	Input dataset.	required
vcov	str or dict	Variance-covariance estimator (vce= is the canonical alias). Besides the pyfixest values ("iid", "HC1", {"CRV1": "firm"}, ...), accepts the extended menu: vce="CR2" / "CR3" / "jackknife" (with cluster=) — clubSandwich glm bias-reduced cluster-robust SEs on the FE-as-dummies design; matches R clubSandwich::vcovCR(glm). vce="wild" (with cluster=) — restricted score wild cluster bootstrap (Kline-Santos 2012) with Stata boottest's exact studentization; bit-exact vs boottest in the enumerated regime.	None
weights	str	Column name for regression weights (not supported with the extended vce= menu).	None
ssc	optional	Small-sample correction.	None
fixef_rm	(singleton, none)	Singleton fixed-effect groups: "singleton" drops the observations in a fixed-effect group of one, as R fixest and pyfixest do by default; "none" keeps them. .. versionchanged:: 1.40.0 The default was "none" (singletons kept) before 1.40.	'singleton'
collin_tol	float	Collinearity tolerance.	1e-6
iwls_tol	float	IWLS convergence tolerance.	1e-8
iwls_maxiter	int	Max IWLS iterations.	25
cluster	str	Cluster id column for the extended vce= menu; also a shorthand for one-way {"CRV1": cluster}.	None
wild_reps	int	Replications for vce="wild". When 2**G <= wild_reps the full Rademacher grid is enumerated (deterministic).	9999
wild_weight_type	str	Wild weight distribution ("rademacher" or "webb").	"rademacher"
seed	int	RNG seed for sampled (non-enumerated) vce="wild" draws.	None
backend	(auto, pyfixest, native)	'auto' fits on StatsPAI's own PPML-HDFE engine (sp.fast.fepois, on the Rust kernels when statspai-hdfe is installed) whenever the call is y ~ x1 + ... | fe1 + ... with numeric regressors, vcov in iid / hetero / {'CRV1': 'a'} / {'CRV1': 'a + b'} and no ssc= or extra pyfixest keyword, and through pyfixest otherwise. Both follow fixest: observations in a fixed-effect group whose outcomes are all zero are dropped, rows with missing values are dropped, singletons are dropped unless fixef_rm='none'. model_info['backend'] says which ran. 'pyfixest' always delegates; 'native' refuses calls the native engine cannot fit.	'auto'
**kwargs	Any	Additional arguments passed to pyfixest.fepois().	{}

Returns:

Type	Description
EconometricResults or list of EconometricResults

Examples:

Poisson regression (PPML-style) with firm fixed effects:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> firm = rng.integers(0, 8, n)
>>> x1 = rng.normal(size=n)
>>> mu = np.exp(0.3 + 0.5 * x1 + 0.1 * firm)
>>> df = pd.DataFrame({"y": rng.poisson(mu), "x1": x1, "firm": firm})
>>> res = sp.fepois("y ~ x1 | firm", data=df)
>>> "x1" in res.params.index
True

Bias-reduced cluster-robust SEs (matches R clubSandwich::vcovCR):

>>> res = sp.fepois("y ~ x1 | firm", data=df, vce="CR2",
...                 cluster="firm")