The package fixest provides a family of functions to perform estimations with multiple fixed-effects. The two main functions are feols for linear models and feglm for generalized linear models. In addition, the function femlm performs direct maximum likelihood estimation, and feNmlm extends the latter to allow the inclusion of non-linear in parameters right-hand-sides. Each of these functions supports any number of fixed-effects and is implemented with full fledged multi-threading in c++. Functions feols and feglm further support variables with varying slopes.

This package is currently (Nov. 2019) the fastest software available to perform fixed-effects estimations (see the project’s homepage for a benchmarking).

The standard-errors of the estimates can be easily and intuitively clustered (up to four-way).

Two specific functions are implemented to seamlessly export the results of multiple estimations into either a data.frame (function esttable), or a Latex table of “article-like” quality (function esttex).

The main features of the package are illustrated in this vignette. The theory used to obtain the fixed-effects is based on Berge (2018), “Efficient estimation of maximum likelihood models with multiple fixed-effects: the R package fixest.” CREA Discussion Papers, 13 (https://wwwen.uni.lu/content/download/110162/1299525/file/2018_13).

1 Simple example using Trade data

This example deals with international trade, which is a setup that usually requires performing estimations with many fixed-effects. We estimate a very simple gravity model in which we are interested in finding out the negative effect of geographic distance on trade. The sample data consists of European trade extracted from Eurostat. Let’s load the data contained in the package:

library(fixest)
data(trade)

This data is a sample of bilateral importations between EU15 countries from 2007 and 2016. The data is further broken down according to 20 product categories. Here is a sample of the data:

Destination Origin Product Year dist_km Euros
LU BE 1 2007 139.5719 2966697
BE LU 1 2007 139.5719 6755030
LU BE 2 2007 139.5719 57078782
BE LU 2 2007 139.5719 7117406
LU BE 3 2007 139.5719 17379821
BE LU 3 2007 139.5719 2622254

The dependent variable of the estimation will be the level of trade between two countries while the independent variable is the geographic distance between the two countries. To obtain the elasticity of geographic distance net of the effects of the four clusters, we estimate the following:

\(E\left(Trade_{i,j,p,t}\right)=\gamma_{i}^{Exporter}\times\gamma_{j}^{Importer}\times\gamma_{p}^{Product}\times\gamma_{t}^{Year}\times Distance_{ij}^{\beta}\),

where the subscripts \(i\), \(j\), \(p\) and \(t\) stand respectively for the exporting country, the importing country, the type of product and the year, and the \(\gamma_{v}^{c}\) are fixed-effects for these groups. Here \(\beta\) is the elasticity of interest.

Note that when you use the Poisson/Negative Binomial families, this relationship is in fact linear because the right hand side is exponentialized to avoid negative values for the Poisson parameter. This leads to the equivalent relation:1

\(E\left(Trade_{i,j,p,t}\right)=\exp\left(\gamma_{i}^{Exporter}+\gamma_{j}^{Importer}+\gamma_{p}^{Product}+\gamma_{t}^{Year}+\beta\times \ln Distance_{ij}\right)\).

1.1 Estimation

The estimation of this model using a Poisson likelihood is as follows:

gravity_results <- feglm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade)

Note that you need not provide the argument family since the Poisson model is the default.

The results can be shown directly with the print method:

print(gravity_results)
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Fixed-effects: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors: Clustered (Origin) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5279   0.115699 -13.206 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -7.025e+11   Adj. Pseudo-R2: 0.76403 
#>            BIC:  1.405e+12     Squared Cor.: 0.61202

The print reports the coefficient estimates and standard-errors as well as some other information. Among the quality of fit information, the squared-correlation corresponds to the correlation between the dependent variable and the expected predictor; it reflects somehow to the idea of R-square in OLS estimations.

1.2 Clustering the standard-errors

To cluster the standard-errors, we can simply use the argument se of the summary method. Let’s say we want to cluster the standard-errors according to the first two clusters (i.e. the Origin and Destination variables). Then we just have to do:

