# From overdispersion to underdispersion: comparing linear regression models with DHARMa

I have been investigating the relationship between the occurence of certain weather phenomena and time. To aid me in evaluating the fit of my (simple linear-regression) models, I have been using the DHARMa-package. So far everything has been working just fine, but today I stumbled on a problem that i cant make sense of.

My data (val) consists of the number of observed events per month, has 540 datapoints and looks like this:

The first model I evaluated was a simple poisson regression:

glm(formula = val ~ time, family = poisson((link = "log")))

Deviance Residuals:
Min      1Q  Median      3Q     Max
-3.567  -3.131  -2.414   1.742   7.021

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.8549351  0.0362014  51.239  < 2e-16 ***
time        -0.0009107  0.0001236  -7.366 1.76e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 4969.4  on 539  degrees of freedom
Residual deviance: 4914.8  on 538  degrees of freedom
AIC: 5979.6

Number of Fisher Scoring iterations: 6


where time<-c(1:length(val)). As you can see the model is obviously overdispersed with residual deviance/ residual deegres of freedom = 9.135 ( this is how Ive been taught to check a model for over/underdispersion). Using testDispersion() in DHARMa to check the dispersion we get a very similar result:

DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 9.2728, p-value < 2.2e-16
alternative hypothesis: two.sided


This leads me to switch to a negative binomial regression using glm.nb() in package MASS:

Call:
glm.nb(formula = val ~ time, init.theta = 0.2718374341, link = log)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.3190  -1.2655  -0.6708   0.3576   1.2432

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.8603669  0.1692592  10.991   <2e-16 ***
time        -0.0009316  0.0005439  -1.713   0.0868 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.2718) family taken to be 1)
Null deviance: 513.73  on 539  degrees of freedom
Residual deviance: 510.86  on 538  degrees of freedom
AIC: 2628
Number of Fisher Scoring iterations: 1
Theta:  0.2718
Std. Err.:  0.0218
2 x log-likelihood:  -2622.0430


Checking the dispersion manually like before gets me a result of 0.949, which seems fine to me.testDipserion() however, detects significant underdispersion:

   DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated

data:  simulationOutput
dispersion = 0.46026, p-value < 2.2e-16
alternative hypothesis: two.sided


Which leaves me to wonder which one is correct. Is my model underdispersed or not ? And does anybody know what causes the difference?

....................

I am including some more tests from DHARMa from my neg-binom model in the case the answer somehow lies there:

• First, I wouldn't rely too much on p values from tests of assumptions. They are conflated with sample size. Second, it sure looks like you have zero-inflation. I would try a zero-inflated Poisson, that might deal with your dispersion issue. Without those 0's it looks like the DV is almost uniform. Jan 6 at 20:17
• I am not sure if a zero-inflated model is appropriate here. For some more context: my original data describes daily air tempratures at one specific station. If the daily temperature exceeds 10 degress C an event occurs and is recorded as a 1 in the dataset. If the threshold is not met, no event happens which is recorded as a 0. This produces a dataset of binary values for each day of the year. I have aggreagated this Data into monthly counts, which I am using in my model. If I understand zero inflation right, there is one process for generating 0 and one for counts. Continued -> Jan 6 at 22:17
• Since my data is only dependent on daily air temperature "false zeros" shouldnt really be possible, since there are no interferences like in the often cited fishing example. BUt im probably misunderstanding something here Jan 6 at 22:24
• ZIP doesn't require false 0's. It's true that there are two processes but they can use the same variables. Jan 6 at 22:28

The most common overdispersion tests in the literature — comparing (resid deviance)/(resid df) to 1, or (sum of [Pearson resids]^2)/(resid df) to 1, or computing the one-tailed p-value of resid deviance or Pearson SSQ against a $$\chi^2_{\textrm{rdf}}$$ distribution) are approximate tests: from Venables and Ripley

A common way to 'discover' over- or underdispersion is to notice that the residual deviance is appreciably different from the residual degrees of freedom, since in the usual theory the expected value of the residual deviance should equal the degrees of freedom. This can be seriously misleading. The theory is asymptotic, and only applies for large $$n_i p_i$$ for a binomial and for large $$\mu_i$$ for a Poisson ...

(They also discuss the Pearson SSQ criterion, which they say is much less biased, but which is also asymptotic.)

In contrast, the DHARMa test is based on simulations from the fitted model — it is computationally intensive, but should be much more reliable. That said, any problem with your model (zero-inflation, outliers, nonlinearity, etc. etc.) can in principle lead to overdispersion. Underdispersion is a little bit harder to create, so it could indeed be an issue with the shape of your conditional distribution.

That said, people often worry less about under- than overdispersion, because the general effect of underdispersion is to make statistical inference conservative ...

Venables, W., and Brian D. Ripley. 2002. Modern Applied Statistics with S. 4th ed. New York: Springer.

What Ben says but also:

1. What I mainly see in your plots is a distributional misfit (wiggle in qq plot) and autocorrelation. The underdispersion is likely a secondary problem.
2. The misfit is not particularly surprising because your data-generating process, as you describe it, seems the be a clear case for a binomial rather than a Poisson distribution (essentially you have a coin flip per day, so you have k/n data per Month, not counts)

-> my recommendation would be to switch to a k/n binomial, check and then address possible overdispersion, 0/1 inflation or autocorrelation as discussed, e.g. in https://theoreticalecology.github.io/AdvancedRegressionModels/