Pyfixest

pyfixest Skill

pyfixest: fast high-dimensional fixed effects estimation for Python. Covers OLS, Poisson, and IV regression with multi-way fixed effects; difference-in-differences estimators (TWFE, did2s, lpdid, Sun-Abraham); clustered standard errors; wild bootstrap; and publication output (etable regression tables, coefplot, iplot event study plots). Use when running fixed effects regressions, difference-in-differences designs, Poisson count models with FE, or producing publication-ready regression tables. For panel random/between effects, use linearmodels; for GLM/time series without FE, use statsmodels.

Comprehensive skill for fixed effects regression, instrumental variables, and difference-in-differences estimation with pyfixest. Use decision trees below to find the right guidance, then load detailed references.

What is pyfixest?

pyfixest is a Python implementation of the R fixest package (Berge, Butts, & McDermott, 2026):

• Fast: Multi-way FE demeaning via alternating projections with numba/JAX/GPU backends
• Concise formula syntax: Fixed effects after |, IV after second |, multiple estimation via sw()/csw()
• Modern DiD: Built-in did2s, local projections DiD (lpdid), and Sun-Abraham saturated estimator
• Flexible inference: Switch SE types post-estimation; wild bootstrap, randomization inference, CCV
• Publication output: etable() for regression tables, coefplot() and iplot() for coefficient visualization

Version Notes

This skill targets pyfixest 0.40.0, the major release aligning with R fixest 0.13. Breaking changes from earlier versions:

• Default standard errors changed from "cluster by first FE" to "iid" — old code silently produces different SEs
• ssc() arguments renamed: adj → k_adj, fixef_k → k_fixef, cluster_adj → G_adj, cluster_df → G_df
• fixef_rm default changed from "none" to "singleton" — singletons now dropped by default
• Multicollinearity tolerance reduced from 1e-10 to 1e-09

How to Use This Skill

Reference File Structure

Each topic in ./references/ contains focused documentation:

Reading Order

1. New to pyfixest? Start with quickstart.md then fixed-effects.md 2. Running DiD? Read quickstart.md, then difference-in-differences.md 3. Need IV? Read quickstart.md, then instrumental-variables.md 4. Making tables? Check tables-and-plots.md 5. Coming from R fixest? Read quickstart.md then gotchas.md

Related Skills

Quick Decision Trees

"I need to run a regression"

What kind of regression?
├─ OLS with fixed effects → ./references/quickstart.md
├─ OLS without fixed effects → ./references/quickstart.md
├─ IV / 2SLS → ./references/instrumental-variables.md
├─ Poisson (count data) → ./references/integration.md
├─ Logit / Probit → ./references/integration.md
├─ Quantile regression → ./references/integration.md
└─ Multiple models at once → ./references/integration.md

"I need difference-in-differences"

DiD design?
├─ Simple 2x2 DiD (one treatment date) → ./references/difference-in-differences.md
├─ Staggered treatment timing → ./references/difference-in-differences.md
│   ├─ did2s (Gardner imputation) → ./references/difference-in-differences.md
│   ├─ Local projections DiD → ./references/difference-in-differences.md
│   └─ Sun-Abraham saturated → ./references/difference-in-differences.md
├─ Event study plot → ./references/difference-in-differences.md
├─ Visualize treatment patterns → ./references/difference-in-differences.md
└─ Parallel trends assessment → ./references/difference-in-differences.md

"I need to choose standard errors"

What inference?
├─ Heteroskedasticity-robust (HC1) → ./references/fixed-effects.md
├─ Clustered (one-way / two-way) → ./references/fixed-effects.md
├─ Few clusters (<20) → ./references/advanced-inference.md
│   └─ Wild cluster bootstrap → ./references/advanced-inference.md
├─ HAC / Newey-West → ./references/fixed-effects.md
├─ Randomization inference → ./references/advanced-inference.md
├─ Multiple hypothesis testing → ./references/advanced-inference.md
└─ Causal cluster variance (CCV) → ./references/advanced-inference.md

"I need to present results"

Presenting results?
├─ Regression table (multiple models) → ./references/tables-and-plots.md
├─ Coefficient plot → ./references/tables-and-plots.md
├─ Event study plot → ./references/tables-and-plots.md
├─ Descriptive statistics table → ./references/tables-and-plots.md
└─ LaTeX output → ./references/tables-and-plots.md

"Something isn't working"

Having issues?
├─ Different results from old code → ./references/gotchas.md
├─ feglm with fixed effects error → ./references/gotchas.md
├─ numba installation problems → ./references/gotchas.md
├─ CRV3 memory issues → ./references/gotchas.md
├─ Poisson convergence → ./references/gotchas.md
├─ Formula parsing errors → ./references/gotchas.md
├─ R fixest vs pyfixest differences → ./references/gotchas.md
└─ Singleton warnings → ./references/gotchas.md

File-First Execution in Research Workflows

Important: In data research pipelines (see CLAUDE.md), pyfixest regressions are executed through script files, not interactively. This ensures auditability and reproducibility.

The pattern: 1. Write regression code to scripts/stage8_analysis/{step}_{task-name}.py 2. Execute via Bash with automatic output capture wrapper script 3. Validation results get automatically embedded in scripts as comments 4. If failed, create versioned copy for fixes

Closely read agent_reference/SCRIPT_EXECUTION_REFERENCE.md for the mandatory file-first execution protocol covering complete code file writing, output capture, and file versioning rules. All regression scripts must follow the Inline Audit Trail (IAT) standard — see agent_reference/INLINE_AUDIT_TRAIL.md. For regression code, document model specification choices (why this estimator, why this clustering level, what identifying assumptions) with # INTENT:, # REASONING:, and # ASSUMES: comments.

See:

• agent_reference/WORKFLOW_PHASE4_ANALYSIS.md — Stage 8 (Analysis & Visualization)
• agent_reference/INLINE_AUDIT_TRAIL.md — IAT documentation standard

The examples below show pyfixest syntax. In research workflows, wrap them in scripts following the file-first pattern.

---

Quick Reference

Essential Import

import pyfixest as pf

Core Estimation Functions

Formula Syntax Quick Reference

Post-Estimation Essentials

fit = pf.feols("Y ~ X1 + X2 | fe", data=df)

fit.summary()                          # Print results
fit.tidy()                             # DataFrame of coefficients
fit.vcov("hetero")                     # Re-estimate with robust SEs (requires arg)
fit.vcov({"CRV1": "state"})            # Re-estimate with clustered SEs
fit.coef()                             # Coefficient values
fit.se()                               # Standard errors
fit.confint()                          # Confidence intervals
fit.predict()                          # Fitted values
fit.resid()                            # Residuals
fit.fixef()                            # Dict of FE name → numpy array (not a DataFrame)

Reporting

pf.etable([fit1, fit2, fit3])          # Regression table
pf.coefplot([fit1, fit2])              # Coefficient plot
pf.iplot(fit)                          # Event study / interaction plot
pf.panelview(data, unit, time, treat)  # Treatment pattern visualization

Topic Index

Citation

When this library is used as a primary analytical tool, include in the report's Software & Tools references:

Berge, L., Butts, K., & McDermott, G. (2026). pyfixest: Fast high-dimensional fixed effects estimation [Computer software]. Based on fixest (R).

Cite when: pyfixest is used for regression estimation (OLS, Poisson, IV) or difference-in-differences analysis. Do not cite when: Only imported but no estimation performed.

For method-specific citations (e.g., individual DiD estimators or inference techniques), consult the reference files in this skill and agent_reference/CITATION_REFERENCE.md.

Advanced Inference

Contents

• Wild Cluster Bootstrap
• Randomization Inference
• Multiple Testing Corrections
• Gelbach Decomposition
• Causal Cluster Variance (CCV)
• Wald Tests

Wild Cluster Bootstrap

When the number of clusters is small (<20), asymptotic cluster-robust standard errors (CRV1) have poor finite-sample properties — rejection rates can far exceed nominal levels. Wild cluster bootstrap provides more reliable inference.

Basic Usage

import pyfixest as pf

fit = pf.feols("Y ~ treatment | entity + year", data=df,
               vcov={"CRV1": "state"})

# Bootstrap test for the treatment coefficient
boot = fit.wildboottest(
    param="treatment",     # Parameter to test
    reps=9999,             # Bootstrap replications (more = more precise p-value)
    cluster="state",       # Cluster variable
    seed=42,               # Reproducibility
)

The result contains:

• p-value: Bootstrap p-value for H0: coefficient = 0
• Confidence interval: Bootstrap confidence interval

Weight Types

Wild cluster bootstrap perturbs cluster-level residuals using random weights:

• Rademacher weights (default): +1 or -1 with equal probability. Standard choice, works well with ≥10 clusters
• Webb weights: 6-point distribution. Better for very few clusters (<10), as Rademacher has only 2^G distinct bootstrap datasets

When to Use

Visualization

# Plot the bootstrap distribution
fit.plot_ritest()  # If using ritest; for wildboottest, inspect the returned object

Package dependency: Wild cluster bootstrap requires the wildboottest package (pip install wildboottest), a Python port of the R fwildclusterboot package.

Randomization Inference

