2
$\begingroup$

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): enter image description here

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

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:

enter image description here

enter image description here

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – PBulls
    Dec 28, 2023 at 17:33
  • $\begingroup$ @PBulls Thank you for the DHARMa Tip! If I understand these residual diagnostics correctly, there shouldnt be any major problems with my nb-model. $\endgroup$
    – Sargnagel
    Dec 28, 2023 at 19:31

1 Answer 1

0
$\begingroup$

Well, my first thought while reading this was "what about the month of the year?" and the plots reinforce that, as there are 12 bands. You note that this is all at one location, and, while you don't say which location, at most places in the world there are rainy seasons and dry seasons. So, I'd add "month" to the model, and maybe change "time" to "year".

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).

$\endgroup$
2
  • $\begingroup$ Thank you for your answer! While your suggestion of adding months to my independent variables doesnt seem to br necessary to fix the issue at hand, i will still to do so as they offer valuable insight into seasonal data trends. $\endgroup$
    – Sargnagel
    Dec 28, 2023 at 19:40
  • $\begingroup$ The bands are due to the outcome ranging from 0-11 occurrences (by chance, and there's no 10 observed so only 11 distinct bands) rather than that being the number of months in a year, but it's still a good suggestion to check if there are any cyclical patterns. $\endgroup$
    – PBulls
    Dec 28, 2023 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.