This function estimates maximum likelihood models with any number of fixed-effects.

## Usage

```
femlm(
fml,
data,
family = c("poisson", "negbin", "logit", "gaussian"),
vcov,
start = 0,
fixef,
fixef.rm = "perfect",
offset,
subset,
split,
fsplit,
split.keep,
split.drop,
cluster,
se,
ssc,
panel.id,
fixef.tol = 1e-05,
fixef.iter = 10000,
nthreads = getFixest_nthreads(),
lean = FALSE,
verbose = 0,
warn = TRUE,
notes = getFixest_notes(),
theta.init,
combine.quick,
mem.clean = FALSE,
only.env = FALSE,
only.coef = FALSE,
data.save = FALSE,
env,
...
)
fenegbin(
fml,
data,
vcov,
theta.init,
start = 0,
fixef,
fixef.rm = "perfect",
offset,
subset,
split,
fsplit,
split.keep,
split.drop,
cluster,
se,
ssc,
panel.id,
fixef.tol = 1e-05,
fixef.iter = 10000,
nthreads = getFixest_nthreads(),
lean = FALSE,
verbose = 0,
warn = TRUE,
notes = getFixest_notes(),
combine.quick,
mem.clean = FALSE,
only.env = FALSE,
only.coef = FALSE,
data.save = FALSE,
env,
...
)
```

## Arguments

- fml
A formula representing the relation to be estimated. For example:

`fml = z~x+y`

. To include fixed-effects, insert them in this formula using a pipe: e.g.`fml = z~x+y|fixef_1+fixef_2`

. Multiple estimations can be performed at once: for multiple dep. vars, wrap them in`c()`

: ex`c(y1, y2)`

. For multiple indep. vars, use the stepwise functions: ex`x1 + csw(x2, x3)`

. The formula`fml = c(y1, y2) ~ x1 + cw0(x2, x3)`

leads to 6 estimation, see details. Square brackets starting with a dot can be used to call global variables:`y.[i] ~ x.[1:2]`

will lead to`y3 ~ x1 + x2`

if`i`

is equal to 3 in the current environment (see details in`xpd`

).- data
A data.frame containing the necessary variables to run the model. The variables of the non-linear right hand side of the formula are identified with this

`data.frame`

names. Can also be a matrix.- family
Character scalar. It should provide the family. The possible values are "poisson" (Poisson model with log-link, the default), "negbin" (Negative Binomial model with log-link), "logit" (LOGIT model with log-link), "gaussian" (Gaussian model).

- vcov
Versatile argument to specify the VCOV. In general, it is either a character scalar equal to a VCOV type, either a formula of the form:

`vcov_type ~ variables`

. The VCOV types implemented are: "iid", "hetero" (or "HC1"), "cluster", "twoway", "NW" (or "newey_west"), "DK" (or "driscoll_kraay"), and "conley". It also accepts object from`vcov_cluster`

,`vcov_NW`

,`NW`

,`vcov_DK`

,`DK`

,`vcov_conley`

and`conley`

. It also accepts covariance matrices computed externally. Finally it accepts functions to compute the covariances. See the`vcov`

documentation in the vignette.- start
Starting values for the coefficients. Can be: i) a numeric of length 1 (e.g.

`start = 0`

, the default), ii) a numeric vector of the exact same length as the number of variables, or iii) a named vector of any length (the names will be used to initialize the appropriate coefficients).- fixef
Character vector. The names of variables to be used as fixed-effects. These variables should contain the identifier of each observation (e.g., think of it as a panel identifier). Note that the recommended way to include fixed-effects is to insert them directly in the formula.

- fixef.rm
Can be equal to "perfect" (default), "singleton", "both" or "none". Controls which observations are to be removed. If "perfect", then observations having a fixed-effect with perfect fit (e.g. only 0 outcomes in Poisson estimations) will be removed. If "singleton", all observations for which a fixed-effect appears only once will be removed. Note, importantly, that singletons are removed in just one pass, there is no recursivity implemented. The meaning of "both" and "none" is direct.

- offset
A formula or a numeric vector. An offset can be added to the estimation. If equal to a formula, it should be of the form (for example)

`~0.5*x**2`

. This offset is linearly added to the elements of the main formula 'fml'.- subset
A vector (logical or numeric) or a one-sided formula. If provided, then the estimation will be performed only on the observations defined by this argument.

- split
A one sided formula representing a variable (eg

`split = ~var`

) or a vector. If provided, the sample is split according to the variable and one estimation is performed for each value of that variable. If you also want to include the estimation for the full sample, use the argument`fsplit`

instead. You can use the special operators`%keep%`

and`%drop%`

to select only a subset of values for which to split the sample. E.g.`split = ~var %keep% c("v1", "v2")`

will split the sample only according to the values`v1`

and`v2`

of the variable`var`

; it is equivalent to supplying the argument`split.keep = c("v1", "v2")`

. By default there is partial matching on each value, you can trigger a regular expression evaluation by adding a`'@'`

first, as in:`~var %drop% "@^v[12]"`

which will drop values starting with`"v1"`

or`"v2"`

(of course you need to know regexes!).- fsplit
A one sided formula representing a variable (eg

`fsplit = ~var`

) or a vector. If provided, the sample is split according to the variable and one estimation is performed for each value of that variable. This argument is the same as`split`

but also includes the full sample as the first estimation. You can use the special operators`%keep%`

and`%drop%`

