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!

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.