summary(gravity_results, se = "twoway")
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Fixed-effects: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors: Two-way (Origin & Destination) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5279   0.132276 -11.551 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -7.025e+11   Adj. Pseudo-R2: 0.76403 
#>            BIC:  1.405e+12     Squared Cor.: 0.61202

The clustering can be done on one (se="cluster"), two (se="twoway"), three (se="threeway") or up to four (se="fourway") variables. If the estimation includes fixed-effects, then by default the clustering will be done using these fixed-effects, in the original order. This is why the Origin and Destination variables were used for the two-way clustering in the previous example. If, instead, you wanted to perform one-way clustering on the Product variable, you need to use the argument cluster:

# Equivalent ways of clustering the SEs:
# One-way clustering is deduced from the arguent 'cluster'
# - using the vector:
summary(gravity_results, cluster = trade$Product)
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Fixed-effects: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors: Clustered 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5279   0.098318  -15.54 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -7.025e+11   Adj. Pseudo-R2: 0.76403 
#>            BIC:  1.405e+12     Squared Cor.: 0.61202
# - by reference:
summary(gravity_results, cluster = "Product")
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Fixed-effects: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors: Clustered (Product) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5279   0.098318  -15.54 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -7.025e+11   Adj. Pseudo-R2: 0.76403 
#>            BIC:  1.405e+12     Squared Cor.: 0.61202
# - with a formula:
summary(gravity_results, cluster = ~Product)
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Fixed-effects: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors: Clustered (Product) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5279   0.098318  -15.54 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -7.025e+11   Adj. Pseudo-R2: 0.76403 
#>            BIC:  1.405e+12     Squared Cor.: 0.61202

Note that you can always cluster the standard-errors, even when the estimation contained no fixed-effect. Buth then you must use the argument cluster:

gravity_simple = feglm(Euros ~ log(dist_km), trade)
# Two way clustering is deduced from the argument 'cluster'
# Using data:
summary(gravity_simple, cluster = trade[, c("Origin", "Destination")])
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Standard-errors: Two-way 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)    24.709    1.15980 21.3040 < 2.2e-16 ***
#> log(dist_km)   -1.029    0.16316 -6.3064  2.86e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -2.426e+12   Adj. Pseudo-R2: 0.18502 
#>            BIC:  4.852e+12     Squared Cor.: 0.05511
# Using a formula (note that the values of the variables are 
#  fetched directly in the original database):
summary(gravity_simple, cluster = ~Origin+Destination)
#> GLM estimation, family = poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Standard-errors: Two-way (Origin & Destination) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)    24.709    1.15980 21.3040 < 2.2e-16 ***
#> log(dist_km)   -1.029    0.16316 -6.3064  2.86e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -2.426e+12   Adj. Pseudo-R2: 0.18502 
#>            BIC:  4.852e+12     Squared Cor.: 0.05511

1.3 Other estimation functions

Now we estimate the same relationship by OLS. We need to put the left hand side in logarithm (since the right-hand-side is not exponentialized):

gravity_results_ols <- feols(log(Euros) ~ log(dist_km)|Origin+Destination+Product+Year, trade)

Of course you can use different families in feglm, exactly as in glm.

To get the estimation for the fixed-effects Negative Binomial:

gravity_results_negbin <- fenegbin(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade)

1.4 Viewing the results in R

Now let’s say that we want a compact overview of the results of several estimations. The best way is to use the function esttable. This function summarizes the results of several fixest estimations into a data.frame. To see the fixed-effects results with the three different likelihoods, we just have to type:

esttable(gravity_results, gravity_results_negbin, gravity_results_ols,
         se = "twoway", titles = c("Poisson", "Negative Binomial", "Gaussian"))
