# Comparing quasi-Poisson and negative binomial fit on panel data

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 =

#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. • it looks like your dataset is small, so bear in mind, it can be really tricky to model the variance - mean relationship. I would expect the significance of the coefficients to be different, but not the coefficient itself Sep 9, 2021 at 5:24
• from what you showed, it looks like coefficients are really different, which is weird.. might be good to check that the fit was done correctly Sep 9, 2021 at 5:25
• @StupidWolf Sorry for the dumb question, but how do I check that the fit was done correctly? Sep 9, 2021 at 19:16
• Could you check the package version? The error when printing the results does not pop for me. Sep 15, 2021 at 11:51

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.