Our package offers a fast algorithm to estimate generalized linear models with high-dimensional fixed effects. The linear predictor of such models takes the following form: \[ \boldsymbol{\eta} = \mathbf{Z} \boldsymbol{\gamma} = \mathbf{D} \boldsymbol{\alpha} + \mathbf{X} \boldsymbol{\beta} = \sum_{j=1}^{k} \mathbf{D}_j \boldsymbol{\alpha}_{j} + \mathbf{X} \boldsymbol{\beta} \, , \] where the matrix \(\mathbf{D}\) arises from dummy encoding of \(k\) high-dimensional categorical variables and \(\mathbf{X}\) contains the variables of interest. We refer to \(\boldsymbol{\beta}\) and \(\boldsymbol{\alpha}\) as structural parameters and fixed effects. The latter are essentially nuisance parameters that are used to control for unobserved heterogeneity.

Brute force estimation of this kind of models is often restricted to
computational limitations (either due to memory or time limitations). We
tackle this problem by providing a fast and memory efficient algorithm
suggested by Stammann (2018) based on the
combination of the Frisch-Waugh-Lovell theorem and the method of
alternating projections. We restrict ourselves to non-linear models
because Gaure (2013) already offers a
great package for linear models. Further, in the case of binary choice
models with only one high-dimensional fixed effects we recommend using
the package `bife`

.

In the following we utilize an example from labor economics to
demonstrate the capabilities of `feglm()`

. More precisely, we
use a balanced micro panel data set from the *Panel Study of Income
Dynamics* to analyze the intertemporal labor force participation of
1,461 married women observed for nine years. A similar empirical
illustration is used in Fernández-Val
(2009).

Before we start, we briefly inspect the data set to get an idea about its structure and potential covariates.

```
data(psid, package = "bife")
head(psid)
```

```
## ID LFP KID1 KID2 KID3 INCH AGE TIME
## 1: 1 1 1 1 1 58807.81 26 1
## 2: 1 1 1 0 2 41741.87 27 2
## 3: 1 1 0 1 2 51320.73 28 3
## 4: 1 1 0 1 2 48958.58 29 4
## 5: 1 1 0 1 2 53634.62 30 5
## 6: 1 1 0 0 3 50983.13 31 6
```

`ID`

and `TIME`

are individual and
time-specific identifiers, `LFP`

is an indicator equal to one
if a woman is in labor force, `KID1`

- `KID3`

are
the number of children in a certain age group, `INCH`

is the
annual income of the husband, and `AGE`

is the age of a
woman.

First, we use a specification similar to Fernández-Val (2009) and estimate a static model of female labor supply where we control for individual and time-specific unobserved heterogeneity (so called individual and time fixed effects).

```
library(alpaca)
<- feglm(
stat ~ KID1 + KID2 + KID3 + log(INCH) | ID + TIME,
LFP data = psid,
family = binomial("probit")
)summary(stat)
```

```
## binomial - probit link
##
## LFP ~ KID1 + KID2 + KID3 + log(INCH) | ID + TIME
##
## Estimates:
## Estimate Std. error z value Pr(> |z|)
## KID1 -0.676898 0.056301 -12.023 < 2e-16 ***
## KID2 -0.344370 0.049897 -6.902 5.14e-12 ***
## KID3 -0.007045 0.035344 -0.199 0.842
## log(INCH) -0.234128 0.054403 -4.304 1.68e-05 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## residual deviance= 6069.65,
## null deviance= 8152.05,
## n= 5976, l= [664, 9]
##
## ( 7173 observation(s) deleted due to perfect classification )
##
## Number of Fisher Scoring Iterations: 7
```

As `glm()`

, the summary statistic of the model provides
detailed information about the coefficients and some information about
the model fit (`residual deviance`

and
`null deviance`

). Furthermore, we report statistics that are
specific to fixed effects models. More precisely, we learn that only
5,976 observations out of 13,149 contribute to the identification of the
structural parameters. This is indicated by the message that 7,173
observations are deleted due to perfect classification. With respect to
binary choice models those are observations that are related to women
who never change their labor force participation status during the nine
years observed. Thus those women were either always employed or
unemployed.^{1} Overall the estimation results are based on
664 women observed for nine years.