Poisson Negative Binomial Gaussian
Dependent Var.: Euros Euros log(Euros)
log(dist_km) -1.5279*** (0.1323) -1.7108*** (0.1797) -2.1699*** (0.1739)
Overdispersion: 0.548774
Fixed-Effects: ——————- ——————- ——————-
Origin Yes Yes Yes
Destination Yes Yes Yes
Product Yes Yes Yes
Year Yes Yes Yes
___________________ ___________________ ___________________ ___________________
Family: Poisson Neg. Bin. OLS
Observations 38,325 38,325 38,325
Squared Corr. 0.612 0.438 0.706
Pseudo R2 0.76403 0.03473 0.2364
BIC 1.405e+12 1,294,419.32 152,589.34

We added the argument se="twoway" to cluster the standard-errors for all estimations. As can be seen this function gives an overview of the estimates and standard-errors, as well as some quality of fit measures. The argument titles is used to add information on each estimation column.

In the previous example, we directly added the estimation results as arguments of the function esttable. But the function also accepts lists of estimations. Let’s give an example. Say you want to see the influence of the introduction of fixed-effects on the estimate of the elasticity of distance. You can do it with the following code where we use the argument fixef to include fixed-effects (instead of inserting them directly in the formula):

gravity_subcluster = list()
all_clusters = c("Year", "Destination", "Origin", "Product")
for(i in 1:4){
    gravity_subcluster[[i]] = feglm(Euros ~ log(dist_km), trade, fixef = all_clusters[1:i])
}

The previous code performs 4 estimations with an increasing number of fixed-effects and store their results into the list named gravity_subcluster. To show the results of all 4 estimations, it’s easy:

esttable(gravity_subcluster, cluster = ~Origin+Destination)
model 1 model 2 model 3 model 4
log(dist_km) -1.0293*** (0.1632) -1.2257*** (0.2084) -1.5176*** (0.1297) -1.5279*** (0.1323)
Fixed-Effects: ——————- ——————- ——————- ——————-
Year Yes Yes Yes Yes
Destination No Yes Yes Yes
Origin No No Yes Yes
Product No No No Yes
___________________ ___________________ ___________________ ___________________ ___________________
Observations 38,325 38,325 38,325 38,325
Squared Corr. 0.057 0.164 0.385 0.612
Pseudo R2 0.18833 0.35826 0.59312 0.76403
BIC 4.833e+12 3.821e+12 2.423e+12 1.405e+12

We have a view of the 4 estimations, all reporting two-way clustered standard-errors thanks to the use of the argument cluster.

1.5 Exporting the results to Latex

So far we have seen how to report the results of multiple estimations on the R console. Now, with the function esttex, we can export the results to high quality Latex tables. The function esttex works exactly as the function esttable, it takes any number of fixest estimations. By default, it reports Latex code on the R console:

# with two-way clustered SEs
esttex(gravity_subcluster, cluster = ~Origin+Destination)
#> \begin{table}[htbp]\centering
#> \caption{no title}
#> \begin{tabular}{lcccc}
#>  & & & & \tabularnewline
#> \hline
#> \hline
#> Dependent Variable:&\multicolumn{4}{c}{Euros}\\
#> Model:&(1)&(2)&(3)&(4)\\
#> \hline
#> \emph{Variables}\tabularnewline
#> log(dist\_km)&-1.0293$^{***}$&-1.2257$^{***}$&-1.5176$^{***}$&-1.5279$^{***}$\\
#>   &(0.1632)&(0.2084)&(0.1297)&(0.1323)\\
#> \hline
#> \emph{Fixed-Effects}&  & & & \\
#> Year&Yes&Yes&Yes&Yes\\
#> Destination&No&Yes&Yes&Yes\\
#> Origin&No&No&Yes&Yes\\
#> Product&No&No&No&Yes\\
#> \hline
#> \emph{Fit statistics}&  & & & \\
#> Observations& 38,325&38,325&38,325&38,325\\
#> Squared Correlation & 0.057&0.164&0.385&0.612\\
#> Pseudo R$^2$ & 0.18833&0.35826&0.59312&0.76403\\
#> BIC & $4.833\times 10^{12}$&$3.821\times 10^{12}$&$2.423\times 10^{12}$&$1.405\times 10^{12}$\\
#> \hline
#> \hline
#> \multicolumn{5}{l}{\emph{Two-way (Origin & Destination) standard-errors in parentheses. Signif Codes: ***: 0.01, **: 0.05, *: 0.1}}\\
#> \end{tabular}
#> \end{table}

