R help - Jul 2013 - Differences between glmmPQL and lmer and AIC calculation

Dear R Community,



I?m relatively new in the field of R and I hope someone of you can
help me to solve my nerv-racking problem.



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). 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.

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)


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

                               Max. No. of iterations reached without
convergence

Warnmeldungen:

                               1: In model.response(mf, "numeric") :

                               using type="numeric" with a factor
response will be ignored

                               2: In glmmML.fit(X, Y, weights,
cluster.weights, start.coef, start.sigma,  :
                3: In glmmML(y ~ activity,  :

                               'vmmin' did not converge. Increase
'maxit'??



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.



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



 3) Is there any other option (library) to analyze my data including an AIC?



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).



Many thanks in advance!


All the best,

Tonio




*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

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