The not-necessarily-obvious thing is that none of the parameter estimates change when you go from poisson, binomial, etc. to their quasi- equivalents -- the variance on all points is inflated by the same amount, so the point estimates don't change. So you can go ahead and use the quasi- variant to get the estimated scale. I'm still struggling with exactly what "sigma" is in the quasi- case (as is Doug Bates) -- it seems NOT equal to the Pearson residuals^2/(resid df) ? It would be worth exploring this some more ... cheers Ben Bolker set.seed(1001) ntot <- 1000 x <- runif(ntot) f <- rep(1:20,each=50) reff <- rnorm(20,mean=0,sd=0.5) y <- rnbinom(ntot,mu=exp(2*x-1+reff[f]),size=1) library(lme4) m1 <- glmer(y~x+(1|f),family=poisson) m2 <- update(m1,family=quasipoisson) lme4:::sigma(m1) ## 1 lme4:::sigma(m2) ## 1.04 (??) ## residuals returned are Pearson residuals ## (uncorrected for overdispersion) myres <- (y-fitted(m1))/sqrt(fitted(m1)) all.equal(myres,residuals(m1)) ## not quite identical, but nearly all(abs(myres-residuals(m1))<1e-4) all.equal(residuals(m1),residuals(m2)) ## TRUE sum(residuals(m1)^2)/ntot ## NOT sigmaML m1 at deviance["sigmaML"] m2 at deviance["sigmaML"]
Amy Wade wrote:
Ben, When I fit a quasi- variant it gives what seems to be a sensible value for the sigma. The trouble is I would like to know if my models are overdispersed without the quasi-variant. I'm using data with poisson distribution from which I know what the estimated scale is and it still gives '1' when it should be 0.97. There was a discussion about this previously and it was suggested to use 'lme4:::sigma(model)', that just gives me the same result as using 'summary(model)@sigma'. Thanks very much for your help. Amy Should I have put this on the public thread? Not really sure how to do that! -----Original Message----- From: Ben Bolker [mailto:bolker at ufl.edu] Sent: 16 February 2009 14:11 To: Amy Wade Subject: Re: [R-sig-ME] No estimated scale value given What happens if you fit the same model with a quasi- variant and then try to access sigma? Ben Bolker Amy Wade wrote:
Hello, I'm trying to run generalized linear mixed models using the lmer() function in the lme4 package. The problem is that the output does not give me a value for the Estimated scale. The rest of the output is as it should be. This makes it very difficult to assess whether the model is overdispersed. My colleague tried the same code on his computer and it did give an estimated scale value. I tried unistalling R and reinstalling the latest edition (2.8.1) with the latest lme4, but this made no difference. I also tried 'summary(model)@sigma' which people suggested on the CRAN forums to extract the scale parameter. This always gave me an answer of '1' even when I used data where I knew this was not the case. When I load up the lme4 library it does warn me about several objects being masked
from 'stats' and 'base'.
Any idea why the output is not giving me the estimate scale? Or any idea how to extract it somehow? Many thanks, Amy [[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- Ben Bolker Associate professor, Biology Dep't, Univ. of Florida bolker at ufl.edu / www.zoology.ufl.edu/bolker GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc
Ben Bolker Associate professor, Biology Dep't, Univ. of Florida bolker at ufl.edu / www.zoology.ufl.edu/bolker GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc