[R] Differences between glmmPQL and lmer and AIC calculation

[Ben Bolker]

Tonio Pieterek <t.pieterek <at> googlemail.com> writes:

>

[snip]

This is really more appropriate for r-sig-mixed-models at r-project.org :
please refer any followup questions there.

> For my Master thesis I collected some behavioral data of fish using
> acoustic telemetry. The aim of the study is to compare two different
> groups of fish (coded as 0 and 1 which should be the dependent
> variable) based on their swimming activity, habitat choice, etc.
> (independent variables). 

  I don't quite understand this part.  You're trying to figure out
whether a particular observation falls into one category or
another (0/1)?  Do individual fish (id in the formula below)
always fall into one category or the other?

> Each fish has several observations over time
> (repeated measurements) which I included as random factor in my models
> using library glmmPQL (package MASS). Because I have a binary data
> structure, I am using generalized linear mixed models.

  For what it's worth, penalized quasi-likelihood (used by glmmPQL)
is generally considered to be a bit dicey with binary responses
(see e.g. Bolker et al 2008 TREE paper).

> Using library glmmPQL the results reflect my descriptive analyses and
> the results are sound. However, we also want to rank several candidate
> models using AIC. And this is where the problems start. Because
> glmmPQL does not provide AIC values or comparable measures, I also
> tried to calculate the same models using function lmer. Against
> expectations, I got completely different results from these two
> libraries (glmmPQL = highly significant; lmer = far away from being
> significant with p = 0.9xx).
> 
> I used the following codes:
> 
> cal1=glmmPQL(y ~ activity, random=~1|id, data=data, family=binomial,
> na.action=na.omit)
> 
> > WORKS FINE
> 
> cal1 = lmer(y ∼ activity + (1 | id ), family = binomial, data=data,
> na.action=na.omit)
> 
> > PRODUCED misleading and totally different results compared to glmmPQL
>  (e.g. sometimes error message
> occurs: In mer_finalize(ans) : false convergence (8); even 
> for very simple models)

  Do you possibly have complete separation, i.e. some individuals
with all-zero responses?

  Have you tried the development version of lme4?

> 
> A glmmML did not work since we got the following failure message, for
> which we were not able to find out the reason and therefore could not
> go on with this model:
> 
> “[glmmml] fail = 1

 [snip]
 
> The questions are:
> 
> 1) Why did glmmPQL and lmer produce completely different results and
> how can I solve this problem? Following Zuur et al. 2009* the models
> should provide very similar results, but they didn`t.

  I strongly suspect that there's something wrong with your setup.
In particular, if the response variable cal1 (0/1) only varies at
the level of individuals (id), and not within id, then you should
probably just calculate the mean activity per individual id and
run a simple logistic regression.  My guess would be that glmmPQL
may have papered over some cracks for you ...

> 2) Can I calculate AIC values (or something comparable) 
> using library glmmPQL?

  No, not without a great deal of difficulty.
> 
>  3) Is there any other option (library) to analyze my data including an AIC?

  You can use JAGS/WinBUGS with data cloning (the dclone package), or
glmmADMB.  

> 
> If something remained unclear or if you have any question about
> details, please let me know.
> 
> I would really appreciate any kind of help referring to my problem(s).
> 
> 
> *Alain F. Zuur,  Elena N. Ieno,  Neil J. Walker, Anatoly A. Saveliev,
> Graham M. Smith. (2009). Mixed Effects Models and Extensions in
> Ecology with R. Springer Science+Business Media, New York, USA.
> 
> ISSN 1431-8776
> ISBN 978-0-387-87457-9
> DOI 10.1007/978-0-387-87458-6


   I suspect