[R-sig-ME] No estimated scale value given

[Douglas Bates]

On Tue, Feb 17, 2009 at 2:23 PM, Ben Bolker <bolker at ufl.edu> wrote:
>  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 (??)

The first one is 1 by definition.  The second is the penalized
weighted residual sum of squares at the parameter estimates divided by
the number of observations.  You did not include the penalty term in
your calculations.

As we have discussed, it is not clear if this is a meaningful quantity
in the quasi-Poisson case.

> ## 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>