[R] Quasi-poisson glm and calculating a qAIC and qAICc...trying to modilfy Bolker et al. 2009 function to work for a glm model

[Gavin Simpson]

On Thu, 2011-01-27 at 08:20 -0500, Jason Nelson wrote:
> Sorry about re-posting this, it never went out to the mailing list when I
> posted this to r-help forum on Nabble and was pending for a few days, now
> that I am subscribe to the mailing list I hope that this goes out:
> 
> I've been a viewer of this forum for a while and it has helped out a lot,
> but this is my first time posting something.
> 
> I am running glm models for richness and abundances.  For example, my beetle
> richness is overdispersed:
> 
> > qcc.overdispersion.test(beetle.richness)
> 
> Overdispersion test Obs.Var/Theor.Var Statistic   p-value
>        poisson data          2.628131  23.65318 0.0048847
> 
> So, I am running a simple glm with my distribution as quasipoisson
> 
> > glm.richness1<-glm(beetle.richness~log.area, family = quasipoisson)
> 
> 
> Now I want to calculate a qAIC and qAICc.  I was trying to modify the
> equation that I found in Bolker et al 2009 supplemental material:
> 
> QAICc <- function(mod, scale, QAICc=TRUE){
>         LL <- logLik(mod)

You are out of luck there then; logLik is not defined for the quasi
families.

>         ll <- as.numeric(LL)
>         df <- attr(LL, "df")
>         n <- length(mod$y)   #used $ to replace @ to make a S3 object
>         if(QAICc)
>                 qaic = as.numeric( -2*ll/scale + 2*df +
> 2*df*(df+1)/(n-df-1))
>         else qaic =as.numeric( -2*ll/scale + 2*df)
>         qaic
> }
> 
> The only problem is that I have no idea how to scale it.  In Bolker at al.
> 2009 it is scaled to "phi":
> 
> phi = lme4:::sigma(model)

phi, IIRC, is the dispersion parameter. You can get the estimated value
from `summary(model)$dispersion` where model is the result of a call to
glm().

> But I am not running a mixed model and I cannot run the qAICc function
> without scaling it.  I am comparing models to each other trying to find the
> best model for both landscape land use land cover data and patch variables.
>  How would I set the scale if I run this function?
> 
> QAICc(glm.richness1, scale = ?)
> 
> Should I set the scale to the square root of the deviance?  phi =
> sqrt(glm.richness1$deviance)
> 
> Your help is much appreciated.

Instead of resorting to these ad-hoc approaches to correct for
overdispersion, you would be better off fitting models that model the
overdispersion. Try a negative binomial (glm.nb() in MASS) or the
zeroinfl() and hurdle() functions in the pscl package. Those all have
proper log likelihoods and you can compute/extract AIC simply using
them.

> Regards,
> Jason

HTH

G

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%