Randomization inference (RI) tests the sharp null hypothesis that treatment had zero effect on every unit. It constructs a reference distribution by repeatedly reassigning treatment and computing the test statistic.

Basic Usage

fit = pf.feols("Y ~ treatment | entity + year", data=df,
               vcov={"CRV1": "state"})

# Randomization inference for treatment
ri = fit.ritest(
    resampvar="treatment",    # Variable to reshuffle
    reps=1000,                # Number of permutations
    cluster="state",          # Reshuffle at cluster level
    type="randomization-c",   # Inference type
)

Visualizing the RI Distribution

# Plot: observed test statistic vs. permutation distribution
fit.plot_ritest()

The plot shows where the actual estimate falls in the distribution of estimates under random reassignment. A p-value close to 0 means the observed effect would be very unlikely under the sharp null.

When to Use

• Randomized experiments: RI is the natural inference framework when treatment is actually randomized
• Testing the sharp null: When you want to test "did the treatment have any effect on anyone?" (rather than "is the average effect nonzero?")
• Small samples: RI does not rely on large-sample asymptotics
• As a complement: Report alongside conventional inference for robustness

Comparison to Standard Inference

Multiple Testing Corrections

When testing the same hypothesis across multiple outcomes or specifications, the probability of at least one false positive increases rapidly. pyfixest provides three correction methods.

Bonferroni Correction

The simplest correction: multiply each p-value by the number of tests.

fit1 = pf.feols("Y1 ~ treatment | fe", data=df, vcov={"CRV1": "cluster"})
fit2 = pf.feols("Y2 ~ treatment | fe", data=df, vcov={"CRV1": "cluster"})
fit3 = pf.feols("Y3 ~ treatment | fe", data=df, vcov={"CRV1": "cluster"})

# Bonferroni-adjusted p-values
pf.bonferroni([fit1, fit2, fit3], param="treatment")

Bonferroni is conservative — it controls the family-wise error rate (FWER) but may reject too few hypotheses.

Romano-Wolf Step-Down

# Resampling-based correction — less conservative than Bonferroni
pf.rwolf(
    [fit1, fit2, fit3],
    param="treatment",
    reps=999,
    seed=42,
)

Romano-Wolf uses a step-down procedure: it tests hypotheses sequentially, dropping rejected ones and re-computing critical values. This yields more power while still controlling FWER.

Westfall-Young

pf.wyoung(
    [fit1, fit2, fit3],
    param="treatment",
    reps=999,
    seed=42,
)

Another resampling-based FWER correction, similar in spirit to Romano-Wolf.

When to Use Multiple Testing Corrections

Apply corrections when:

• Testing the same treatment on multiple outcomes (e.g., does a school reform affect test scores, attendance, AND graduation?)
• Running the same specification on multiple subgroups
• Presenting a family of related hypothesis tests

Do NOT routinely apply corrections to:

• Different specifications of the same outcome (these are robustness checks, not independent tests)
• Exploratory analysis (corrections are for confirmatory testing)

Gelbach Decomposition

Gelbach (2016) decomposes the change in a coefficient when additional controls are added. This answers: "Which specific controls explain most of the change in the coefficient of interest?"

Basic Usage

# Full model with all controls
fit = pf.feols("wage ~ gender + education + experience + industry | state",
               data=df, vcov={"CRV1": "state"})

# Decompose: how much does the gender coefficient change when controls are added?
gb = fit.decompose(
    decomp_var="gender[T.male]",      # Coefficient to decompose
    combine_covariates={              # Group related covariates
        "education": re.compile("education"),
        "experience": re.compile("experience"),
    },
)

Viewing Results

import re

# Table of decomposition results
gb.etable(panels="levels")           # By variable levels
gb.etable(panels="all")              # Including normalized contributions

# DataFrame output
gb.tidy()

# Visualization
gb.coefplot()

Interpretation

The decomposition shows how much each control variable (or group of controls) explains of the gap between the bivariate coefficient (gender only) and the multivariate coefficient (gender + controls). This is more informative than simply noting "the coefficient changed when I added controls" — it pinpoints which controls matter and by how much.

Example interpretation: "The raw gender wage gap is 15%. Adding education controls reduces it by 4 percentage points, experience by 3 points, and industry by 6 points, leaving an unexplained gap of 2%."

Causal Cluster Variance (CCV)

Following Abadie, Athey, Imbens, and Wooldridge (2023), CCV provides design-based variance estimation that accounts for the specific randomization or sampling design.

fit = pf.feols("Y ~ treatment | entity + year", data=df)

ccv = fit.ccv(
    treatment="treatment",     # Treatment variable
    cluster="state",           # Cluster variable
    pk=0.05,                   # Proportion of treated clusters
    qk=1.0,                    # Proportion of units sampled per cluster
    seed=42,
    n_splits=8,                # Number of sample splits
)

When CCV Applies

CCV is appropriate when:

• Treatment is assigned at the cluster level (e.g., state-level policy)
• Clusters are sampled from a larger population
• You want to account for both sampling uncertainty and treatment effect heterogeneity

CCV standard errors can be substantially smaller than conventional cluster-robust SEs when the treatment effect heterogeneity across clusters is small relative to sampling uncertainty.

Wald Tests

Test linear hypotheses about estimated coefficients.

Using wald_test()

The wald_test() method takes a restriction matrix R and optional vector q for the hypothesis H0: Rβ = q:

import numpy as np

fit = pf.feols("Y ~ X1 + X2 + X3 | fe", data=df)

# Test: X1 = X2 (i.e., X1 - X2 = 0)
# R matrix has one row: [1, -1, 0] for [X1, X2, X3]
fit.wald_test(R=np.array([[1, -1, 0]]))

# Test: X1 = 0 AND X2 = 0 (joint significance)
fit.wald_test(R=np.eye(2, 3))  # First two coefficients

The distribution parameter controls whether the test uses an F-distribution (default) or chi-squared.

With marginaleffects (Recommended for Complex Hypotheses)

For string-based hypotheses and nonlinear tests (which require the delta method), use the marginaleffects package:

from marginaleffects import hypotheses

# Linear hypothesis
hypotheses(fit, "X1 - X2 = 0")

# Nonlinear hypothesis (delta method for standard errors)
hypotheses(fit, "(X1 / Intercept - 1) * 100 = 0")

The marginaleffects package automates gradient computation for the delta method, making nonlinear hypothesis tests straightforward.

References and Further Reading

• MacKinnon, J.G., Nielsen, M.Ø., and Webb, M.D. (2023). "Cluster-Robust Inference: A Guide to Empirical Practice." Journal of Econometrics, 232(2), 272-299
• Romano, J.P. and Wolf, M. (2005). "Stepwise Multiple Testing as Formalized Data Snooping." Econometrica, 73(4), 1237-1282
• Westfall, P.H. and Young, S.S. (1993). Resampling-Based Multiple Testing. Wiley
• Gelbach, J.B. (2016). "When Do Covariates Matter? And Which Ones, and How Much?" Journal of Labor Economics, 34(2), 509-543
• Abadie, A., Athey, S., Imbens, G.W., and Wooldridge, J.M. (2023). "When Should You Adjust Standard Errors for Clustering?" Quarterly Journal of Economics, 138(1), 1-35
• Fisher, R.A. (1935). The Design of Experiments. Oliver and Boyd. (Original randomization inference framework)
• Young, A. (2019). "Channeling Fisher: Randomization Tests and the Statistical Insignificance of Seemingly Significant Experimental Results." Quarterly Journal of Economics, 134(2), 557-598
• pyfixest documentation — Inference: https://pyfixest.org

Difference-in-Differences

Contents

• TWFE (Traditional Two-Way Fixed Effects)
• Event Study Specification
• Modern DiD Estimators
• panelview for Treatment Visualization
• Event Study Plotting
• Parallel Trends Assessment

TWFE (Traditional Two-Way Fixed Effects)

Basic TWFE DiD

import pyfixest as pf

# Classic 2x2 DiD with entity + time FE
fit = pf.feols("Y ~ treatment | entity + year", data=df,
               vcov={"CRV1": "entity"})
fit.summary()

When TWFE Works

TWFE DiD produces unbiased estimates when:

• Single treatment date: All treated units adopt treatment simultaneously
• Homogeneous effects: Treatment effect is the same across all units and time periods
• No anticipation: Units don't change behavior before treatment

When TWFE Fails

With staggered treatment timing (units adopt at different times) and heterogeneous treatment effects, TWFE can produce severely biased estimates — including sign reversals. This occurs because TWFE implicitly uses already-treated units as controls for newly-treated units, creating "forbidden comparisons" with negative weights.

When you have staggered treatment: Use one of the modern estimators below (did2s, lpdid, or Sun-Abraham saturated).

For methodological details on the TWFE problem, load the data-scientist skill's causal inference reference.

Event Study Specification

Manual Event Study with i()

# Create relative-time variable: periods since treatment
df["rel_year"] = df["year"] - df["treatment_year"]

# Event study: i() creates dummies for each relative year, omitting -1 as reference
fit = pf.feols("Y ~ i(rel_year, ref=-1) | entity + year", data=df,
               vcov={"CRV1": "entity"})

# Plot the event study
fit.iplot()

The ref=-1 normalizes to the period immediately before treatment, which is the standard convention.

