Comparing quasi-Poisson and negative binomial fit on panel data

Question

I am using the R package fixest to model over-dispersed count data with fixed effects.

While I have read that quasi-poisson and negative binomial regressions typically return similar results, this is not true in my case. When I use a quasi-poisson distribution with fixed effects, the model returns a statistically significant coefficient with a positive sign. But when I use a negative binomial distribution with fixed effects, the model does not produce a statistically significant result and the sign of the coefficient is negative (it reverses).

My first questions is, why do I get such different results?

    library(tidyverse)
    library(fixest)
    df_wall_deaths =   
 read.csv('https://gist.githubusercontent.com/Alex-Devoid/ecf9a727dc885430737f40b38ceed95f/raw/2bec498a5b9ad1520a75f206cb612671639366a0/panel_data.csv')

    #The percent_unwalled_border variable is a percentage 
    # expressed as a decimal
    #CORNO are corridors
    res_quasipoisson = feglm(death_counts ~ 
        percent_unwalled_border | CORNO, data =df_wall_deaths,  
          family = "quasipoisson")
    esttable(res_quasipoisson,order = "Dist")

    res_negbin = femlm(death_counts ~ percent_unwalled_border | 
      CORNO, df_wall_deaths,family = 'negbin')
      esttable(res_negbin,order = "Dist")
    # when I try to print res_negbin without esttable, it also 
    # returns an error:
    # Error in coeftable[1, , FALSE] : incorrect number of 
    # dimensions

My second question is, how can I compare the quasi-Poisson and negative binomial fit on this panel data?

Since I am working with panel data and using fixed effects, how can I compare the observed relationship between the mean and the variance in my dataset to the expected relationship from the models?

Below I have compared this relationship as if the models were not using fixed effects. But since fixed effects compare each corridor (CORNO) to itself, is it as simple as finding the mean and variance of arbitrary groups within each corridor? I have tried several things, including using a facet chart to run the comparison in each corodior, but I'm just not sure what the correct approach is.

    k <- res_negbin$theta
    
    THETA <- (sum(residuals(res_quasipoisson, 
   type='pearson')^2)/degrees_freedom(res_quasipoisson_clustered, 
   type = "resid"))
    
    ##### 
    #Arbitrarily divide data into 10 groups
    df_wall_deaths<-df_wall_deaths %>%    
    mutate(groups=cut_number(percent_unwalled_border, n=10))
    #Calculate mean and variance for each group
    mean_variance <- df_wall_deaths %>%
      group_by(groups)%>%
      summarize(percenrt_unwalled_mean = 
        mean(percent_unwalled_border),
                death_counts_mean = mean(death_counts),
                variance = var(death_counts))     
    THETA <- (sum(residuals(res_quasipoisson, 
   type='pearson')^2)/degrees_freedom(res_quasipoisson_clustered,
   type = "resid"))
    
    #plot the predicted mean-variance curves
    ggplot(mean_variance, aes(death_counts_mean,variance)) +
      geom_point() +
      #quasi-poisson
      geom_abline(slope=THETA) +
      #Negbin
      stat_function(fun=function(x)x+x^2/k,color="red")

EDIT: I'm not sure if this is within the scope of this website, but I wanted to share a bit more:

Since I posted my questions, made a judgment call to use to the quasi-poisson distribution because it gives more weight to higher counts, while the negative binomial distribution weights counts more evenly, giving comparatively more weight to smaller counts within the dataset.

I think it's more appropriate to let the quasi-poisson distribution give more importance or weight to higher counts because we want to know more about the corridors where most migrant deaths are found. Do corridors with high counts also have a large share of Arizona’s unwalled border in a given year, or not? I asked if this was an acceptable reason to choose the quasi-poisson distribution in this post: Comparing fixed effects quasi-Poisson and negative binomial fix

This data starts in 1990 when things were very different at the border. While I don't consider this dataset to be to be zero-inflated, for reasons I won't get into here, in the early 1990s more corridors had zero and low counts in these earlier years before it was Border Patrol’s strategy to control urban areas and redirect traffic to more remote areas. I think a negative binomial model that weighs these low counts the same as higher counts in later years could skew the model.

Answer 1

On the first point: the assumption that the quasi-poisson leads to the same estimates as the negative binomial is wrong. They are two different models.

I guess what the OP has in mind is unbiasedness. The two estimators are indeed unbiased ($E[\hat{\beta}]=\beta$) but this does not mean that they lead to the same estimates. Especially that efficiency differ between the two methods.

---

On the second point: Even discounting for the FE structure, isn't your graph already sufficient to decide between the two methods? It really looks like the Poisson estimate is driven up by the large values (although, again, on expectation this estimator is unbiased -- but the estimate may be far off the true one for this particular sample).