This function has many optional arguments. The user can export the Latex table directly into a file (argument file), add a title (arg. title) and a label to the table (arg. label).

The coefficients can be renamed easily (arg. dict), some can be dropped (arg. drop) and they can be easily reordered with regular expressions (arg. order).

The significance codes can easily be changed (arg. signifCode) and all quality of fit information can be customized. Among others, the number of fixed-effect per cluster can also be displayed using the argument showClusterSize.

1.5.1 An elaborate example

Consider the following example of the exportation of two tables:

# we set the dictionary once and for all
myDict = c("log(dist_km)" = "$\\ln (Distance)$", "(Intercept)" = "Constant")
# 1st export: we change the signif code and drop the intercept
esttex(gravity_subcluster, signifCode = c("a" = 0.01, "b" = 0.05),
       drop = "Int", dict = myDict, file = "Estimation Table.tex", 
       replace = TRUE, title = "First export -- normal Standard-errors")
# 2nd export: clustered S-E + distance as the first coefficient
esttex(gravity_subcluster, se = "cluster", cluster = ~Product, order = "dist", 
       dict = myDict, file = "Estimation Table.tex", 
       title = "Second export -- clustered standard-errors (on Product variable)")

In this example, two tables containing the results of the 5 estimations are directly exported in the file “Estimation Table.tex”. The file is re-created in the first exportation thanks to the argument replace = TRUE.

To change the variable names in the Latex table, we use the argument dict. The variable myDict is the dictionary we use to rename the variables, it is simply a named vector. The original name of the variables correspond to the names of myDict while the new names of the variables are the values of this vector. Any variable that matches the names of myDict will be replaced by its value. Thus we do not care of the order of appearance of the variables in the estimation results.

In the first export, the coefficient of the intercept is dropped by using drop = "Int" (could be anything such that grepl(drop[1], "(Intercept)") is TRUE). In the second, the coefficient of the distance is put before the intercept (which is kept). Note that the actions performed by the arguments drop or order are performed before the renaming takes place with the argument dict.

1.6 Extracting the fixed-effects

To obtain the fixed-effects of the estimation, the function fixef must be performed on the results. This function returns a list containing the fixed-effects coefficients for each dimension. The summary method helps to have a quick overview:

fixedEffects <- fixef(gravity_results)
summary(fixedEffects)
#> Fixed_effects coefficients
#>                         Origin Destination Product   Year
#> Number of fixed-effects     15          15      20     10
#> Number of references         0           1       1      1
#> Mean                      23.3        3.09  0.0129  0.157
#> Variance                  1.63        1.23    1.86 0.0129
#> 
#> COEFFICIENTS:
#>   Origin:    SE    PT    NL    LU    IT                 
#>           23.25 22.44 24.43 20.23 24.33 ... 10 remaining
#> -----
#>   Destination:   SE    PT    NL LU    IT                 
#>                3.57 2.552 3.231  0 4.218 ... 10 remaining
#> -----
#>   Product: 1     2      3     4      5                 
#>            0 1.414 0.6562 1.449 -1.521 ... 15 remaining
#> -----
#>   Year: 2007    2008     2009    2010  2011                
#>            0 0.06912 0.005225 0.07331 0.163 ... 5 remaining

We can see that the fixed-effects are balanced across the dimensions. Indeed, apart from the first dimension, only one coefficient per fixed-effect needs to be set as reference (i.e. fixed to 0) to avoid collinearity across the fixed-effects of the different clusters. This ensures that the fixed-effects coefficients can be compared within cluster. Had there be strictly more than one reference per cluster, their interpretation would have not been possible at all. If this was the case though, a warning message would have been prompted. Note that the mean values are meaningless per se, but give a reference points to which compare the fixed-effects within a cluster. Let’s look specifically at the Year fixed-effects:

fixedEffects$Year
#>        2007        2008        2009        2010        2011        2012 
#> 0.000000000 0.069122284 0.005225473 0.073308208 0.163013386 0.192605170 
#>        2013        2014        2015        2016 
#> 0.230629376 0.242605404 0.282800683 0.310325692

Finally, the plot method helps to distinguish the most notable fixed-effects:

plot(fixedEffects)

For each cluster, the fixed-effects are first centered, then sorted, and finally the most notable (i.e. highest and lowest) are reported. The exponential of the coefficient is reported in the right hand side to simplify the interpretation for models with log-link (as the Poisson model). As we can see from the country of destination cluster, trade involving France (FR), Italy (IT) and Germany (DE) as destination countries is more than 2.7 times higher than the EU15 average. Further, the highest heterogeneity come from the product category, where trade in product 4 (dairy products) is roughly 2.7 times the average while product 14 (vegetable plaiting materials) represents a negligible fraction of the average.

Note however that the interpretation of the fixed-effects must be taken with extra care. In particular, here the fixed-effects can be interpreted only because they are perfectly balanced.

2 Additional features

Now we present some other features of the package. First the possibility to add variables with varying slopes. Second how to combine several fixed-effects. Third, in the case of difference-in-difference analysis, the estimation and graph of the yearly average treatment effects. Fourth the lag.formula utility to lag variables easily. Fifth the possibility for non-linear in parameter estimation. Finally the use of parallelism to accelerate the estimation.

2.1 Varying slopes

You can introduce variables with varying slopes directly in the fixed-effects part of the formula using square brackets. Let’s go through a simple example using iris data:

base_vs = iris
names(base_vs) = c(paste0("x", 1:4), "species")

We want to estimate x1 as a function of x2 and the variable x3 with slopes varying according to species. We also want the species fixed-effect. We just have to do:

est_vs = feols(x1 ~ x2 | species[x3], base_vs)
est_vs
#> OLS estimation, Dep. Var.: x1
#> Observations: 150 
#> Fixed-effects: species: 3
#> Varying slopes: x3 (species: 3)
#> Standard-errors: Clustered (species) 
#>    Estimate Std. Error z value Pr(>|z|)    
#> x2 0.450006   0.157823  2.8513 0.004354 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -31.59   Adj. R2: 0.86351 
#>                        R2-Within: 0.17894

If you want to see the slopes for x3, just use the function fixef:

summary(fixef(est_vs))
#> Fixed-effects/Slope coefficients
#>                                species x3 (slopes: species)
#> Number of fixed-effects/slopes       3                    3
#> Number of references                 0                    0
#> Mean                               1.7                0.639
#> Variance                          1.74                0.069
#> 
#> COEFFICIENTS:
#>   species: versicolor virginica setosa
#>                 1.879    0.3036  2.927
#> -----
#>   x3 (slopes: species): versicolor virginica setosa
#>                             0.6599    0.8909 0.3667

2.2 Combining several fixed-effects

Let’s use the data we created in the previous section, and add a new variable:

# we create another "fixed-effect"
base_vs$fe = rep(1:5, 30)
head(base_vs)
#>    x1  x2  x3  x4 species fe
#> 1 5.1 3.5 1.4 0.2  setosa  1
#> 2 4.9 3.0 1.4 0.2  setosa  2
#> 3 4.7 3.2 1.3 0.2  setosa  3
#> 4 4.6 3.1 1.5 0.2  setosa  4
#> 5 5.0 3.6 1.4 0.2  setosa  5
#> 6 5.4 3.9 1.7 0.4  setosa  1

Say we want to “combine” the variable species with the variable fe and create a brand new fixed-effect variable. We can do it simply using ^:

est_comb = feols(x1 ~ x2 | species^fe, base_vs)
est_comb
#> OLS estimation, Dep. Var.: x1
#> Observations: 150 
#> Fixed-effects: species^fe: 15
#> Standard-errors: Clustered (species^fe) 
#>    Estimate Std. Error z value Pr(>|z|)    
#> x2 0.782815   0.125551   6.235 4.52e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -77.92   Adj. R2: 0.72986 
#>                        R2-Within: 0.28023

The function ^ does the same as paste0(species, "_", fe) but is more convenient (and faster for large data sets). You can still extract the fixed-effects the same way:

fixef(est_comb)[[1]]
#>     setosa_1     setosa_2     setosa_3     setosa_4     setosa_5 
#>     2.443630     2.384084     2.164943     2.296256     2.323630 
#> versicolor_1 versicolor_2 versicolor_3 versicolor_4 versicolor_5 
#>     3.713320     3.800694     4.003367     3.745539     3.575086 
#>  virginica_1  virginica_2  virginica_3  virginica_4  virginica_5 
#>     4.513272     3.986351     4.423725     4.216804     4.159382

2.3 Yearly treatment effect

In some difference-in-difference analyses, it is often useful not only to have the total treatment effect but to trace the evolution of the treatment. Package fixest offers a simple tool to do just that. Let’s take an example:

# Sample data illustrating the DiD
data(base_did)
head(base_did)
#>             y         x1 id period post treat
#> 1  2.87530627  0.5365377  1      1    0     1
#> 2  1.86065272 -3.0431894  1      2    0     1
#> 3  0.09416524  5.5768439  1      3    0     1
#> 4  3.78147485 -2.8300587  1      4    0     1
#> 5 -2.55819959 -5.0443544  1      5    0     1
#> 6  1.72873240 -0.6363849  1      6    1     1
# Estimation of yearly effect (they are automatically added)
# We also add individual/time fixed-effects:
est_did = did_estimate_yearly_effects(y ~ x1 | id + period, base_did,
                                      treat_time = ~treat+period, reference = 5)
est_did
#> OLS estimation, Dep. Var.: y
#> Observations: 1,080 
#> Fixed-effects: id: 108,  period: 10
#> Standard-errors: Clustered (id) 
#>           Estimate Std. Error   z value  Pr(>|z|)    
#> x1        0.973490   0.048174 20.208000 < 2.2e-16 ***
#> treat_1  -1.403000   1.170900 -1.198200   0.23083    
#> treat_2  -1.247500   1.152900 -1.082100  0.279216    
#> treat_3  -0.273206   1.167400 -0.234025  0.814966    
#> treat_4  -1.795700   1.147400 -1.565000  0.117583    
#> treat_6   0.784452   1.084600  0.723274  0.469512    
#> treat_7   3.598900   1.161800  3.097800   0.00195 ** 
#> treat_8   3.811800   1.315700  2.897200  0.003765 ** 
#> treat_9   4.731400   1.157100  4.089200   4.3e-05 ***
#> treat_10  6.606200   1.181700  5.590300  2.27e-08 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -2,984.58   Adj. R2: NA 
#>                           R2-Within: NA

In the example above, we must provide the treatment and time identifiers in argument treat_time, we also must provide a reference period (in this case 5). A treatment variable is added for each period (but the reference), and then a regular OLS estimation is performed with feols. You can change the estimation method with the argument estfun.

Now to display the yearly treatment effects on a graph:

did_plot_yearly_effects(est_did)

2.4 Lagging variables

To lag variables in a panel setting, a simple and fast method has been implemented: lag.formula. Let’s give an example with the previous data set:

base_lag = base_did
# we create a lagged value of the variable x1
base_lag$x1.l1 = lag(x1~id+period, 1, base_lag)
head(base_lag)
#>             y         x1 id period post treat      x1.l1
#> 1  2.87530627  0.5365377  1      1    0     1         NA
#> 2  1.86065272 -3.0431894  1      2    0     1  0.5365377
#> 3  0.09416524  5.5768439  1      3    0     1 -3.0431894
#> 4  3.78147485 -2.8300587  1      4    0     1  5.5768439
#> 5 -2.55819959 -5.0443544  1      5    0     1 -2.8300587
#> 6  1.72873240 -0.6363849  1      6    1     1 -5.0443544