Unified Event Study Interface

# pf.event_study() provides a clean unified interface
fit = pf.event_study(
    data=df,
    yname="Y",            # Outcome variable
    idname="entity",      # Unit identifier
    tname="year",         # Time period
    gname="treatment_year",  # Cohort (year of treatment adoption)
    estimator="twfe",     # Estimator: "twfe", "did2s", or "saturated"
    att=True,             # True = pooled ATT; False = dynamic event study
    cluster="entity",     # Clustering variable
)
fit.summary()

Set att=False for dynamic (period-by-period) estimates, att=True for a single pooled treatment effect.

Modern DiD Estimators

did2s — Gardner (2022) Two-Stage Imputation

fit = pf.did2s(
    data=df,
    yname="Y",
    first_stage="~ 0 | entity + year",     # FE to estimate from untreated obs
    second_stage="~ i(rel_year, ref=-1)",   # Treatment effect specification
    treatment="treated",                     # Binary treatment indicator
    cluster="entity",                        # Clustering variable
)
fit.summary()
fit.iplot()  # Event study plot

Data requirement: The dataset must include units that are never treated (treatment indicator = 0 for all periods). The did2s estimator uses these never-treated units in Stage 1 to estimate time fixed effects. Datasets where all units are eventually treated will cause estimation failure (shape mismatches or singular matrix errors).

How it works: 1. Stage 1: Estimate entity and time FE using only untreated (and not-yet-treated) observations 2. Stage 2: Impute the counterfactual for treated observations, then regress the residual on treatment indicators

Advantages: Avoids the negative-weighting problem of TWFE. Consistent under staggered adoption with heterogeneous effects.

For a pooled ATT (single treatment effect number):

fit = pf.did2s(
    data=df,
    yname="Y",
    first_stage="~ 0 | entity + year",
    second_stage="~ treated",               # Single treatment dummy → pooled ATT
    treatment="treated",
    cluster="entity",
)

lpdid — Local Projections DiD (Dube, Girardi, Jorda, & Taylor, 2023)

result = pf.lpdid(
    data=df,
    yname="Y",
    idname="entity",
    tname="year",
    gname="treatment_year",     # Cohort variable
    att=True,                   # True = pooled ATT, False = period-specific
    vcov={"CRV1": "entity"},   # Defaults to CRV1 by idname
    pre_window=5,               # Number of pre-treatment periods
    post_window=10,             # Number of post-treatment periods
    never_treated=0,            # Value of gname for never-treated units (default: 0)
)

Important: lpdid() returns a DataFrame (not a Feols object) — its API differs from feols() / did2s().

Advantages:

• Flexible dynamics: does not assume a specific functional form for treatment effects over time
• Allows non-absorbing treatment (treatment can turn off)
• Robust to misspecification of the outcome model

Sun-Abraham Saturated Estimator

# Via event_study() with estimator="saturated"
fit = pf.event_study(
    data=df,
    yname="Y",
    idname="entity",
    tname="year",
    gname="treatment_year",
    estimator="saturated",
    att=False,                  # Dynamic event study
    cluster="entity",
)

# Full interaction-weighted estimates
fit.summary()

# Aggregate to overall treatment effect
agg = fit.aggregate(weighting="shares")

How it works: Fully saturates the model with cohort-by-period indicators, then aggregates using appropriate weights. Properly handles staggered timing by never using already-treated units as controls.

Visualization:

# Event study plot of saturated estimates
fit.iplot()

Note: The R fixest package provides additional Sun-Abraham utilities (aggregate.fixest(), treatment heterogeneity tests) that may not all be ported to pyfixest yet. Check the pyfixest changelog and documentation for the latest available methods on the saturated estimator result object.

Choosing Among Modern Estimators

All three are consistent under staggered treatment with heterogeneous effects. Choice often depends on what you want to test and how you want to present results.

panelview for Treatment Visualization

Before running any DiD model, visualize the treatment assignment pattern:

# Heatmap of treatment status across units and time
pf.panelview(
    data=df,
    unit="entity",
    time="year",
    treat="treated",       # Binary treatment indicator
)

This produces a heatmap showing which units are treated in which periods — essential for understanding:

• How many units are treated vs. control
• Whether treatment timing is staggered
• Whether there are gaps or reversals in treatment
• How many never-treated units exist

Event Study Plotting

iplot() for Models with i() Terms

fit = pf.feols("Y ~ i(rel_year, ref=-1) | entity + year", data=df,
               vcov={"CRV1": "entity"})

# Basic event study plot
fit.iplot()

# With joint confidence bands (Bonferroni + Scheffe)
fit.iplot(joint="both")

# With customization
fit.iplot(
    alpha=0.05,              # Significance level
    figsize=(10, 6),         # Figure size
    joint="both",            # Both Bonferroni and Scheffe bands
    yintercept=0,            # Reference line at zero
    coord_flip=False,        # Horizontal orientation
)

Joint confidence bands are wider than pointwise CIs but account for multiple testing across periods — they provide a valid simultaneous test of no pre-trends and no treatment effect.

Comparing Estimators Visually

# Run multiple estimators
fit_twfe = pf.event_study(data=df, yname="Y", idname="entity",
                          tname="year", gname="g", estimator="twfe", att=False)
fit_did2s = pf.event_study(data=df, yname="Y", idname="entity",
                           tname="year", gname="g", estimator="did2s", att=False)

# Compare with coefplot
pf.coefplot([fit_twfe, fit_did2s])

Clustering in DiD Designs

The choice of cluster level is especially important in DiD. Following Cameron & Miller (2015): cluster at the level of treatment assignment. If a policy varies at the state level, cluster at the state level — not the entity level, even if entities are the panel units.

# State-level policy → cluster at state
fit = pf.feols("Y ~ treatment | entity + year", data=df,
               vcov={"CRV1": "state"})          # NOT "entity"

# did2s with state-level clustering
fit = pf.did2s(data=df, yname="Y",
               first_stage="~ 0 | entity + year",
               second_stage="~ i(rel_year, ref=-1)",
               treatment="treated",
               cluster="state")                   # NOT "entity"

See fixed-effects.md for full guidance on choosing cluster levels and handling few-cluster inference.

Parallel Trends Assessment

The parallel trends assumption — that treated and control groups would have followed the same trajectory absent treatment — is fundamentally untestable. Pre-treatment event study coefficients provide suggestive evidence but cannot confirm the assumption.

Visual Assessment

# Pre-treatment coefficients close to zero suggest (but don't prove) parallel trends
fit = pf.feols("Y ~ i(rel_year, ref=-1) | entity + year", data=df,
               vcov={"CRV1": "entity"})
fit.iplot(joint="both")  # Joint bands make the test more conservative

Joint F-Test for Pre-Trends

# Test that all pre-treatment coefficients are jointly zero
# Use joint confidence bands in iplot() as a visual joint test
fit.iplot(joint="both")  # Bonferroni + Scheffe simultaneous bands

# For a formal Wald test, construct the restriction matrix
# targeting the pre-period coefficient indices, or use
# marginaleffects.hypotheses() for string-based tests

Important Caveat

Failing to reject the null of no pre-trends does not confirm parallel trends. It may simply reflect low statistical power. Roth (2022) shows that pre-tests have low power against violations that would meaningfully bias treatment effect estimates. When parallel trends are critical to your identification:

• Present pre-treatment coefficients transparently
• Discuss the plausibility of the assumption based on institutional knowledge
• Consider robustness to violations (e.g., Rambachan & Roth, 2023, HonestDiD)

References and Further Reading

• Gardner, J. (2022). "Two-Stage Differences in Differences." arXiv:2207.05943
• Dube, A., Girardi, D., Jorda, O., and Taylor, A.M. (2023). "A Local Projections Approach to Difference-in-Differences." NBER Working Paper 31184
• Sun, L. and Abraham, S. (2021). "Estimating Dynamic Treatment Effects in Event Studies with Heterogeneous Treatment Effects." Journal of Econometrics, 225(2), 175-199
• Callaway, B. and Sant'Anna, P.H.C. (2021). "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, 225(2), 200-230
• Roth, J. (2022). "Pretest with Caution: Event-Study Estimates after Testing for Parallel Trends." American Economic Review: Insights, 4(3), 305-322
• Roth, J., Sant'Anna, P.H.C., Bilinski, A., and Poe, J. (2023). "What's Trending in Difference-in-Differences? A Synthesis of the Recent Econometrics Literature." Journal of Econometrics, 235(2), 2218-2244
• de Chaisemartin, C. and D'Haultfoeuille, X. (2020). "Two-Way Fixed Effects Estimators with Heterogeneous Treatment Effects." American Economic Review, 110(9), 2964-2996
• Goodman-Bacon, A. (2021). "Difference-in-Differences with Variation in Treatment Timing." Journal of Econometrics, 225(2), 254-277
• pyfixest documentation — DiD Estimation: https://pyfixest.org

Fixed Effects and Standard Errors

Contents

• Multi-Way Fixed Effects
• Standard Error Types
• Clustered Standard Errors
• Wild Cluster Bootstrap
• HAC Standard Errors
• Causal Cluster Variance
• Small Sample Corrections
• Backend Options for Performance