to select only a subset of values for which to split the sample. E.g.`fsplit = ~var %keep% c("v1", "v2")`

will split the sample only according to the values`v1`

and`v2`

of the variable`var`

; it is equivalent to supplying the argument`split.keep = c("v1", "v2")`

. By default there is partial matching on each value, you can trigger a regular expression evaluation by adding an`'@'`

first, as in:`~var %drop% "@^v[12]"`

which will drop values starting with`"v1"`

or`"v2"`

(of course you need to know regexes!).- split.keep
A character vector. Only used when

`split`

, or`fsplit`

, is supplied. If provided, then the sample will be split only on the values of`split.keep`

. The values in`split.keep`

will be partially matched to the values of`split`

. To enable regular expressions, you need to add an`'@'`

first. For example`split.keep = c("v1", "@other|var")`

will keep only the value in`split`

partially matched by`"v1"`

or the values containing`"other"`

or`"var"`

.- split.drop
A character vector. Only used when

`split`

, or`fsplit`

, is supplied. If provided, then the sample will be split only on the values that are not in`split.drop`

. The values in`split.drop`

will be partially matched to the values of`split`

. To enable regular expressions, you need to add an`'@'`

first. For example`split.drop = c("v1", "@other|var")`

will drop only the value in`split`

partially matched by`"v1"`

or the values containing`"other"`

or`"var"`

.- cluster
Tells how to cluster the standard-errors (if clustering is requested). Can be either a list of vectors, a character vector of variable names, a formula or an integer vector. Assume we want to perform 2-way clustering over

`var1`

and`var2`

contained in the data.frame`base`

used for the estimation. All the following`cluster`

arguments are valid and do the same thing:`cluster = base[, c("var1", "var2")]`

,`cluster = c("var1", "var2")`

,`cluster = ~var1+var2`

. If the two variables were used as fixed-effects in the estimation, you can leave it blank with`vcov = "twoway"`

(assuming`var1`

[resp.`var2`

] was the 1st [resp. 2nd] fixed-effect). You can interact two variables using`^`

with the following syntax:`cluster = ~var1^var2`

or`cluster = "var1^var2"`

.- se
Character scalar. Which kind of standard error should be computed: “standard”, “hetero”, “cluster”, “twoway”, “threeway” or “fourway”? By default if there are clusters in the estimation:

`se = "cluster"`

, otherwise`se = "iid"`

. Note that this argument is deprecated, you should use`vcov`

instead.- ssc
An object of class

`ssc.type`

obtained with the function`ssc`

. Represents how the degree of freedom correction should be done.You must use the function`ssc`

for this argument. The arguments and defaults of the function`ssc`

are:`adj = TRUE`

,`fixef.K="nested"`

,`cluster.adj = TRUE`

,`cluster.df = "min"`

,`t.df = "min"`

,`fixef.force_exact=FALSE)`

. See the help of the function`ssc`

for details.- panel.id
The panel identifiers. Can either be: i) a one sided formula (e.g.

`panel.id = ~id+time`

), ii) a character vector of length 2 (e.g.`panel.id=c('id', 'time')`

, or iii) a character scalar of two variables separated by a comma (e.g.`panel.id='id,time'`

). Note that you can combine variables with`^`

only inside formulas (see the dedicated section in`feols`

).- fixef.tol
Precision used to obtain the fixed-effects. Defaults to

`1e-5`

. It corresponds to the maximum absolute difference allowed between two coefficients of successive iterations. Argument`fixef.tol`

cannot be lower than`10000*.Machine$double.eps`

. Note that this parameter is dynamically controlled by the algorithm.- fixef.iter
Maximum number of iterations in fixed-effects algorithm (only in use for 2+ fixed-effects). Default is 10000.

- nthreads
The number of threads. Can be: a) an integer lower than, or equal to, the maximum number of threads; b) 0: meaning all available threads will be used; c) a number strictly between 0 and 1 which represents the fraction of all threads to use. The default is to use 50% of all threads. You can set permanently the number of threads used within this package using the function

`setFixest_nthreads`

.- lean
Logical, default is

`FALSE`

. If`TRUE`

then all large objects are removed from the returned result: this will save memory but will block the possibility to use many methods. It is recommended to use the arguments`se`

or`cluster`

to obtain the appropriate standard-errors at estimation time, since obtaining different SEs won't be possible afterwards.- verbose
Integer, default is 0. It represents the level of information that should be reported during the optimisation process. If

`verbose=0`

: nothing is reported. If`verbose=1`

: the value of the coefficients and the likelihood are reported. If`verbose=2`

:`1`

+ information on the computing time of the null model, the fixed-effects coefficients and the hessian are reported.- warn
Logical, default is

`TRUE`

. Whether warnings should be displayed (concerns warnings relating to convergence state).- notes
Logical. By default, two notes are displayed: when NAs are removed (to show additional information) and when some observations are removed because of only 0 (or 0/1) outcomes in a fixed-effect setup (in Poisson/Neg. Bin./Logit models). To avoid displaying these messages, you can set

`notes = FALSE`

. You can remove these messages permanently by using`setFixest_notes(FALSE)`

.- theta.init
Positive numeric scalar. The starting value of the dispersion parameter if

`family="negbin"`

. By default, the algorithm uses as a starting value the theta obtained from the model with only the intercept.- combine.quick
Logical. When you combine different variables to transform them into a single fixed-effects you can do e.g.

`y ~ x | paste(var1, var2)`

