PyFixest: High-Dimensional Fixed Effects Regression

PyFixest is a high-performance Python library designed for regression analysis with high-dimensional fixed effects, serving as a faithful port of the R fixest package.

Core Estimation Routines:

• OLS: via feols()
• IV (Instrumental Variables): via feols()
• Poisson Regression: via fepois()
• GLM (Generalized Linear Models): via feglm()
• Quantile Regression

Advanced Features:

• Event Studies: Supports linear projections, two-stage difference-in-differences imputation, and fully-saturated event study estimators.
• Inference: Includes cluster robust variance estimators (CRV1, CRV3), wild cluster bootstrap, randomization inference, causal cluster variance estimators, and Romano-Wolf correction for family-wise error rate control.
• Diagnostics & Decomposition: Supports Gelbach's regression decomposition, Lal's event study specification test, and weak instrument diagnostics.

Design Philosophy: PyFixest aims for syntax parity with R's fixest. Users can transition between R and Python seamlessly, as the input arguments and results (point estimates, standard errors, $R^2$) are designed to be identical.

Pyfixest uses a specific vocabulary to distinguish between different stages of data transformation during estimation. This prevents confusion between raw user data, transformed data, and ephemeral solver data.

The performance of fixed-effects solvers in PyFixest depends on the connectivity of the underlying bipartite (or tripartite) graph formed by your fixed effects.

• Use MapDemeaner(backend="rust") (MAP): When your data has a well-connected graph. This occurs with high worker mobility, low assortative matching (random sorting), and cross-cutting fixed effects (e.g., occupations that are not nested within firms or workers). In these cases, MAP is often faster, though LsmrDemeaner might be slightly faster on very simple problems where its setup overhead is negligible.

• Use LsmrDemeaner() (within LSMR): When your data has a sparse/poorly-connected graph. This occurs with low worker mobility, strong assortative matching (high sorting), or nested fixed effects (e.g., occupations that are highly concentrated within specific firms or tied to specific workers). LsmrDemeaner is significantly more robust and maintains stable runtimes in these difficult scenarios where vanilla MAP convergence degrades sharply.

# Use MapDemeaner for well-connected/dense graphs
MapDemeaner(backend="rust")

# Use LsmrDemeaner for sparse/difficult/nested graphs
LsmrDemeaner()

When estimating models with multiple fixed effects, a simple 'within transformation' (subtracting group means) is insufficient because demeaning one set of effects (e.g., workers) leaves traces of the other (e.g., firms) in the residuals. PyFixest provides several algorithms to solve this demeaning problem:

1. Method of Alternating Projections (MAP):

   • Accessed via pf.MapDemeaner(backend="rust").
   • Sequentially sweeps through each fixed effect to demean the target variable.
   • In PyFixest, this implementation does not currently use accelerations (unlike R's fixest), so it may be slower in certain scenarios.

2. LSMR:

   • Accessed via pf.LsmrDemeaner(...).
   • A solver based on the approach used in FixedEffectsModels.jl.

3. Modified LSMR with Schwarz preconditioning (via within):

   • Accessed via pf.LsmrDemeaner(backend="within").
   • This is the PyFixest default.
   • It builds and exploits the block structure of the normal equations to form a high-quality preconditioner, making it highly efficient for complex structures.

Pyfixest implements several frozen, slotted state values to ensure data integrity during the estimation lifecycle.

Key State Domains

• ModelMatrix: Formula-materialized pandas tables that keep dependent, independent, fixed-effect, IV, weight, and offset roles separate.
• ObservationWeights: Canonical user-scale weights (aweights or fweights). Setting values=None enables an allocation-free unweighted path.
• WithinLinearData: Unpremultiplied within-scale response and design arrays.
• WithinIvData: Extends WithinLinearData with instrument and endogenous arrays (specific to IV models).
• GlmWorkingState: Final within-scale working response and design, IRLS working weights, predictors, means, and response/working-residual domains.
• DemeanedData: Array-native cache entries using ordered column names as metadata.

Weighting Semantics

• Analytic weights: The retained row count is the effective sample size.
• Frequency weights: The effective sample size is the sum of the user-scale values.
• Probability weights (pweights): Currently unsupported.
• Fixed-effect projection: May use observation or IRLS weights, but returned arrays remain in the original scale.
• Solver Temporaries: OLS, IV, and IRLS fit primitives create square-root-weighted design and response arrays only as local solver temporaries. They do not persist solver-scale copies of canonical data.

The Selection Matrix defines which tests are required based on the type of change made:

PyFixest's API is heavily inspired by the R package fixest. Users familiar with fixest will find familiar syntax and conventions. Key borrowed features include:

• Formula Syntax: Uses feols(), fepois(), and feglm() for estimation. Supports the i() interaction operator and multiple estimation syntax via sw(), sw0(), csw(), and csw0().
• Small-sample corrections: Controlled via the ssc() function. Default parameters like adj, fixef_K, cluster_adj, and cluster_df mirror fixest exactly.
• Post-estimation: Functions like etable(), coefplot(), iplot(), and coef() mirror fixest output and plotting.
• Variance-Covariance Adjustments: You can adjust the vcov post-estimation by calling the .vcov() method on the results object (e.g., an Feols object).
• Fixed Effects & Predictions: The .predict() and .fixef() methods provide functionality for obtaining fitted values, fixed effects, and linear predictions, mirroring fixest.

The Method of Alternating Projections (MAP) is an iterative algorithm used to solve for fixed effects (FE) coefficients. It works by repeatedly sweeping through each fixed effect and subtracting group means from the current residual.

In terms of the Gramian matrix, this is equivalent to block Gauss-Seidel, where each sweep solves one diagonal block at a time. For a model with worker ($W$), firm ($F$), and year ($Y$) effects, the process is:

1. Start with $r = y$
2. Subtract worker means: $r \leftarrow r - D_W G_{WW}^{-1} D_W^ op r$
3. Subtract firm means: $r \leftarrow r - D_F G_{FF}^{-1} D_F^ op r$
4. Subtract year means: $r \leftarrow r - D_Y G_{YY}^{-1} D_Y^ op r$
5. Repeat steps 2-4 until convergence.

Key Characteristics:

• Efficiency: Each step is individually cheap because the diagonal blocks (e.g., $G_{WW}$) are easy to invert.
• Connectivity Dependency: Convergence speed depends heavily on the 'graph' of connections. In dense graphs (many movers, high connectivity), MAP converges very quickly. In sparse graphs (few movers, workers staying in one firm), MAP convergence stalls because information about fixed effects only 'trickles' through residuals from one sweep to the next.

Pyfixest uses a shared estimation pipeline to ensure consistency across different model types. The process follows a structured flow from configuration to result generation:

1. Configuration: An EstimationConfig is built.
2. Parsing: parse_formula expands multiple-estimation syntax.
3. Preparation: runner.run_estimation and fit_one prepare the model matrices.
4. Fitting: The process dispatches to estimator-specific fit paths:
   • Linear/IV (feols/feiv): Constructs within-scale arrays, drops collinear columns, and solves OLS/2SLS using local weight transforms.
   • GLM/Poisson (fepois/feglm): Uses Iteratively Reweighted Least Squares (IRLS) with explicit observation and working weights.
   • Quantile (quantreg): Uses a Frisch-Newton solve (note: absorbed fixed effects are currently unsupported).
5. Inference: The pipeline proceeds to vcov calculation and get_inference.
6. Output: Returns a fitted result or a FixestMulti container for multiple estimations.

shared estimator API
  -> EstimationConfig
  -> parse_formula
  -> prepare_model_matrix
  -> estimator-specific get_fit
       -> feols / feiv
            -> construct within-scale arrays
            -> drop collinear columns
            -> solve OLS / 2SLS with local weight transforms
       -> fepois / feglm
            -> IRLS with explicit observation and working weights
       -> quantreg
            -> Frisch-Newton solve (absorbed fixed effects are unsupported)
  -> vcov
  -> get_inference
  -> fitted result / FixestMulti