# Transfering the approach of GLM/GLM.NB to FEGLM in R to find the best dispersion parameter

I would like to analyze my dataset with around 1 million observations and 10 thousand fixed effects with a negative binomial regression model. Due to the high number of fixed effects I cannot apply 'glm' / 'glm.nb' from 'MASS' package.

To handle high dimensional fixed effects I found 'feglm' ('alpaca' package). It is working, but I have to provide a value for theta as in 'glm'. To find the optimal theta you can use 'glm.nb' (if you don't have that much variables or fixed effects). A similar function is missing for 'feglm'.

I tried to transfer the approach of 'glm.nb'. This is my result:

i=1
theta=c(1.2345)

repeat {

feglm_negbin <- feglm(y ~ x1 + x2 + x3 + x4 | fe1 + fe2 , data=dataset, family = negative.binomial(theta[i]))

i=i+1
theta.value=theta.ml(data\$Aenderungen_Technik,fitted(feglm_negbin), limit=1000)
theta[i]=theta.value
print(theta[i])

if (theta[i]-theta[i-1]<=0.00001) {break} }


As stopping criterion I chose a small value (0.00001), because I do not know the statistical correct definition of it.

• Does this procedure make sense? Is it correct?
• What is the statistical stopping criterion in 'glm.nb'? Can I transfer it to my approach?

your idea was fine. Find attached an example that shows that glm.nb() and feglm() produce the same estimates.

library(alpaca)
library(MASS)
library(pglm)

data("PatentsRD", package = "pglm")

### benchmark MASS glm.nb() ###
nb1 <- glm.nb(patent ~ rdexp + spil + factor(firm) + factor(year) + 0, data = PatentsRD)
summary(nb1)

### using alpaca feglm() ####
i <- 1
theta <- 1
repeat {
nb.alpaca <- feglm(patent ~ rdexp + spil | firm + year, data = PatentsRD,
family = negative.binomial(theta))
theta.old <- theta
theta <- theta.ml(PatentsRD$$patent, fitted(nb.alpaca), limit = 1000) print(theta) if (abs(theta - theta.old) <= sqrt(.Machine$$double.eps)) break
i <- i + 1
}

### Final estimation with converged \theta ###
nb.alpaca <- feglm(patent ~ rdexp + spil | firm + year,
data = PatentsRD, family = negative.binomial(theta))

### Compare ###
coef(summary(nb.alpaca))
coef(summary(nb1))[1:2,]
$$$$
`