Because coefficients itself are not very meaningful, econometricians
are usually interested in so called partial effects (also known as
marginal or ceteris paribus effects). A commonly used statistic is the
average partial effect. `alpaca`

offers a post-estimation
routine to estimate average partial effects and their corresponding
standard errors.^{2}

```
<- getAPEs(stat)
apes.stat summary(apes.stat)
```

```
## Estimates:
## Estimate Std. error z value Pr(> |z|)
## KID1 -0.0880155 0.0077937 -11.293 < 2e-16 ***
## KID2 -0.0447776 0.0068196 -6.566 5.17e-11 ***
## KID3 -0.0009161 0.0050062 -0.183 0.855
## log(INCH) -0.0304432 0.0077090 -3.949 7.85e-05 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

A widespread reason that prevents the use of non-linear fixed effects
models in practice is the so-called incidental parameter bias problem
(*IPP*) first mentioned by Neyman and
Scott (1948). Fortunately, for classical panel data sets, like in
this example, there already exist several asymptotic bias corrections
tackling the *IPP* (see Fernández-Val and
Weidner (2018) for an overview).^{3} Our package provides a
post-estimation routine that applies the analytical bias correction
derived by Fernández-Val and Weidner
(2016).^{4} Technical details on how the bias
correction is accelerated using the method of alternating projections
can be found in Czarnowske and Stammann
(2020).

```
<- biasCorr(stat)
stat.bc summary(stat.bc)
```

```
## binomial - probit link
##
## LFP ~ KID1 + KID2 + KID3 + log(INCH) | ID + TIME
##
## Estimates:
## Estimate Std. error z value Pr(> |z|)
## KID1 -0.596285 0.055528 -10.738 < 2e-16 ***
## KID2 -0.303346 0.049517 -6.126 9e-10 ***
## KID3 -0.006117 0.035211 -0.174 0.862081
## log(INCH) -0.207061 0.053928 -3.840 0.000123 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## residual deviance= 6069.65,
## null deviance= 8152.05,
## n= 5976, l= [664, 9]
##
## ( 7173 observation(s) deleted due to perfect classification )
##
## Number of Fisher Scoring Iterations: 7
```

```
<- getAPEs(stat.bc)
apes.stat.bc summary(apes.stat.bc)
```

```
## Estimates:
## Estimate Std. error z value Pr(> |z|)
## KID1 -0.096501 0.007620 -12.664 < 2e-16 ***
## KID2 -0.049093 0.006766 -7.255 4.01e-13 ***
## KID3 -0.000990 0.004987 -0.198 0.843
## log(INCH) -0.033510 0.007588 -4.416 1.00e-05 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

Whereas analytical bias corrections for static models get more and more attention in applied work, it is not well known that they can also be used for dynamic models with fixed effects.

Before we can adjust our static to a dynamic specification, we first have to generate a lagged dependent variable.

```
library(data.table)
setDT(psid)
:= shift(LFP), by = ID] psid[, LLFP
```

Contrary to the bias correction for the static models, we need to
additionally provide a bandwidth parameter (`L`

) that is
required for the estimation of spectral densities (see Hahn and Kuersteiner (2011)). Fernández-Val and Weidner (2016) suggest to do a
sensitivity analysis and try different values for `L`

but not
larger than four. Note that in this case the order of factors to be
concentrated out, specified in the second part of the formula, is
important (cross-sectional identifier first and time identifier
second).

```
<- feglm(
dyn ~ LLFP + KID1 + KID2 + KID3 + log(INCH) | ID + TIME,
LFP data = psid,
family = binomial("probit")
)<- biasCorr(dyn, L = 1L)
dyn.bc summary(dyn.bc)
```

