2
$\begingroup$

I am running Quasi-Poisson generalized linear models (overdispersed data) on my count data, and plotting data and the fitted model with ggplot.

My data with DPUT:

structure(list(YR = c(1960, 1961, 1962, 1963, 1964, 1965, 1966, 
1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 
1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 
1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 
2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 
2011, 2012, 2013, 2014, 2015, 2016), count = c(1L, 3L, 2L, 4L, 
3L, 5L, 3L, 10L, 9L, 6L, 9L, 7L, 7L, 10L, 7L, 11L, 13L, 7L, 10L, 
12L, 15L, 12L, 11L, 9L, 10L, 16L, 21L, 9L, 10L, 9L, 18L, 19L, 
14L, 14L, 13L, 21L, 20L, 18L, 28L, 21L, 27L, 26L, 25L, 24L, 25L, 
34L, 37L, 31L, 59L, 45L, 49L, 42L, 48L, 65L, 52L, 62L, 49L), 
LEAD = c(1L, 3L, 2L, 4L, 3L, 5L, 3L, 10L, 9L, 6L, 9L, 7L, 
7L, 10L, 7L, 11L, 13L, 7L, 10L, 12L, 15L, 12L, 11L, 9L, 10L, 
16L, 21L, 9L, 10L, 9L, 18L, 19L, 14L, 14L, 13L, 21L, 20L, 
18L, 28L, 21L, 27L, 26L, 25L, 24L, 25L, 34L, 37L, 31L, 59L, 
45L, 49L, 42L, 48L, 65L, 52L, 62L, 49L)), .Names = c("YR", 
"count", "LEAD"), row.names = 27:83, class = "data.frame")

Here's my script:

mod <- glm(count~YR, data = table60, family = "poisson")
pchisq(mod$deviance,df=mod$df.residual,lower.tail = F)  # P value for GOF
mod$deviance/mod$df.residual   # overdispersion/underdispersion
mod.q <- glm(count~YR, data = table60, family = "quasipoisson")

predq.df <- data.frame(YR = seq(min(table60$YR), max(table60$YR), length.out = 100))
predq <- predict(mod.q, newdata = predq.df, se.fit = TRUE,type="response")
predq.df$count <- predq$fit
predq.df$countmin <- predq$fit - 2 * predq$se.fit
predq.df$countmax <- predq$fit + 2 * predq$se.fit

ggplot(table60,aes(x=YR,y=count)) + 
  scale_y_continuous(limits=c(-1,70),breaks=c(0,10,20,30,40,50,60,70,80),expand=c(0,0)) +
  scale_x_continuous(limits=c(1960,2018),breaks=c(1960,1965,1970,
                                              1975,1980,1985,1990,1995,
                                              2000,2005,2010,2015)) +
  geom_point() +
  geom_ribbon(data = predq.df, aes(y=count, ymin = countmin, ymax = countmax), alpha=0.3) +
  geom_line(data = predq.df,aes(x=YR,y=count)) +
  xlab(" ") + ylab("Count") +
  ggtitle("ALL DATA - Quasi-Poisson Generalized Linear Model")

As I want to publish/report this with some p-value indicating the goodness-of-fit for the model I was advised to include the synthax in the 2nd line of the script:

pchisq(mod$deviance,df=mod$df.residual,lower.tail = F)  # P value for GOF

that would give a p-value for the GOF of the model.

I have now done this for many subsets of my data, but I can't make sense out of these p-values. The above script produces the plot below, with a p-value of 0.22, which seems to me to be too high. And many of the subsets of data treated with the same procedure come p with p-values for GOF that seems dead wrong:

enter image description here

Are anyone able to see what could be wrong, if anything?

$\endgroup$
2
$\begingroup$

Here is the analysis of deviance table from R for the Poisson regression:

> anova(mod, test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: count

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                    56     687.16              
YR    1   624.55        55      62.61 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R does not set the table out very well, so the table needs some explanation. The total deviance for this data is 687.16 on 56 df. That represents the total variability of the counts about their mean. The trend with Year is very strong. It explains most of the deviance by itself. The deviance for the trend is 624.55, which is 91% of the total deviance, on just one df. The trend with Year is very highly significant (P < 2.2e-6). Actually the correct p-value for the trend is 7.7e-138, which is very, very small:

> pchisq(624.55, df=1, lower.tail=FALSE)
[1] 7.65869e-138

The trend with Year is so strong that it explains all except for 62.61 of the deviance. In other words, the residual deviance is 62.61 on 55 df. The residual deviance is so small that it is compatible with Poisson variability. This means that the regression on Year explains so much that the rest just seems to be random variation. Trying to put any more predictors in the regression would likely be over-fitting.

The goodness of fit test simply compares the residual deviance to a chi-square distribution. This is a reasonable thing to do as long as the counts are not too small (less than or equal to 2). The goodness of fit test is not significant:

> pchisq(62.61, df=55, lower.tail=FALSE)
[1] 0.2243739

This tells you that (i) there is no evidence that you need to include any more predictors in the model and (ii) the Poisson model seems reasonable. In other words, your regression on Year is as good as it is possible to be.

Note that the Poisson regression and the quasi-Poisson regression give identical fitted values and identical trends with time. The only difference is that the p-values from the quasi-Poisson regression will be slightly higher:

> fit.q <- glm(count ~ YR, family = "quasipoisson", data=table60)
> anova(fit.q, test="F")
Analysis of Deviance Table

Model: quasipoisson, link: log

Response: count

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
NULL                    56     687.16                     
YR    1   624.55        55      62.61 537.48 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The exact p-value now is 4.6e-30:

> pf(537.48, df1=1, df2=55, lower.tail=FALSE)
[1] 4.587874e-30

It's still very small.

$\endgroup$
  • $\begingroup$ I have many subsets from these data, where many of the annual counts are zero, and where the p-value for GOF is a lot higher even though the fit looks ok. I don't understand the p-value of 0.22 when the fit seems to be good. What does it tell me? What estimates and/or p-values do I need to report when explaining it in text in a results chapter? Also, for the subsets the mod$deviance/mod$df.residual would vary between 0.5 and 1.5, which would indicate under- and overdispersion, in which case a Quasi-Poisson would be more correct than Poisson, am I right? $\endgroup$ – Dag Feb 21 '17 at 8:58
  • $\begingroup$ No, I don't think so. I'm having trouble following what you don't understand. The p-value obviously tells you that the null hypothesis of the no over-dispersion is accepted. The p-value helps justify the Poisson regression and appears to be what you have been asked for. There seems no logic to considering subsets of the data. It seems that your agenda is to justify quasi-Poisson rather than to assess data objectively. You can estimate an over-dispersion coefficient if you want, it does little harm, but any reader will see from the goodness of fit test that it was not really needed. $\endgroup$ – Gordon Smyth Feb 21 '17 at 20:57
  • $\begingroup$ Thank you again Gordon. I can be accused for not understanding enough statistics, but not in trying not to assess my data objectively, because that is exactly what I want. I just don't understand it well enough. What I don't understand is simple. The plot clearly states that there is an increase in the count with YEAR. But as I understand the p-value of the GOF is that the H0 must be accepted i.e. that there is no increase in the count with YEAR? I must be wrong about something, but I am not entirely certain where I am wrong. $\endgroup$ – Dag Feb 21 '17 at 21:08
  • $\begingroup$ Oh, I see now. Yes, there is a very strong trend with Year. That's not what the goodness of fit is testing. I'll edit my answer to try to make it clear for you. $\endgroup$ – Gordon Smyth Feb 21 '17 at 22:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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