The first two arguments are mandatory. The formula informs on the variable to be lagged (on the left hand side), and the two panel identifiers. Note that the time index must appear second. The second argument tells how much lags we want. Using negative values gives leads. Finally the last argument informs on where to find the variables.

In case you use the popular package data.table, you can create lagged variables very simply:

library(data.table)
base_lag_dt = as.data.table(base_did)
# we create a lagged value of the variable x1
base_lag_dt[, x1.l1 := lag(x1~id+period, 1)]

2.5 Non-linear in parameters example

The function feNmlm is similar to femlm but allows to have non-linear in parameters right-hand-sides (RHS). First an example without fixed-effects, the one with fixed-effects is given later. Let’s say we want to estimate the following relation with a Poisson model:

\(E\left(z_i\right) = a\times x_i + b\times y_i\).

In fact, this type of model is non-linear in the context of a Poisson model because the sum is embedded within the log:

\(E\left(z_i\right) = \exp\left(\log\left(a\times x_i + b\times y_i\right)\right)\).

So let’s estimate such a relation. (Note that you can estimate this relation with GLM and identity link, but I carry on for the example.) First we generate the data:

# Generating data:
n = 1000
# x and y: two positive random variables
x = rnorm(n, 1, 5)**2
y = rnorm(n, -1, 5)**2
# E(z) = 2*x + 3*y and some noise
z = rpois(n, 2*x + 3*y) + rpois(n, 1)
base = data.frame(x, y, z)

To estimate the non-linear relationship, we need to use the argument NL.fml where we put the non-linear part. We also have to provide starting values with the argument NL.start. Finally, to ensure the RHS can be evaluated in any situation, we add lower bounds for the parameters with the argument lower.

result_NL = feNmlm(z~0, base, NL.fml = ~ log(a*x + b*y), NL.start = list(a=1, b=1), lower = list(a=0, b=0))

Note that the arguments NL.start and lower are named lists. Setting lower = list(a=0, b=0) means that the optimization algorithm will never explore parameters for \(a\) and \(b\) that are lower than 0. The results obtained can be interpreted similarly to results with linear RHS. We can see them with a print:

print(result_NL)
#> Non-linear ML estimation, family = Poisson, Dep. Var.: z
#> Observations: 1,000 
#> Standard-errors: Standard 
#>   Estimate Std. Error z value  Pr(>|z|)    
#> a   2.0319   0.011298  179.84 < 2.2e-16 ***
#> b   3.0157   0.012959  232.72 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -3,626.84   Adj. Pseudo-R2: 0.9367 
#>            BIC:  7,281.31     Squared Cor.: 0.99159

We can see that we obtain coefficients close to the generating values.

2.5.1 Adding fixed-effects to non-linear in parameters models

Adding fixed-effects is identical to the linear case. The user must only be well aware of the functional form. Indeed, the fixed-effects must enter the estimation linearly. This means that the previous equation with one set of fixed-effects writes:

\(E\left(z_i\right) = \gamma_{id_i} \left( a\times x_i + b\times y_i \right)\),

where \(id_i\) is the class of observation \(i\) and \(\gamma\) is the vector of fixed-effects. Here the fixed-effects are in fact linear because in the context of the Poisson model we estimate:

\(E\left(z_i\right) = \exp\left(\gamma_{id_i}+\log\left(a\times x_i + b\times y_i\right)\right)\).

Further, remark that there exists an infinity of values of \(\gamma^{\prime}\), \(a^{\prime}\) and \(b^{\prime}\) such that:

\(\gamma_{k} \left( a\times x_i + b\times y_i \right) = \gamma_{k}^{\prime} \left(a^{\prime}\times x_i + b^{\prime}\times y_i \right),\forall i,k\).

