[R-sig-ME] generalized linear mixed models: large differences when using glmmPQL or lmer with laplace approximation
Douglas Bates
bates at stat.wisc.edu
Tue Oct 7 23:18:42 CEST 2008
On Tue, Oct 7, 2008 at 4:07 PM, Ken Beath <ken at kjbeath.com.au> wrote:
> There is a large difference between the estimated std of the random effect,
> usually a sign that the glmmPQL approximation isn't working.
Or that there is a mistake in the calculation of the standard errors
for the random effects, which is more likely in this case.
The actual optimization is with respect to the relative standard
deviation of the random effects (relative to the scale parameter in
the conditional standard deviation of the response). For the Poisson
family or the binomial family that scale parameter is fixed at 1 (you
could also consider the situation to be that there isn't a scale
parameter in those cases). For the quasipoisson and quasibinomial
families you maybe estimate a value there or maybe not. I don't know.
I believe Ben's simulations showed that I was doing the wrong thing
there.
> I would try using adaptive Gauss-Hermite. Set the nAGQ parameter to
> increasing values until the results stabilise.
>
> Ken
>
> On 08/10/2008, at 12:21 AM, Martijn Vandegehuchte wrote:
>
>> 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 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
>>
>
> _______________________________________________
> 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