[R-sig-ME] generalized linear mixed models: large differences whenusing glmmPQL or lmer with laplace approximation
Doran, Harold
HDoran at air.org
Tue Oct 7 16:17:12 CEST 2008
Do you have access to Generalized, Linear and Mixed Models (McCullough
and Searle)? There is a nice little discussion comparing PQL and laplace
and discussing why PQL may result in estimates that are biased and
inconsistent.
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf
> Of Martijn Vandegehuchte
> Sent: Tuesday, October 07, 2008 9:29 AM
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] generalized linear mixed models: large
> differences whenusing glmmPQL or lmer with laplace approximation
>
> Dear list,
>
> First of all, I am a mere ecologist, trying to get the truth
> out of his data, and not a statistician, so forgive me my
> lack of statistical background and possible conceptual
> misunderstandings.
>
> I am currently comparing generalized linear mixed models in
> glmmPQL and lmer, with a quasipoisson family, and have found
> out that parameter estimates are quite different for both
> methods. I read some of the discussions on the R-forum and it
> seems that the Laplace approximation used in the current
> version of lmer is generally preferred to the PQL method. I
> am an ex-SAS user, and in proc glimmix in SAS the default is
> PQL, and the estimates and p-values are almost exact the same
> as with glmmPQL in R. But lmer gives quite different results,
> and now I am wondering what would be the best option for me.
>
> First of all, parameter estimates of a same model can be
> somewhat different in lmer or glmmPQL. Second of all, in
> lmer, I only get t-values but no associated p-values
> (apparently they are omitted because of the uncertainty about
> the df). But if I compare the t-values generated by glmmPQL
> with those of a same model in lmer, the differences are
> substantial. My dataset consists of 120 observations, so
> basically you could guess the order of magnitude of the
> p-values in lmer based on the t-value and a "large" df.
>
> First example:
> In lmer:
>
> >
> model<-lmer(schirufu~diameter+leafvit+densroot+cover+nemcm+(1|site),fa
> > mily=quasipoisson)
> > summary (model)
> Generalized linear mixed model fit by the Laplace approximation
> Formula: schirufu ~ diameter + leafvit + densroot + cover +
> nemcm + (1 | site)
> AIC BIC logLik deviance
> 2045 2068 -1015 2029
> Random effects:
> Groups Name Variance Std.Dev.
> site (Intercept) 12.700 3.5638
> Residual 15.182 3.8964
> Number of obs: 120, groups: site, 6
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 1.31017 1.47249 0.890
> diameter -0.24799 0.29180 -0.850
> leafvit 1.29007 0.21041 6.131
> densroot 0.31024 0.04939 6.281
> cover -0.24544 0.22179 -1.107
> nemcm 0.24817 0.12028 2.063
>
> Correlation of Fixed Effects:
> (Intr) diamtr leafvt densrt cover
> diameter 0.031
> leafvit -0.083 0.321
> densroot 0.011 -0.017 -0.202
> cover 0.021 -0.448 0.016 0.214
> nemcm -0.014 0.114 0.114 0.310 -0.017
> >
>
> Although no p-values are given, it suggests that fixed
> effects leafvit, densroot and nemcm would be significant.
> In glmmPQL:
>
> >
> model<-glmmPQL(schirufu~diameter+leafvit+densroot+cover+nemcm,random=~
> > 1|site,family=quasipoisson)
> iteration 1
> iteration 2
> iteration 3
> iteration 4
> iteration 5
> > summary(model)
> Linear mixed-effects model fit by maximum likelihood
> Data: NULL
> AIC BIC logLik
> NA NA NA
>
> Random effects:
> Formula: ~1 | site
> (Intercept) Residual
> StdDev: 0.7864989 4.63591
>
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Fixed effects: schirufu ~ diameter + leafvit + densroot +
> cover + nemcm
> Value Std.Error DF t-value p-value
> (Intercept) 1.4486735 0.4174843 109 3.470007 0.0007
> diameter -0.2600504 0.3477017 109 -0.747913 0.4561
> leafvit 1.2236406 0.2489291 109 4.915619 0.0000
> densroot 0.3236446 0.0596342 109 5.427164 0.0000
> cover -0.2523163 0.2698555 109 -0.935005 0.3519
> nemcm 0.2336305 0.1451751 109 1.609301 0.1104
> Correlation:
> (Intr) diamtr leafvt densrt cover
> diameter 0.130
> leafvit -0.335 0.313
> densroot 0.027 -0.022 -0.203
> cover 0.090 -0.463 0.015 0.214
> nemcm -0.056 0.097 0.107 0.301 -0.014
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -2.4956188 -0.4154369 -0.1333850 0.1724601 4.7355928
>
> Number of Observations: 120
> Number of Groups: 6
> >
>
> Note the difference in parameter estimates. Also, the fixed
> effect nemcm now is not significant any more.
>
> Second example,now with an offset:
> In lmer:
>
> >
> model<-lmer(nemcm~diameter+leafvit+densroot+rootvit+cover+schirufu+(1|
> > site), offset= loglength, family=quasipoisson) summary (model)
> Generalized linear mixed model fit by the Laplace approximation
> Formula: nemcm ~ diameter + leafvit + densroot + rootvit +
> cover + schirufu + (1 | site)
> AIC BIC logLik deviance
> 1593 1618 -787.4 1575
> Random effects:
> Groups Name Variance Std.Dev.
> site (Intercept) 21.522 4.6392
> Residual 173.888 13.1867
> Number of obs: 120, groups: site, 6
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 0.06733 1.92761 0.0349
> diameter 0.14665 0.60693 0.2416
> leafvit -0.19902 0.48802 -0.4078
> densroot -0.49178 0.64221 -0.7658
> rootvit 0.37699 0.46810 0.8054
> cover -0.23545 0.57896 -0.4067
> schirufu 0.23226 0.46866 0.4956
>
> Correlation of Fixed Effects:
> (Intr) diamtr leafvt densrt rootvt cover
> diameter -0.016
> leafvit 0.015 0.396
> densroot 0.055 -0.233 -0.291
> rootvit -0.038 -0.251 -0.629 0.277
> cover 0.024 -0.796 -0.133 0.253 0.117
> schirufu -0.032 0.137 -0.029 -0.505 -0.078 -0.121
> >
>
> This suggests no significant effects at all.
> In glmmPQL:
>
> >
> model<-glmmPQL(nemcm~diameter+leafvit+densroot+rootvit+cover+schirufu+
> > offset(loglength),random=~1|site, family=quasipoisson)
> iteration 1
> iteration 2
> iteration 3
> > summary (model)
> Linear mixed-effects model fit by maximum likelihood
> Data: NULL
> AIC BIC logLik
> NA NA NA
>
> Random effects:
> Formula: ~1 | site
> (Intercept) Residual
> StdDev: 0.2684477 4.507758
>
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Fixed effects: nemcm ~ diameter + leafvit + densroot +
> rootvit + cover + schirufu + offset(loglength)
> Value Std.Error DF t-value p-value
> (Intercept) 0.1131898 0.1656949 108 0.6831220 0.4960
> diameter 0.1225231 0.1976568 108 0.6198779 0.5366
> leafvit -0.2191361 0.1697784 108 -1.2907181 0.1996
> densroot -0.4733839 0.2221562 108 -2.1308604 0.0354
> rootvit 0.3858120 0.1615706 108 2.3878846 0.0187
> cover -0.2075038 0.1922054 108 -1.0795940 0.2827
> schirufu 0.2028444 0.1633954 108 1.2414323 0.2171
> Correlation:
> (Intr) diamtr leafvt densrt rootvt cover
> diameter -0.050
> leafvit 0.077 0.360
> densroot 0.217 -0.168 -0.262
> rootvit -0.163 -0.202 -0.632 0.257
> cover 0.084 -0.772 -0.098 0.200 0.073
> schirufu -0.103 0.099 -0.050 -0.483 -0.068 -0.075
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -1.1146287 -0.5208003 -0.1927005 0.2462878 7.9755368
>
> Number of Observations: 120
> Number of Groups: 6
> >
>
> Again some differences in parameter estimates, but now the
> two fixed effects densroot and rootvit turn out to be significant.
> So my questions are:
> - what would you recommend me to use? lmer or glmmPQL
> (laplace approximation or penalized quasi-likelihood)?
> - if lmer is the better option, is there a way to get a
> reliable p-value for the fixed effects?
> I have experienced that deleting a term and comparing models
> using anova() always overestimates the significance of that
> term, probably because the quasipoisson correction for
> overdispersion is not taken into account.
>
> Thank you very much beforehand,
>
> Martijn.
>
> --
> Martijn Vandegehuchte
> Ghent University
> Department Biology
> Terrestrial Ecology Unit
> K.L.Ledeganckstraat 35
> B-9000 Ghent
> telephone: +32 (0)9/264 50 84
> e-mail: martijn.vandegehuchte at ugent.be
>
> website TEREC: www.ecology.ugent.be/terec
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list