On 3/31/06, Marco Geraci <[EMAIL PROTECTED]> wrote: > Dear R users, > I am estimating Poisson mixed models using glmmPQL > (MASS) and lmer (lme4). We know that glmmPQL do not > provide the correct loglikelihood for such models (it > gives the loglike of a 'pseudo' or working linear > mixed model). I would like to know how the loglike is > calculated by lmer.
Good point. The person who wrote the documentation (i.e. I) should have mentioned that. With lmer, one can fit a generalized linear mixed model using PQL or by optimizing the Laplacian approximation to the deviance directly. Even when you use PQL, which is the default method, the Laplace approximation is evaluated at the converged parameter estimates. This is the value of the loglikelihood that is reported. I am reconsidering having PQL as the default method for fitting generalized linear mixed models with lmer. I would appreciate it if you and others who do fit such models could tell me if it is noticeably slower or less reliable to do direct optimization of the Laplace approximaiton. That is, if you use the combination of optional arguments method = "Laplace" and control = list(usePQL = FALSE) does the fit take substantially longer? On your example I get > system.time(fit1.lmer <- lmer(z ~ X-1 + (1|id), family=poisson)) [1] 0.48 0.01 0.54 0.00 0.00 > system.time(fit2.lmer <- lmer(z ~ X-1 + (1|id), family=poisson, method = > "Laplace", control = list(usePQL = FALSE))) [1] 0.61 0.00 0.62 0.00 0.00 > fit1.lmer Generalized linear mixed model fit using PQL Formula: z ~ X - 1 + (1 | id) Family: poisson(log link) AIC BIC logLik deviance 922.2406 934.8844 -458.1203 916.2406 Random effects: Groups Name Variance Std.Dev. id (Intercept) 0.82446 0.908 number of obs: 500, groups: id, 100 Estimated scale (compare to 1) 1.021129 Fixed effects: Estimate Std. Error z value Pr(>|z|) X1 1.003639 0.098373 10.202 < 2.2e-16 *** X2 2.075037 0.052099 39.829 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: X1 X2 -0.337 > fit2.lmer Generalized linear mixed model fit using Laplace Formula: z ~ X - 1 + (1 | id) Family: poisson(log link) AIC BIC logLik deviance 921.958 934.6019 -457.979 915.958 Random effects: Groups Name Variance Std.Dev. id (Intercept) 0.8895 0.94313 number of obs: 500, groups: id, 100 Estimated scale (compare to 1) 0.04472136 Fixed effects: Estimate Std. Error z value Pr(>|z|) X1 0.985721 0.101645 9.698 < 2.2e-16 *** X2 2.075060 0.052114 39.818 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: X1 X2 -0.326 The only unexpected part of that output is the estimated scale which is wrong (well, it is never calculated in this case and consequently should not be displayed). > A minor question is: why do glmmPQL and lmer give > different degrees-of-freedom for the same estimated > model? Does glmmPQL consider the scale parameter 'phi' > as a degree of freedom? I believe that is the reason. The class of an object fit by glmmPQL is > class(fm1) [1] "glmmPQL" "lme" and the only method I could find for the "glmmPQL" class is > methods(class = "glmmPQL") [1] predict.glmmPQL* Non-visible functions are asterisked Generic functions other than predict will choose the method for the lme class of linear mixed effects models (or, if there isn't an lme method, the default method). The lme methods defined in the nlme package are appropriate for linear mixed effects models (which is what Jose and I wrote them for) and typically are not appropriate for a generalized linear mixed model. > A toy example > > set.seed(100) > m <- 5 > n <- 100 > N <- n*m > X <- cbind(1,runif(N)) > Z <- kronecker(diag(n),rep(1,m)) > z <- rpois(N, exp(X%*%matrix(c(1,2)) + > Z%*%matrix(rnorm(n)))) > id <- rep(1:n,each=m) > fit.glmm <- glmmPQL(z ~ X-1, random = ~1|id, > family="poisson",verbose=F) > fit.lmer <- lmer(z ~ X-1 + (1|id), > family="poisson",verbose=F) > > > logLik(fit.glmm) > 'log Lik.' -386.4373 (df=4) > > logLik(fit.lmer) > 'log Lik.' -458.1203 (df=3) > > Thanks, > Marco > > > > > sessionInfo() > R version 2.2.1, 2005-12-20, i386-pc-mingw32 > > attached base packages: > [1] "methods" "stats" "graphics" "grDevices" > "utils" > [6] "datasets" "base" > > other attached packages: > mvtnorm SemiPar cluster lme4 lattice > Matrix > "0.7-2" "1.0-1" "1.10.4" "0.995-2" "0.12-11" > "0.995-5" > nlme MASS > "3.1-68.1" "7.2-24" > > > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html