. The algorithm provides a shorthand to do the same operation:`y ~ x | var1^var2`

. Because pasting variables is a costly operation, the internal algorithm may use a numerical trick to hasten the process. The cost of doing so is that you lose the labels. If you are interested in getting the value of the fixed-effects coefficients after the estimation, you should use`combine.quick = FALSE`

. By default it is equal to`FALSE`

if the number of observations is lower than 50,000, and to`TRUE`

otherwise.- mem.clean
Logical, default is

`FALSE`

. Only to be used if the data set is large compared to the available RAM. If`TRUE`

then intermediary objects are removed as much as possible and`gc`

is run before each substantial C++ section in the internal code to avoid memory issues.- only.env
(Advanced users.) Logical, default is

`FALSE`

. If`TRUE`

, then only the environment used to make the estimation is returned.- only.coef
Logical, default is

`FALSE`

. If`TRUE`

, then only the estimated coefficients are returned. Note that the length of the vector returned is always the length of the number of coefficients to be estimated: this means that the variables found to be collinear are returned with an NA value.- data.save
Logical scalar, default is

`FALSE`

. If`TRUE`

, the data used for the estimation is saved within the returned object. Hence later calls to predict(), vcov(), etc..., will be consistent even if the original data has been modified in the meantime. This is especially useful for estimations within loops, where the data changes at each iteration, such that postprocessing can be done outside the loop without issue.- env
(Advanced users.) A

`fixest`

environment created by a`fixest`

estimation with`only.env = TRUE`

. Default is missing. If provided, the data from this environment will be used to perform the estimation.- ...
Not currently used.

## Value

A `fixest`

object. Note that `fixest`

objects contain many elements and most of
them are for internal use, they are presented here only for information.
To access them, it is safer to use the user-level methods
(e.g. `vcov.fixest`

, `resid.fixest`

, etc) or functions (like for instance
`fitstat`

to access any fit statistic).

- nobs
The number of observations.

- fml
The linear formula of the call.

- call
The call of the function.

- method
The method used to estimate the model.

- family
The family used to estimate the model.

- data
The original data set used when calling the function. Only available when the estimation was called with

`data.save = TRUE`

- fml_all
A list containing different parts of the formula. Always contain the linear formula. Then, if relevant:

`fixef`

: the fixed-effects;`NL`

: the non linear part of the formula.- nparams
The number of parameters of the model.

- fixef_vars
The names of each fixed-effect dimension.

- fixef_id
The list (of length the number of fixed-effects) of the fixed-effects identifiers for each observation.

- fixef_sizes
The size of each fixed-effect (i.e. the number of unique identifier for each fixed-effect dimension).

- convStatus
Logical, convergence status.

- message
The convergence message from the optimization procedures.

- obs_selection
(When relevant.) List containing vectors of integers. It represents the sequential selection of observation vis a vis the original data set.

- fixef_removed
(When relevant.) In the case there were fixed-effects and some observations were removed because of only 0/1 outcome within a fixed-effect, it gives the list (for each fixed-effect dimension) of the fixed-effect identifiers that were removed.

- coefficients
The named vector of estimated coefficients.

- coeftable
The table of the coefficients with their standard errors, z-values and p-values.

- loglik
The log-likelihood.

- iterations
Number of iterations of the algorithm.

- ll_null
Log-likelihood of the null model (i.e. with the intercept only).

- ll_fe_only
Log-likelihood of the model with only the fixed-effects.

- ssr_null
Sum of the squared residuals of the null model (containing only with the intercept).

- pseudo_r2
The adjusted pseudo R2.

- fitted.values
The fitted values are the expected value of the dependent variable for the fitted model: that is \(E(Y|X)\).

- residuals
The residuals (y minus the fitted values).

- sq.cor
Squared correlation between the dependent variable and the expected predictor (i.e. fitted.values) obtained by the estimation.

- hessian
The Hessian of the parameters.

- cov.iid
The variance-covariance matrix of the parameters.

- se
The standard-error of the parameters.

- scores
The matrix of the scores (first derivative for each observation).

- residuals
The difference between the dependent variable and the expected predictor.

- sumFE
The sum of the fixed-effects coefficients for each observation.

- offset
(When relevant.) The offset formula.

- weights
(When relevant.) The weights formula.

## Details

Note that the functions `feglm`

and `femlm`

provide the same results when using
the same families but differ in that the latter is a direct maximum likelihood
optimization (so the two can really have different convergence rates).

## Combining the fixed-effects

You can combine two variables to make it a new fixed-effect using `^`

.
The syntax is as follows: `fe_1^fe_2`

. Here you created a new variable which is the combination
of the two variables fe_1 and fe_2. This is identical to doing `paste0(fe_1, "_", fe_2)`

but more convenient.

Note that pasting is a costly operation, especially for large data sets.
Thus, the internal algorithm uses a numerical trick which is fast, but the drawback is
that the identity of each observation is lost (i.e. they are now equal to a meaningless
number instead of being equal to `paste0(fe_1, "_", fe_2)`

). These “identities”
are useful only if you're interested in the value of the fixed-effects (that you can
extract with `fixef.fixest`

). If you're only interested in coefficients of the variables,
it doesn't matter. Anyway, you can use `combine.quick = FALSE`

to tell the internal
algorithm to use `paste`

instead of the numerical trick. By default, the numerical
trick is performed only for large data sets.

## Lagging variables

