# Estemating trends for extreme weather using Poisson/NegativeBinomial Regression

I have a large set of weather data containing historic information on how many days a month a certain threshholf of daily rainfall at one location is reached. My goal is to investigate if a statistcally significant trend is present in my data (not to forecast future extreme weather events ) and if it is positive or negative. Since my dataset consists of the number of events in a month (i.e. 0 = no events, 1 = 1 event a month etc.), and the month as a date (i.e "1900-01-01" for January of 1900), I believe a Poisson- or Negative Binomial Model to be appropriate.

I have used both to try and find a trend and was pretty happy with my results, but plotting the residuals left me with some questions.

Since I want to find out if the number of events increases or decreases with time, I am using the number of events per month (data) as the dependent variable and a vector representing the number of months since the first month (time <-(c(1:length(data))) as the independent variable.

Here is my code with my data sample at the bottom of this post:

library(MASS)

hist(data)
time<-c(1:length(data))
poisreg<-glm(data~time,family = poisson(link = "log")) # Trying Poisson Model to calculate Trend
cd1<-cooks.distance(poisreg)
which(cd1>=1)           # No outliers
summary(poisreg)
dispersiontest(poisreg) # Test reveals that model is overdispersed. Dispersion = 1.717436

nbreg<-glm.nb(data~time) # Switching to Negative Binomial Model, to see if Dispersion is improved
cd2<-cooks.distance(nbreg)
which(cd2>=1)          # No outliers
summary(nbreg)        # Dispersion is a lot better than it was

plot(data~time,type="h", xlab= "Months since Start", ylab = "Number of occurences per month")
abline(h=mean(data),lwd=3,col="green") # Mean of data set
lines(fitted(nbreg),lwd=2,col="red")   # Plotting Trend of nbreg with fitted values. Trend of both models is basically the same.


This gets me a significant, positive Trend

Call:
glm.nb(formula = data ~ time, init.theta = 2.780757975, link = log)
Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.8847  -0.7779  -0.1715   0.4344   2.6846
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.5035181  0.0778014   6.472 9.68e-11 ***
time        0.0006797  0.0002180   3.118  0.00182 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Negative Binomial(2.7808) family taken to be 1)
Null deviance: 684.01  on 599  degrees of freedom
Residual deviance: 674.41  on 598  degrees of freedom
AIC: 2244
Number of Fisher Scoring iterations: 1
Theta:  2.781
Std. Err.:  0.412
2 x log-likelihood:  -2238.008


with a plausible graph (red line is the trend, green line is the mean):

The residuals however are pretty bad, with clear trends present when there shouldnt be any:

Ive tried to remedy this by introducing quadratic terms to my model using poly(), which didnt help. I am guessing the problem lies somewhere with the big number of months where data=0 (since the threshold was never reached). What can I do to improve the residual plots? Is it even necessary if im only trying to find a trend in historic data (and dont want to forecast)?

Are there any models that would be more appriopriate to use for finding trends in this dataset?

Is my approach here generally ok or are there any big mistakes Ive made or things Ive missed?

Thank you!

Here is the data for the number of extreme weather events per month:

data<-c(0, 5, 1, 4, 1, 0, 1, 0, 0, 4, 0, 0, 6, 1, 0, 3, 1, 1, 2, 1,
1, 3, 2, 1, 1, 2, 2, 2, 1, 0, 1, 0, 4, 1, 1, 0, 2, 1, 0, 2, 1,
2, 4, 0, 3, 2, 2, 2, 1, 3, 1, 0, 2, 0, 2, 0, 2, 0, 2, 4, 0, 2,
1, 5, 1, 4, 3, 0, 3, 3, 1, 2, 3, 2, 0, 0, 1, 1, 3, 0, 1, 0, 3,
1, 0, 1, 2, 3, 2, 3, 1, 0, 4, 1, 2, 0, 5, 3, 2, 0, 0, 0, 1, 2,
1, 1, 2, 0, 2, 1, 2, 0, 1, 0, 2, 3, 1, 6, 2, 1, 1, 3, 2, 4, 1,
0, 1, 1, 4, 0, 1, 5, 0, 2, 4, 0, 2, 2, 2, 2, 2, 2, 0, 1, 1, 2,
4, 1, 1, 0, 5, 1, 0, 2, 1, 2, 0, 1, 0, 1, 0, 2, 1, 3, 2, 2, 3,
0, 3, 1, 2, 6, 5, 4, 3, 2, 6, 1, 2, 2, 4, 1, 4, 1, 2, 0, 6, 2,
2, 0, 9, 2, 1, 2, 4, 2, 5, 0, 1, 2, 6, 0, 2, 2, 1, 1, 2, 1, 1,
1, 4, 9, 3, 4, 0, 5, 0, 1, 2, 2, 1, 1, 2, 4, 1, 0, 1, 2, 2, 2,
2, 4, 2, 6, 3, 0, 2, 2, 0, 0, 0, 1, 2, 1, 4, 1, 1, 9, 1, 0, 3,
3, 0, 3, 3, 0, 0, 2, 6, 9, 6, 2, 0, 1, 2, 0, 4, 0, 5, 0, 1, 1,
4, 0, 0, 0, 3, 2, 0, 2, 1, 0, 2, 2, 4, 4, 1, 4, 1, 2, 3, 1, 3,
1, 6, 5, 1, 3, 1, 1, 1, 2, 1, 1, 3, 6, 3, 5, 2, 1, 1, 4, 7, 1,
6, 4, 5, 4, 2, 1, 1, 0, 0, 3, 1, 0, 1, 1, 1, 1, 2, 5, 1, 2, 2,
0, 2, 3, 2, 0, 0, 4, 0, 2, 3, 1, 1, 6, 0, 1, 2, 5, 3, 1, 2, 0,
1, 3, 0, 1, 0, 2, 0, 0, 1, 2, 6, 2, 0, 0, 0, 2, 2, 4, 0, 1, 1,
2, 0, 4, 0, 0, 2, 0, 6, 0, 5, 2, 0, 8, 4, 0, 3, 3, 0, 3, 5, 2,
4, 5, 2, 2, 1, 3, 4, 0, 1, 1, 3, 1, 2, 1, 0, 5, 3, 0, 0, 2, 4,
1, 4, 4, 1, 1, 0, 2, 2, 0, 3, 0, 1, 0, 0, 2, 3, 2, 1, 3, 0, 3,
2, 11, 2, 8, 1, 1, 0, 3, 5, 3, 2, 2, 4, 1, 0, 2, 0, 2, 2, 2,
6, 5, 2, 4, 3, 9, 0, 7, 0, 0, 2, 0, 1, 0, 0, 0, 5, 2, 1, 0, 0,
1, 4, 5, 2, 1, 0, 2, 6, 0, 7, 2, 1, 3, 1, 1, 2, 6, 6, 1, 1, 1,
4, 3, 0, 3, 3, 4, 5, 1, 5, 6, 5, 1, 3, 0, 1, 3, 1, 3, 1, 5, 1,
4, 4, 0, 2, 2, 1, 1, 5, 0, 2, 3, 1, 2, 6, 0, 1, 2, 1, 5, 4, 2,
1, 0, 3, 3, 7, 5, 1, 0, 1, 2, 4, 4, 4, 3, 9, 1, 2, 2, 5, 2, 1,
2, 1, 0, 1, 0, 6, 6, 6, 3, 2, 2, 0, 3, 5, 2, 0, 0, 1, 1, 2, 3,
2, 1, 1, 1, 0, 1, 1, 0, 2, 3, 3, 7, 0, 2)


EDIT:

Afte applying various DHARMa residual diagnostics, the neg-binom model seems to be ok. Ive added them here if anyone has a similiar problem:

• Welcome to CV. I can't comment on your analysis as a whole, but you really can't expect GLM residual plots to look like (general) LM ones. For one, the band's in the plots are because only discrete response levels are possible in these models. I'd recommend you have a look at DHARMa for GL(M)M residual diagnostics. If that still shows problems you could consider a zero-inflated model. Dec 28, 2023 at 17:33
• @PBulls Thank you for the DHARMa Tip! If I understand these residual diagnostics correctly, there shouldnt be any major problems with my nb-model. Dec 28, 2023 at 19:31

The other idea I had was to use time series methods that allow you to disaggregate a time series into trend, cyclical, and random components. There is an R package tempdisagg and I think there are other packages that do this, as well. tempdisagg has a vignette and an article associated with it (see the package documentation).