How can I have both "firm' and 'year' fixed effect using 'bife' fixed-effect logistic regression in R? I am running a fixed-effect logistic regression using 'bife' command from the 'bife' R-package. Here, I am trying to have both FIRM and YEAR fixed-effects simultaneously.
However, the examples I can find only has one fixed-effect not two or more.
Therefore, I set up the function as below and wonder this is correct if I am trying to fix both FIRM and YEAR at the same time by using "FIRM&YEAR" since I still get the answers.
 bife(Y~ X1+ X2+X3+X4+X5 | FIRM&YEAR, data=DATA, model="logit")

May I know whether the form above is correct? 
Otherwise, may I know what the correct form for this purpose is?
 A: The R package bife does not allow for more than one fixed effect. The first impression for this is given by the help page ?bife as it only mentions the individual fixed effect; this is corroborated by looking at the description of the formula parameter:

formula
  an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. formula must be of type y ~ x | id where the id refers to an individual identifier (fixed effects).

Further, by looking at the code itself, one can conclude the package only supports one fixed effect. This line in file bife.R extracts only one variable (see, e.g. here):
id <- model.part(formula, data = mf, rhs = 2L)[[1L]]

where the last part of that line (extraction by [[1L]]) ensures only one column of the data is bound to variable id which serves for holding the fixed effects if you look further down the code in the file.
Also, the paper describing the model and method does mentions only one fixed effect (see the package vignette for the reference).
A: It's correct that bife is only designed for one-way fixed effects.
However you can include other sets of fixed effects as dummies using factor() on the RHS of | in the formula interface.
In the following I provide an example how to estimate a one-way model (individual fixed effects) and a two-way model (individual and time fixed effects).
library("bife")

# Load 'psid' dataset
dataset <- psid
head(dataset)

# One-way fixed effects logit model w/o bias-correction
mod_1 <- bife(LFP ~ AGE + I(INCH / 1000) + KID1 + KID2 + KID3 | ID, data = dataset, bias_corr = "no")


# Two-way fixed effects model w/o bias-correction
mod_2 <-  bife(LFP ~ AGE + I(INCH / 1000) + KID1 + KID2 + KID3 + factor(TIME) | ID, data = dataset, bias_corr = "no")

Please note, if you specify the bias-correction argument, the bias-correction will be only performed with respect to ID.
Also note, that the two-way specification only makes sense, if the second set of fixed effects does not have too many levels.
Otherwise this will quickly become infeasible or very time consuming, but there is a quick and feasible software routine available:
the package alpaca is designed for high-dimensional k-way fixed effects GLMs.