To use leads/lags of variables in the estimation, you can: i) either provide the argument
`panel.id`

, ii) either set your data set as a panel with the function
`panel`

, `f`

and `d`

.

You can provide several leads/lags/differences at once: e.g. if your formula is equal to
`f(y) ~ l(x, -1:1)`

, it means that the dependent variable is equal to the lead of `y`

,
and you will have as explanatory variables the lead of `x1`

, `x1`

and the lag of `x1`

.
See the examples in function `l`

for more details.

## Interactions

You can interact a numeric variable with a "factor-like" variable by using
`i(factor_var, continuous_var, ref)`

, where `continuous_var`

will be interacted with
each value of `factor_var`

and the argument `ref`

is a value of `factor_var`

taken as a reference (optional).

Using this specific way to create interactions leads to a different display of the
interacted values in `etable`

. See examples.

It is important to note that *if you do not care about the standard-errors of
the interactions*, then you can add interactions in the fixed-effects part of the formula,
it will be incomparably faster (using the syntax `factor_var[continuous_var]`

, as explained
in the section “Varying slopes”).

The function `i`

has in fact more arguments, please see details in its associated help page.

## On standard-errors

Standard-errors can be computed in different ways, you can use the arguments `se`

and `ssc`

in `summary.fixest`

to define how to compute them. By default, in the presence
of fixed-effects, standard-errors are automatically clustered.

The following vignette: On standard-errors describes in details how the standard-errors are computed in
`fixest`

and how you can replicate standard-errors from other software.

You can use the functions `setFixest_vcov`

and `setFixest_ssc`

to
permanently set the way the standard-errors are computed.

## Multiple estimations

Multiple estimations can be performed at once, they just have to be specified in the formula.
Multiple estimations yield a `fixest_multi`

object which is ‘kind of’ a list of
all the results but includes specific methods to access the results in a handy way.
Please have a look at the dedicated vignette:
Multiple estimations.

To include multiple dependent variables, wrap them in `c()`

(`list()`

also works).
For instance `fml = c(y1, y2) ~ x1`

would estimate the model `fml = y1 ~ x1`

and
then the model `fml = y2 ~ x1`

.

To include multiple independent variables, you need to use the stepwise functions.
There are 4 stepwise functions: `sw`

, `sw0`

, `csw`

, `csw0`

, and `mvsw`

. Of course `sw`

stands for stepwise, and `csw`

for cumulative stepwise. Finally `mvsw`

is a bit special,
it stands for multiverse stepwise. Let's explain that.
Assume you have the following formula: `fml = y ~ x1 + sw(x2, x3)`

.
The stepwise function `sw`

will estimate the following two models: `y ~ x1 + x2`

and
`y ~ x1 + x3`

. That is, each element in `sw()`

is sequentially, and separately,
added to the formula. Would have you used `sw0`

in lieu of `sw`

, then the model
`y ~ x1`

would also have been estimated. The `0`

in the name means that the model
without any stepwise element also needs to be estimated.
The prefix `c`

means cumulative: each stepwise element is added to the next. That is,
`fml = y ~ x1 + csw(x2, x3)`

would lead to the following models `y ~ x1 + x2`

and
`y ~ x1 + x2 + x3`

. The `0`

has the same meaning and would also lead to the model without
the stepwise elements to be estimated: in other words, `fml = y ~ x1 + csw0(x2, x3)`

leads to the following three models: `y ~ x1`

, `y ~ x1 + x2`

and `y ~ x1 + x2 + x3`

.
Finally `mvsw`

will add, in a stepwise fashion all possible combinations of the variables
in its arguments. For example `mvsw(x1, x2, x3)`

is equivalent to
`sw0(x1, x2, x3, x1 + x2, x1 + x3, x2 + x3, x1 + x2 + x3)`

. The number of models
to estimate grows at a factorial rate: so be cautious!

Multiple independent variables can be combined with multiple dependent variables, as in
`fml = c(y1, y2) ~ cw(x1, x2, x3)`

which would lead to 6 estimations. Multiple
estimations can also be combined to split samples (with the arguments `split`

, `fsplit`

).

You can also add fixed-effects in a stepwise fashion. Note that you cannot perform
stepwise estimations on the IV part of the formula (`feols`

only).

If NAs are present in the sample, to avoid too many messages, only NA removal concerning the variables common to all estimations is reported.

A note on performance. The feature of multiple estimations has been highly optimized for
`feols`

, in particular in the presence of fixed-effects. It is faster to estimate
multiple models using the formula rather than with a loop. For non-`feols`

models using
the formula is roughly similar to using a loop performance-wise.

## Argument sliding

When the data set has been set up globally using
`setFixest_estimation`

`(data = data_set)`

, the argument `vcov`

can be used implicitly.
This means that calls such as `feols(y ~ x, "HC1")`

, or `feols(y ~ x, ~id)`

, are valid:
i) the data is automatically deduced from the global settings, and ii) the `vcov`

is deduced to be the second argument.

## Piping

Although the argument 'data' is placed in second position, the data can be piped to the
estimation functions. For example, with R >= 4.1, `mtcars |> feols(mpg ~ cyl)`

works as
`feols(mpg ~ cyl, mtcars)`

.

## Tricks to estimate multiple LHS

To use multiple dependent variables in `fixest`

estimations, you need to include them
in a vector: like in `c(y1, y2, y3)`

.

