[R-sig-ME] problem with lme4

[Steve Walker]

Mathilde,

My first thought is complete/quasi-complete separation, but I can't be 
sure without looking at the data.  If I'm right, as the tolerance is 
turned down, glmer tries to reach larger and larger values of the 
as.factor(AP2)1 coefficient.  Once you get too large, numerical 
instabilities take over and you get error messages.

For example, when tolPwrss = 1e-5, the linear predictor for observations 
in the as.factor(AP2)1 category is about 25, which is an extremely large 
number on the logit scale (corresponding to a probability that is pretty 
darn close to one).

If I'm right, you might want to try Vince Dorie's blme package, which 
can put prior distributions on fixed effect coefficients to keep them 
from blowing up.

Cheers,
Steve

On 3/11/2014, 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 :
>> 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
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>