generalized linear mixed models: large differences when using glmmPQL or lmer with laplace approximation

Tue, Oct 7, 2008 9:44 AM

Martijn,

Absolute truth is the realm of pure mathematics, philosophy, and religion.  R deals with statistics where we are more concerned with usefulness than absolute correctness.

One famous quote among statisticians is: "All models are wrong, some models are useful".

Given that, I would suggest looking at the quality of the predictions from your models.  Get predicted values from the 2 different fits and compare them to each other and the original data.  If one set of predictions is clearly superior, then you have your answer.  If the predictions are not that different, then that can tell you that possibly there are multiple models that are equally useful for this data.  Then you can look at the predictions from models fit with fewer predictors to see how that affects the predictions.

Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111

-----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 7:21 AM
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] generalized linear mixed models: large differences
when using 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),family=qua
sipoisson)

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,f
amily=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(l
oglength),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

generalized linear mixed models: large differences when using glmmPQL or lmer with laplace approximation

Thread (8 messages)