[R-sig-ME] problem with lme4

[Ben Bolker]

On 14-03-11 11:43 AM, Saussac Mathilde wrote:
> Dear all,
> 
>  
> 
> I am a French student in BioStatistic and I have a question about the
> package lme4.
> 
>  
> 
> I am using the glmer function but it doesn’t work … I have these errors :

  Can you say what version you are using?

  From your results below it looks very much like you have issues with
complete separation (i.e., some categories in your data have
all-success or all-failure outcomes).  In general a binomial model
with extreme parameters (absolute value > 10) suggests complete or
near-complete separation [on the probability scale a change of 10
log-odds units corresponds to a change from p=0.5 to p=0.99995, or
from 0.006 to 0.993 ...]  More recent versions of lme4 are somewhat
better at dealing with the complete-separation issue; these are
available from Github or from http://lme4.r-forge.r-project.org/repos,
but we are also hoping that a new version will arrive on CRAN within
the next few days (week?).

  More generally, adjusting tolerances is (as you suggest) a little
bit dangerous, as it may get you silly results.

  The most principled thing to do with complete-separation results is
to use some form of penalization or regularization or prior (closely
related concepts for ways of squashing extreme values in toward zero);
at present you can use the blme package or MCMCglmm to impose such
priors.

  Can you give a reproducible example
<http://tinyurl.com/reproducible-000>, or perhaps the results of
with(d,table(AP2,Motif)) ?  I would expect that some categories
of AP2 have all successes or none ...

> 
> */>model <- glmer(Motif ~ as.factor(AP2) + (1|ID), data=d,
> family=binomial(link="logit"))/*
> 
> /Erreur dans pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac,
> verbose) : /
> 
> /  Downdated VtV is not positive definite /
> 
> Or sometimes this one :
> 
> /Erreur : pwrssUpdate did not converge in 30 iterations/
> 
> I found on the web, that we can change the tolerance (Tolpwrss) and
> indeed it works.
> 
> But, for different values, my results are completely different (see the
> examples for tolPwrss = 1e-5 and 1e-2) :
> 
> */> model <- glmer(Motif ~ as.factor(AP2) + (1|ID), data=d,
> family=binomial(link="logit"),control= glmerControl(tolPwrss=1e-5,
> optimizer="bobyqa"))/*
> 
> /Generalized linear mixed model fit by maximum likelihood ['glmerMod']/
> 
> /Family: binomial ( logit )/
> 
> /Formula: Motif ~ as.factor(AP2) + (1 | ID) /
> 
> /   Data: d /
> 
> / /
> 
> /     AIC      BIC   logLik deviance /
> 
> / 78.9974  93.2092 -35.4987  70.9974 /
> 
> / /
> 
> /Random effects:/
> 
> /Groups Name        Variance Std.Dev./
> 
> /ID     (Intercept) 578.5    24.05   /
> 
> /Number of obs: 258, groups: ID, 218/
> 
> / /
> 
> /Fixed effects:/
> 
> /                 Estimate Std. Error z value Pr(>|z|)  /
> 
> /(Intercept)        10.205      6.149   1.660    0.097 ./
> 
> /as.factor(AP2)1    15.057  46157.700   0.000    1.000  /
> 
> /as.factor(AP2)2     9.356     10.219   0.916    0.360  /
> 
> /---/
> 
> /Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ //0.1 ‘ ’ 1 /
> 
> / /
> 
> /Correlation of Fixed Effects:/
> 
> /            (Intr) a.(AP2)1/
> 
> /as.fc(AP2)1  0.000         /
> 
> /as.fc(AP2)2 -0.303  0.000  /
> 
> / /
> 
> */> model <- glmer(Motif ~ as.factor(AP2) + (1|ID), data=d,
> family=binomial(link="logit"),control= glmerControl(tolPwrss=1e-2,
> optimizer="bobyqa"))/*
> 
> /Generalized linear mixed model fit by maximum likelihood ['glmerMod']/
> 
> /Family: binomial ( logit )/
> 
> /Formula: Motif ~ as.factor(AP2) + (1 | ID) /
> 
> /   Data: d /
> 
> / /
> 
> /     AIC      BIC   logLik deviance /
> 
> /112.9273 127.1392 -52.4637 104.9273 /
> 
> / /
> 
> /Random effects:/
> 
> /Groups Name        Variance Std.Dev./
> 
> /ID     (Intercept) 3.794    1.948   /
> 
> /Number of obs: 258, groups: ID, 218/
> 
> / /
> 
> /Fixed effects:/
> 
> /                Estimate Std. Error z value Pr(>|z|)    /
> 
> /(Intercept)       2.8039     0.4435   6.322 2.58e-10 ***/
> 
> /as.factor(AP2)1  10.5258   117.9609   0.089   0.9289    /
> 
> /as.factor(AP2)2   2.1162     1.0125   2.090   0.0366 *  /
> 
> /---/
> 
> /Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 /
> 
> / /
> 
> /Correlation of Fixed Effects:/
> 
> /            (Intr) a.(AP2)1/
> 
> /as.fc(AP2)1 -0.004         /
> 
> /as.fc(AP2)2 -0.415  0.002  /
> 
> / /
> 
> So I want to know, until which threshold of tolPwrss are results reliable ??
> 
> Is it correct to change this value for all my model ? (Sometimes it
> doesn’t work until I put tolPwrss=0.01 or 0.1 ..)
> 
>  
> 
> I hope you will understand what I explained and you will be able to help
> me ..
>  
> If there is few things I haven't been clear, don't hesitate to ask me
> again.
>  
> Thank you for your attention to these matters.
> 
>  
> 
> Best regards,
> 
>  
> 
> Mathilde
> 
>  
>