generalized linear mixed models: large differences when using glmmPQL or lmer with laplace approximation
Due to some travel and the need to attend to other projects, I haven't been keeping up as closely with this list as I normally do. Regarding the comparison between the PQL and Laplace methods for fitting generalized linear mixed models, I believe that the estimates produced by the Laplace method are more reliable than those from the PQL method. The objective function optimized by the Laplace method is a direct approxmation, and generally a very good approximation, to the log-likelihood for the model being fit. The PQL method is indirect (the "QL" part of the name stands for "quasi-likelihood") and, because it involves alternating conditional optimization, can alternate back-and-forth between two potential solutions, neither of which is optimal. (To be fair, such alternating occurs more frequently in the analogous method for nonlinear mixed-models, in which I was one of the co-conspirators, than in the PQL method for GLMMs.) It may be that the problem you are encountering has more to do with the use of the quasipoisson family than with the Laplace approximation. I am not sure that the derivation of the standard errors in lmer when using the quasipoisson family is correct, in part because I don't really understand the quasipoisson and quasibinomial families. As far as I know, they don't correspond to probability distributions so the theory is a bit iffy. Do you need to use the quasipoisson family or could you use the poisson family? Generally the motivation for the quasipoisson familiy is to accomodate overdispersion. Often in a generalized linear mixed model the problem is underdispersion rather than overdispersion. In one of Ben's replies in this thread he discusses the degrees of freedom attributed to certain t-statistics. Regular readers of this list are aware that degrees of freedom is one of my least favorite topics. If one has a reasonably large number of observations and a reasonably large number of groups then the issue is unimportant. (Uncertainty in degrees of freedom is important only when the value of the degrees of freedom is small. In fact, when I first started studying statistics we used the standard normal in place of the t-distribution whenever the degrees of freedom exceeded 30). Considering that the quasi-Poisson doesn't correspond to a probability distribution in the first place, (readers should feel free to correct me if I am wrong about this) I find the issue of the number of degrees of freedom that should be attributed to a distribution of a quantity calculated from a non-existent distribution to be somewhat off the point. I think the problem is more likely that the standard errors are not being calculated correctly. Is that what you concluded from your simulations, Ben? On Tue, Oct 7, 2008 at 8:21 AM, Martijn Vandegehuchte
<martijn.vandegehuchte at ugent.be> 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