I started what seemed like a straightforward analysis, but I've gotten stuck with overdispersion in my negative binomial model.
I would like to know which sites are different from each other in terms of number of calls. Can anyone tell me please how they would approach this? (we don't have any covariates).
Calls have been counted at each site for several nights over four years. I have aggregated counts of bird calls so that there is one mean per Site and Year; 'n' is the sample size of each mean.
Here is some example data:
site <- as.factor(rep(letters[1:11], each=4)) year <- as.factor(rep(c("2017","2018","2019","2020"),11)) calls <- c(222, 3778,11472,3642,2251,3008,41924,1718,284,29,2508,1610, 16,5,128,8,130,108,75,78,32,54,40,23,4,13,67,11,60,20,26,3,99,26,82,13, 2325,3487,12696,2849,48929,18309,34645,34625) n <- c(10,8,7,8,12,8,7,8,4,6,7,7,9,6,7,7,9,5,7,8,8,5,7,8,6, 7,7,8,8,7,7,7,8,7,7,7,9,9,7,7,8,10,7,9) birds <- data.frame(site,year,calls,n)
And the nb model:
require(MASS) m1 <- glm.nb(calls ~ site, weights=n,link='log',data=birds) summary(m1)
Which is overdispersed:
df_resid <- nrow(model.frame(m_nb1)) - length(coef(m_nb1)+1) pearson_resid <- residuals(m_nb1, type = "pearson") pearson_sq <- sum(pearson_resid^2) pearson_sq / df_resid
Any thoughts will be greatly appeciated!