Fast and user-friendly econometrics estimations: The R package fixest

fixest: A fast and feature-rich framework for econometric estimations in R

Introduction

Fixed-effects estimation is a cornerstone of modern applied econometrics. By absorbing unobserved heterogeneity through categorical controls, fixed-effects models allow researchers to isolate causal relationships in panel data, firm-employee records, trade flows, and many other settings. The computational challenge is that fixed-effects can have thousands or millions of levels, making standard estimation routines impractical. This has motivated a body of work on efficient algorithms for “demeaning” the fixed-effects out of the data before estimating the coefficients of interest, grounded in the canonical Frisch-Waugh-Lovell (FWL) theorem frisch1933PartialTimeRegressions, lovell1963SeasonalAdjustmentEconomic.

A variety of software packages have been developed to tackle the demeaning problem. In R rproglang2025, lfe gaure2013lfe introduced efficient fixed-effects estimation using the method of alternating projections gaure2013ols, with support for ordinary least squares (OLS), instrumental varaibles (IV), and clustered standard errors. alpaca stammann2018FastFeasibleEstimation extended the same approach to generalized linear models (GLMs). In Julia bezanson2017julia, \textbf{\texttt{FixedEffectModels}} gomez2024FixedEffectModelsJlFast uses the LSMR algorithm fong2011LSMRIterativeAlgorithm for OLS, while \textbf{\texttt{GLFixedEffectModels}} boehm2025GLFixedEffectModelsJl handles GLMs. In \textsc{Python} python2025, \textbf{\texttt{PyFixest}} fischer2024Pyfixest offers similar functionality.\footnote{In this paper, our comparison focuses on free and open source software. Note that in the not free and open source language \textsc{Stata}, correia2017_reghdfe and correia2020fast use a variant of the method of alternating projections to handle fixed-effects in OLS and GLM.}

These packages have generally made fixed-effects estimation accessible and fast. However, they each address only part of the typical workflow for applied researchers. For example, lfe does not support GLMs and is incompatible with the sandwich zeileis2020VariousVersatileVariances ecosystem for flexible variance-covariance (VCOV) estimation. alpaca handles GLMs but not OLS or IV. The core alternating projections algorithm of these two packages is also not well suited to “sparse” network structures, which is common in many real-life datasets (e.g., employee-firm settings). Similarly, the \textbf{\texttt{PyFixest}} library, developed by independent authors, implements the same algorithm and hence suffers from the same limitations. It does, however, provide a set of features close to \textbf{\texttt{fixest}}, since its original objective was to port \textbf{\texttt{fixest}} syntax and features to \textsc{Python}. In contrast, the two Julia packages offer strong performance for their respective use-cases, but with limited post-estimation methods and options (e.g., conducting inference on different VCOVs). In sum, none of the existing packages provides a unified interface that spans estimation, inference, and presentation within a single, highly performant framework.

It is against this background that we introduce the R package fixest. The package name reflects the fact that fixest originated as a specialized library for fixed-effects regression. However, it has since evolved into a general purpose tool that covers a wide range of econometric methods and use cases. Its defining ambition is to be both fast and coherent, and thus to fill a gap left by existing software options. For many applied research projects and researchers, \texttt{fixest} should suffice as a standalone package able to cover all of their econometric needs.

fixest offers best-in-class performance for many frequentist estimators, supported by a unified framework that enables estimation, inference, and reporting in a single package. It supports OLS, IV, GLM, and maximum likelihood (MLE) estimation through a common syntax that is fully compatible with standard R stats methods. For fixed-effects estimation, specifically, it introduces a novel algorithm based on fixed-point acceleration using the method of irons1969version, which achieves fast convergence across many different data contexts. The performance-critical parts of the codebase are implemented in parallel C++ via \textbf{\texttt{Rcpp}} eddelbuettel2011jss_rcpp, further ensuring optimal efficiency.

We illustrate fixest performance in a wide array of benchmarks. Using simulated data, we generate scenarios with simple or difficult fixed-effects convergence and compare fixest to leading alternatives in R, Python and \textsc{Julia}. \textbf{\texttt{fixest}} is consistently the fastest in all situations, closely followed by the \textsc{Julia} software. A second set of benchmarks uses a large real-life data set to illustrate \textbf{\texttt{fixest}} specific features useful in applied econometrics (detailed below). In those cases, \textbf{\texttt{fixest}} speed is unmatched.

In addition to speed, fixest contributes several features that distinguish it from existing software. For example, one of the core design principles is that inference is separate from estimation. This means that users can swap out a variety of VCOVs---including heteroskedasticity-robust, clustered, Newey-West newey1987simple, Driscoll-Kraay driscoll1998consistent, Conley conley1999gmm, and other standard errors---at virtually zero computational cost, since the underlying model does not have to be re-estimated. An expressive formula interface supports programmatic variable interpolation, factor manipulation, and multiple estimations in a single call, where shared computations are pooled across models for additional speed gains. fixest also provides built-in tools for publication-ready regression tables and coefficient plots, and invests heavily in informative error messages that help users diagnose problems quickly.