```
## binomial - probit link
##
## LFP ~ LLFP + KID1 + KID2 + KID3 + log(INCH) | ID + TIME
##
## Estimates:
## Estimate Std. error z value Pr(> |z|)
## LLFP 1.01607 0.04759 21.350 < 2e-16 ***
## KID1 -0.45387 0.06811 -6.664 2.67e-11 ***
## KID2 -0.15736 0.06116 -2.573 0.01008 *
## KID3 0.01561 0.04406 0.354 0.72315
## log(INCH) -0.18833 0.06231 -3.023 0.00251 **
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## residual deviance= 4777.58,
## null deviance= 6549.14,
## n= 4792, l= [599, 8]
##
## ( 1461 observation(s) deleted due to missingness )
## ( 6896 observation(s) deleted due to perfect classification )
##
## Number of Fisher Scoring Iterations: 6
```

```
<- getAPEs(dyn.bc)
apes.dyn.bc summary(apes.dyn.bc)
```

```
## Estimates:
## Estimate Std. error z value Pr(> |z|)
## LLFP 0.186310 0.006686 27.864 < 2e-16 ***
## KID1 -0.072321 0.007832 -9.235 < 2e-16 ***
## KID2 -0.025074 0.007003 -3.580 0.000343 ***
## KID3 0.002487 0.005008 0.497 0.619447
## log(INCH) -0.030009 0.007002 -4.286 1.82e-05 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

`alpaca`

is also compatible with
`linearHypothesis()`

of the `car`

package and
offers the possibility to compute robust and multi-way clustered
standard errors. Further it supports the two- and three-way bias
corrections suggested by Hinz, Stammann, and
Wanner (2020). For examples, see vignette “Estimating the
intensive and extensive margin of trade”.

Cruz-Gonzalez, Mario, Iván Fernández-Val, and Martin Weidner. 2017.
“Bias Corrections for Probit and Logit Models with Two-Way Fixed
Effects.” *The Stata Journal* 17 (3): 517–45.

Czarnowske, Daniel, and Amrei Stammann. 2020. “Fixed Effects
Binary Choice Models: Estimation and Inference with Long Panels.”
*ArXiv e-Prints*.

Fernández-Val, Iván. 2009. “Fixed Effects Estimation of Structural
Parameters and Marginal Effects in Panel Probit Models.”
*Journal of Econometrics* 150 (1): 71–85.

Fernández-Val, Iván, and Martin Weidner. 2016. “Individual and
Time Effects in Nonlinear Panel Models with Large n, t.”
*Journal of Econometrics* 192 (1): 291–312.

———. 2018. “Fixed Effects Estimation of Large-t Panel Data
Models.” *Annual Review of Economics* 10 (1): 109–38.

Gaure, Simen. 2013. “lfe: Linear Group
Fixed Effects.” *The R Journal* 5 (2): 104–17.

Hahn, Jinyong, and Guido Kuersteiner. 2011. “Bias Reduction for
Dynamic Nonlinear Panel Models with Fixed Effects.”
*Econometric Theory* 27 (6): 1152–91.

Hinz, Julian, Amrei Stammann, and Joschka Wanner. 2020. “State
Dependence and Unobserved Heterogeneity in the Extensive Margin of
Trade.” *ArXiv e-Prints*.

Neyman, Jerzy, and Elizabeth L. Scott. 1948. “Consistent Estimates
Based on Partially Consistent Observations.”
*Econometrica* 16 (1): 1–32.

Stammann, Amrei. 2018. “Fast and Feasible Estimation of
Generalized Linear Models with High-Dimensional k-Way Fixed
Effects.” *ArXiv e-Prints*.

Note that in this specification (with individual and time fixed effects) also observations related to a specific time period where all women are either in labor force or not can be dropped. However this is very unlikely in practice.↩︎

The routine is currently restricted to binary choice models but will be extended in the future.↩︎

Cruz-Gonzalez, Fernández-Val, and Weidner (2017) apply the same bias correction to a pseudo panel of bilateral trade flows.↩︎

See footnote 2.↩︎