# Modelling the difference of two zero-inflated count variables

I'm trying to model count data on house sparrows numbers before and after mitigation measures and want to find out which of these measures affect their numbers (if any). I'm however not entirely sure how to tackle this issue.

My first idea was to combine the housesparrow numbers before and after the mitigation into one 'difference' variable. Because it is count data I was going to use a Poisson distribution, but the 'difference' variable of course contains negative values, so I can't. There's also the problem that a lot of migation was conducted without housesparrows actually being present, so the dataset has a lot of zeros and is thus probably zero-inflated. It is at the very least not normally distributed, so a model assuming a normal distribution is likely out of the question.

In other threads with similar issues I have seen people suggest using a negative binomial model, but as far as I know, the dependent variable in that model can't be negative either. Another thing I've seen suggested is keeping both counts of housesparrows (before and after mitigation) in the model, but I'm not sure how I would do that without making one of the counts the dependant variable and the other an independent variable alongside the already present independent variable of different mitigation measures.

Any help at all would be greatly appreciated. Thanks in advance!

• Conditional normality in the observations is often nice to have, but surprisingly often not all that necessary (what is important for statistical inference is the normal distribution of parameters, which is a different thing). I don't think there are standard tools for your kind of problem - even the difference of two non-zero-inflated Poissons, a skellam-distribution, is rarely seen in a regression context. You may want to look at nonparametric alternatives, like permutation tests. (And if you redo your analysis with OLS, you may find there is little difference.) Commented Jan 4, 2022 at 14:53
• So as I understand it, you have a count variable (number of house sparrows) and one explanatory variable (before migration, after migration). Why not just run some sort of GLM, like a negative binomial regression as you mentioned, adjusting for time period as a predictor? Commented Jan 4, 2022 at 15:03
• @StephanKolassa A Skellam distribution was indeed another thing I've seen suggested for such an issue. The 'skellam' package actually features a regression option, but the problem is I'm only used to rather simple glm's, so I'm not sure how to use and interpret it. Using nonparametric alternatives is a good idea, so I'll try and look into using permutations tests. Commented Jan 5, 2022 at 9:48
• @DemetriPananos I have two count variables before and after a treatment and want to see which treatment yields the most positive (or negative) result by making a regression of the difference of these two count variables. Time isn't really an independent variable, but more just to differentiate between the two count variables (so before and after the treatment). And like I said in my post, I can't use a negative binomial because the difference of these two count measures contain negative values because of it being a difference measure. Does that make it more clear? Commented Jan 5, 2022 at 9:52

You do not have to take the difference in measurements in order to compare before/after treatment. Use a binary indicator variable to represent before (0) and after(1) and then simply run a regression on those data.

Here is an example of what your data might look like


period count
<int> <dbl>
1      0     7
2      1     4
3      0     4
4      1     4
5      0     2
6      1     5
7      0     4
8      1     3
9      0     5
10      1     3


Again, 0 represents before the treatment and 1 represents after the treatment. Conditional on period, all the counts are non-negative, so we can still apply negative binomial regression. In R...

fit = MASS::glm.nb(count~period, data = d, control = glm.control(maxit=100))
summary(fit)
Call:
MASS::glm.nb(formula = count ~ period, data = d, control = glm.control(maxit = 100),
init.theta = 0.9889630259, link = log)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.7268  -0.8516  -0.4340   0.3434   2.8197

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.96241    0.03732  25.785  < 2e-16
period       0.28376    0.05188   5.469 4.52e-08

(Intercept) ***
period      ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.989) family taken to be 1)

Null deviance: 2247.0  on 1999  degrees of freedom
Residual deviance: 2217.1  on 1998  degrees of freedom
AIC: 9027.4

Number of Fisher Scoring iterations: 1

Theta:  0.9890
Std. Err.:  0.0459

2 x log-likelihood:  -9021.4490


We see that post treatment, the average countk changed by a factor of $$\exp(0.28)$$ (the expoenntial of the period coefficient). The estimate is statistically significant as well.

There is no need to take differences if your intent is to compare means and you are willing to make assumptions about the likelihood of the data. Alternatively, if you have lots of data, a t-test would also be applicable despite the discrete nature of the data. If neither of these are interesting to you, you could always bootstrap the differences.