First, if names are stored in a vector, they can readily be inserted in a formula to
perform multiple estimations using the dot square bracket operator. For instance if
`my_lhs = c("y1", "y2")`

, calling `fixest`

with, say `feols(.[my_lhs] ~ x1, etc)`

is
equivalent to using `feols(c(y1, y2) ~ x1, etc)`

. Beware that this is a special feature
unique to the *left-hand-side* of `fixest`

estimations (the default behavior of the DSB
operator is to aggregate with sums, see `xpd`

).

Second, you can use a regular expression to grep the left-hand-sides on the fly. When the
`..("regex")`

feature is used naked on the LHS, the variables grepped are inserted into
`c()`

. For example `..("Pe") ~ Sepal.Length, iris`

is equivalent to
`c(Petal.Length, Petal.Width) ~ Sepal.Length, iris`

. Beware that this is a
special feature unique to the *left-hand-side* of `fixest`

estimations
(the default behavior of `..("regex")`

is to aggregate with sums, see `xpd`

).

## Dot square bracket operator in formulas

In a formula, the dot square bracket (DSB) operator can: i) create manifold variables at once, or ii) capture values from the current environment and put them verbatim in the formula.

Say you want to include the variables `x1`

to `x3`

in your formula. You can use
`xpd(y ~ x.[1:3])`

and you'll get `y ~ x1 + x2 + x3`

.

To summon values from the environment, simply put the variable in square brackets. For example:
`for(i in 1:3) xpd(y.[i] ~ x)`

will create the formulas `y1 ~ x`

to `y3 ~ x`

depending on the
value of `i`

.

You can include a full variable from the environment in the same way:
`for(y in c("a", "b")) xpd(.[y] ~ x)`

will create the two formulas `a ~ x`

and `b ~ x`

.

The DSB can even be used within variable names, but then the variable must be nested in
character form. For example `y ~ .["x.[1:2]_sq"]`

will create `y ~ x1_sq + x2_sq`

. Using the
character form is important to avoid a formula parsing error. Double quotes must be used. Note
that the character string that is nested will be parsed with the function `dsb`

, and thus it
will return a vector.

By default, the DSB operator expands vectors into sums. You can add a comma, like in `.[, x]`

,
to expand with commas--the content can then be used within functions. For instance:
`c(x.[, 1:2])`

will create `c(x1, x2)`

(and *not* `c(x1 + x2)`

).

In all `fixest`

estimations, this special parsing is enabled, so you don't need to use `xpd`

.

One-sided formulas can be expanded with the DSB operator: let `x = ~sepal + petal`

, then
`xpd(y ~ .[x])`

leads to `color ~ sepal + petal`

.

You can even use multiple square brackets within a single variable, but then the use of nesting
is required. For example, the following `xpd(y ~ .[".[letters[1:2]]_.[1:2]"])`

will create
`y ~ a_1 + b_2`

. Remember that the nested character string is parsed with `dsb`

,
which explains this behavior.

When the element to be expanded i) is equal to the empty string or, ii) is of length 0, it is
replaced with a neutral element, namely `1`

. For example, `x = "" ; xpd(y ~ .[x])`

leads to
`y ~ 1`

.

## References

Berge, Laurent, 2018, "Efficient estimation of maximum likelihood models with multiple fixed-effects: the R package FENmlm." CREA Discussion Papers, 13 ().

For models with multiple fixed-effects:

Gaure, Simen, 2013, "OLS with multiple high dimensional category variables", Computational Statistics & Data Analysis 66 pp. 8--18

On the unconditionnal Negative Binomial model:

Allison, Paul D and Waterman, Richard P, 2002, "Fixed-Effects Negative Binomial Regression Models", Sociological Methodology 32(1) pp. 247--265

## See also

See also `summary.fixest`

to see the results with the appropriate standard-errors,
`fixef.fixest`

to extract the fixed-effects coefficients, and the function
`etable`

to visualize the results of multiple estimations.
And other estimation methods: `feols`

, `feglm`

, `fepois`

, `feNmlm`

.

## Examples