Multi-Way Fixed Effects

Syntax

import pyfixest as pf

# One-way FE
fit = pf.feols("Y ~ X1 | entity", data=df)

# Two-way FE (entity + time)
fit = pf.feols("Y ~ X1 | entity + year", data=df)

# Three-way FE
fit = pf.feols("Y ~ X1 | entity + year + industry", data=df)

# Interacted FE (entity-by-year pairs)
fit = pf.feols("Y ~ X1 | entity ^ year", data=df)

How FE Demeaning Works

pyfixest absorbs fixed effects via the Frisch-Waugh-Lovell theorem using alternating projections (iterative demeaning). This is fast and memory-efficient — it avoids creating dummy variable matrices.

Key parameters controlling demeaning:

• fixef_tol=1e-08: Convergence tolerance (decrease for more precision)
• fixef_maxiter=100000: Maximum demeaning iterations
• fixef_rm="singleton": Remove singleton FE groups (default since v0.40)

Extracting Fixed Effects

fit = pf.feols("Y ~ X1 | entity + year", data=df)

# Get estimated fixed effects as a dict of numpy arrays
fe_dict = fit.fixef()
# fe_dict["entity"] → numpy array of entity FE estimates
# fe_dict["year"] → numpy array of year FE estimates
# To inspect: {name: vals[:5] for name, vals in fe_dict.items()}

The fixef() method returns a `dict` mapping FE names to numpy arrays (not a DataFrame). It recovers the absorbed intercepts via the algorithm in Berge (2018). Parameters atol and btol control recovery precision.

When to Use Fixed Effects

Fixed effects are appropriate when:

• Entity FE: Control for time-invariant unobserved heterogeneity across units (states, firms, individuals)
• Time FE: Control for common shocks affecting all units in a period
• Two-way FE: Entity + time together for panel data — the standard panel regression specification
• Interacted FE (entity ^ year): When you need entity-specific time trends or very flexible controls

For methodology guidance on when FE identification is credible, load the data-scientist skill's causal inference references.

Standard Error Types

Complete Reference Table

Syntax Examples

fit = pf.feols("Y ~ X1 | fe", data=df)

# IID (default)
fit.vcov("iid")

# Heteroskedasticity-robust
fit.vcov("hetero")
fit.vcov("HC1")  # same as "hetero"

# HC2, HC3 (no FE or IV allowed)
fit_no_fe = pf.feols("Y ~ X1", data=df)
fit_no_fe.vcov("HC2")
fit_no_fe.vcov("HC3")

Clustered Standard Errors

One-Way Clustering

# At estimation time
fit = pf.feols("Y ~ X1 | entity + year", data=df, vcov={"CRV1": "state"})

# Or switch post-estimation
fit = pf.feols("Y ~ X1 | entity + year", data=df)
fit.vcov({"CRV1": "state"})

Two-Way Clustering

# Cluster by state AND year
fit.vcov({"CRV1": "state+year"})

Two-way clustering accounts for correlation within states (across years) AND within years (across states).

CRV3 for Few Clusters

# Jackknife cluster variance — more conservative, appropriate with few clusters
fit.vcov({"CRV3": "state"})

CRV3 (cluster jackknife) is more reliable than CRV1 when the number of clusters is small (roughly 10-30). For fewer than ~10-15 clusters, consider wild cluster bootstrap instead (see advanced-inference.md).

Choosing the Cluster Level

Rule of thumb from Cameron & Miller (2015): Cluster at the level of treatment assignment. If a policy varies at the state level, cluster at the state level. If treatment is at the individual level in a clustered sample, cluster at the sampling unit level.

When uncertain, clustering at a coarser level is generally conservative (wider CIs).

Wild Cluster Bootstrap

For models with very few clusters (<20), asymptotic cluster-robust SEs are unreliable. Wild cluster bootstrap provides better finite-sample inference. See advanced-inference.md for full syntax, weight type options (Rademacher vs Webb), and guidance on when to use bootstrap vs. CRV3.

HAC Standard Errors

Newey-West (Time Series)

fit = pf.feols("Y ~ X1 | entity", data=df,
               vcov="NW",
               vcov_kwargs={"time_id": "year"})

Driscoll-Kraay (Panel with Cross-Sectional Dependence)

fit = pf.feols("Y ~ X1 | entity", data=df,
               vcov="DK",
               vcov_kwargs={"time_id": "year"})

Driscoll-Kraay SEs are robust to both serial correlation and cross-sectional dependence, making them appropriate for macro panels where shocks are correlated across units.

Causal Cluster Variance

Following Abadie, Athey, Imbens, and Wooldridge (2023), CCV provides design-based inference when treatment assignment is clustered and clusters are sampled from a larger population. See advanced-inference.md for full CCV syntax, parameters, and when-to-use guidance.

Small Sample Corrections

The ssc() function controls degrees-of-freedom adjustments:

# Default behavior
fit = pf.feols("Y ~ X1 | fe", data=df,
               ssc=pf.ssc(k_adj=True, k_fixef="none", G_adj=True, G_df="min"))

# Match Stata's conventional two-way clustering
fit = pf.feols("Y ~ X1 | fe", data=df,
               vcov={"CRV1": "state+year"},
               ssc=pf.ssc(G_df="conventional"))

ssc() Parameters

v0.40 breaking change: These parameter names were all renamed. See gotchas.md for the mapping.

Backend Options for Performance

The demeaner_backend parameter controls the FE demeaning algorithm:

# Use JAX backend for GPU acceleration
fit = pf.feols("Y ~ X1 | f1 + f2", data=df, demeaner_backend="jax")

# Use scipy fallback if numba is problematic
fit = pf.feols("Y ~ X1 | f1 + f2", data=df, demeaner_backend="scipy")

GPU acceleration targets the iterative alternating-projections demeaning step, which dominates computation time for models with many FE levels. For small datasets (<100K observations) the overhead of GPU transfer may exceed the speedup — numba on CPU is typically fastest.

The solver parameter separately controls the linear algebra solver for the regression itself:

• "scipy.linalg.solve" (default)
• "numpy.linalg.solve"
• "jax" (for GPU)

References and Further Reading

• Berge, L., Butts, K., and McDermott, G. (2026). "Fast and User-Friendly Econometrics Estimations: The R Package fixest." arXiv:2601.21749
• Cameron, A.C. and Miller, D.L. (2015). "A Practitioner's Guide to Cluster-Robust Inference." Journal of Human Resources, 50(2), 317-372
• Abadie, A., Athey, S., Imbens, G.W., and Wooldridge, J.M. (2023). "When Should You Adjust Standard Errors for Clustering?" Quarterly Journal of Economics, 138(1), 1-35
• MacKinnon, J.G., Nielsen, M.Ø., and Webb, M.D. (2023). "Cluster-Robust Inference: A Guide to Empirical Practice." Journal of Econometrics, 232(2), 272-299
• pyfixest documentation — Standard Errors: https://pyfixest.org

Common Gotchas and Troubleshooting

Contents

• v0.40 Breaking Changes
• feglm Does NOT Support Fixed Effects
• numba Dependency
• Formula Parsing
• CRV3 Memory Usage
• Convergence in fepois
• Singleton Fixed Effects
• lpdid Returns a DataFrame
• HC2/HC3 Restrictions
• fixest (R) vs pyfixest Differences
• Matching Stata Results

v0.40 Breaking Changes

Version 0.40.0 aligned pyfixest with R fixest 0.13, introducing several breaking changes that silently change results:

Default Standard Errors Changed

Before v0.40: Default SE was cluster-robust by the first fixed effect variable. After v0.40: Default SE is "iid".

# Old behavior: fit.vcov was auto-set to {"CRV1": "f1"}
# New behavior: fit.vcov is "iid"

# If you want the old behavior, specify explicitly:
fit = pf.feols("Y ~ X | f1", data=df, vcov={"CRV1": "f1"})

Impact: Code that relied on the old default will produce different standard errors, t-statistics, and p-values without any error or warning. Always specify vcov explicitly to avoid ambiguity.

ssc() Arguments Renamed

# Old (will error in v0.40+)
pf.ssc(adj=True, fixef_k="nested", cluster_adj=True)

# New
pf.ssc(k_adj=True, k_fixef="nonnested", G_adj=True)

Note: the option value "nested" was also renamed to "nonnested".

Singleton Removal Default Changed

Before v0.40: fixef_rm="none" — singletons kept by default. After v0.40: fixef_rm="singleton" — singletons dropped by default.

Singleton fixed effects are groups with only one observation. Keeping them can inflate degrees of freedom and produce misleading inference.

# To preserve old behavior (not recommended):
fit = pf.feols("Y ~ X | fe", data=df, fixef_rm="none")

Multicollinearity Tolerance

Default collin_tol changed from 1e-10 to 1e-09. This may cause some near-collinear variables to be dropped that were previously kept.

feglm Does NOT Support Fixed Effects

This is the most common source of confusion. feglm() (logit, probit, Gaussian GLM) does not currently support FE demeaning:

# This RAISES NotImplementedError:
pf.feglm("binary_Y ~ X | entity", data=df, family="logit")

Workarounds

The linear probability model with heteroskedasticity-robust or clustered SEs is the most common approach in applied economics when FE are needed with a binary outcome.

