[R-sig-ME] quasi-binomial family in lme4

Jarrod Hadfield j.hadfield at ed.ac.uk
Wed Nov 10 10:53:32 CET 2010


Hi,

It seems worrying to me that the Herd:Id estimate and the Id estimate  
are exactly the same. I get a (slightly) different estimate with  
version 0.999375-35 on Linux, but again the estimates are identical.  
MCMCglmm fitting the same model gives a very different answer.

Cheers,

Jarrod



> summary(gm2)
Generalized linear mixed model fit by the Laplace approximation
Formula: incidence ~ period + (1 | Id/herd)
    Data: cbpp
    AIC   BIC logLik deviance
  102.2 114.4 -45.11    90.21
Random effects:
  Groups  Name        Variance Std.Dev.
  herd:Id (Intercept) 0.29654  0.54456
  Id      (Intercept) 0.29654  0.54456
Number of obs: 56, groups: herd:Id, 56; Id, 56

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   1.1129     0.2477   4.492 7.04e-06 ***
period2      -1.1969     0.4185  -2.860 0.004233 **
period3      -1.4223     0.4381  -3.246 0.001169 **
period4      -2.0049     0.5292  -3.788 0.000152 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Correlation of Fixed Effects:
         (Intr) perid2 perid3
period2 -0.592
period3 -0.565  0.335
period4 -0.468  0.277  0.265

Using a parameter expanded prior in order to improve mixing as the  
herd variance approaches zero:

prior<-list(R=list(V=1, nu=0), G=list(G1=list(V=1, nu=1, alpha.mu=0,  
alpha.V=1000)))

mcmc2<-MCMCglmm(incidence~period, random=~herd, family="poisson", data=cbpp)

> summary(mcmc2)

  Iterations = 12991
  Thinning interval  = 3001
  Sample size  = 1000

  DIC: 175.9590

  G-structure:  ~herd

      post.mean  l-95% CI u-95% CI eff.samp
herd   0.01433 1.067e-16   0.0713    61.58

  R-structure:  ~units

       post.mean l-95% CI u-95% CI eff.samp
units    0.7347   0.1564    1.421    222.9

  Location effects: incidence ~ period

             post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)    1.0749   0.4584   1.5987    844.5  0.002 **
period2       -1.2054  -2.0862  -0.3556    771.8  0.004 **
period3       -1.4568  -2.4362  -0.4839    674.8  0.004 **
period4       -2.0584  -3.1851  -0.9430    446.5 <0.001 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1


Quoting Gosselin Frederic <frederic.gosselin at cemagref.fr>:

> Dear Colleague,
>
> good question:
>
> the commands are under R2.5.1:
> ************************************************************************
> library(lme4)
> cbppbis<-cbind.data.frame(cbpp,Id=as.factor(1:dim(cbpp)[1]))
> gm1 <- lmer(incidence ~ period + (1 | herd), family = quasipoisson,  
> data = cbpp)
> summary(gm1)
>
> gm2 <- lmer(incidence ~ period + (1 | Id/herd), family = poisson,  
> data = cbppbis)
> summary(gm2)
> **************************************************************************
> (fope it is declared in the good order for the random effects in gm2)
>
> The results do show a slight discrepancy between both methods:
>
> *************************************************************************
>> summary(gm1)
> Generalized linear mixed model fit using Laplace
> Formula: incidence ~ period + (1 | herd)
>    Data: cbpp
>  Family: quasipoisson(log link)
>    AIC   BIC logLik deviance
>  112.2 122.3 -51.11    102.2
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  herd     (Intercept) 0.35085  0.59233
>  Residual             1.40470  1.18520
> number of obs: 56, groups: herd, 15
>
> Fixed effects:
>             Estimate Std. Error t value
> (Intercept)   1.2812     0.2200   5.824
> period2      -1.1240     0.3315  -3.391
> period3      -1.3203     0.3579  -3.689
> period4      -1.9477     0.4808  -4.051
>
> Correlation of Fixed Effects:
>         (Intr) perid2 perid3
> period2 -0.339
> period3 -0.314  0.219
> period4 -0.233  0.163  0.151
>
>
>
>
>> summary(gm2)
> Generalized linear mixed model fit using Laplace
> Formula: incidence ~ period + (1 | Id/herd)
>    Data: cbppbis
>  Family: poisson(log link)
>    AIC   BIC logLik deviance
>  102.2 114.4 -45.11    90.21
> Random effects:
>  Groups  Name        Variance Std.Dev.
>  herd:Id (Intercept) 0.29608  0.54413
>  Id      (Intercept) 0.29608  0.54413
> number of obs: 56, groups: herd:Id, 56; Id, 56
>
> Estimated scale (compare to  1 )  0.9249959
>
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept)   1.1149     0.2476   4.503 6.69e-06 ***
> period2      -1.2013     0.4184  -2.871 0.004089 **
> period3      -1.4224     0.4378  -3.249 0.001159 **
> period4      -2.0089     0.5294  -3.795 0.000148 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Correlation of Fixed Effects:
>         (Intr) perid2 perid3
> period2 -0.592
> period3 -0.565  0.335
> period4 -0.468  0.277  0.264
> ************************************************************************
>
> The estimates and the standard errors are not exactly the same -  
> which might not be unlogical given that the relationship between  
> variance and mean is not the same in both models. They however are  
> not very far one from the other.
>
> Of course, this needs further investigation.
>
> I remember of a paper that motivated the use of the quasi-poisson  
> method based on empirical relationships between the residual  
> variance and the mean:
>
> Ver Hoef J.M. et Boveng P.L., 2007, Quasi-poisson vs. negative  
> binomial regression: How should we model overdispersed count data?,  
> Ecology, 88, 11, p. 2766-2772.
>
> Sincerely,
>
> Frédéric
>
> -----Message d'origine-----
> De : John Maindonald [mailto:john.maindonald at anu.edu.au]
> Envoyé : mercredi 10 novembre 2010 09:51
> À : Gosselin Frederic
> Cc : r-sig-mixed-models at r-project.org; tiflo at csli.stanford.edu
> Objet : Re: [R-sig-ME] quasi-binomial family in lme4
>
> I wonder if you have compared the results that you quote with the  
> result you get with observation level random effects in a poisson  
> model.
>
> As I see it, use of observation level random effects should, unless  
> there is evidence that a multiplicative effect on the scale of the  
> response is a better fit, replace use of the quasi- models in glm()  
> as well as in generalised linear mixed models.
>
> John Maindonald             email: john.maindonald at anu.edu.au
> phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
> Centre for Mathematics & Its Applications, Room 1194, John Dedman  
> Mathematical Sciences Building (Building 27) Australian National  
> University, Canberra ACT 0200.
> http://www.maths.anu.edu.au/~johnm
>
> On 10/11/2010, at 6:39 PM, Gosselin Frederic wrote:
>
>> Hi Florian,
>>
>> a different perpsective on the quasi-likelihood debate - that comes  
>> out sporadically on this list:
>>
>> (i) I globally agree with the previous repliers that a fully
>> probabilistic solution looks better - at least aesthetically - than a
>> quasi-likelihood;
>>
>> (ii) however, as I have already mentioned on the list (cf. below),  
>> earlier versions of lme4 give much more sensible results than the  
>> latest versions:
>> http://markmail.org/message/s4abxhhdacqjkunm
>>
>> This is why in the following papers:
>> Elek Z., Dauffy-Richard & Gosselin F., 2010, Carabid species  
>> responses to hybrid poplar plantation in floodplains in France,  
>> Forest Ecology and Management, 260, 9, p. 1446-1455.
>>
>> and
>>
>> Vuidot A., Paillet Y., Archaux F. & Gosselin F. (In Press)  
>> Influence of tree characteristics and forest management on tree  
>> microhabitats in France, Biological Conservation.
>>
>> we used version the R version 2.5.1 and the associated lme4 version  
>> (here with quasi-poisson, not quasi-bionomial).
>>
>> Hope this helps.
>>
>> Sincerely,
>>
>> Frédéric Gosselin
>> Engineer & Researcher (PhD) in Forest Ecology Cemagref Domaine des
>> Barres F-45290 Nogent sur Vernisson France
>>
>> http://www.cemagref.fr/les-contacts/les-pages-personnelles-professionn
>> elles/gosselin-frederic/english-short-scientific-cv
>>
>>
>> 	[[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.




More information about the R-sig-mixed-models mailing list