On 3/31/06, Marco Geraci <marcodoc75 at yahoo.com>
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 at 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
>