I’m working with a poisson model to look at relationships between fish and habitat. I have multiple-pass count data (the same pool snorkeled multiple times), and I’m using a program called “unmarked”, which uses hierarchical models to estimate both abundance (poisson) and detection probability (binomial). My models don’t seem to have great fit, and I’d like to be able to consider and account for overdispersion. The program doesn’t allow me to use a quasipoisson or have a function to get deviance residuals, so I’m looking to calculate deviance residuals, find a dispersion parameter, and adjust my AIC and SE values by “hand”.

I’ve put some code together based on my statistics book (The Statistical Sleuth), and what I’ve read from several sources online, but I’d very much appreciate it if some folks with more experience could take a look at the code and tell me if it makes sense and is appropriate. After running this code, I do get a dispersion parameter that’s similar to what I get when I run a quasipoisson glm on just the abundance model, and everything in the code seems to function just fine, but I’d like a second opinion. Also, I do have a couple of observed values of 0, which give me a negative infinity value for the individual deviance residual, and I’m not sure how to deal with this.

Thanks very much for your thoughts, I really appreciate it! Please let me know if you have any questions.

Code:

observed <- getY(full.model@data)

expected <- fitted(full.model)

# Individual deviance residuals

dev.i<-(sign(observed-expected)*sqrt(2*(observed*log(observed/expected)-(observed-expected))))

dev.i.sq<-dev.i^2

#sum of squared residuals

sumsqdev<-sum(dev.i.sq,na.rm=T)

df<-(length(observed)-sum(is.na(observed)))-k

dispersion<-sumsqdev/df

chat<-dispersion

# Deviance of full model

dev.part<-(observed*log(observed/expected))

sum.dev.part<-sum(dev.part,na.rm=T)

dev.sum<-2*sum.dev.part

#Quasi-AIC

q.aic.m1<-(-2*logLik(full.model)/chat)+(2*k*(k+1)/(n-k-1))

# Adjusted SE and P-values

se.full.model<-SE(full.model)

se.adj<-((sqrt(chat))*se.full.model)

z.adj<-(coef(full.model)/se.adj)

p.adj<-2*pnorm(-abs(z.adj))