[R] loglikelihood and lmer

[Douglas Bates]

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
>