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.