numba Dependency

pyfixest uses numba for the default FE demeaning backend. numba can be tricky to install:

Common Issues

Problem: numba fails to install (especially on M-series Macs or minimal environments).

# Error: ModuleNotFoundError: No module named 'numba'
# Or: numba compilation errors on first use

Fix: Use the scipy fallback backend:

fit = pf.feols("Y ~ X | fe", data=df, demeaner_backend="scipy")

Problem: First call is very slow (numba JIT compilation).

Fix: This is expected — numba compiles on first use, then caches. Subsequent calls are fast. For scripts, this is a one-time cost.

Formula Parsing

pyfixest uses formulaic (not patsy) for formula parsing. This produces some syntax differences from statsmodels:

Categorical Variables

# pyfixest (formulaic)
"Y ~ C(state)"              # Basic categorical
"Y ~ i(state, ref='CA')"    # With reference level (preferred)

# statsmodels (patsy)
"Y ~ C(state, Treatment('CA'))"   # Reference via Treatment() — NOT supported in pyfixest

Interactions

# Both pyfixest and statsmodels
"Y ~ X1 * X2"      # Main effects + interaction
"Y ~ X1 : X2"      # Interaction only (no main effects)

# pyfixest-specific: i() for categorical interactions
"Y ~ i(group, X1, ref='control')"  # Group-specific slopes

Common Parsing Errors

# Error: variable names with spaces or special characters
# Fix: rename columns before estimation
df = df.rename(columns={"my variable": "my_variable"})

# Error: transformations not recognized
# formulaic supports: C(), np.log(), np.sqrt(), etc.
# Use numpy explicitly:
import numpy as np
"Y ~ np.log(X1) + X2"

CRV3 Memory Usage

CRV3 (cluster jackknife) standard errors require storing a G × k matrix where G = number of clusters and k = number of parameters.

Problem: With many clusters and many parameters, this can exhaust memory.

# Example: 1000 clusters × 500 parameters = large matrix
fit = pf.feols("Y ~ X1 + ... + X500 | fe", data=df)
fit.vcov({"CRV3": "cluster"})  # May run out of memory

Fix: Use CRV1 or wild bootstrap instead:

fit.vcov({"CRV1": "cluster"})  # Much less memory
# Or for few clusters:
fit.wildboottest(param="X1", cluster="cluster", reps=9999)

Convergence in fepois

Symptoms

# Warning: Maximum number of iterations reached
# Warning: Separation detected

Causes and Fixes

Slow convergence (many FE with sparse data):

fit = pf.fepois("Y ~ X | f1 + f2 + f3", data=df,
                iwls_maxiter=100,     # Increase from default 25
                iwls_tol=1e-06,       # Relax tolerance slightly
                )

Separation (FE levels that perfectly predict zero):

Some combinations of FE levels may have zero counts in all observations. These separated observations have infinite likelihood and must be removed. pyfixest can detect separation, but you may need to investigate which FE levels are problematic.

# Check for zero-count FE groups
print(df.groupby(["f1", "f2"])["Y"].sum().value_counts())

Singleton Fixed Effects

Singletons are FE groups with exactly one observation. Since v0.40, pyfixest drops them by default (fixef_rm="singleton").

Why Singletons Are Dropped

A singleton FE perfectly fits that observation's residual, contributing nothing to parameter estimation while consuming a degree of freedom. Keeping singletons inflates R² and can bias standard errors.

Warning Messages

# "X singleton observations removed"
# This is expected and correct behavior

If many singletons are removed, investigate whether your panel is very unbalanced or whether your FE specification is too fine-grained.

lpdid Returns a DataFrame

Unlike feols(), did2s(), and event_study(), the lpdid() function returns a pandas DataFrame, not a Feols object:

result = pf.lpdid(data=df, yname="Y", idname="entity",
                  tname="year", gname="treatment_year")

# result is a DataFrame with columns like:
# period, estimate, std_error, ci_lower, ci_upper, etc.

# This does NOT work:
# result.summary()     # AttributeError
# result.iplot()       # AttributeError
# pf.etable([result])  # TypeError

To visualize lpdid() results, use the returned DataFrame directly with matplotlib or plotnine.

HC2/HC3 Restrictions

HC2 and HC3 standard errors are not supported with fixed effects or instrumental variables:

# These will error:
fit = pf.feols("Y ~ X | fe", data=df)
fit.vcov("HC2")  # Error: HC2 not supported with FE

fit = pf.feols("Y ~ 1 | 0 | X ~ Z", data=df)
fit.vcov("HC3")  # Error: HC3 not supported with IV

Why: HC2 and HC3 require the hat matrix, which is expensive to compute when FE are absorbed via demeaning. Use HC1 ("hetero") or cluster-robust SEs instead.

fixest (R) vs pyfixest Differences

Features in R fixest Not Yet in pyfixest

Check the pyfixest GitHub issues and changelog for current status:

• Some etable() customization options
• Some specialized FE features
• Certain post-estimation utilities

When a feature is missing in pyfixest, consider whether statsmodels, linearmodels, or manual implementation can fill the gap.

Matching Stata Results

Clustered Standard Errors

Stata and pyfixest use slightly different default small-sample corrections:

# To match Stata's one-way clustering:
fit = pf.feols("Y ~ X | fe", data=df,
               vcov={"CRV1": "cluster"},
               ssc=pf.ssc(k_adj=True, k_fixef="none", G_adj=True))

# To match Stata's two-way clustering:
fit = pf.feols("Y ~ X | fe", data=df,
               vcov={"CRV1": "cluster1+cluster2"},
               ssc=pf.ssc(G_df="conventional"))

HC3 Standard Errors

# To match Stata's HC3 (robust, small):
fit_no_fe = pf.feols("Y ~ X", data=df)
fit_no_fe.vcov("HC3")
# Use ssc=pf.ssc(k_adj=False) if results don't match

OLS Precision

With IID standard errors, pyfixest and R fixest match to ~10^-18 precision. Poisson matches to ~10^-8 to 10^-9. Differences beyond these thresholds suggest a specification mismatch, not a numerical issue.

Polars DataFrame Input

pyfixest expects a pandas DataFrame as input. If your data pipeline uses Polars (as DAAF recommends), convert before passing to estimation functions:

# Convert Polars → pandas before estimation
df = df_polars.to_pandas()
fit = pf.feols("Y ~ X1 | fe", data=df)

Passing a Polars DataFrame directly may raise a TypeError or produce unexpected behavior. Always convert explicitly. See quickstart.md for the full conversion pattern.

Quick Diagnostic Table

References and Further Reading

• pyfixest changelog: https://py-econometrics.github.io/pyfixest/changelog.html
• pyfixest GitHub issues: https://github.com/py-econometrics/pyfixest/issues
• Berge, L., Butts, K., and McDermott, G. (2026). "Fast and User-Friendly Econometrics Estimations: The R Package fixest." arXiv:2601.21749
• Cameron, A.C. and Miller, D.L. (2015). "A Practitioner's Guide to Cluster-Robust Inference." Journal of Human Resources, 50(2), 317-372

Instrumental Variables

Contents

• IV Formula Syntax
• First Stage Diagnostics
• Weak Instrument Tests
• IV Diagnostics Summary
• Common IV Designs

IV Formula Syntax

pyfixest uses a three-part formula for IV estimation. The third part (after the second |) specifies endogenous ~ instruments:

Basic IV (No Fixed Effects)

import pyfixest as pf

# Y ~ exogenous | 0 (no FE) | endogenous ~ instruments
fit = pf.feols("Y ~ X_exog | 0 | X_endog ~ Z_instrument", data=df)
fit.summary()

Use 0 for the FE slot when you have no fixed effects but need IV.

IV with Fixed Effects

# Y ~ exogenous | FE | endogenous ~ instruments
fit = pf.feols("Y ~ X_exog | entity + year | X_endog ~ Z_instrument", data=df)

Multiple Instruments (Over-Identification)

# Two instruments for one endogenous variable
fit = pf.feols("Y ~ 1 | fe | X_endog ~ Z1 + Z2", data=df)

Over-identification (more instruments than endogenous variables) allows for Sargan/Hansen tests of instrument validity (though pyfixest does not currently implement these directly).

No Exogenous Regressors

When the only non-FE regressor is the endogenous variable, use 1 for the exogenous part:

# Only endogenous variable (plus FE)
fit = pf.feols("Y ~ 1 | entity + year | X_endog ~ Z1", data=df)

IV with Clustered Standard Errors

fit = pf.feols("Y ~ X_exog | entity + year | X_endog ~ Z1", data=df,
               vcov={"CRV1": "entity"})

Or switch post-estimation:

fit = pf.feols("Y ~ X_exog | entity + year | X_endog ~ Z1", data=df)
fit.vcov({"CRV1": "entity"}).summary()

First Stage Diagnostics

The first stage regression estimates the relationship between the instrument(s) and the endogenous variable. A strong first stage is essential for reliable IV estimates.

Accessing First Stage Results

fit = pf.feols("Y ~ 1 | fe | X_endog ~ Z1 + Z2", data=df)

# Access the first-stage Feols object (internal attribute)
first_stage = fit._model_1st_stage
first_stage.summary()

