I have estimated the following negative binomial regression model with group fixed effects in Stata. The data are time series cross sectional. The panelvar is group and the timevar is time.
tsset group time
xtnbreg y x1 x2 x3 + x4 + x5, fe
I want to replicate these findings in R. To do this, I have tried these 4 models:
nb1 <- femlm(y ~ x1 + x2 + x3 + x4 + x5 | group, panel.id = ~group + time, family = "negbin", mydata)
nb2 <- fenegbin(y ~ x1 + x2 + x3 + x4 + x5 | group, panel.id = ~group + time, mydata)
nb3 <- glm.nb(y ~ x1 + x2 + x3 + x4 + x5 + factor(group), data=mydata)
nb4 <- glmmadmb(y ~ x1 + x2 + x3 + x4 + x5 + factor(group), data = mydata, family="nbinom")
The results produced by nb1-4 are all identical, but different from the results produced by xtnbreg in Stata. The coefficients, standard errors, and p-values are all substantively different.
I have tried replicating a standard negative binomial regression in Stata and R and have been able to do so successfully.
Does anyone have any idea what's going on here? I have reviewed related posts on this forum (such as this one: is there an R function for Stata's xtnbreg?) and have not found any answers.
SOLVED (mostly): The R code that reproduces the results generated by xtnbreg, fe in Stata:
nb5 <- pglm(y ~ x1 + x2 + x3 + x4 + x5 ,family = negbin, data = mydata, effect = "individual", model="within", index = "group")
I found the solution on RPubs: https://rpubs.com/cuborican/xtpoisson. I still do not know why this works, only that it does. I suspect that Ben is correct and it has something to do with estimating conditional vs unconditional ML. If anyone knows for sure, please share.