[R-sig-ME] Passing lme4 optimizer arguments to pbmodcomp and afex/mixed
Henrik Singmann
henrik.singmann at psychologie.uni-freiburg.de
Wed Sep 25 19:45:17 CEST 2013
Hello again,
I can confirm Tobias' observation that simple changes in the values of the predictors lead to the pwrssUpdate error message. I can demonstrate it with data from a recent but unrelated Cross-Validated question (that is desperately awaiting answers): http://stats.stackexchange.com/q/70783/442
require(lme4)
options(contrasts=c('contr.sum', 'contr.poly'))
dat <- read.table("http://pastebin.com/raw.php?i=vRy66Bif")
dat$V1 <- factor(dat$V1)
m1 <- glmer(true ~ distance*consequent*direction*dist + (consequent*direction*dist|V1), dat, family = binomial)
# runs but takes quite some time, so better stop it.
# transform centered consequent variable by squaring it (retaining the sign):
dat$consequent2 <- with(dat, ifelse(consequent < 0, -consequent^2, consequent^2))
m2 <- glmer(true ~ distance*consequent2*direction*dist + (consequent2*direction*dist|V1), dat, family = binomial)
# produces the error:
# Error: pwrssUpdate did not converge in 30 iterations
I hope this helps.
Best,
Henrik
PS: My session Info:
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
[5] LC_TIME=German_Germany.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_1.1-0 Matrix_1.0-12 lattice_0.20-15
loaded via a namespace (and not attached):
[1] grid_3.0.1 MASS_7.3-26 minqa_1.2.1 nlme_3.1-109 Rcpp_0.10.3
[6] splines_3.0.1
Am 25.09.2013 18:31, schrieb Tobias Heed:
> Hi,
>
> I get the same error about pwrssUpdate in a binomial GLMM in lme4 1.1.0 (downloaded about a week ago). I don't know whether my error is related to the ones reported here. Maybe it points into some direction about where the error might come from.
> I am using several categorical variables as fixed effects. One of my variables is a time difference between two stimuli (-110, -50, 50, 110 ms). To reflect the values, I used a contrast with these values, which resulted in the pwrssUpdate error. I then divided the values by 100 (all other factors unchanged) and ran the same glmer command. Then, the GLMM works normal.
> If needed, I can provide the data.
>
> Best, Tobias
>
> --
> --------------------------------------------------------------------------------------------------------------
> Tobias Heed, PhD
> Biological Psychology and Neuropsychology | University of Hamburg
> Von-Melle-Park 11, Room 206 | D-20146 Hamburg, Germany
> Phone: (49) 40 - 42838 5831 | Fax: (49) 40 - 42838 6591
> tobias.heed at uni-hamburg.de | Website | Google Scholar | ResearcherID
> --------------------------------------------------------------------------------------------------------------
>
>> On 25.09.2013, at 17:48, Tom Wenseleers <Tom.Wenseleers at bio.kuleuven.be> wrote:
>>
>> Hi Ben & Henrik,
>> Many thanks for the feedback!
>> (@Henrik: using your lme4 version with 10000 pwrssUpdate steps also didn't help unfortunately)
>>
>> In case you would have the time to look into it, here is my example:
>>
>> #library(pbkrtest) # latest dev version from http://people.math.aau.dk/~sorenh/software/pbkrtest/devel/
>> library(devtools)
>> install_github("lme4",user="lme4")
>> install.packages("afex",repos="http://r-forge.r-project.org",type="source")
>> library(afex)
>> data=read.csv("http://www.kuleuven.be/bio/ento/temp/alates.csv",header=T,fill=T)
>> summary(data)
>> fit=glmer(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|REPLICATE),family=poisson,data=data)
>> summary(fit) # this works OK
>> # Generalized linear mixed model fit by maximum likelihood ['glmerMod']
>> # Family: poisson ( log )
>> # Formula: ALATES ~ ANT_TENDED * MELEZITOSE + NR_START + (1 | CLONES) + (1 | PLANT) + (1 | REPLICATE)
>> # Data: data_alates
>> #
>> # AIC BIC logLik deviance
>> # 741.1707 761.6854 -362.5853 725.1707
>> #
>> # Random effects:
>> # Groups Name Variance Std.Dev.
>> # PLANT (Intercept) 0.7423 0.8616
>> # REPLICATE (Intercept) 0.2388 0.4886
>> # CLONES (Intercept) 0.6739 0.8209
>> # Number of obs: 96, groups: PLANT, 96; REPLICATE, 12; CLONES, 8
>> #
>> # Fixed effects:
>> # Estimate Std. Error z value Pr(>|z|)
>> # (Intercept) 0.69149 1.61200 0.429 0.668
>> # ANT_TENDED1 0.15985 0.09792 1.632 0.103
>> # MELEZITOSE1 -0.28012 0.30861 -0.908 0.364
>> # NR_START 0.15771 0.16280 0.969 0.333
>> # ANT_TENDED1:MELEZITOSE1 0.08052 0.09706 0.830 0.407
>> #
>> # Correlation of Fixed Effects:
>> # (Intr) ANT_TENDED1 MELEZI NR_STA
>> # ANT_TENDED1 -0.134
>> # MELEZITOSE1 0.057 -0.010
>> # NR_START -0.978 0.136 -0.060
>> # ANT_TENDED1: 0.028 0.023 -0.005 -0.029
>>
>> # LRTs also work:
>> mixed(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|REPLICATE),family=poisson,type=3,method="LRT",data=data_alates)
>> #Fitting 6 (g)lmer() models:[......]
>> # Effect df.large df.small chisq df p
>> #1 (Intercept) 8 7 0.18 1 0.6744
>> #2 ANT_TENDED 8 7 2.59 1 0.1074
>> #3 MELEZITOSE 8 7 0.78 1 0.3783
>> #4 NR_START 8 7 0.89 1 0.3444
>> #5 ANT_TENDED:MELEZITOSE 8 7 0.68 1 0.4109
>>
>> But if I try parametric bootstrapping it throws convergence errors:
>> cl <- makeCluster(rep("localhost", 4))
>> mixed(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|REPLICATE),family=poisson,method="PB",args.test = list(nsim = 10,cl=cl),data=data_alates)
>> # Error in checkForRemoteErrors(lapply(cl, recvResult)) :
>> # 4 nodes produced errors; first error: pwrssUpdate did not converge in 10000 iterations
>>
>> Same if one tries this from within lme4:
>> mm1=glmer(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|REPLICATE),family=poisson,data=data)
>> mm <- model.matrix(mm1)
>> mm1 <- glmer(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|REPLICATE),family=poisson,data=data)
>> form <- reformulate(c("-1",colnames(mm)[-1],"(1|CLONES)","(1|PLANT)","(1|REPLICATE)"),response="ALATES")
>> newdat <- data.frame(list(model.matrix(mm1)[,-1],getME(mm1,"flist"),model.frame(mm1)["ALATES"]),check.names=FALSE)
>> mm2 <- update(mm1,formula=form,data=newdat)
>> confint(mm2,method="boot",nsim=20)
>> #Computing bootstrap confidence intervals ...
>> # 2.5 % 97.5 %
>> #sd_(Intercept)|PLANT 0.69800722 1.0135147
>> #sd_(Intercept)|REPLICATE 0.19014163 0.6788588
>> #sd_(Intercept)|CLONES 0.56855996 1.0145768
>> #ANT_TENDED1 0.01637193 0.3186532
>> #MELEZITOSE1 -0.81867836 0.1685135
>> #NR_START 0.13921792 0.2738027
>> #ANT_TENDED1:MELEZITOSE1 -0.02224318 0.1300941
>> There were 27 warnings (use warnings() to see them)
>>> warnings()
>> Warning messages:
>> 1: In optwrap(object at optinfo$optimizer, ff, x0, lower = lower, ... :
>> convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations exceeded
>> 2: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 3: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 4: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 5: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 6: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 7: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 8: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 9: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 10: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 11: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 12: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 13: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 14: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 15: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 16: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 17: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 18: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 19: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, ... :
>> Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431
>> 20: In bootMer(object, bootFun, nsim = nsim, ...) :
>> some bootstrap runs failed (10/20)
>> 21: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>> 22: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>> 23: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>> 24: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>> 25: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>> 26: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>> 27: In norm.inter(t, alpha) : extreme order statistics used as endpoints
>>
>> I presume it has to do with some simulated datasets not allowing parameters to be accurately estimated. I wonder if it wouldn't be possible perhaps to let PBmodcomp just ignore the ones for which no convergence could be achieved (or else to conservatively assume that the simulated LRT value in those cases was more extreme than the observed ones - which is probably better and less likely to lead to strong biases - at least you wouldn't have strong evidence that in those cases the LRT value for the simulated dataset was smaller than the observed one), and throw a warning in which it would be reported in how many cases no convergence could be achieved (if this number of small, the result might still be more or less usable). Or is there something that could be done on the lme4 side with respect to optimizer settings to get over this? Also, is there perhaps some provisions to calculate p-values for all factors from the lme4 likelihood profiles?
>>
>> Cheers,
>> Tom
>>
>>
>> -----Original Message-----
>> From: Henrik Singmann [mailto:henrik.singmann at psychologie.uni-freiburg.de]
>> Sent: 25 September 2013 15:36
>> To: Ben Bolker
>> Cc: Tom Wenseleers; SÞren HÞjsgaard; r-sig-mixed-models at r-project.org
>> Subject: Re: Passing lme4 optimizer arguments to pbmodcomp and afex/mixed
>>
>> I am also currently working on a binomial GLMM on which I encounter this problem regularly (depending on which model I fit). But as my data is somewhat large (12000 observations) and fitting takes a while I wouldn't consider it a nice demonstration of the issue.
>>
>> However, I did try to increase the number of steps in pwrssUpdate by altering the source of lme4 and, as predicted by Ben, this does not help at all (it just reaches the new boundary).
>>
>> But, if nevertheless someone is interested to try it out for him or herself, you can install a version of lme4 1.1.0 that is identical to the development version as of yesterday that only stops after 10000 pwrssUpdate steps from my github:
>> library("devtools"); install_github("lme4",user="singmann")
>>
>> Best,
>> Henrik
>>
>> -----Original Message-----
>> From: Ben Bolker [mailto:bbolker at gmail.com]
>> Sent: 25 September 2013 15:21
>> To: Tom Wenseleers; SÞren HÞjsgaard; henrik.singmann at psychologie.uni-freiburg.de; r-sig-mixed-models at r-project.org
>> Subject: Re: Passing lme4 optimizer arguments to pbmodcomp and afex/mixed
>>
>>
>> At the moment it's not possible to increase the number of iterations; that's on the list to allow, but my guess is that that's not going to help. Unfortunately, our test cases for lme4 did not include a lot of binomial GLMMs with complete separation and/or extreme values of the beta parameters (say |beta|>5). I suspect that what's going on here is similar to the examples listing in these issues ..
>>
>> https://github.com/lme4/lme4/issues/134
>> https://github.com/lme4/lme4/issues/124
>>
>> Is it in fact the case that if you fit this model with lme4.0 (or an older version of lme4) you get large values of beta?
>>
>> Looks like we're going to have to put higher priority on fixing this ...
>>
>> Ben Bolker
>>
>>
>>> On 13-09-25 09:13 AM, Tom Wenseleers wrote:
>>> Dear all, Just a quick question re pbmodcomp (I'm using version
>>> 0.4-0 of pbkrtest, and latest devel version of lme4). I am calling
>>> pbmodcomp via afex, but I am getting the dreaded error "3 nodes
>>> produced errors; first error: pwrssUpdate did not converge in 30
>>> iterations" (my call is
>>> mixed(ALATES~ANT_TENDED*MELEZITOSE+NR_START+(1|CLONES)+(1|PLANT)+(1|RE
>>> PLICATE),type=3,family=poisson,method="PB",args.test
>>> = list(nsim = 10,cl=cl),data=data_alates) ).
>>>
>>> I would now like to try passing some optimizer arguments to the lmer
>>> calls. To do this I first called options(glmerControl=list(optimizer
>>> = c("bobyqa", "bobyqa"),tolPwrss = 1e-03
>>> ,optCtrl=list(maxfun=1000))) But then I still get " pwrssUpdate did
>>> not converge in 30 iterations". Does one of you happen to know how I
>>> could perhaps change the nr of maximum iterations to >30? And also if
>>> it is possible to pass such arguments in the call to mixed and/or
>>> pbmodcomp?
>>>
>>> Cheers, Tom
>>>
>>>
>>> ______________________________________________________________________
>>> _________________
>>>
>>> Prof. Tom Wenseleers P Lab. of Socioecology and Social
>>> Evolution Dept. of Biology Zoological Institute K.U.Leuven
>>> Naamsestraat 59, box 2466 B-3000 Leuven Belgium +32 (0)16 32 39 64 /
>>> +32 (0)472 40 45 96 //tom.wenseleers at bio.kuleuven.be
>>> http://bio.kuleuven.be/ento/wenseleers/twenseleers.htm
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> [[alternative HTML version deleted]]
>
>
>
--
Dipl. Psych. Henrik Singmann
PhD Student
Albert-Ludwigs-Universität Freiburg, Germany
http://www.psychologie.uni-freiburg.de/Members/singmann
More information about the R-sig-mixed-models
mailing list