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=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@ugent.be 

website TEREC: www.ecology.ugent.be/terec

	[[alternative HTML version deleted]]