# Check first-stage F-statistic
# Rule of thumb: F > 10 suggests instruments are not weak
# (Staiger & Stock, 1997)

Note: The first-stage model is accessed via the _model_1st_stage attribute. While this is a private attribute (underscore prefix), it is the documented access pattern. The comprehensive IV_Diag() method (below) is the preferred entry point for all IV diagnostics.

Interpreting First Stage

The first stage estimates: X_endog = π₀ + π₁·Z1 + π₂·Z2 + FE + error

Key checks:

• Sign and significance of π: Instruments should predict the endogenous variable in the expected direction
• F-statistic: Joint significance of excluded instruments
• Partial R²: How much variation in X_endog the instruments explain (beyond other covariates and FE)

Weak Instrument Tests

Weak instruments produce unreliable IV estimates — biased toward OLS, with severely distorted inference.

Effective F-Statistic and Weak Instrument Tests

The preferred approach is to use the comprehensive IV_Diag() method, which reports all relevant diagnostics:

fit = pf.feols("Y ~ 1 | fe | X_endog ~ Z1 + Z2", data=df)

# Comprehensive diagnostic output — includes:
# - Effective F-statistic (Olea & Pflueger 2013)
# - Cragg-Donald F
# - Kleibergen-Paap rk Wald F
fit.IV_Diag()

• Effective F-statistic: Generalizes Stock-Yogo to non-homoskedastic settings. Critical values depend on the desired maximal bias/size distortion — consult Olea & Pflueger (2013) tables.
• Cragg-Donald F: Assumes iid errors — compare to Stock-Yogo critical values
• Kleibergen-Paap rk Wald F: Robust to heteroskedasticity/clustering — preferred with non-iid errors

Anderson-Rubin Confidence Intervals

For weak-instrument-robust inference, Anderson-Rubin (AR) confidence intervals remain valid regardless of instrument strength. These are wider than standard IV CIs but have correct coverage even with weak instruments.

IV Diagnostics Summary

fit = pf.feols("Y ~ 1 | fe | X_endog ~ Z1 + Z2", data=df)

# Comprehensive IV diagnostic output
fit.IV_Diag()

IV_Diag() combines first-stage statistics, weak instrument tests, and other diagnostic information into a single summary.

Common IV Designs

Brief descriptions of common instrument strategies. For methodology guidance on identification assumptions, load the data-scientist skill's causal inference references.

Lottery / Randomization Instruments

An experimental or quasi-experimental lottery determines treatment eligibility, but not everyone complies.

# Charter school lottery: lottery_win instruments for charter_attendance
fit = pf.feols("test_score ~ demographics | district | charter_attend ~ lottery_win",
               data=df, vcov={"CRV1": "school"})

Identification: Random assignment ensures the instrument is independent of potential outcomes. Estimates a Local Average Treatment Effect (LATE) for compliers.

Geographic / Distance Instruments

Proximity to a facility instruments for use of that facility.

# Distance to college instruments for years of education
fit = pf.feols("log_wage ~ experience | 0 | education ~ college_proximity",
               data=df, vcov="hetero")

Key assumption: Distance affects the outcome only through the endogenous variable (exclusion restriction). Threats: geographic sorting, distance correlated with local labor markets.

Policy / Regulatory Instruments

Exogenous policy variation instruments for the behavior the policy targets.

# Compulsory schooling laws instrument for education
fit = pf.feols("log_wage ~ 1 | birth_cohort + state | education ~ compulsory_years",
               data=df, vcov={"CRV1": "state"})

Reference: Angrist & Krueger (1991) quarter-of-birth design; Acemoglu & Angrist (2001) compulsory schooling.

Shift-Share / Bartik Instruments

Combines local industry shares (exposure) with national industry trends (shifts):

# Bartik instrument: local exposure to national industry shocks
# Construct the instrument as: sum(local_share_j * national_growth_j)
df["bartik"] = compute_bartik(df)  # user-constructed
fit = pf.feols("Y ~ controls | region + year | employment_change ~ bartik",
               data=df, vcov={"CRV1": "region"})

Modern guidance: Borusyak, Hull, & Jaravel (2022) and Goldsmith-Pinkham, Sorkin, & Swift (2020) provide complementary identification frameworks — one based on exogeneity of shares, the other on exogeneity of shifts.

Judge / Examiner Fixed Effects

Random assignment of cases to judges with varying leniency instruments for the decision:

# Judge leniency instruments for incarceration
# Construct leave-out judge leniency (excluding the current case)
df["judge_leniency"] = compute_leave_out_mean(df)
fit = pf.feols("Y ~ controls | court + year | incarcerated ~ judge_leniency",
               data=df, vcov={"CRV1": "judge"})

Reference: Kling (2006); Dobbie, Goldin, & Yang (2018); Stevenson (2018).

References and Further Reading

• Staiger, D. and Stock, J.H. (1997). "Instrumental Variables Regression with Weak Instruments." Econometrica, 65(3), 557-586
• Olea, J.L.M. and Pflueger, C. (2013). "A Robust Test for Weak Instruments." Journal of Business & Economic Statistics, 31(3), 358-369
• Cunningham, S. (2021). Causal Inference: The Mixtape. Yale University Press. Ch. 7: Instrumental Variables. https://mixtape.scunning.com/
• Angrist, J.D. and Pischke, J.-S. (2009). Mostly Harmless Econometrics. Princeton University Press. Ch. 4: Instrumental Variables in Action
• Borusyak, K., Hull, P., and Jaravel, X. (2022). "Quasi-Experimental Shift-Share Research Designs." Review of Economic Studies, 89(1), 181-213
• Goldsmith-Pinkham, P., Sorkin, I., and Swift, H. (2020). "Bartik Instruments: What, When, Why, and How." American Economic Review, 110(8), 2586-2624
• pyfixest documentation — IV Estimation: https://pyfixest.org

Integration and Advanced Features

Contents

• Multiple Estimation
• Poisson Regression
• GLM (Logit/Probit)
• Quantile Regression
• marginaleffects Integration
• Online Learning / Streaming
• Compressed Regression
• Performance Tuning
• Related Packages

Multiple Estimation

pyfixest can estimate many related models from a single function call using stepwise operators. This is efficient and produces well-organized output.

Stepwise Operators

Examples

import pyfixest as pf

data = pf.get_data()

# Build up controls cumulatively — classic "robustness table" pattern
fits = pf.feols("Y ~ csw0(X1, X2) | f1", data=data)
pf.etable(fits)

# Multiple dependent variables
fits = pf.feols("Y + Y2 ~ X1 | f1", data=data)
pf.etable(fits)

# Combine: multiple outcomes × cumulative controls
fits = pf.feols("Y + Y2 ~ csw0(X1, X2) | f1", data=data)
pf.etable(fits)  # Produces 2 outcomes × 3 control sets = 6 models

Sample Splitting

# Estimate by subgroup only
fits = pf.feols("Y ~ X1 | f1", data=data, split="f2")

# Full sample + each subgroup
fits = pf.feols("Y ~ X1 | f1", data=data, fsplit="f2")
pf.etable(fits)

split runs the regression separately for each level of the splitting variable. fsplit adds the full-sample estimate as the first column.

Cartesian Product of Operators

Operators combine multiplicatively:

# csw on controls × sw on FE specifications
fits = pf.feols("Y ~ csw(X1, X2) | sw(f1, f1 + f2)", data=data)
# Produces: 2 control specs × 2 FE specs = 4 models

Poisson Regression

fepois() estimates Poisson pseudo-maximum likelihood (PPML) regression with multi-way FE, following the ppmlhdfe algorithm (Correia, Guimarães, & Zylkin, 2020).

Basic Usage

# Count outcome with fixed effects
fit = pf.fepois("count_Y ~ X1 + X2 | entity + year", data=df,
                vcov={"CRV1": "entity"})
fit.summary()

When to Use Poisson

• Count data: Outcomes that are non-negative integers (patents, publications, trade flows)
• Log-linear models: Poisson PPML is consistent for E[Y|X] = exp(Xβ) even when Y is not a count — making it appropriate for gravity models in trade, for example
• Zeros in the outcome: Unlike log-OLS, Poisson handles zeros naturally without requiring log(Y+1) transformations

Convergence and Separation

# Increase iterations if convergence is slow
fit = pf.fepois("Y ~ X1 | fe1 + fe2", data=df,
                iwls_maxiter=50,       # Max IWLS iterations (default 25)
                iwls_tol=1e-08,        # Convergence tolerance
                separation_check=True, # Check for separated observations
                )

Separation occurs when some FE levels perfectly predict zero counts. Separated observations have infinite likelihood and must be detected and handled. pyfixest can check for separation but the user should be aware of this possibility with sparse count data.

Poisson with Multiple Estimation

# Stepwise controls with Poisson
fits = pf.fepois("Y ~ csw0(X1, X2, X3) | entity + year", data=df)
pf.etable(fits)

GLM (Logit/Probit)

feglm() estimates generalized linear models.

Basic Usage

# Logit model
fit = pf.feglm("binary_Y ~ X1 + X2", data=df, family="logit",
               vcov="hetero")
fit.summary()

# Probit model
fit = pf.feglm("binary_Y ~ X1 + X2", data=df, family="probit",
               vcov="hetero")

# Gaussian GLM
fit = pf.feglm("Y ~ X1 + X2", data=df, family="gaussian")

Fixed Effects Limitation

`feglm()` does NOT currently support fixed effects demeaning. This is a work in progress. Attempting pf.feglm("Y ~ X | fe", ...) raises NotImplementedError.

Workarounds:

• Linear probability model: Use pf.feols("binary_Y ~ X | fe", data=df) — interprets coefficients as percentage point changes in probability
• Manual dummies: Include FE as explicit dummy variables (slow for many levels)
• statsmodels: For logit/probit with moderately many FE levels, use statsmodels with dummy variables
• Conditional logit: For binary outcomes with entity FE, Chamberlain's conditional logit eliminates the incidental parameters problem

Quantile Regression

Estimate conditional quantile functions (experimental feature).

# Single quantile (median regression)
fit = pf.quantreg("Y ~ X1 + X2", data=df, quantile=0.5)
fit.summary()

# Multiple quantiles
fits = pf.quantreg("Y ~ X1 + X2", data=df,
                   quantile=[0.1, 0.25, 0.5, 0.75, 0.9])

# Visualize coefficient estimates across quantiles
pf.qplot(fits, nrow=2)

Standard Errors for Quantile Regression

fit = pf.quantreg("Y ~ X1 + X2", data=df, quantile=0.5,
                  vcov="iid")           # IID
fit = pf.quantreg("Y ~ X1 + X2", data=df, quantile=0.5,
                  vcov="hetero")        # Heteroskedasticity-robust
fit = pf.quantreg("Y ~ X1 + X2", data=df, quantile=0.5,
                  vcov={"CRV1": "g"})   # Clustered

Solver Options

Note: Quantile regression is marked experimental in pyfixest. Check the changelog for stability updates.

marginaleffects Integration

The marginaleffects Python package provides post-estimation interpretation for pyfixest models: average marginal effects, predictions, comparisons, and hypothesis tests.

Installation

pip install marginaleffects

Average Marginal Effects

from marginaleffects import avg_slopes, predictions, comparisons

fit = pf.feols("Y ~ X1 + X2 + X1:X2 | fe", data=df)

# Average marginal effect of X1 (accounting for interaction)
avg_slopes(fit, variables="X1")

Predictions

# Predicted values at specific covariate values
predictions(fit, newdata=datagrid(X1=[0, 1, 2], X2=df["X2"].mean()))

Hypothesis Testing (Delta Method)

from marginaleffects import hypotheses

# Linear hypothesis
hypotheses(fit, "X1 - X2 = 0")

# Nonlinear hypothesis (delta method computes gradient automatically)
hypotheses(fit, "(X1 / Intercept - 1) * 100 = 0")

Returns a DataFrame with estimate, standard error, z-statistic, p-value, and confidence interval.

Compatibility and SE Limitation

marginaleffects works directly with Feols, Fepois, and Feglm objects — no conversion needed.

Important: When using marginaleffects with pyfixest models that include fixed effects, standard errors are computed with vcov=False (no uncertainty quantification for FE parameters). This means marginaleffects SEs do not account for uncertainty in the absorbed fixed effects. For most applied work this is acceptable (FE are nuisance parameters), but be aware of this limitation when interpreting confidence intervals from avg_slopes() or predictions().

Online Learning / Streaming

Update regression coefficients for new observations without re-fitting the entire model, using the Sherman-Morrison formula.

fit = pf.feols("Y ~ X1 + X2", data=df_initial)

# Update with new data
fit_updated = fit.update(X_new, y_new, inplace=False)

• inplace=False (default): Returns a new object, original unchanged
• inplace=True: Modifies the fitted object in place

Useful for very large datasets or sequential data processing where re-estimation is expensive.

Compressed Regression

For large datasets, compressed regression reduces memory usage by compressing the data before estimation.

fit = pf.feols("Y ~ X1 + X2 | fe", data=df,
               use_compression=True)

Returns a FeolsCompressed object. Point estimates are identical to standard feols(); inference adjusts for the compression.

Performance Tuning

Backend Selection

When to Switch Backends

# GPU-accelerated estimation
fit = pf.feols("Y ~ X1 | f1 + f2 + f3", data=large_df,
               demeaner_backend="jax", solver="jax")

Other Performance Parameters

fit = pf.feols("Y ~ X1 | f1 + f2", data=df,
               lean=True,              # Reduce memory footprint of result
               copy_data=False,        # Don't copy input data (careful: may modify in place)
               store_data=False,       # Don't store data in result object
               fixef_rm="singleton",   # Drop singleton FE (default, improves speed)
               )

Related Packages

pyfixest sits within a broader Python econometrics ecosystem. These packages complement its functionality:

References and Further Reading

• Correia, S., Guimarães, P., and Zylkin, T. (2020). "Fast Poisson Estimation with High-Dimensional Fixed Effects." Stata Journal, 20(1), 95-115
• Arel-Bundock, V. (2024). "marginaleffects: Predictions, Comparisons, Slopes, Marginal Means, and Hypothesis Tests." https://marginaleffects.com/
• Berge, L., Butts, K., and McDermott, G. (2026). "Fast and User-Friendly Econometrics Estimations: The R Package fixest." arXiv:2601.21749
• Koenker, R. (2005). Quantile Regression. Cambridge University Press
• pyfixest documentation: https://pyfixest.org

pyfixest Quickstart

Contents

• Installation
• Your First Regression
• Switching Standard Errors After Estimation
• Formula Syntax Overview
• Quick Comparison: pyfixest vs statsmodels vs R fixest
• Using pyfixest with Polars DataFrames

Installation

Basic Install

pip install pyfixest
# or
uv add pyfixest

Optional Dependencies

pip install pyfixest[plots]        # lets-plot for visualization
pip install pyfixest[gt]           # great_tables for etable output
pip install pyfixest[jax]          # JAX backend for GPU demeaning
pip install wildboottest           # Wild cluster bootstrap
pip install marginaleffects        # Post-estimation interpretation

Verify Installation

import pyfixest as pf
print(pf.__version__)  # Should be 0.40.0+

Your First Regression

Basic OLS (No Fixed Effects)

import pyfixest as pf
import pandas as pd

# Load example data
data = pf.get_data()

# Simple OLS
fit = pf.feols("Y ~ X1 + X2", data=data)
fit.summary()

OLS with Fixed Effects

# One-way FE: absorb entity-level intercepts
fit = pf.feols("Y ~ X1 | f1", data=data)
fit.summary()

# Two-way FE: entity + time fixed effects
fit = pf.feols("Y ~ X1 | f1 + f2", data=data)
fit.summary()

Reading the Summary Output

The .summary() output displays:

• Dep. var.: The outcome variable
• Observations: Sample size (after dropping singletons if applicable)
• S.E. type: Standard error type used (e.g., iid, hetero, CRV1)
• Coefficient table: Estimate, Std. Error, t value, Pr(>|t|), confidence interval
• R2 / R2 Within / R2 Adj.: Model fit statistics
• Fixed effects info: Number of levels absorbed for each FE

Switching Standard Errors After Estimation

A key workflow pattern: estimate once, try different SE assumptions without re-estimating.

fit = pf.feols("Y ~ X1 | f1", data=data)

# IID (default in v0.40+)
fit.summary()

# Heteroskedasticity-robust
fit.vcov("hetero").summary()

# Clustered by f1
fit.vcov({"CRV1": "f1"}).summary()

# Two-way clustered
fit.vcov({"CRV1": "f1+f2"}).summary()

This works because SE computation is independent of point estimation under OLS — changing the variance estimator only affects standard errors, not coefficients.

Formula Syntax Overview

pyfixest uses a three-part formula separated by |:

depvar ~ exogenous_vars | fixed_effects | endogenous ~ instruments

Part 1: Dependent Variable and Exogenous Regressors

# Basic regressors
"Y ~ X1 + X2"

# Interaction with main effects
"Y ~ X1 * X2"          # equivalent to X1 + X2 + X1:X2

# Interaction only (no main effects)
"Y ~ X1:X2"

# Categorical variable
"Y ~ C(state)"

# Transformations (via formulaic)
"Y ~ np.log(X1) + X2"

Part 2: Fixed Effects (After First |)

# One-way FE
"Y ~ X1 | entity"

# Two-way FE
"Y ~ X1 | entity + year"

# Three-way FE
"Y ~ X1 | entity + year + industry"

# Interacted FE (entity-by-year)
"Y ~ X1 | entity ^ year"

Part 3: Instrumental Variables (After Second |)

# IV: endogenous ~ instrument
"Y ~ X_exog | fe | X_endog ~ Z_instrument"

# IV with no FE (use 0 for empty FE)
"Y ~ X_exog | 0 | X_endog ~ Z1 + Z2"

# Multiple instruments
"Y ~ 1 | fe | X_endog ~ Z1 + Z2"

The i() Operator for Interactions and Categoricals

# Categorical with reference level
"Y ~ i(year, ref=2000) | entity"

# Numeric interaction with categorical
"Y ~ i(group, X1, ref='control')"

# Two categoricals interacted
"Y ~ i(race, gender, ref='white', ref2='male')"

