[R-sig-ME] problem with lme4

[Saussac Mathilde]

Thank you very much for your answers ! 

I was using the version 2.15.3 of R and 1.0-4 of lme4. 

You are right, there a complete separation in the case I told you: 
> with(d,table(AP2,Motif))
      Motif
AP2     0       1
  0         14    84
  1         0       44
  2         3      113 

I downloaded the latest versions of R and lme4 and I don't have mistakes in this case anymore ! 

But for others variables, like this case, I still have the mistake  : 
> modelM31R <- glmer(Motif~as.factor(recours) + (1|ID), data=d31, family=binomial(link="logit"))
Error: pwrssUpdate did not converge in 30 iterations
> with(d31,table(Motif,recours))
        recours
Motif     1        2     3     4
    0          7        8     1     1
    1          185   41   11   4
Or I have some mistakes when I want to use the function interaction in the model. 

As you suggest, I am going try to use the blme package or MCMCglmm.

Thank you very much for your help and your advices,  
Mathilde

-----Message d'origine-----
De : Ben Bolker [mailto:bbolker at gmail.com] 
Envoyé : mardi 11 mars 2014 17:03
À : Saussac Mathilde; r-sig-mixed-models at r-project.org
Cc : Pereira Bruno
Objet : Re: problem with lme4

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