An example is \(\gamma^{\prime}_{k} = 2\times \gamma_k\), \(a^{\prime} = a/2\) and \(b^{\prime} = b/2\). Thus estimating this relation directly will lead to a problem to uniquely identify the coefficients. To circumvent this problem, we just have to fix one of the coefficient, this will ensure that we uniquely identify them.

Let’s generate this relation:

# the class of each observation
id = sample(20, n, replace = TRUE)
base$id = id
# the vector of fixed-effects
gamma = rnorm(20)**2
# the new vector z_bis
z_bis = rpois(n, gamma[id] * (2*x + 3*y)) + rpois(n, 1)
base$z_bis = z_bis

Now we estimate it with the fixed-effects while fixing one of the coefficients (we fix \(a\) to its true value but it could be any value):

# we add the fixed-effect in the formula
result_NL_fe = feNmlm(z_bis~0|id, base, NL.fml = ~ log(2*x + b*y), NL.start = list(b=1), lower = list(b=0))
# The coef should be around 3
coef(result_NL_fe)
#>       b 
#> 3.00565
# the gamma and the exponential of the fixed-effects should be similar
rbind(gamma, exp(fixef(result_NL_fe)$id))
#>                2         17           20         4           7          5
#> gamma 2.79492534 0.08872863 0.0008413015 0.2259703 0.002133664 1.00874978
#>       0.09334733 0.75513419 1.2657439756 0.2345213 0.947952628 0.01072864
#>               6         12        13           3        16         14
#> gamma 0.9428807 0.05975179 0.7352481 1.746031537 0.1058210 0.02670944
#>       1.0118936 0.03349028 1.4687853 0.008358586 0.3603113 1.61192121
#>                8         9         18        11        10         1
#> gamma 1.46179126 1.6094976 0.03813667 0.3518205 0.7450913 0.7069153
#>       0.06603475 0.7411054 0.71242931 0.1140938 1.7626529 2.8267866
#>               19         15
#> gamma 0.02223765 1.26242843
#>       0.03169395 0.04452463

As we can see, we obtain the “right” estimates.

2.6 Multi-threading

The package fixest integrates multi-platform parallelism to hasten the estimation process. By default it makes use of all the available threads minus 2. To change the number of threads used, just use the argument nthreads:

# Sample of results:
# 1 nthreads: 3.13s
system.time(fenegbin(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, nthreads = 1))
# 2 nthreads: 1.82s
system.time(fenegbin(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, nthreads = 2))
# 4 nthreads: 1.17s
system.time(fenegbin(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, nthreads = 4))

As you can see, the efficiency of increasing the number of threads is not 1 to 1. Two threads do not divide the computing time by 2, nor four threads by 4. However it still reduces significantly the computing time, which might be valuable for large sample estimations.

You can permanently set the number of threads used by fixest using setFixest_nthreads(nthreads).

2.7 Collineartiy

The user ought to estimate the coefficient of variables that are not collinear: neither among each other, neither with the fixed-effects. Estimation with collinear variables leads to a non invertible Hessian (leading to the absence of Variance-Covariance matrix for the coefficients). In such cases, the estimating functions will raise a warning and suggest to use the function collinearity to spot the problem.

Let’s take an example in which we want to make a fixed-effects estimation with a variable which is constant. Of course it makes no sense (this variable is perfectly collinear with the fixed-effects), so a warning will be raised suggesting to use the function collinearity to figure out what is wrong.

base_coll = trade
base_coll$constant_variable = 1
res <- femlm(Euros ~ log(dist_km) + constant_variable|Origin+Destination+Product+Year, base_coll)
#> Warning: [femlm]: The optimization algorithm did not converge, the
#> results are not reliable. The information matrix is singular: presence of
#> collinearity. Use function collinearity() to pinpoint the problems.
collinearity(res)
#> [1] "Variable constant_variable is constant, thus collinear with the fixed-effects."

As we can see, the function collinearity spots the collinear variables and name them. Even in elaborate cases of collinearity, the algorithm tries to find out the culprit and informs the user accordingly.


  1. Since the \(\gamma\) are parameters, I omit to put them in logarithmic form.