# Binning levels
"Y ~ i(age, bin={'young': [18,19,20], 'old': [60,61,62]})"

i() is especially important for event study specifications — see difference-in-differences.md.

Quick Comparison: pyfixest vs statsmodels vs R fixest

pyfixest syntax is nearly identical to R fixest, making cross-language work straightforward.

Using pyfixest with Polars DataFrames

DAAF pipelines use Polars for data processing (Stages 5-7), but pyfixest expects a pandas DataFrame as input. Convert before estimation:

import polars as pl
import pyfixest as pf

# Load processed data (Polars DataFrame from earlier pipeline stages)
df_polars = pl.read_parquet("data/processed/analysis_data.parquet")

# Convert to pandas for pyfixest
df = df_polars.to_pandas()

# Now estimate
fit = pf.feols("Y ~ X1 + X2 | entity + year", data=df)

Post-estimation results (.tidy(), .coef(), etc.) return pandas objects. If you need to rejoin results with Polars data downstream, convert back with pl.from_pandas().

pyfixest also accepts PyArrow-backed pandas DataFrames, but explicit .to_pandas() conversion is the most reliable approach.

Next Steps

• Learn about fixed effects and standard error types → fixed-effects.md
• Set up instrumental variables → instrumental-variables.md
• Run difference-in-differences designs → difference-in-differences.md
• Create publication tables and plots → tables-and-plots.md

References and Further Reading

• Berge, L., Butts, K., and McDermott, G. (2026). "Fast and User-Friendly Econometrics Estimations: The R Package fixest." arXiv:2601.21749. https://arxiv.org/abs/2601.21749
• pyfixest documentation: https://pyfixest.org
• pyfixest GitHub: https://github.com/py-econometrics/pyfixest
• R fixest documentation: https://lrberge.github.io/fixest/

Tables and Plots

Contents

• etable: Publication-Quality Regression Tables
• coefplot: Coefficient Plots
• iplot: Interaction / Event Study Plots
• dtable: Descriptive Statistics Tables
• panelview: Treatment Visualization

etable: Publication-Quality Regression Tables

etable() produces formatted regression tables comparing multiple models side by side.

Basic Usage

import pyfixest as pf

data = pf.get_data()

fit1 = pf.feols("Y ~ X1", data=data)
fit2 = pf.feols("Y ~ X1 + X2", data=data)
fit3 = pf.feols("Y ~ X1 + X2 | f1", data=data)
fit4 = pf.feols("Y ~ X1 + X2 | f1 + f2", data=data)

# Compare all models
pf.etable([fit1, fit2, fit3, fit4])

Output Formats

# Interactive HTML table (default, requires great_tables)
pf.etable([fit1, fit2], type="gt")

# Markdown table
pf.etable([fit1, fit2], type="md")

# LaTeX (requires booktabs, threeparttable, makecell)
pf.etable([fit1, fit2], type="tex")

# DataFrame for further processing
df_table = pf.etable([fit1, fit2], type="df")

# Save to file
pf.etable([fit1, fit2], type="tex", file_name="regression_table.tex")

Customization

pf.etable(
    [fit1, fit2, fit3],
    # Variable selection
    keep=["X1", "X2"],           # Show only these variables (regex supported)
    drop=["Intercept"],          # Hide these variables
    exact_match=False,           # False = regex matching (default)

    # Formatting
    coef_fmt="b \n (se)",        # Format: b=coef, se=SE, p=p-val, ci=CI
    digits=3,                    # Decimal places
    signif_code=[0.001, 0.01, 0.05],  # Significance stars

    # Labels
    labels={"X1": "Education", "X2": "Experience"},
    felabels={"f1": "Entity FE", "f2": "Year FE"},
    cat_template="{value}",      # For i() categoricals: show value only

    # Structure
    caption="Main Regression Results",
    model_heads=["(1)", "(2)", "(3)"],
)

Coefficient Format Templates

The coef_fmt parameter controls how coefficients are displayed:

Fixed Effects Rows

etable() automatically adds rows showing which fixed effects are included (checkmarks) in each model. Use felabels to provide readable names.

Mixing Estimators in One Table

feols(), did2s(), and event_study() all return Feols objects, so they can be combined in a single etable() call:

fit_twfe = pf.feols("Y ~ treatment | entity + year", data=df, vcov={"CRV1": "state"})
fit_did2s = pf.did2s(data=df, yname="Y", first_stage="~ 0 | entity + year",
                      second_stage="~ treated", treatment="treated", cluster="state")

# Compare TWFE vs did2s in one table
pf.etable([fit_twfe, fit_did2s])

lpdid() returns a DataFrame (not a Feols object) and cannot be included in etable(). Present lpdid results separately or extract coefficients manually.

Working with Multiple Estimation

When using sw(), csw(), or multiple dependent variables, feols() returns a FixestMulti object. Pass it directly to etable() — do not wrap it in a list:

# Multiple models from single estimation call
fits = pf.feols("Y ~ csw0(X1, X2, X3) | f1", data=data)
pf.etable(fits)              # Correct: pass FixestMulti directly
# pf.etable([fits])          # Wrong: wrapping in list causes TypeError

Individual Feols models from separate calls can still be combined as pf.etable([model1, model2]).

coefplot: Coefficient Plots

coefplot() visualizes estimated coefficients with confidence intervals.

Basic Usage

# Single model
pf.coefplot(fit1)

# Compare coefficients across models
pf.coefplot([fit1, fit2, fit3])

Customization

pf.coefplot(
    [fit1, fit2],
    keep=["X1", "X2"],          # Variables to include
    drop=["Intercept"],         # Variables to exclude
    coord_flip=True,            # Horizontal layout (default)
    title="Treatment Effects",
    figsize=(8, 5),
)

Use Cases

• Comparing the same coefficient across different specifications
• Visualizing the effect of adding controls (from csw0() specifications)
• Presenting results for non-technical audiences

iplot: Interaction / Event Study Plots

iplot() is specifically designed for models with i() interaction terms, particularly event studies.

Basic Event Study Plot

fit = pf.feols("Y ~ i(rel_year, ref=-1) | entity + year", data=df,
               vcov={"CRV1": "entity"})

# Default event study plot
fit.iplot()

Joint Confidence Bands

# Pointwise CIs only
fit.iplot(joint=None)

# Both Bonferroni and Scheffe simultaneous bands
fit.iplot(joint="both")

# Only Bonferroni
fit.iplot(joint="bonferroni")

# Only Scheffe
fit.iplot(joint="scheffe")

Joint bands account for multiple testing across time periods. They answer: "Can we reject that ALL pre-treatment coefficients are zero simultaneously?"

Customization

fit.iplot(
    alpha=0.05,              # Significance level for CIs
    figsize=(12, 6),         # Figure dimensions
    yintercept=0,            # Horizontal reference line
    coord_flip=False,        # Vertical (standard) orientation
    title="Event Study: Treatment Effect Over Time",
)

Comparing Estimators

fit_twfe = pf.feols("Y ~ i(rel_year, ref=-1) | entity + year", data=df)
fit_did2s = pf.did2s(data=df, yname="Y",
                      first_stage="~ 0 | entity + year",
                      second_stage="~ i(rel_year, ref=-1)",
                      treatment="treated", cluster="entity")

# Side-by-side comparison
pf.iplot([fit_twfe, fit_did2s])

dtable: Descriptive Statistics Tables

dtable() produces summary statistics tables. Note: this function is being migrated to the maketables package (maketables.DTable()).

Basic Usage

# Summary statistics for selected variables
pf.dtable(data, vars=["Y", "X1", "X2"])

By Group

# Summary by treatment group
pf.dtable(data, vars=["Y", "X1", "X2"], by="treated")

Available Statistics

Default statistics include mean, standard deviation, min, max, and count. Custom statistics can be specified.

panelview: Treatment Visualization

panelview() creates heatmap-style visualizations of treatment assignment patterns across units and time.

pf.panelview(
    data=df,
    unit="entity",        # Unit identifier column
    time="year",          # Time period column
    treat="treated",      # Treatment indicator column
)

This is essential for DiD analysis — visualize the treatment pattern before running any estimator. It shows:

• Which units are treated and when
• Whether adoption is staggered
• The size of the never-treated comparison group
• Any treatment reversals or gaps

Saving Plots

All pyfixest plots use matplotlib or lets-plot as backends. To save:

import matplotlib.pyplot as plt

# Method 1: Use matplotlib's savefig after iplot/coefplot
fit.iplot()
plt.savefig("event_study.png", dpi=300, bbox_inches="tight")
plt.close()

# Method 2: For more control, access the figure object
fig = fit.iplot()
fig.savefig("event_study.png", dpi=300, bbox_inches="tight")

For research pipeline scripts, save all figures to output/figures/ following the file naming conventions in CLAUDE.md.

References and Further Reading

• Berge, L., Butts, K., and McDermott, G. (2026). "Fast and User-Friendly Econometrics Estimations: The R Package fixest." arXiv:2601.21749
• pyfixest documentation — Visualization: https://pyfixest.org
• great_tables Python package: https://posit-dev.github.io/great-tables/
• maketables package (successor to dtable): check pyfixest changelog for migration guidance