Calculating deviance residuals and dispersion for a count model

[kirkby]

Hello,

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))

[ledzep]

My models don’t seem to have great fit, and I’d like to be able to consider and account for overdispersion.

Did you formally test for over dispersion first?
A rule of thumb is there is an overdispersion if your residual mean deviance are much larger than 1 [because of the chi-squared proporty that expected value of a chi-squared variable is its degree of freedom].

You have used the most complex model that you can fit to the data to estimate your dispersion parameter. So, your approach to estimate the dispersion parameter to scale the standard error seems correct.

What software you using? It looks like R to me. There is already an option in R to fit a quasi-Poisson model.

model1<-glm(formula,family=quasipoisson)

[kirkby]

Thanks for your response. I'm using a program within R called "unmarked". This program let's me run a hierarchical model linking fish abundance and detection probability, the second of which I couldn't account for using a standard quasi-/poisson glm for abundance. The program uses a poisson regression, but there's no built in method for a quasi-poisson, so I've been attempting to make those adjustments by hand. I don't have enough stats experience to be entirely comfortable with my own methods, so I wanted to get some feedback and/or corrections.

Thanks again for reading and responding!

Kristen