```
# Load trade data
data(trade)
# We estimate the effect of distance on trade => we account for 3 fixed-effects
# 1) Poisson estimation
est_pois = femlm(Euros ~ log(dist_km) | Origin + Destination + Product, trade)
# 2) Log-Log Gaussian estimation (with same FEs)
est_gaus = update(est_pois, log(Euros+1) ~ ., family = "gaussian")
# Comparison of the results using the function etable
etable(est_pois, est_gaus)
#> est_pois est_gaus
#> Dependent Var.: Euros log(Euros+1)
#>
#> log(dist_km) -1.528*** (0.1156) -2.168*** (0.1542)
#> Fixed-Effects: ------------------ ------------------
#> Origin Yes Yes
#> Destination Yes Yes
#> Product Yes Yes
#> _______________ __________________ __________________
#> Family Poisson Gaussian
#> S.E.: Clustered by: Origin by: Origin
#> Observations 38,325 38,325
#> Squared Cor. 0.60377 0.70434
#> Pseudo R2 0.76039 0.23563
#> BIC 1.43e+12 152,008.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Now using two way clustered standard-errors
etable(est_pois, est_gaus, se = "twoway")
#> est_pois est_gaus
#> Dependent Var.: Euros log(Euros+1)
#>
#> log(dist_km) -1.528*** (0.1307) -2.168*** (0.1712)
#> Fixed-Effects: ------------------ ------------------
#> Origin Yes Yes
#> Destination Yes Yes
#> Product Yes Yes
#> _______________ __________________ __________________
#> Family Poisson Gaussian
#> S.E.: Clustered by: Orig. & Dest. by: Orig. & Dest.
#> Observations 38,325 38,325
#> Squared Cor. 0.60377 0.70434
#> Pseudo R2 0.76039 0.23563
#> BIC 1.43e+12 152,008.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comparing different types of standard errors
sum_hetero = summary(est_pois, se = "hetero")
sum_oneway = summary(est_pois, se = "cluster")
sum_twoway = summary(est_pois, se = "twoway")
sum_threeway = summary(est_pois, se = "threeway")
etable(sum_hetero, sum_oneway, sum_twoway, sum_threeway)
#> sum_hetero sum_oneway sum_twoway
#> Dependent Var.: Euros Euros Euros
#>
#> log(dist_km) -1.528*** (0.0220) -1.528*** (0.1156) -1.528*** (0.1307)
#> Fixed-Effects: ------------------ ------------------ ------------------
#> Origin Yes Yes Yes
#> Destination Yes Yes Yes
#> Product Yes Yes Yes
#> _______________ __________________ __________________ __________________
#> S.E. type Heteroskedas.-rob. by: Origin by: Orig. & Dest.
#> Observations 38,325 38,325 38,325
#> Squared Cor. 0.60377 0.60377 0.60377
#> Pseudo R2 0.76039 0.76039 0.76039
#> BIC 1.43e+12 1.43e+12 1.43e+12
#>
#> sum_threeway
#> Dependent Var.: Euros
#>
#> log(dist_km) -1.528*** (0.1340)
#> Fixed-Effects: ------------------
#> Origin Yes
#> Destination Yes
#> Product Yes
#> _______________ __________________
#> S.E. type by: Orig. & Dest. & Prod.
#> Observations 38,325
#> Squared Cor. 0.60377
#> Pseudo R2 0.76039
#> BIC 1.43e+12
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Multiple estimations:
#
# 6 estimations
est_mult = femlm(c(Ozone, Solar.R) ~ Wind + Temp + csw0(Wind:Temp, Day), airquality)
# We can display the results for the first lhs:
etable(est_mult[lhs = 1])
#> est_mult[lhs = 1].1 est_mult[lhs = ..2 est_mult[lhs = ..3
#> Dependent Var.: Ozone Ozone Ozone
#>
#> Constant 0.5334** (0.1943) 0.8890* (0.4430) 0.8028. (0.4441)
#> Wind -0.0761*** (0.0052) -0.1151** (0.0439) -0.1152** (0.0440)
#> Temp 0.0483*** (0.0020) 0.0441*** (0.0051) 0.0444*** (0.0051)
#> Wind x Temp 0.0005 (0.0005) 0.0005 (0.0005)
#> Day 0.0051*** (0.0015)
#> _______________ ___________________ __________________ __________________
#> S.E. type IID IID IID
#> Observations 116 116 116
#> Squared Cor. 0.62615 0.63129 0.63379
#> Pseudo R2 0.53377 0.53402 0.53759
#> BIC 1,564.7 1,568.7 1,561.5
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# And now the second (access can be made by name)
etable(est_mult[lhs = "Solar.R"])
#> est_mult[lhs = ..1 est_mult[lhs = ..2 est_mult[lhs = "..3
#> Dependent Var.: Solar.R Solar.R Solar.R
#>
#> Constant 3.777*** (0.0714) 4.343*** (0.1689) 4.529*** (0.1708)
#> Wind 0.0119*** (0.0020) -0.0415** (0.0146) -0.0433** (0.0147)
#> Temp 0.0169*** (0.0008) 0.0098*** (0.0021) 0.0088*** (0.0021)
#> Wind x Temp 0.0007*** (0.0002) 0.0007*** (0.0002)
#> Day -0.0061*** (0.0007)
#> _______________ __________________ __________________ ___________________
#> S.E. type IID IID IID
#> Observations 146 146 146
#> Squared Cor. 0.08076 0.08166 0.08983
#> Pseudo R2 0.06045 0.06207 0.07100
#> BIC 8,192.6 8,183.5 8,110.7
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Now we focus on the two last right hand sides
# (note that .N can be used to specify the last item)
etable(est_mult[rhs = 2:.N])
#> est_mult[rhs = ..1 est_mult[rhs = ..2 est_mult[rhs = ..3
#> Dependent Var.: Ozone Solar.R Ozone
#>
#> Constant 0.8890* (0.4430) 4.343*** (0.1689) 0.8028. (0.4441)
#> Wind -0.1151** (0.0439) -0.0415** (0.0146) -0.1152** (0.0440)
#> Temp 0.0441*** (0.0051) 0.0098*** (0.0021) 0.0444*** (0.0051)
#> Wind x Temp 0.0005 (0.0005) 0.0007*** (0.0002) 0.0005 (0.0005)
#> Day 0.0051*** (0.0015)
#> _______________ __________________ __________________ __________________
#> S.E. type IID IID IID
#> Observations 116 146 116
#> Squared Cor. 0.63129 0.08166 0.63379
#> Pseudo R2 0.53402 0.06207 0.53759
#> BIC 1,568.7 8,183.5 1,561.5
#>
#> est_mult[rhs = 2..4
#> Dependent Var.: Solar.R
#>
#> Constant 4.529*** (0.1708)
#> Wind -0.0433** (0.0147)
#> Temp 0.0088*** (0.0021)
#> Wind x Temp 0.0007*** (0.0002)
#> Day -0.0061*** (0.0007)
#> _______________ ___________________
#> S.E. type IID
#> Observations 146
#> Squared Cor. 0.08983
#> Pseudo R2 0.07100
#> BIC 8,110.7
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Combining with split
est_split = fepois(c(Ozone, Solar.R) ~ sw(poly(Wind, 2), poly(Temp, 2)),
airquality, split = ~ Month)
# You can display everything at once with the print method
est_split
#> Standard-errors: IID
#>
#>
#> # SAMPLE: 5
#>
#>
#> ### Dep. var.: Ozone
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.30403 0.042415 77.89732 < 2.2e-16 ***
#> poly(Wind, 2)1 -6.41466 0.617300 -10.39149 < 2.2e-16 ***
#> poly(Wind, 2)2 3.07584 0.570194 5.39437 6.8765e-08 ***
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 4.05816 0.109394 37.09674 < 2.2e-16 ***
#> poly(Temp, 2)1 11.28467 1.448123 7.79261 6.5637e-15 ***
#> poly(Temp, 2)2 2.50568 1.171805 2.13830 3.2492e-02 *
#>
#> ### Dep. var.: Solar.R
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.207784 0.018194 286.23787 < 2.2e-16 ***
#> poly(Wind, 2)1 -0.569626 0.271696 -2.09656 3.6033e-02 *
#> poly(Wind, 2)2 -1.834732 0.250895 -7.31274 2.6176e-13 ***
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.10262 0.052341 97.48802 < 2.2e-16 ***
#> poly(Temp, 2)1 -2.63000 0.635659 -4.13744 3.512e-05 ***
#> poly(Temp, 2)2 -5.59769 0.432368 -12.94657 < 2.2e-16 ***
#>
#>
#> # SAMPLE: 6
#>
#>
#> ### Dep. var.: Ozone
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.07100 0.124756 24.61599 < 2.2e-16 ***
#> poly(Wind, 2)1 6.09901 1.922777 3.17198 0.001514 **
#> poly(Wind, 2)2 -2.28442 1.031893 -2.21381 0.026842 *
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.41843 0.080824 42.29472 < 2.2e-16 ***
#> poly(Temp, 2)1 3.26499 1.293450 2.52425 0.01159456 *
#> poly(Temp, 2)2 5.66897 1.562592 3.62792 0.00028571 ***
#>
#> ### Dep. var.: Solar.R
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.218537 0.014191 367.74415 < 2.2e-16 ***
#> poly(Wind, 2)1 2.220448 0.177688 12.49633 < 2.2e-16 ***
#> poly(Wind, 2)2 -0.624761 0.139873 -4.46664 7.9458e-06 ***
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.200307 0.019697 264.014030 < 2.2e-16 ***
#> poly(Temp, 2)1 3.286120 0.364969 9.003834 < 2.2e-16 ***
#> poly(Temp, 2)2 0.155763 0.322767 0.482586 0.62939
#>
#>
#> # SAMPLE: 7
#>
#>
#> ### Dep. var.: Ozone
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.778200 0.043180 87.49820 < 2.2e-16 ***
#> poly(Wind, 2)1 -6.722474 0.703358 -9.55768 < 2.2e-16 ***
#> poly(Wind, 2)2 -0.959048 0.681576 -1.40710 0.1594
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.59501 0.190413 13.62828 < 2.2e-16 ***
#> poly(Temp, 2)1 24.89459 3.062021 8.13011 4.2888e-16 ***
#> poly(Temp, 2)2 -8.57801 1.768706 -4.84988 1.2354e-06 ***
#>
#> ### Dep. var.: Solar.R
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.360516 0.014178 378.099115 < 2.2e-16 ***
#> poly(Wind, 2)1 -0.682752 0.216716 -3.150439 0.0016303 **
#> poly(Wind, 2)2 0.013888 0.277451 0.050054 0.9600792
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.10937 0.056186 90.93638 < 2.2e-16 ***
#> poly(Temp, 2)1 4.82809 0.928155 5.20181 1.9735e-07 ***
#> poly(Temp, 2)2 -1.11475 0.615860 -1.81007 7.0285e-02 .
#>
#>
#> # SAMPLE: 8
#>
#>
#> ### Dep. var.: Ozone
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.755844 0.041522 90.453699 < 2.2e-16 ***
#> poly(Wind, 2)1 -6.956477 0.697914 -9.967528 < 2.2e-16 ***
#> poly(Wind, 2)2 -0.344133 0.486153 -0.707869 0.47903
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.14340 0.097610 32.20355 < 2.2e-16 ***
#> poly(Temp, 2)1 17.29925 1.746264 9.90643 < 2.2e-16 ***
#> poly(Temp, 2)2 -5.78006 0.896709 -6.44586 1.1495e-10 ***
#>
#> ### Dep. var.: Solar.R
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.134474 0.016224 316.48225 < 2.2e-16 ***
#> poly(Wind, 2)1 -0.459093 0.272152 -1.68690 0.0916227 .
#> poly(Wind, 2)2 0.683835 0.246613 2.77291 0.0055558 **
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 4.926595 0.041382 119.05115 < 2.2e-16 ***
#> poly(Temp, 2)1 4.222008 0.816701 5.16959 2.3461e-07 ***
#> poly(Temp, 2)2 -0.794082 0.462440 -1.71716 8.5950e-02 .
#>
#>
#> # SAMPLE: 9
#>
#>
#> ### Dep. var.: Ozone
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.36418 0.036953 91.04039 < 2.2e-16 ***
#> poly(Wind, 2)1 -4.04407 0.492839 -8.20566 2.2933e-16 ***
#> poly(Wind, 2)2 2.98908 0.470324 6.35536 2.0794e-10 ***
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.36946 0.038292 87.99344 < 2.2e-16 ***
#> poly(Temp, 2)1 6.10018 0.690949 8.82870 < 2.2e-16 ***
#> poly(Temp, 2)2 2.25482 0.627045 3.59595 0.00032321 ***
#>
#> ### Dep. var.: Solar.R
#>
#> Expl. vars.: poly(Wind, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.122995 0.014639 349.966703 < 2.2e-16 ***
#> poly(Wind, 2)1 -0.564257 0.193672 -2.913473 0.0035743 **
#> poly(Wind, 2)2 0.144413 0.224957 0.641955 0.5209022
#> ---
#> Expl. vars.: poly(Temp, 2)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.111358 0.015302 334.02914 < 2.2e-16 ***
#> poly(Temp, 2)1 1.286054 0.251169 5.12026 3.0511e-07 ***
#> poly(Temp, 2)2 -0.928769 0.260074 -3.57117 3.5539e-04 ***
# Different way of displaying the results with "compact"
summary(est_split, "compact")
#> sample lhs rhs (Intercept) poly(Wind, 2)1
#> 1 5 Ozone poly(Wind, 2) 3.30*** (0.042) -6.41*** (0.617)
#> 2 5 Ozone poly(Temp, 2) 4.06*** (0.109)
#> 3 5 Solar.R poly(Wind, 2) 5.21*** (0.018) -0.570* (0.272)
#> 4 5 Solar.R poly(Temp, 2) 5.10*** (0.052)
#> 5 6 Ozone poly(Wind, 2) 3.07*** (0.125) 6.10** (1.92)
#> 6 6 Ozone poly(Temp, 2) 3.42*** (0.081)
#> 7 6 Solar.R poly(Wind, 2) 5.22*** (0.014) 2.22*** (0.178)
#> 8 6 Solar.R poly(Temp, 2) 5.20*** (0.020)
#> 9 7 Ozone poly(Wind, 2) 3.78*** (0.043) -6.72*** (0.703)
#> 10 7 Ozone poly(Temp, 2) 2.60*** (0.190)
#> 11 7 Solar.R poly(Wind, 2) 5.36*** (0.014) -0.683** (0.217)
#> 12 7 Solar.R poly(Temp, 2) 5.11*** (0.056)
#> 13 8 Ozone poly(Wind, 2) 3.76*** (0.042) -6.96*** (0.698)
#> 14 8 Ozone poly(Temp, 2) 3.14*** (0.098)
#> 15 8 Solar.R poly(Wind, 2) 5.13*** (0.016) -0.459. (0.272)
#> 16 8 Solar.R poly(Temp, 2) 4.93*** (0.041)
#> 17 9 Ozone poly(Wind, 2) 3.36*** (0.037) -4.04*** (0.493)
#> 18 9 Ozone poly(Temp, 2) 3.37*** (0.038)
#> 19 9 Solar.R poly(Wind, 2) 5.12*** (0.015) -0.564** (0.194)
#> 20 9 Solar.R poly(Temp, 2) 5.11*** (0.015)
#> poly(Wind, 2)2 poly(Temp, 2)1 poly(Temp, 2)2
#> 1 3.08*** (0.570)
#> 2 11.3*** (1.45) 2.51* (1.17)
#> 3 -1.83*** (0.251)
#> 4 -2.63*** (0.636) -5.60*** (0.432)
#> 5 -2.28* (1.03)
#> 6 3.26* (1.29) 5.67*** (1.56)
#> 7 -0.625*** (0.140)
#> 8 3.29*** (0.365) 0.156 (0.323)
#> 9 -0.959 (0.682)
#> 10 24.9*** (3.06) -8.58*** (1.77)
#> 11 0.014 (0.277)
#> 12 4.83*** (0.928) -1.11. (0.616)
#> 13 -0.344 (0.486)
#> 14 17.3*** (1.75) -5.78*** (0.897)
#> 15 0.684** (0.247)
#> 16 4.22*** (0.817) -0.794. (0.462)
#> 17 2.99*** (0.470)
#> 18 6.10*** (0.691) 2.25*** (0.627)
#> 19 0.144 (0.225)
#> 20 1.29*** (0.251) -0.929*** (0.260)
# You can still select which sample/LHS/RHS to display
est_split[sample = 1:2, lhs = 1, rhs = 1]
#> Standard-errors: IID
#> Dep. var.: Ozone
#> Expl. vars.: poly(Wind, 2)
#> Sample: 5
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.30403 0.042415 77.89732 < 2.2e-16 ***
#> poly(Wind, 2)1 -6.41466 0.617300 -10.39149 < 2.2e-16 ***
#> poly(Wind, 2)2 3.07584 0.570194 5.39437 6.8765e-08 ***
#> ---
#> Sample: 6
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.07100 0.124756 24.61599 < 2.2e-16 ***
#> poly(Wind, 2)1 6.09901 1.922777 3.17198 0.001514 **
#> poly(Wind, 2)2 -2.28442 1.031893 -2.21381 0.026842 *
```