I am fitting a poisson regression model in R to count time series data to perform an Interrupted Time Series Analysis, the aim of my analysis is to see if an intervention affected the counts. I am basing my analysis on this tutorial: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5407170/
My dependent variable is counts. I have used harmonic terms specifiying sine & cosine pairs to adjust for seasonality.
Having fit a poisson model, I observed overdispersion:
Residual deviance: 2205.4 on 57 degrees of freedom
And so instead I fit a quasipoisson model.
model2 <- glm(count ~ intervention + time + harmonic(month,2,12), family=quasipoisson, data= data)
When modelling checking using the below code for the deviance goodness of fit test:
pchisq(model$deviance, df=model$df.residual, lower.tail=FALSE)
All of my models were showing p< 0.05 which if i am correct means I fail to reject the null hypothesis and so none of the models are correctly specified. Furthermore when checking for autocorrelation using the acf and pacf plots, all of my proposed models residuals were showing autocorrelation (ACF plot attached below).
I have tried fitting a negative binomial model including harmonic terms specifiying sine & cosine pairs to adjust for seasonality
model3 <- glm.nb(count ~ intervention + time + harmonic(month,2,12), data= data)
Using the same code as above for the deviance goodness of fit test my p value was 0.242 (P>0.05 meaning we accept null hypothesis, model is correctly specified.)
However when checking for autocorrelation using acf & pacf plots this model's residuals also showed autocorrelation.
Does this autocorrelation mean that my models are specified wrong? Even though my deviance goodness of fit test is suggesting the model is specified correctly?
How can I rectify this autocorrelation? I do not have any prior statistics knowledge, I have seen a lot of literature suggesting an ARIMA model, but not sure I have the experience for this.
