0
$\begingroup$

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[1]
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?
$\endgroup$
1
$\begingroup$

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)

### Load data set ###
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,]
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.