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 BIClogLik deviance
922.2406 934.8844 -458.1203 916.2406
Random effects:
Groups NameVariance 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 NameVariance 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:
mvtnormSemiParcluster lme4lattice
Matrix
0.7-21.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