This article describes fixest version 0.14.0, available from the Comprehensive R Archive Network at \url{https://cran.r-project.org/package=fixest}. We begin with a self-contained motivating example in Section (ref) that introduces the core syntax through an end-to-end analysis. Section (ref) presents the econometric theory and associated computational considerations for estimating fixed-effects models, as well detailing fixest's novel fixed-point acceleration solution. Sections (ref) through (ref) describe the package's estimation functions, formula interface, VCOV options, and presentation tools. Section (ref) reports benchmarks against leading alternatives. Section (ref) discusses fixest's integration with the broader \textsc{R} ecosystem and pedagogical context, before concluding with the package's development philosophy.

Motivating example: Estimating the returns to science funding

Assume that we have data on scientific researchers in the EU and US. We are able to track their individual productivity levels (as measured by the number of articles produced every year), as well as their access to funding (in \$US'000). Our goal is to estimate the effect that funding has on scientific production, whilst controlling for differences across individuals, years and regions.

For demonstration purposes, let us create a dummy dataset called scipubs with six variables: number of articles (articles), funding received (funding), region (eu_us), a dummy representing a policy shock (policy), scientist identifier (indiv), and the year (\texttt{year}).

Ordinary least squares estimation

We start with a simple ordinary least squares (OLS) model, where we regress science production (articles) on funding (funding) and a set of control variables. OLS estimations in fixest are performed with the function feols. In its simplest guise, feols can be seen as a drop-in replacement for the base \textsc{R} \texttt{stats::lm} function. For example, we can use simple indicator (factor) variables to control for individual and year fixed-effects:

The above feols results are identical to those from an equivalent stats::lm call.\footnote{Specifically, lm(articles\ \textasciitilde\ funding\ +\ as.factor(indiv)\ +\ as.factor(year),\ scipubs).} However, the output is slightly different insofar as feols automatically displays the results in summary form (including standard errors, R-squared, etc.) and avoids displaying all the coefficients unless explicitly requested by the user. In contrast, the default lm print method would only display the point estimates and would do so for all 118 coefficients.

This interchangeability between feols and lm is convenient, since only the function call changes. However, it is not the preferred way to estimate a fixed-effects model with feols. Rather, we can take advantage of fixest's specialized architecture for fixed-effects by invoking the syntax:

\[ \mathtt{dep\_var} \sim \mathtt{indep\_vars} \mathtt{\ |\ } \mathtt{fixed\_effects} \]

In a fixest formula, any variable after the pipe (\textbar) is treated as a fixed-effect and passed to a highly optimized internal routine. Applying this principle to our example, we have:

Observe that our estimate of the key funding coefficient remains identical to our previous result of 0.096365. Moreover, because we specified indiv and year as fixed-effects, the return object correctly recognizes these variables as controls (“nuisance” parameters) and avoids displaying the 100+ coefficients on them. At the same time, it is important to note that the benefits extend beyond aesthetics. This second estimation completes quicker than the first one, exactly because the fixed-effects are handled by a specialized and highly optimized routine. We demonstrate the speed advantages in a benchmarking section later in the paper.

To interpret the coefficients, we need to know whether they are statistically significant. Hence, we need to report their standard-errors. By default feols assumes that the errors are homoskedastic and non-correlated (labelled for convenience as “IID”, or “independent and identically distributed” in our return object above).

But fixest makes it easy to specify a wide variety of alternate standard-error types. For example, because our data is a panel (i.e., repeated cross-sections of the same individuals over time), we may be worried that the errors are correlated within individuals. In econometric parlance, this suggests the need to cluster the standard-errors at the individual level. We can implement this directly at estimation time by using the argument vcov\ =\ \textasciitildeindiv:

Of course, other standard-error formulations are also possible. While one may be tempted to re-run multiple versions of the same model with slight variations to the standard-errors, this quickly becomes inefficient. A core design principle of fixest is that computation should be separate from inference. As such, fixest makes it easy to adjust standard errors for existing models “on the fly”. Almost all of the package's post-estimation routines and methods allow for this kind of post hoc adjustment. Here we briefly demonstrate with a simple side-by-side comparison of two standard types, using the etable function to display the results in regression table form:

The key argument above is vcov\ =\ list("iid",\ \textasciitildeindiv)), which prompts the display (and computation of) of both IID and clustered standard errors alongside each other. Note that this type of functionality is easily extended and that fixest supports many other standard-error adjustments. We defer details until later in the paper and leave this as a taster.

Let us return to our broader research question. Suppose we want to know whether the effect of science funding is similar across the EU and US? While there are multiple approaches to tackling this question, a potential starting point is to estimate separate models for each region. The split argument provides a convenient way to do this in a single call. Again, we can pass the resulting object to etable for convenient side-by-side display.\footnote{In truth, passing to etable is redundant since fixest multiple estimation objects are printed with etable by default. But we do so here to avoid assuming too much familiarity with the package at this introductory stage.}

Thus far, we have made minimal assumptions about the causal direction of our results. An obvious complication is the likely endogeneity of scientific funding (e.g., better researchers with more publications get more funding). We may therefore wish to instrument funding with an unexpected policy shock that changed funding availability in an exogenous way. For our toy example here, we can use the policy variable---corresponding to a dummy on certain years---as an instrument. To perform an instrumental variable (IV) estimation in fixest, we use the following syntax:

\[ \mathtt{dep\_var} \sim \mathtt{exogenous\_vars} \mathtt{\ |\ } \mathtt{fixed\_effects} \mathtt{\ |\ } \mathtt{endo\_vars} \sim \mathtt{instruments} \]

Applying this syntax to our research productivity question, we have:

The results of the second stage estimation are reported, as well as some fit statistics, like the (weak instrument) F-test for the first stage and the Wu-Hausman endogeneity test.

Generalized linear model estimation

All of our estimations thus far have modeled articles as the outcome of interest. But number of articles is a positive count (integer). Hence it is arguably more natural to model the outcome as a Poisson process. In fixest this is easily done via the fepois function:

As the reader will note, the syntax is identical to the feols case; we simply switch out the main function call. The same is true for a variety of other generalized model families that fixest supports through its feglm and femlm functions. Again, however, we leave this a taster and defer details until later in the paper.

Presentation

Having iterated through various models, let us conclude this introductory example by considering how we might present our findings. Assume we have settled on the following variants: Poisson model with the EU/US region split (and full sample), and the OLS-IV model.

The same etable function from before can be used to compile all of our preferred results into a single \LaTeX table. We need simply specify a valid file path to write \LaTeX output instead of printing the table to screen. We also invoke some additional arguments:

The resulting table is shown in Table (ref).

In addition to a regression table, we might also wish to present our results in the form of a coefficient plot via the dedicated coefplot function. Below we demonstrate with a simple example, limited to the split-sample Poisson estimation. The resulting plot is displayed in Figure (ref).

Thus concludes our introductory example. Having established the core intuition and syntax, we now examine the econometric theory and computational algorithms that make fixest estimations fast and reliable. We begin a theory review, before turning to a detailed review of the package's main functions and features.

Fixed-effects: Theory and the fixest implementation

fixest originated as a package specifically optimized for tackling fixed-effects estimations. However, since the definition of the term `fixed-effects' can differ greatly across fields, we first explain the meaning of a fixed-effects estimation from our perspective. We also show how it is tightly linked to variable demeaning. We then document existing solutions to demeaning problems and detail fixest's specific approach.

We offer two expositional notes before proceeding. First, we focus on OLS estimation. The extension to generalized linear models is discussed in Section (ref). Second, we focus on a very simple setup; namely, a balanced panel with only two fixed-effects dimensions (individual and time). The logic readily extends to imbalanced panels and any number of fixed-effects. However, this level of generality would necessitate ugly notations that can obscure the key insights without corresponding benefit. Hence, we refrain from attempting this here. See Berge2018_FENmlm for a treatment of the general case in the context of maximum likelihood estimation.

Fixed-effects and the FWL theorem: A demeaning problem

For illustration, consider the setup of a balanced panel in which we follow individuals, indexed by \(i\), across years, indexed by \(t\). There are \(N_I\) individuals and \(N_T\) periods, leading to \(N = N_I\times N_T\) observations. Let \(y\) be the \(N\)-vector representing the dependent variable and \(X\) a \(N \times K\) matrix of explanatory variables. We are interested in the following data generating process:

with \(X_{it}\) the row-vector of explanatory variables for individual \(i\) and time \(t\), and \(\gamma\) the \(K\)-vector of the coefficients of interest. In the context of the previous equation, the coefficients \(\alpha_i\) capture any confounder (i.e. explanatory variables that are not in the equation) that is individual-specific and time-invariant. Similarly, \(\beta_t\) captures any time-specific confounder. Finally, \(\epsilon_{it}\) is a random variable representing the error term.

It is important to stress that we are not directly interested in the values of \(\alpha_i\) and \(\beta_t\). Rather, they are effectively nuisance parameters, insofar as they provide a means of controlling for a set of confounders. Indeed, \(\alpha_i\) and \(\beta_t\) here represent the fixed effects of the estimation, since they capture invariant (fixed) confounders across the unit and time dimensions, respectively. Fixed effects are often written as \(\alpha_i\), but this is really short-hand notation for a set of \(N_I\) indicator variables for each unit with a coefficient \(\alpha_i\) for each (equivalent to factor(unit) in a regression). Extending this idea more generally, we define “fixed-effects estimation” in this article simply as an estimation that includes categorical variables as controls.

Now let us link fixed-effects with demeaning. Consider the estimation of the coefficients of Equation (ref) by OLS. With a hat referring to the estimated coefficients and \(r_{it}\) to the residual, we have:

Let us define \(\tilde{y}\) as the residual of the OLS estimation of \(y\) on the fixed-effects only:

\[ \tilde{y}_{it} \overset{\text{def}}{=} y_{it} - \hat{\alpha}_i^y - \hat{\beta}_t^y, \]

with \(\hat{\alpha}_i^y\) the OLS estimator of the fixed-effect coefficient.

In a similar fashion, let us define matrix \(\tilde{X}\) whose \(k^{th}\) column is defined as the residual of the estimation of the \(k^{th}\) column of \(X\) on the fixed-effects only. If we estimate \(\tilde{y}\) on \(\tilde{X}\) by OLS, we have:

The Frisch-Waugh-Lovell frisch1933PartialTimeRegressions, lovell1963SeasonalAdjustmentEconomic theorem states that the coefficients \(\hat{\gamma}\) and the residuals \(r_{it}\) of Equation (ref) are the same as the ones in Equation (ref). That is \(\hat{\gamma} = \hat{\gamma}_{bis}\) and \(r_{it} = r_{it}^{bis}\). It follows that the variance-covariance matrices of the coefficients are also similar: \(V\left(\hat{\gamma}\right)=s\times V\left(\hat{\gamma}_{bis}\right)\), with \(s \in \left[0;1\right]\) a scaling factor.\footnote{See ding2021FrischWaughLovell for intuitive demonstrations.}

The equivalence between Equation (ref) and Equation (ref) means that instead of estimating \(N_I + N_T + K\) coefficients, we can estimate only \(K\) coefficients. Given that solving an OLS problem is of complexity \(\mathcal{O}\left(n_{coef}^3\right)\) with \(n_{coef}\) the number of coefficients to estimate, avoiding the estimation of \(N_I + N_T\) coefficients greatly reduces the dimensionality of the problem.\footnote{For large $N_I$, or $N_T$, the direct estimation of Equation (ref) is simply intractable.}

However, to estimate Equation (ref), one still needs to compute \(\tilde{y}\) and \(\tilde{X}\) first. The good news is that computing these values is generally fast, thanks the simple structure of the fixed-effects. Specifically, the computing step amounts to a demeaning exercise, which we now illustrate.

To compute \(\tilde{y}\), we need to find the coefficients \(\hat{\alpha}_i^y\) and \(\hat{\beta}_t^y\) that solve the following OLS minimization problem:

We cannot solve for \(\alpha_i\) and \(\beta_t\) directly, since we would be confronted by the aforementioned \(\mathcal{O}\left(n_{coef}^3\right)\) complexity of the OLS solution. Instead, we can solve it iteratively by minimizing over \(\alpha_i\) and \(\beta_t\) sequentially until we obtain the optimal coefficients. Let \(\alpha_i^{(m)}\) be the coefficients minimizing Equation (ref) given \(\beta_t^{(m-1)}\) at iteration \(m\). We define \(\beta_t^{(m)}\) similarly, with \(\beta_t^{(0)} = 0, \forall t\). The first order conditions (FOC) of Equation (ref) lead to:

As we can see, the coefficients satisfying the FOCs are simple group means; hence our describing it as a demeaning problem. The complexity of obtaining all \(N_I\) coefficients \(\alpha_i^{m}\) is \(\mathcal{O}\left(n_{obs}\right)\). Let \(n_{iter}\) be the number of iterations required for convergence, such that \(|\alpha_i^{(n_{iter}+1)}-\alpha_i^{(n_{iter})}| \approx 0, \forall i\) and \(|\beta_t^{(n_{iter}+1)}-\beta_t^{(n_{iter})}| \approx 0, \forall t\). Then for large \(N_I\), \(N_T\) and for small \(n_{iter}\), obtaining \(\hat{\gamma}\) from Equation (ref) is much faster than directly estimating Equation (ref).

To estimate \(\hat{\gamma}\) we end up with a total complexity of \(\mathcal{O}\left(K \times n_{iter} \times n_{obs} \times n_{fixef}^2\right)\). We can see that the number of iterations, \(n_{iter}\), plays a critical role in reducing complexity. Defining an algorithmic routine that ensures a low value of \(n_{iter}\) is thus essential to reducing computing time and cost. Before detailing fixest's novel approach, we first describe existing solutions offered by other software.

Fixed-effects solutions in other software

As we have shown, the key to a fast OLS estimation with fixed-effects is finding a computationally efficient solution to the fixed-effects-only regression described by Equation (ref). In fact, this problem corresponds to a sparse linear regression, and various solutions have been offered to solve it. The most widely known of which is the method of alternating projects.

Alternating projections

The method of alternating projections gaure2013ols is almost equivalent to the demeaning problem defined in the previous section. The only difference is that, instead of finding the solution coefficients, only the residuals of the OLS regression is computed (\(\tilde{y}\)).

To illustrate this method, let us rewrite the problem with matrices. We drop the double indexing \(it\) and use a single index \(o\) to refer to the observations. Let \(D_I\) be the \(N \times N_I\) matrix such that \((D_I)_{o,i}\) is equal to 1 if observation \(o\) refers to individual \(i\) and 0 otherwise. We define \(D_T\) similarly.

Using these matrix notations, the regression of \(y\) on the individual fixed-effects leads to the classic OLS solution \(\hat\alpha = (D_I^\top D_I)^{-1}D_I^\top y\). Let the \(N \times N\) matrix \(R_I\) be the residual maker, defined as:

\[ R_I \overset{\text{def}}{=} I_N - D_I(D_I^\top D_I)^{-1}D_I^\top. \]

The product \(R_I y\) gives the residual of \(y\) such that \((R_I y)_{o} = y_o - \hat\alpha_i\), with observation \(o\) referring to individual \(i\). You can see a direct equivalence between the value of this residual and the right hand numerators in Equation (ref).

With \(R_T\) defined similarly, we can reach the residual of the OLS estimation on the two fixed-effects by sequentially applying \(R_I\) and \(R_T\), such that:

\[ \tilde{y} = R_T R_I ... R_T R_I y. \]

One algorithm to obtain \(\tilde y\) is: set \(y^{(0)} = y\), then compute iteratively \(y^{m} = R_T R_I y^{(m-1)}\) until \(||y^{m} - y^{(m-1)}||^\infty \approx 0\). Each update of \(y^{(m)}\) is known as an Halperin projection. This approach avoids to compute explicitly the coefficients of the OLS estimations at each step and is then slightly faster than the demeaning algorithm of the previous section, although they have similar convergence properties.\footnote{The astute reader may wonder: why not simply applying a quick exponentiation to \((R_T R_I)\) before right multiplying by \(y\)? The problem is that while \(R_T\) and \(R_I\) are sparse matrices, their product is dense, which quickly becomes problematic for \(N \times N\) matrices.}

This method is used in the R packages lfe gaure2013lfe and alpaca stammann2018FastFeasibleEstimation, and in the Python package \textbf{\texttt{PyFixest}} fischer2024Pyfixest.

Other methods

In the language Stata, reghdfe correia2016feasible uses the method of kelner2013SimpleCombinatorialAlgorithm to solve the linear system of equations given by Equation (ref).

In the language Julia, the software FixedEffectModels gomez2024FixedEffectModelsJlFast solves the linear system using the LSMR method fong2011LSMRIterativeAlgorithm.

The fixest solution: A fixed-point problem with acceleration

The package fixest uses a demeaning algorithm to compute \(\tilde{y}\) and \(\tilde{X}\). But it does not use the method of alternating projections shown above. Instead it estimates the solution coefficients directly. Depending on the structure of the fixed-effects, how they are distributed across observations, the number of iterations required to reach convergence can be arbitrarily high. To ensure good convergence properties in all circumstances, fixest uses a fixed-point acceleration algorithm whose implementation we now detail.

Before starting, we must say that the current example of a balanced panel converges in exactly one iteration, making the demeaning algorithm of Equation (ref) the fastest available for this situation and any acceleration unnecessary. Even though we proceed with this balanced example for illustration purposes, it may be assumed that our method is applied to more difficult contexts where convergence requires a high number of iterations.\footnote{A canonical example of difficult convergence is the context of an estimation with employees and firms fixed-effects. Typically the majority of employees will be employed in only a handful of firms across their lifetime and the probability that two employees are employed in the same company is low. From a network perspective, where employees are connected to other employees through companies, this means that the average network distance between all employees is very high; making shocks very slow to propagate. This contrasts sharply against the balanced panel case, where the average network distance between each individual is exactly one, and shocks propagate instantly. In Section (ref), we provide a simple benchmark example where changing the structure of the fixed-effects increases the computational time by several orders of magnitude.}

Let us rewrite the demeaning algorithm as a fixed-point problem. Having the demeaning algorithm in this form proves useful as it unlocks the use of fixed-point acceleration algorithms.

We define the function \(f_I:\mathbb{R}^{N_T} \rightarrow \mathbb{R}^{N_I}\) such that \(\alpha^{(m + 1)} = f_I (\beta^{(m)})\), with \(\alpha^{(m + 1)}\) and \(\beta^{(m)}\) the vector versions of the values in Equation (ref). Defining \(f_T: \mathbb{R}^{N_I} \rightarrow \mathbb{R}^{N_T}\) similarly, we can rewrite Equation (ref) using these functions:

\[

\]

Replacing \(\alpha^{(m + 1)}\) by this definition leads to:

\[ \beta^{(m + 1)} = f_T (f_I (\beta^{(m)})). \]

We clearly see that \(\beta^{(m + 1)}\) is only a function of \(\beta^{(m)}\) and does not depend on \(\alpha^{(m + 1)}\) in any way. Defining \(F = f_T \circ f_I\), \(F: \mathbb{R}^{N_T} \rightarrow \mathbb{R}^{N_T}\) our problem reduces to finding \(\beta^*\) such that:

\[ \beta^* = F(\beta^*). \]

This is a fixed-point problem with \(N_T\) parameters and, hence, a canonical fixed-point algorithm can be leveraged to hasten the convergence process ramiere2015iterative. In the fixest case, we use algorithm proposed by irons1969version.

Starting at \(\beta^0_t=0, \forall t\), the Irons and Tuck algorithm replaces \(\beta^{(m + 1)} = F(\beta^{(m)})\) with the following iteration:

\[ \beta^{(m+1)} = F(F(\beta^{(m)})) - \dfrac{\Delta F(\beta^{(m)}) \cdot \Delta^2 \beta^{(m)}}{||\Delta^2 \beta^{(m)}||^2}\Delta F(\beta^{(m)}), \]

where \(\Delta \beta^{(m)} = F(\beta^{(m)}) - \beta^{(m)}\), \(\Delta F(\beta^{(m)}) = F(F(\beta^{(m)})) - F(\beta^{(m)})\), and \(\Delta^2 \beta^{(m)} = \Delta F(\beta^{(m)}) - \Delta \beta^{(m)}\).

Once the coefficients \(\beta^*\) are found, we can reconstruct the other coefficients \(\alpha^* = f_I (\beta^*)\), compute the demeaned variables \(\tilde{y}\) and \(\tilde{X}\), and obtain our coefficients of interest \(\hat{\gamma}\) as the OLS coefficients of estimating \(\tilde{y}\) on \(\tilde{X}\).

Varying slopes

On top of regular fixed-effects estimations, a unique strength of fixest is that it efficiently handles variables with varying slopes. The concept of varying slopes parallels that of random slopes in the mixed-effects literature, where the coefficient on one variable is modulated by the value of another (categorical) variable. For instance, we might wish to model individual-specific time trends in our panel, under the following data generating process:

\[ y_{it} = \alpha_{i} + t \times \beta_i + X_{it} \gamma + \epsilon_{it}. \]

Each individual in this setting is modeled as having their own intercept, \(\alpha_i\), as well as their own linear time trend coefficient \(\beta_i\). Hence there would be \(N_I\) trends to estimate. fixest is able to estimate the set of coefficients \(\beta_i\) all at once. We briefly expose the idea of this algorithm.

Generalizing the above example, consider a model with individual fixed-effects and a varying slopes variable \(z_{it}\):

\[ y_{it} = \alpha_{i} + z_{it} \times \beta_i + X_{it} \gamma + \epsilon_{it}. \]

The first order conditions for the OLS problem leads to the following optimal values for \(\alpha_i\) and \(\beta_i\):

\[ \alpha_{i} = \dfrac{1}{n} \sum_{t=0}^{N_T}\left( y_{it} - z_{it} \beta_i - X_{it} \gamma \right) \quad \beta_{i} = \dfrac{1}{\sum_{t=0}^{N_T} z_{it}^2} \sum_{t=0}^{N_T}\left( y_{it} - \alpha_i - X_{it} \gamma \right) \]

Taking all other values as given, we can solve for \(\alpha_i\) and \(\beta_i\) at once, since this is a system of two equations with two unknowns. There is a unique solution for them.\footnote{Internally, when \(\sum_{t=0}^{N_T} z_{it}^2 = 0\), we assume the absence of the coefficient \(\beta_i\).} Hence, to find all the coefficients \(\alpha_i\) and \(\beta_i\), we need only solve \(N_I\) tiny systems of equations; no iteration is needed. Internally, fixest uses Gaussian elimination to solve these systems.

Finally, to obtain the coefficients \(\gamma\), we proceed exactly as in the previous section. We regress \(y\) on the fixed-effects and the variables with varying slopes to obtain the residual \(\tilde{y}\), and obtain \(\tilde{X}\) similarly. Then we regress \(\tilde{y}\) on \(\tilde{X}\). Using the technique mentioned above, the regression to obtain \(\tilde{y}\) is direct and does not require any iteration.

Now let us generalize the previous example. Instead of having one fixed-effect and one variable with varying slopes, we have a set of \(L_I\) variables with varying slopes contained in the \(N \times L_I\) matrix \(Z_I\), with the coefficients varying across the individual dimension. For an individual \(i\), we have to estimate \(L_I\) coefficients, i.e. the vector \(\vec{\alpha}_i\). We can contrast this against the simple fixed-effects case, where we had only had to estimate one coefficient \(\alpha_i\). Similarly for the time dimension, we have a set of \(L_T\) variables with varying slopes in the matrix \(Z_T\). The model we estimate is:

\[ y_{it} = Z_I\vec{\alpha}_{i} + Z_T \vec{\beta}_t + X_{it} \gamma + \epsilon_{it}. \]

Note that the individual fixed-effects, not attached to a variable, can be represented by a column of 1s in the matrix \(Z_I\).

As always, to obtain the coefficient of interest \(\gamma\), we first need to obtain the residual of the OLS regression of \(y\) on all the coefficients with varying slopes (in \(Z_I\) and \(Z_T\)). To obtain the \(L_I\) coefficients of \(\vec{\alpha}_{i}\) and the \(L_T\) coefficients of \(\vec{\beta}_{t}\), we apply the same fixed-point algorithm detailed in the previous section, where we iterate across these two dimensions.

The fact that all coefficients in a given dimension (e.g., \(\vec{\alpha}_{i}\)) are found at once---irrespective of the number of variable with varying slopes---greatly reduces the complexity of the algorithm and increases numerical stability. To the best of our knowledge, this feature is unique to fixest.

Extending OLS to GLM and MLE

Our exposition thus far has been limited to the OLS framework. However, fixest also handles fixed-effects efficiently in the context of generalized linear models (GLM) estimations and maximum likelihood estimation (MLE). Here we briefly outline fixest's GLM and MLE implementations, and discuss the pros and cons of each approach for similar likelihood families.

In OLS, the conditional expectation of the dependent variable is a linear combination of the predictor variables and the coefficients: \(E\left(y \ |\ X\right) = X\gamma\). GLM/MLE offers more flexibility as they model the conditional expectation of \(y\) as a function of the linear combination of the predictors and their coefficients: \(E\left(y \ |\ X\right) = f\left(X\gamma\right)\). Further, in both GLM and MLE the dependent variable is modeled with a certain likelihood. The objective of these algorithms is to maximize that likelihood with respect to the coefficients \(\gamma\). For example, if \(y\) is a count variable (i.e. only zero or positive integers), one natural way to model it is with a Poisson model, such that its expectation is modeled as: \(E\left(y \ |\ X\right) = \exp\left(X\gamma\right)\).

In fixest, the maximization algorithms of GLM and MLE differ: GLM uses the iteratively re-weighted least squares (IRLS), while MLE uses the direct maximum likelihood estimation with a Newton-Raphson algorithm. We describe the key differences between these two methods below.

fixest's GLM estimation

Given that GLM's IRLS is a succession of OLS estimations, the OLS algorithm detailed in the previous section can be directly leveraged. Hence fixed-effects can be estimated directly, as well as variables with varying slopes. Further, the method works for all likelihood families that can be estimated via GLM (Gaussian, Poisson, binomial, Gamma, etc).

fixest's MLE estimation

fixest implements the MLE method for four likelihoods: Gaussian, Poisson, Logit and Negative Binomial. Fixed-effects are implemented via a concentration of the likelihood, now briefly detailed. See Berge2018_FENmlm for a full description.

Let \(\alpha\) represent all the fixed-effect coefficients and \(\gamma\) the coefficients of interest. The log-likelihood can be written as \(\ln L(y, X, \gamma, \alpha)\) and we want to find the maximizing coefficients \(\hat\alpha\) and \(\hat\gamma\). Let \(\ln \tilde L(y, X, \gamma)\) be the concentrated likelihood such that \(\ln \tilde L(y, X, \gamma) = \max_{\alpha}\ln L(y, X, \gamma, \alpha)\).

In the case of the Gaussian likelihood, the coefficients to concentrate the likelihood are obtained in a similar fashion as the OLS problem, given that the first order conditions are similar. For the Poisson likelihood, closed forms for the coefficients also exist while this is not true for the other two likelihoods. For those last two, various tricks are leveraged to hasten the computation of the concentrated likelihood, detailed in Berge2018_FENmlm.

Differences between GLM and MLE procedures in fixest

Contrary to GLM, fixest's MLE method suffers from the following limitations: i) variables with varying slopes cannot be used, ii) analytical weights cannot be used, iii) the families are limited to the four hard coded ones, iv) the algorithm is more sensitive to starting values.

However, contrary to GLM, the MLE estimation has two main advantages. First, it includes the Negative Binomial. This likelihood family has an ancillary parameter governing overdispersion, which cannot be directly estimated via GLM. Second, the functional form of the right hand side is not restricted to linear in parameter predictors. Indeed, the IRLS algorithm of GLM estimations cannot be used directly when the predictors are non-linearly related.

Econometric estimations with fixest

Core estimation functions

Table (ref) summarizes the core fixest functions for performing econometric estimations:

The user-facing API for these six core functions is virtually identical, lowering the switching cost of moving between different estimation routines. For example, the estimation functions share many of the same arguments, of which the main ones are summarized in Table (ref). Among the listed arguments, fml and vcov deserve special mention, since they govern much of fixest's special behavior. We cover these two arguments in separate deep dives in Sections (ref) and (ref), respectively.

In addition to the high-level estimation functions, fixest also offers low-level routines for OLS (feols.fit) and GLM (feglm.fit), respectively. Similar to base stats::(g)lm.fit, these low-level functions offer performance boosts at the cost of placing greater strictures on the inputs; i.e. they require a matrix of predictors alongside the dependent variable in vector form.

\makesavenoteenv{tabular} \makesavenoteenv{table}

The fixest object class and associated methods

Consider a simple estimation on the base R airquality dataset:

The core fixest estimation functions all return an object of class fixest.\footnote{An exception is the case of multiple estimation, which returns an object of (sister) class fixest_multi. However, the distinction is not important here. We revisit multiple estimation in Section (ref).} This object is a list containing many elements that are useful for post-estimation routines, including the model scores (first derivatives), the unscaled VCOV, the coefficients, etc. While it is possible to access these list elements directly---e.g, est\$nobs---that is not the recommended approach. Rather, \texttt{fixest} implements a wide array of dedicated methods and functions that can be applied to \texttt{fixest} objects. Users should thus invoke these methods---e.g., \texttt{nobs(est)}---since they provide additional convenience facilities and safeguards.

Methods from the stats package

The base R stats package implements a variety of generic statistical methods. Table (ref) lists the stats methods supported by fixest objects, sometimes with bespoke modifications.

Here is a simple example, where we use the update.fixest method to add another explanatory variable (Wind) to our previous model. We'll also save the resulting estimation for re-use:

As noted, fixest supports bespoke customizations of stats methods. Here is an illustration where we apply the predict.fixest method to the est2 object that we just created. The first two \texttt{predict} calls adopt standard arguments to generate the same output as from, say, an equivalent \texttt{update.lm} call. However, the third call demonstrates a bespoke \textbf{\texttt{fixest}} customization, where we request predicted values for the full original data sample (i.e, before any \texttt{NA} values have been removed).

fixest-specific methods and functions

In addition to the generic stats methods, fixest introduces a variety of bespoke methods and functions for dedicated use with fixest objects. These are summarized in Table (ref), with a further set of worked examples provided in the main text below.

\oldparagraph{Coefficients}

The goal of the fixest-specific function and methods is to supplement and, in many cases, expand upon the generic methods defined by base R. Consider the difference between the (generic) coef method and the more detailed coeftable output:

Readers will note the special attribute "vcov_type" associated with the coefficient table. Like many fixest functions, coeftable accepts a vcov argument for on-the-fly VCOV adjustment. We return to this idea in the VCOV deep-dive in Section (ref).

As detailed by our theory review in Section (ref), fixest views the fixed-effects as a set of nuisance parameters that need to be controlled for---as efficiently as possible---in order to derive valid results for the main parameters of interest. However, it can still be useful to interrogate the fixed-effects coefficients themselves (e.g., to see if they accord with intuition). These coefficients can be retrieved with the fixef method. Here we show that fixef yields identical coefficients to an equivalent, but less computationally efficient, estimation with dummy variables:

\oldparagraph{Data accessors}

The data-accessing methods of fixest deserve special attention because they enable important internal functionality, as well as provide building blocks for downstream packages that depend on fixest. Let us consider some brief examples, continuing with our est2 estimation object.

The airquality dataset contains missing values and, as a result, not all of the previous observations were used in the preceding estimations. We can extract the indices of the used observations with obs:

(Notice that the fifth observation is missing.)

Similarly, we can access the full data set with fixest_data, or the left-hand-side and right-hand-sides with model.matrix:

\oldparagraph{Fit statistics}

While fixest supports several functions for assessing model fit, the single most important one to remember is fitstat. The fitstat function provides a general interface for computing and accessing fit statistics for fixest estimations. It is tightly coupled with other \textbf{\texttt{fixest}} methods and functions like \texttt{print.fixest}, where it serves as the internal engine for retrieving relevant model fit information.

Consider an example where we want to report the following statistics from our est2 object: (a) the number of observations, (b) R2 (R-squared), (c) adjusted-R2, and (d) the Wald test of joint coefficient nullity.

fitstat returns an object of class fixest_fitstat which has a dedicated print method for pretty printing, as illustrated above. Since this object is a list, each of the statistics can be easily accessed programmatically. We can also pass the desired fit arguments as a one-sided formula instead of a character vector. In the example below, we again compute the R2 and Wald test, but this time assign it to an object. This facilitates convenient post-computation access:

As already mentioned, fitstat is tightly linked with other fixest methods and functions like print.fixest and etable. Indeed, many of these downstream functions accept a fitstat argument for customizing their return output. Consider a simple example, where we override the default \texttt{print.fixest} return and instead retrieve the (raw) R2 and Wald test as the required statistics:

Note that it is also possible to enable these customizations globally via the family of setFixest_* functions.\footnote{We have elided over them for brevity, but the setFixest_* functions (and their getFixest_* counterparts) provide convenient APIs for enabling and customizing many global settings. More information can be found in the package vignettes and relevant help files: ?setFixest_estimation, ?setFixest_etable, etc.} For example with the following command, the \texttt{print} methods always display the R2 and the Wald test:

This can be easily reverted to the factory default with:

Although we do not show it here, the exact same idea can be applied to etable, which also accepts a fitstat argument. We make two additional observations, which apply to both print.fixest and etable. First, when using a formula, you can add to the default fit statistics---rather than override them---by using a point (e.g., \textasciitilde.\ +\ wald). Second, when a particular fit statistic contains sub-elements, like the Wald test, you can refer to specific component directly (e.g, wald.p for the p-value, \texttt{wald.stat} for the statistic, etc.)

Finally, please note that many other fit statistics are available beyond those that we have used here; see fitstat(show_types = TRUE) for a full list. It is also very easy to add your own custom fit statistics with the function fitstat_register. Any custom statistics that are defined in this way will automatically be available for use with downstream methods and functions like print and etable. Again, we must omit examples for brevity, but readers may consult the ?fitstat_register documentation for more information.

\oldparagraph{Presentation}

We refrain from a discussion of the fixest presentation methods here, deferring instead to the deep-dive in Section (ref).

The fixest formula

Any fixest estimation must start with a formula describing the relationship between variables. This formula plays a critical role in fixest. Indeed, many features are implemented directly within the formula and cannot be otherwise accessed. In this section we describe how to: i) add fixed-effects, ii) estimate instrumental variable models, iii) implement multiple estimations, iv) interpolate variables, v) creating categorical variables and interactions, and vi) use fixest's panel-data features. All these features are enabled via the formula.

Fixed-effects

Table (ref) summarizes the various ways by which fixed-effects can be included in a formula. The table also provides equivalences to non-fixest syntax. While we provide this reference as a convenient “cheat sheet” guide, we expand on specific fixed-effects use cases in more detail below.

Individual fixed effects

To add fixed-effects in an estimation, you must add them after a pipe, \textbar, immediately after the independent variables. The basic syntax is:

\[ \mathtt{dep\_var} \sim \mathtt{indep\_vars} \mathtt{\ |\ } \mathtt{fixed\_effects} \]

Let us consider an example from the trade dataset, which is included in fixest.

Below we estimate a so-called “gravity” trade model, where the trade volume between two countries (in \texteuro) is modeled as a function of their geographic distance. Alongside controlling for product and year fixed-effects, note that we also include origin and destination fixed-effects to control for bilateral trade relationships. As shown in the third line of the output, there are 15 countries of origin and destination, and 20 different products.

Combined fixed-effects

Suppose that we believe the confounding effect of any particular product on the gravity equation depends on the original and destination country. fixest provides a special \^ (caret) operator to quickly combine fixest effects. Here we combine product with the origin and destination country, respectively:

Now there are 300 country-product coefficients (the number of origin countries times the number of products). The important increase in the number of regressors lead to an improvement of the fit statistics, as shown with the RMSE going down from 1.7 to 1.4.

Varying slopes

As described in Section (ref), fixest is able to handle varying slopes very efficiently. Again, we provide a special formula syntax to denote variables with varying slopes. Say fe is the fixed-effect and z1 and z2 are variables whose slopes vary according to fe. Then we specify the model as:

\[ \mathtt{dep\_var} \sim \mathtt{indep\_vars} \mathtt{\ |\ } \mathtt{fe[z_1, z_2]} \]

Continuing with our trade example, we can add a unique year trend for each country of origin as follows:

\[ \ln(Euros)_{ijt} = \alpha_i + \beta_j + t \times \theta_i + \gamma \ln(dist\_km) + \epsilon, \]

with \(i\) the index of the country of origin, \(j\) the country of destination, and \(t\) the year.\footnote{This estimation is analogous, albeit much faster and with an exact solution, to a mixed-effects model of form: lme4::lmer(log(Euros) \textasciitilde log(dist_km) + Year + (Year \textbar Origin) + (1 \textbar Destination),\ trade).}

Instrumental variables

While it is not available for all model families, fixest supports two stage least squares (2SLS) instrumental variables (IV) estimations for linear models in feols. Table (ref) provides a summary cheat sheet for how to perform 2SLS estimations, including without fixed-effects.

Similar to the fixed-effects component, the IV part of the formula is specified in its own slot after a \textbar. The syntax is:

\[ \mathtt{dep\_var} \sim \mathtt{exogenous\_vars} \mathtt{\ |\ } \mathtt{fixed\_effects} \mathtt{\ |\ } \mathtt{endo\_vars} \sim \mathtt{instruments} \]

Again, we emphasize the order importance of our multipart formula. The fixed-effect part is optional (and comes second). The IV part always comes always last. If there is no exogenous regressor, the part containing the exogenous regressors should be replaced with a single 1.s

In 2SLS estimations, by default only the second stage is reported. To access the first stage(s), use the argument stage of the summary method (Section (ref)). Let us use data from the Fulton Fish Market the introductory example from Section (ref), where we instrument science funding with a policy shock to estimate the causal effect on productivity. We display both the second stage and the first stage:

The function fitstat can access the following fit statistics when applied to instrumental variables estimations:

Categorical variables and interactions: i()

To tell an estimation procedure to handle numeric variables as categorical variables, as.factor is the common way to do it in R. fixest introduces a simple function i, to handle categorical variables and provides easy access to: i) reordering, ii) binning, iii) interactions.

Use the function i to treat any variable as a categorical variable. This specific use can be seen as a replacement for as.factor. When i is applied, the names of the categories is as follows: "var_name::value", whereas when using as.factor the new names are "as.factor(var_name)value". Further, the function \texttt{i} allows for richer integration with some \textbf{\texttt{fixest}} functions like \texttt{etable} or \texttt{iplot}.

In the example below, let us treat the numeric variable year as categorical, either using as.factor, either using i. With coef we then show i's specific formatting of the coefficient names. We end by showing the interaction with etable and how we end up with a display of greater clarity.

The function i offers two useful arguments to manipulate the categorical variable it creates:

Using the example above, let us first take Month\ =\ 9 as a reference (instead of 5), then run the same initial estimation but grouping the months 8 and 9, and finally display the results with etable.

As we can see, the first two estimations are the same, only the reference of the categorical variable has changed, leading to different coefficient estimates (but same R2, also note that \(76.90 - 11.35 = 65.55\)). In the last estimation where we grouped the coefficients for the two last months, we can see that the estimates for the other coefficients stay the same.

Now say we want to estimate yearly time trends for the EU and the US separately, we can do so with i(cat,\ numerical):

If we wanted to include indicators for region by year fixed effects with i(cat,\ i.cat2). The i. specifies that the second variable is intended to be treated as a categorical as well (borrowing the syntax from Stata). We could alternatively specify this in the fixed-effect portion with \textbar\ eu_us\^year.

Multiple estimations

A core design feature of fixest is that it should be easy to estimate multiple models (or functional forms) all at once in a single function call. In this section we focus on: i) multiple dependent variables and ii) stepwise regressions. Again, we provide a cheat sheet for fixest's multiple estimation syntax in Table (ref).

Multiple dependent variables

To estimate a model for multiple dependent variables, you can nest them within the c() function on the formula LHS:

\[ \mathtt{c(y_1, y_2)} \sim \mathtt{indep\_vars} \mathtt{\ |\ } \mathtt{fixed\_effects}. \]

This will run two regressions: one with \(\mathtt{y_1}\) as dependent variable and one with \(\mathtt{y_2}\). Here is a simple example using the built-in airquality dataset. Note that the default print method automatically reverts to etable, conveniently displaying both regressions alongside each other.

Beyond syntactic convenience, estimating multiple dependent outcomes in this way can yield dramatic performance benefits. The bulk of the computing time in a fixed-effects estimation is spent dealing with the fixed-effects themselves. In a multiple estimation, fixest is able to pool this computationally-expensive step across the various models, thus avoiding the effort duplication that would arise from estimating the models separately. The practical upside is that estimating, say, ten models will typically take the same amount of time as estimating a single model on its own.\footnote{An important caveat is that this pooling acceleration shortcut only applies to OLS. For GLM models, the more complex model structure means that the same shortcut cannot readily be used.} We benchmark these gains with a large dataset in Section (ref).

Stepwise regressions

Multiple estimation is also supported on the formula RHS. For example, all fixest estimations support a family of special functions for stepwise regressions: sw (stepwise), sw0 (stepwise, starting from the empty element), csw (cumulative stepwise), csw0 (cumulative stepwise starting with the empty element), and \texttt{mvsw} (multiverse stepwise, i.e. crossing). These functions accept any number of variables as argument, and can be used directly in a formula as regular variables. Moreover, they can be used in both the regular part of the formula RHS (containing the explanatory variables), or in the fixed-effects part. However, note that only stepwise function is allowed per part.

Here is an example where we run six estimations in a single feols call. Notice how expressive the syntax is, despite its compactness. Once you know that sw stands for stepwise and csw for cumulative \textbf{s}tep\textbf{w}ise, and that a \texttt{0} suffix enforces an empty starting element, the code is very easy to understand.

Formula interpolation

Another key design feature of fixest formulas is the ease of interpolating variables. Programmatic formula manipulation in R is a long-standing issue, for which users often resort to ad hoc solutions.\footnote{For example, as.formula(paste(...)). While stats::reformulate provides a mechanism for programmatic formula construction, it is limited to simple (single part) formulas and does not support interpolation.} In contrast, \textbf{\texttt{fixest}} support three primary ways to interpolate variables directly in the formula itself: i) with formula macros, ii) with regular expressions, and iii) with the dot-square-bracket operator.

A cheat sheet of all possible interpolations in Table (ref). We also work through some explicit examples below, and explain the role of the xpd function, which powers all of fixest's interpolation features under the hood.

Formula macros

Users can define (global) formula macros in fixest with the function setFixest_fml. While the function accepts any number of arguments, each argument name must start with two dots and the arguments themselves must be either a character vector or a one-sided formula.\footnote{The ..\textlessmacro\textgreater naming restriction ensures that macros can easily be distinguished from regular variables.} For example setFixest_fml(..controls\ =\ c("Wind",\ "Temp")) creates a macro variable named ..controls. This variable can now be used in any \textbf{\texttt{fixest}} estimation as shown below:

Users may create as many macros as they wish through setFixest_fml. Macros can also contain fixed-effects if they are created with a formula, e.g. ..controls\ =\ \textasciitildeWind\ \textbar\ Month. Macros can prove especially useful in analysis scripts, with clear definitions of recurring sets of variables. One caveat, though is that we don't recommend using macros in loops. For this use case, interpolation with the dot-square-bracket operator, later described, is better suited.

Regex interpolation

Another way to interpolate variables in formulas is through regular expressions. Use ..("regex") or regex("regex") to include in the formula any variable in the data set that matches the regular expression. When a variable name in the formula ends with two dots, any variable in the data set starting with the same letters are included. Here is an example of these two cases:

Dot-square-bracket

Finally, a powerful way to interpolate variables is with the dot-square-bracket operator. When expressions of the form .{[}x{]} are found in a formula, the value of the variable x is directly inserted in the formula. When x is a vector, the values are combined with a +.

For example, let us reimplement our formula macro via the dot-square-bracket operator:

This interpolation makes it very easy to work with formulas programmatically within loops.

As illustrated above, when the interpolation is on several variables, they are combined with a plus (i.e. if x\ =\ c("a",\ "b"), then .{[}x{]} becomes a\ +\ b in the formula). However, sometimes we may want to combine them with a comma, especially when inserting the variables as several arguments of a function. You can override the default behavior but adding a comma right after the opening bracket, as in .{[},x{]}. Let us estimate multiple models that way and access the coefficients:

xpd: generalized formula interpolation and expansion

Under the hood, all interpolation features in fixest are powered by the function xpd. This function ensures that interpolated formulas are expanded to their correct full form. For example, we could have used xpd to verify that the expanded formula in our previous example, prior to running the actual estimation:

As the reader may have noticed, this example also demonstrates that xpd can be used on “regular” formulas outside of the main fixest estimation functions. The only difference is that these native fixest functions automatically use xpd to expand their formulas internally and thus do not need an explicit \texttt{xpd} call.

The ?xpd documentation covers a variety of use-cases and features, which we do not recapitulate here. Instead, we content ourselves with a final example where the desired outcome is interpolating over variables with numerical suffixes. In this particular case, assume that we wish include the variables x1--x5 in a regression. When the dot-square-bracket operator is attached to a variable, this variable can be used as the interpolation prefix, as follows:

Formula features for panel data

The last features we will highlight in this section are functions that are designed to work with time-series or panel data. In panel data settings, where you have unit-time identifiers, it is often useful to quickly compute the leads, lags, and/or differences of variables. The fixest formula provides a family of convenience functions for computing these transformations on the fly: l, f and d. However, to access these functions \texttt{fixest} must first know what the unit and time identifiers are. This information can be passed either via the run-time \texttt{panel.id} argument, or by explicitly setting the dataset as a panel via the \texttt{panel} function.

Let us demonstrate with a quick example of the former approach. Here we specify the panel identifiers as formula to panel.id, which in turn allows us to concisely estimate the lead of the dependent variable (f(y)) on the dependent variable and its three first lags (l(x1,\ 0:3)).

It should be emphasized that fixest's panel functionality offers rich support for handling ambiguous and irregular timesteps, as well as missing observations. For example, users can pass a time.step argument to the panel function to delineate between daily data that excludes weekends. Again, we refer the interested reader to the package documentation (e.g., ?panel) for more examples.

Differences-in-differences

fixest also has a set of tools for estimating difference-in-differences (DiD) models, which are a particularly popular research design in applied econometrics settings huntington2021effect. Using the base_stagg data set provided in fixest, let us illustrate a DiD analysis. In this example we interact the variable period, used as categorical, with the variable \texttt{treat} and set the reference to the period 5, the period right before the treatment takes place. We also include unit and time fixed-effects.

Alternatively, there has been new proposals for robust DID estimators that are valid when treatment-timing is staggered. The sunab() function implements the method from sun2021estimating.

It is very easy to plot both of these estimates thanks to the iplot function.

The graph is represented in Figure (ref).

Adjusting the standard-errors

Next to the coefficients themselves, the variance-covariance matrix (VCOV) is arguably the most important element of any estimation. This is the arbiter of statistical significance and the determinant for any inferential claims about our results. Given its prominence, fixest aims to streamline the way VCOVs are selected through a simple and intuitive interface. Specifically, all of the core estimating functions, as well as many of the methods, accept a named vcov argument. This interface is tightly integrated with the rest of package and enables many synergies, such as accessing multiple VCOVs at estimation time or when exporting the results. Moreover, a core design principle of fixest is that it should be easy to change a model's VCOV post estimation, which in turn enables on-the-fly inference adjustments to standard-errors.

In this section, we first present fixest's built-in VCOVs and how to parameterize them. We then show how to integrate the package with other software to compute additional VCOVs.

Built-in VCOVs

fixest provides native support for a variety of VCOVs. These built-in VCOVs cover some of the most common adjustments that are required in applied econometrics research and also offer best-in-class performance. The full set of these built-in VCOVs are documented in Table (ref), including specification options and required arguments. For the sake of brevity, we list the keyword shortcuts for these built-in VCOVs below:

To invoke any of these built-in VCOVs, users need simply pass the relevant keyword to the vcov argument. As we have tried to emphasize, this can be done either at estimation time or as part of a post-estimation method. Many fixest functions accept the vcov argument and in the below example we demonstrate using the se function, which retrieves standard errors. First, we specify heteroskedastic-robust errors at estimation, before reverting to spherical errors for comparison in a subsequent call.

While this example is very simple, it demonstrates the ease with which fixest allows users to adjust VCOVs and, thus, standard errors. Moreover, the on-the-fly adjustment underscores a key fixest design principle: estimation is separate from inference. This principle enables user-level flexibility in addition to facilitating computational efficiency, since we only need to recompute the model scores instead of re-estimating the entire model.

Beyond keyword strings, the vcov argument also accepts various other inputs. Again, these are detailed in full in Table (ref). To quickly highlight one particularly common use-case in panel-based econometrics, users can specify clustered standard errors by passing a one-sided formula denoting which variable(s) to cluster on. For example:

Continuing with non-keyword inputs to vcov, note that some VCOVs may require additional information. For example, Driscoll-Kraay VCOVs require unit-time identifiers, while Conley VCOVs require geographic information. In general, fixest tries to deduce any additional VCOV arguments using simple but robust rules of thumb.\footnote{For example, feols(depth\ \textasciitilde\ mag,\ quakes,\ vcov\ =\ "conley") specifies Conley standard errors in the base R quakes dataset to allow for arbitrary spatial auto correlation. Normally this requires explicit latitude and longitude information. But \texttt{fixest} automatically deduces that this information is furnished by the \texttt{lat} and \texttt{long} columns of the dataset.} But users can also pass additional arguments explicitly in one of two ways. First, by specifying a special helper function, whose name is the VCOV keyword prefixed with \texttt{vcov_}. For example, \texttt{vcov_cluster} has the argument \texttt{cluster} accepting a one-sided formula or a character scalar specifying the cluster variable. Second, by invoking a formula with the following syntax:

\[ \mathtt{vcov} = \mathtt{VCOV\_name} \sim \mathtt{variables}. \]

Here, VCOV_name is a function that takes additional arguments. For example, NW and DK accept the argument lag, and conley accept the arguments cutoff, pixel and \texttt{distance}.

Re-using the same estimation from above, let us report the Newey-West and clustered standard-errors using these two approaches. This time, we'll pass the key vcov argument a list of VCOVs from an etable call.

As a final remark about the built-in VCOVs, fixest applies a small sample correction of \(N / (N - K)\), where \(N\) is the number of observations and \(K\) the number of variables. However, this small sample correction can be disabled or adjusted in various ways via the ssc function. Small sample corrections can play a surprisingly prominent role in lining up results from different econometrics software, particularly in the context of clustered VCOVs (where multiple corrections are potentially valid from a theory perspective). We omit further discussion for the sake of brevity here. But the interested reader is referred to the dedicated “\href{https://lrberge.github.io/fixest/articles/standard_errors.html}{On standard errors}” vignette for an exhaustive discussion and comparison between fixest and other popular panel econometrics packages.

Access to external VCOVs

Beyond the built-in VCOVs, the vcov argument also accepts custom functions and VCOVs from packages. For example, fixest has a relatively good compatibility with the well-known sandwich package zeileis2006ObjectorientedComputationSandwich, zeileis2020VariousVersatileVariances, which specializes in computing the variance-covariance matrices. However, there are a few important caveats. Some functions from \texttt{sandwich} do not work well with \textbf{\texttt{fixest}} estimations when fixed-effects are present, due to the way that these fixed-effects are handled. Similarly, while most \textbf{\texttt{sandwich}} functions have a \texttt{fixest} method, there is not a drop-in replacement for \texttt{vcovBS()}, \texttt{vcovJK()}, and \texttt{vcovOPG()}. For these reasons, we recommend that users rely on native \textbf{\texttt{fixest}} VCOVs where possible.

Users can also define their custom VCOVs and pass them to summary, etable, coeftable, etc. In this next example, we compute a bootstrapped VCOV and display the associated results.\footnote{As an aside, fixest provides a dedicated estimation environment that, among other things, allows for extremely efficient Bayesian bootstrap computation. See the ?est_env documentation for an example.} Note that we pass our custom VCOV through a named list for pretty printing and display.

Presentation

In this section, we discuss in turn how to display the results in a table or in a graph.

Formatting the results into a table

As we have already seen from the preceding examples, the etable function allows us to export the results from (multiple) fixest estimations to a nicely formatted regression table. The table can either be displayed directly in the R console or exported as a \LaTeX table.

It has the following main arguments:

The factory defaults of etable are usually good enough for displaying results in the R console and there is little need for further customization. However, the situation is more complicated when creating \LaTeX publication table, since these often require fine-grained control and detailing. The default \LaTeX table produced by etable is quite plain, yet contains all of the necessary information to uniquely identify the context of an estimation, as well as interpret it. Specifically, we report: the sample (if the sample was split), the type of standard errors, the family when several estimations have different families (e.g., OLS or Poisson), convergence state when one model did not converge, etc. We can therefore think of the default etable output as appropriate for an exploratory table, where we want as much information as possible without worrying too much about aesthetics. Table (ref), from our motivating example back in Section (ref), provides an example of such a table.

For published tables, the ability to style, as well as precisely add pieces of information, take on heightened importance. etable caters to the needs of publication-focused users by offering over 60 arguments to customize their table. Striking the balance for usability here requires nuance, since the number of functions can be overwhelming. However, the basic philosophy is to simply set your table preferences once and then forget about them. Thanks to the function setFixest_etable() it is possible to set the styling arguments globally. Further, the function setFixest_dict() sets the labels of the arguments globally.

To illustrate, let us revisit our earlier Table (ref) example. We re-use the saved est_split and est_iv estimations from before and our etable call also remains unchanged (we pass the exact same arguments). Our only changes are i) adding labels to the variables with setFixest_dict, and ii) changing the default \LaTeX style with setFixest_etable. The resulting table is displayed in Table (ref) and reveal a number of aesthetic improvements.

There are many etable options beyond the ones presented here, including special integrations for R Markdown and Quarto documents, as well as support for coefficient highlighting, three-part tables, etc. The fixest website provides further details on how to customize regression tables in the form of standalone documentation (e.g. \href{https://lrberge.github.io/fixest/reference/etable.html}{?etable}, \href{https://lrberge.github.io/fixest/reference/style.tex.html}{?style.tex}) and two vignettes (“\href{https://lrberge.github.io/fixest/articles/exporting_tables.html}{Exporting estimation tables}” and “\href{https://lrberge.github.io/fixest/articles/etable_new_features.html}{etable: new features in \texttt{fixest} 0.10.2}”).

Plotting the coefficients

Alongside tabling support, fixest also provides two dedicated plotting functions for plotting the coefficients and confidence intervals: coefplot and iplot. coefplot plots the point estimates and confidence intervals of coefficients from one or several \texttt{fixest} estimations. \texttt{iplot} does the same, but limited to categorical variables or interaction terms created with the function \texttt{i()}; see Section (ref).

Some of the main arguments are:

Paralleling etable, users can set any of these parameters globally via the setFixest_coefplot() function, and thus define a preferred aesthetic for all of their plots.

In Figure (ref) below, we combine several examples to demonstrate fixest's plotting facilities. We start by performing a single estimation (with "iid" errors) and then obtaining a second version with clustered standard errors.

Our first graph (top left) is a simple coefplot call, demonstrating the default behavior and the ease with which multiple estimations are passed.

The second graph (top right) is a tweaked version of the first, flipping the plot axes and passing various other customizations arguments.

The third graph (bottom left) is the same data plotted with iplot and its default values, which then focuses only on the categorical variable that was included with i(). Moreover, instead of having two estimation objects with different standard-errors, here we demonstrate that users can also have only one estimation object and pass several VCOVs using the argument vcov:

The last graph (bottom right) reports the coefficients of an event study plotted with iplot and a custom style displaying a ribbon:

Even though coefplot/iplot accommodate many options to customize the graphs, it may not be enough to create very precise publication graphs. In that case, the user can access the dataset used to build the graph (without plotting) by invoking the only.params\ =\ TRUE argument. Similarly, these fixest plotting functions yield base R graphics, which may or may not be to the user's taste. The \texttt{ggfixest} package ggfixest offers equivalent functionality in the form of dedicated \textbf{\texttt{ggplot2}} wickham2016Ggplot2ElegantGraphics methods for \textbf{\texttt{fixest}} objects, as well as a set of ancillary plotting features.

Benchmarks

In this section, we present a set of benchmarks that demonstrate fixest's performance against alternative state of the art software, across a broad range of settings. For the sake of parsimony and reproducibility, we restrict our comparisons to free and open-source software (FOSS) which specialize in fixed-effects estimations.\footnote{The alternatives that we do not benchmark here, including base R routines such as (g)lm and proprietary software such as Stata, are generally slower than our comparison FOSS group for this set of tasks.} Specifically, we benchmark \texttt{fixest} against \textbf{\texttt{alpaca}} stammann2018FastFeasibleEstimation and \textbf{\texttt{lfe}} gaure2013lfe in \textsc{R}, \textbf{\texttt{FixedEffectModels}} gomez2024FixedEffectModelsJlFast and \textbf{\texttt{GLFixedEffectModels}} boehm2025GLFixedEffectModelsJl in \textsc{Julia}, and \textbf{\texttt{PyFixest}} fischer2024Pyfixest in \textsc{Python}.\footnote{For the sake of clarification, please note that \textbf{\texttt{PyFixest}} is developed by independent authors. As the name suggests, it is a port of \textbf{\texttt{fixest}}'s syntax and features into \textsc{Python}, although some of the core internal routines are different; as we documented in Section (ref). The \textbf{\texttt{PyFixest}} team recently began work to replace their existing alternating projections routine with \textbf{\texttt{fixest}}'s fixed-point algorithm. This means that the performance of these two libraries should become closer in the near future, as they will both rely on the same method and implementation.} We use the latest available versions of these software, as of 15 January 2026. The benchmarks presented here are run on a Windows PC with 128GB of RAM and an Intel\textsuperscript{\tiny\textregistered} Core\textsuperscript{TM} i7-10850H CPU @ 2.70GHz processor. Separate benchmarks run on an Apple Mac M4 Pro laptop with 48 GB of RAM yield outcomes that are proportionally very similar.

Before continuing, we issue the standard disclaimer that benchmarks should always be interpreted with a degree of caution, given the inherent sensitivity to arbitrary data and estimation choices. Nonetheless, we believe that we provide a fair and representative overview, which seeks to minimize bias by running multiple benchmarks over many use-cases, using both simulated and real-life data. In particular, we present two broad sets of benchmarks. The first uses simulated datasets with different convergence properties, providing insight into how the performance of our different software scales with the number of observations. The second evaluates performance on a demanding real-life dataset, where we vary the type of model estimation and focus on specific advantages provided by fixest's feature set. All benchmarking code is provided with the replication packet for this paper.

Simulated data

Our simulated dataset draws inspiration from labor economics by generating a panel of employees and firms over time. Our specific regression task is estimating wages as a function of covariates (like training), whilst controlling for individual, firm and year fixed-effects. Importantly, we consider two data generating processes (DGPs) for our simulated data: one for which the convergence of the estimated fixed-effects is “simple” and another for which it is “difficult”. The next two paragraphs provide intuition for why convergence can be simple (fast) or difficult (slow). Readers uninterested in the intuition can safely skip these paragraphs and jump straight to the results in Figures (ref) and (ref).

Consider a balanced panel of employees. The assignation of firms to employees is the key determinant of convergence speed, which in turn elevates the importance of algorithmic routines for estimating the individual and firm fixed-effects. In the “simple” case, we assign firms randomly to employees, resulting in each firm being connected to many different employees. From a network perspective, this creates a dense network where firms share common employees, which in turn means that any changes to a single firm coefficient will propagate quickly through the system, enabling fast convergence. In contrast, for the “difficult” case, we assign firms sequentially to employees, resulting in firms being connected to very few employees. Now our network is extremely sparse with few connections between firms, which means that changes to firm coefficients are slow to propagate, yielding slow convergence.

Our DGP is as follows. For \(N\) observations, we generate a panel of \(N_I = N/N_T\) individuals over \(N_T = 10\) years. These employees work across \(N_F = N_I / 23\) companies. In the simple case, we assign companies randomly to employees. In the difficult case, we assign them sequentially. All the fixed-effect coefficients are generated with a standard normal distribution. There are two covariates \(x_1\), generated from a standard normal distribution, and \(x_2 = x_1^2\). Finally the outcome, \(y\) is obtained as follows: \[ y_{ift} = x_{1, ift} + 0.05 x_{2, ift} + \mu_i + \lambda_t + \nu_f + \varepsilon_{ift}, \]

where \(\varepsilon_{ift}\) is drawn from a standard normal distribution. We report below the code to generate the data:

Our simulated benchmarks evaluate both the OLS and Poisson GLM estimators. Specifically, the benchmarks estimate a set of models where: i) the dependent variable is \(y\) for OLS and \(exp(y)\) for Poisson, ii) the explanatory variables are \(x_1\) and \(x_2\), iii) there are two or three sets of fixed-effects, and iv) the convergence of the fixed-effects can be “simple” or “difficult”. Two fixed-effects correspond to individual and firms, while three fixed-effects adds the year. Finally we also vary the number of observations from 10,000 to 10M.

The results of the OLS benchmark are represented in Figure (ref). It represents the average computing time across three replications, in logarithmic scale. The top row represents the results for the “simple” case. There we see that fixest::feols is consistently faster than the alternative software, with a computing time of about 1s for 10M observations and two fixed-effects, while on the other end of the scale lfe::felm takes an order of magnitude longer at 10s.

An interesting observation about this “simple” case is that there are effectively no scaling effects; the relative differences between software do not meaningfully change with the number of observations. However, this is not true for our second, “difficult” set of results, which we report in the bottom row of Figure (ref). Both fixest and Julia's FixedEffectModels continue to scale almost linearly as the dataset grows. Somewhat remarkably, these two implementations are effectively identical for our most difficult case of three-fixed effects and 10M observations at just over 60s. In contrast, the convergence times of \texttt{lfe} and \textbf{\texttt{PyFixest}} explode in this more difficult setting. We do not report benchmark times for these two libraries above 100,000 observations, since neither could fit the model within our constraints.\footnote{Either erroring due to absence of convergence, or exceeding our 1hr timeout.}

Our Poisson benchmarks are reported in Figure (ref) and follow much the same pattern as the OLS results. Again, we observe that fixest times are consistently faster across the board, with good scaling properties for both the “simple” (top) and “difficult” (bottom) DGP cases. GLFixedEffectModels is nearest in performance and also scales well across our different settings. In contrast, alpaca and \textbf{\texttt{PyFixest}} struggle in the difficult convergence case, and were unable to complete within our 1hr timeout constraint for 1M or more observations.

NYC taxi data

Having evaluated the impact that different convergence properties can have on simulated data, we now illustrate fixest's performance on a demanding real-life example. Specifically, we use the well-known NYC Taxi and Limousine Commission dataset.\footnote{See: \url{https://www.nyc.gov/site/tlc/about/tlc-trip-record-data.page}.} Our sample covers three months of taxi trip data from Jan--Mar 2012 and is rather large at over 42M observations.

The results, which are the average timings across three runs, are reported in Table (ref), where we run four basic specifications that are intended to highlight fixest-specific features:

Starting with the baseline estimation in the top row of Table (ref), fixest is fastest at only 9.9s, followed by FixedEffectModels at 18s and then PyFixest at 23s.

Adding varying slopes to one explanatory variable reduces the estimation time for fixest and FixedEffectModels to 8.9s and 15s respectively, while this feature is unavailable in PyFixest.

Estimating two outcomes simultaneously adds less than 4s in fixest, thanks to its highly efficient handling of multiple estimations. Since this feature in unavailable in FixedEffectModels, the estimation time is effectively doubled to 35s. PyFixest also handles multiple estimations, although their computation time increases by 80% to 43s.

Finally, as we have emphasized, fixest's design separates estimation from inference. This means that estimating two models with different VCOVs (standard-errors) incurs virtually no overhead. Conversely, for FixedEffectModels the time is doubled. PyFixest also separates inference from estimation, albeit not as efficiently, thus its timing increases by 60%.

Conclusion and discussion

The fixest package represents a mature econometric toolkit that we hope serves diverse user communities effectively. While this article has focused on the core user-facing features, we also wish to highlight two additional user groups: package developers and educators.

For package developers, fixest provides a robust foundation with numerous convenience features. Many econometric methods are implemented via a regression model and, by relying on fixest for estimation, package developers automatically gain access to all of these features. This allows developers to focus on methodological innovations, rather than worrying about speed or re-implementing standard features like fixed-effects handling, robust standard errors, plotting, tabling, and basic method support. Moreover, the fixest ecosystem extends seamlessly through integrations with \textbf{\texttt{broom}} robinson2025BroomConvertStatistical, \textbf{\texttt{parameters}} parameters, \textbf{\texttt{marginaleffects}} arel-bundock2024marginaleffects, \textbf{\texttt{gtsummary}} gtsummary, and \textbf{\texttt{modelsummary}} arel-bundock2022modelsummary, providing even more functionality that package developers can leverage.

Second, we believe that fixest can prove useful for educators. The range of fixest features are intended to cover large swathes of the applied econometrics toolkit. More to the point, incorporating all of these features into a single package with a unified syntax allows instructors to focus on methodology, rather than software proliferation. As a leading example, consider The Effect textbook huntington2021effect, which aims to teach causal inference at the undergraduate level. Code covering basic regression models, robust inference, instrumental variables, fixed-effects, DiD, and regression discontinuity all are written using the single function feols.

We would like to conclude with a few words on fixest's development philosophy. This software is built on three cornerstones: speed, robustness, and stability.

On speed. All data intensive operations in fixest have been implemented with as many optimizations as possible, involving methodological innovations when needed. This is highlighted in the benchmarks for the core estimating functions, but the same philosophy extends to many other functions.

On robustness. Error-handling is central to fixest since it is critical to a good user experience. We have tried our best to anticipate any potential problems from the user side, even unlikely ones, and the fixest codebase contains hundreds, possibly thousands, of custom error messages and exceptions. For the remaining bugs that are inevitably left over, we could count on the numerous feedback from a fantastic community that has helped iron them out. fixest is by now a mature package, forged by battle-hardened testing over several years and across an array of contexts.

On Stability. As core developers and users of fixest, we value stability very highly. Any design choices made in consideration of the codebase and user-facing API are done with that desideratum in mind. We would make three observations in support of this claim. First, fixest has no hard dependencies that are at risk themselves. The core codebase relies only on three non-base R packages, namely \texttt{Rcpp} rcpp2024software, \textbf{\texttt{stringmagic}} berge2025stringmagic and \textbf{\texttt{dreamerr}} berge2025dreamerr. The former already acts as the solid foundation of thousands of \textsc{R} packages, while the latter two are spinoffs from \textbf{\texttt{fixest}} itself. Second, only a handful of breaking changes to non-core \textbf{\texttt{fixest}} functions have been introduced over its seven years of existence. This, despite the introduction of many new features and considerable growth of the codebase. Third, we strive to avoid making rushed design decisions. Rather, we introspect carefully and then seek feedback to ensure that any new features are as stable and future proof as they can.

By continuing to build on these foundations, we hope that fixest will prove to be a valuable and reliable companion for empirical work in the years to come.

Acknowledgments

Laurent Bergé wishes to thank Grant McDermott whose early publicity and genuine enthusiasm led to unforeseen levels of development, Dirk Eddelbuettel for creating and maintaining Rcpp (without this key technology fixest would have never existed), the community which made this software infinitely more robust thanks to bug reports, and the R authors and the CRAN team for making such an unbelievably reliable, and fun, programming language.