Gavin,
Thanks for your help, but I did find in Booker et al. 2009 supplemental
material where they did use logLik for quasi- distributed data.
Even though it was for mixed models effect. In section 6 in the
supplemental material, this is the model:
mp1 <- lmer(total.fruits ~ nutrient * amd + rack + status + (amd * nutrient
| popu) + (amd * nutrient | gen), data = dat.tf, fammily = poisson
mq1 <- update(mp1, family = quasipoisson)
(phi = lme4:::sigma(mq1))
Then they build the function that I have in my original message.
Then they compute the qAICc:
QAICc(mq1, scale = phi)
I was just following there example and try to fit my data in for a general
linear model instead of a mixed model.
You thoughts were very much appreciated!
Regards,
Jason
On Fri, Jan 28, 2011 at 5:43 AM, Gavin Simpson <gavin.simp...@ucl.ac.uk>wrote:
> 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/<http://www.ucl.ac.uk/%7Eucfagls/>
> UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
>
--
Jason M. Nelson
Master Candidate
Department of Zoology
Miami University
PSN 167F (Lab): 513.529.3391
PSN 149 (office)
Cell: 616.901.5923
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
