[R-sig-ME] Problem with computing gr and false convergence
Ben Bolker
bbolker at gmail.com
Wed Jun 15 23:07:30 CEST 2011
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On 06/15/2011 04:26 PM, Iker Vaquero Alba wrote:
>
> On 06/15/2011 01:14 PM, Iker Vaquero Alba wrote:
>
>> Thank you very much once again. The code for the graphs is really
>> useful. However, I have to be a bit annoying again. When I have
> this model:
>
>> > summary(s.g1.15)
>> Generalized linear mixed model fit by the Laplace approximation
>> Formula: nhatch ~ mod + briventral + brithr + year + (1 | site/pair)
>
> [snip]
>
>
>
>> and I try to simplify, again, the variable "briventral", I get the
>> same error message as usual:
>
>> > s.g1.21<-update(s.g1.15,~.-briventral)
>> Mensajes de aviso perdidos
>> 1: In mer_finalize(ans) :
>> Cholmod warning 'not positive definite' at
>> file:../Cholesky/t_cholmod_rowfac.c, line 432
>> 2: In mer_finalize(ans) : false convergence (8)
>
>> I cannot understand why, I mean, if the models with more parameters
>> in it worked, why this one, which is much less suspicious of being
>> overfitted, doesn't? Obviously it looks like there's a big problem
> with
>> the variable itself. Any ideas about what could it be, or posible
>> solutions?
>> Sorry again. Thank you very much
>
> I agree that this is a big pain.
> I was able to get it to work by dropping the 'site' random effect
> (which has zero variance anyway, so doesn't affect the results much).
>
> So, I guess "(1|site/pair)" means "site+site:pair". Is that right?
>
> I completely agree that from a distance none of this stuff seems to
> have much logic to it. I believe (but this belief is admittedly
> faith-based, and not evidence-based in this case) that *if* one were to
> spend a lot of time digging in the gory details of the shapes of the
> objective functions and the precise workings of the optimization
> algorithm that one could understand why this particular problem is being
> problematic in this particular way. It will be a *little* bit easier
> if/when lme4a has the capability to do likelihood profiles on all
> parameters of GLMMs (not yet). In the meantime, it's back to trial and
> error based on reasonable ways of simplifying the problem.
>
> I have a bigger question about your modeling strategy, though, which
> is: why are you trying to remove this parameter? In order to test its
> significance against the full model? In general, if you're going to
> follow the approach of estimating significance of individual
> predictors/parameters by LRT, you should *not* reduce the model from the
> fullest version that works first (i.e., don't throw out non-significant
> predictors before testing for significance among the remaining
> predictors).
>
> Actually, I normally test the significance of all the individual
> parameters and then drop the least significant one each time, but
> the piece of code I wrote in my post was only part of the whole
> process (the only one which was giving problems)
Fine, but I was questioning the stepwise procedure in general. I
think you should not drop non-significant predictors, due to the issues
that Frank Harrell (for example) eloquently and repeatedly brings up.
But I agree that this procedural point is orthogonal to your specific
technical question.
>
> Something like drop1(g1,test="Chisq") or
> drop1(g2,test="Chisq") [see below for def. of g2] or drop1(g2) [which
> gives AIC rather than Chisq/LRT significance values] is the right way to
> proceed if you want to test hypotheses about the predictors ...
>
> I´ve tried that and, again, I get an error message:
>
> > drop1(g2,test="Chisq")
> Error in x$terms : $ operator not defined for this S4 class
>
> I´ve tried to learn about this problem and, if I got it right,
> this error message appears when you try to use some feature that has
> not been implemented in lme4 but worked in S (and other R
> packages???). Am I right or am I doing something wrong?
I believe this functionality was enabled only in late April, so you
might want to try re-installing the R-forge version of lme4 ...
cheers
Ben
>
> Thank you so much for all the effort you've taken with my "case".
>
> Best wishes. Iker
>
>
> g2 <-
> glmer(nhatch~sex+mod+briventral+brithr+tlength+cond+year+(1|site:pair),
> data=X,family=poisson)
> s.g2.15<-update(g2,~.-sex-tlength-cond)
> update(s.g2.15,~.-briventral,verbose=TRUE)
>
> good luck
> Ben Bolker
>
>
>> Iker
>
>
> ------------------------------------------------------------------------
>> *De:* Ben Bolker <bbolker at gmail.com
> </mc/compose?to=bbolker at gmail.com>>
>> *Para:* Iker Vaquero Alba <karraspito at yahoo.es
> </mc/compose?to=karraspito at yahoo.es>>
>> *CC:* r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>
>> *Enviado:* mié,15 junio, 2011 16:14
>> *Asunto:* Re: [R-sig-ME] Problem with computing gr and false
> convergence
>
>> On 06/14/2011 06:45 PM, Iker Vaquero Alba wrote:
>
>>> Thank you so much for that, Ben.
>
>>> I´m trying to do a split plot simplification with that model,
> but I
>>> still get this message with briventral:
>
>>> s.g1.5<-update(g1,~.-briventral)
>>> Mensajes de aviso perdidos
>>> 1: In mer_finalize(ans) :
>>> Cholmod warning 'not positive definite' at
>>> file:../Cholesky/t_cholmod_rowfac.c, line 432
>>> 2: In mer_finalize(ans) : false convergence (8)
>
>
>>> Do you think I should try to centre and scale that variable as
> well?
>
>>> Thank you.
>
>>> Iker.
>
>> Oddly enough I am able to fit the model better when the continuous
>> variables are *not* centered and scaled. I can't explain this
>> particularly well, but weird things sometimes happen when you have a
>> model that is on the edge of being overfitted (which this one still is
>> -- 12 parameters plus two nested random effects for 94 data points
> in a
>> fairly unbalanced design). Some code is also included below to take a
>> quick look at the distribution of the data. There are some patterns,
>> but the continuous variables don't seem to be doing too much ...
>
>> X <- read.table("nfledge-nhatch insects 2009-2010.txt",header=TRUE)
>> X <-
> transform(X,site=factor(site),pair=factor(pair),year=factor(year))
>> X <- transform(X,ctlength=scale(tlength),
>> cbriventral=scale(briventral))
>> f1 <- nhatch~(sex+mod+briventral+brithr+tlength+cond)^2-sex:mod
>> ncol(model.matrix(f1,data=X)) ## 41 parameters!
>> f2 <- nhatch~sex+mod+briventral+brithr+tlength+cond+year
>> ncol(model.matrix(f2,data=X)) ## 12 parameters (still dicey)
>> with(X,table(site,pair,year))
>> library(lme4)
>> X <-
>
> na.omit(subset(X,select=c(site,pair,year,sex,mod,nhatch,briventral,brithr,tlength,cond)))
>> g1 <-
>
> glmer(nhatch~sex+mod+briventral+brithr+tlength+cond+year+(1|site/pair),
>> data=X,family=poisson)
>> s.g1.5<-update(g1,~.-briventral)
>
>> library(ggplot2)
>> X2 <- melt(X,id.var=1:6)
>> ggplot(X2,
>> aes(x=value,y=nhatch,colour=mod,shape=sex))+
>> geom_point()+facet_grid(year~variable,scale="free_x")
>> ggplot(X,
>
> aes(x=mod,y=nhatch,colour=mod))+geom_boxplot()+facet_grid(year~sex)
>
>> dd <- drop1(g1)
>
>
>
>>>
> ------------------------------------------------------------------------
>>> *De:* Ben Bolker <bbolker at gmail.com
> </mc/compose?to=bbolker at gmail.com> <mailto:bbolker at gmail.com
> </mc/compose?to=bbolker at gmail.com>>>
>>> *Para:* Iker Vaquero Alba <karraspito at yahoo.es
> </mc/compose?to=karraspito at yahoo.es>
>> <mailto:karraspito at yahoo.es </mc/compose?to=karraspito at yahoo.es>>>
>>> *CC:* r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>
>> <mailto:r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>>
>>> *Enviado:* mar,14 junio, 2011 23:02
>>> *Asunto:* Re: [R-sig-ME] Problem with computing gr and false
> convergence
>
>>> On 06/14/2011 04:12 PM, Iker Vaquero Alba wrote:
>>>> Thank you very much. Data is attached.
>
>>> Unfortunately, looking at your data makes it very clear that you
> will
>>> have a lot of trouble fitting this model: maybe this isn't the
> complete
>>> data set ... ?
>
>>> * there are only two years, which makes it nearly impossible to
> handle
>>> year as a random effect
>
>>> * you have a total of 94 observations in the data set, and your
> model
>>> involves 41 (!!) fixed effect parameter. There is just no way
> you can
>>> fit this many parameters (even neglecting the random effects).
>
>>> I had some success with this model by dropping the interactions and
>>> fitting only the main effects. The site-level variance is
> estimated as
>>> zero, but there is some pair*site variance (and a significant
> difference
>>> between years, at least as inferred from a Wald Z test)
>
>>> X <- read.table("nfledge-nhatch insects 2009-2010.txt",header=TRUE)
>>> X <-
> transform(X,site=factor(site),pair=factor(pair),year=factor(year))
>>> X <- transform(X,ctlength=scale(tlength,center=TRUE))
>>> f1 <- nhatch~(sex+mod+briventral+brithr+tlength+cond)^2-sex:mod
>>> ncol(model.matrix(f1,data=X)) ## 41 parameters!
>>> f2 <- nhatch~sex+mod+briventral+brithr+tlength+cond+year
>>> ncol(model.matrix(f2,data=X)) ## 12 parameters
>>> with(X,table(site,pair,year))
>>> library(lme4)
>>> g1 <-
>>>
> glmer(nhatch~sex+mod+briventral+brithr+ctlength+cond+year+(1|site/pair),
>>> data=X,family=poisson,na.action=na.omit)
>
>
>>>>
>>>> To center and scale continuous variables, I've tried
> standardizing
>>>> them, as someone suggested in some post, and I have succesfully done
>>>> before. But when fitting the model with standardized variables,
> I get:
>>>>
>>>> "Error in `contrasts<-`(`*tmp*`, value = "contr.treatment") :
>>>> contrasts can be applied only to factors with 2 or more levels"
>>>>
>>>> Any ideas? Thank you very much!
>>>>
>>>>
>>>>
> ------------------------------------------------------------------------
>>>> *De:* Ben Bolker <bbolker at gmail.com
> </mc/compose?to=bbolker at gmail.com> <mailto:bbolker at gmail.com
> </mc/compose?to=bbolker at gmail.com>>
>> <mailto:bbolker at gmail.com </mc/compose?to=bbolker at gmail.com>
> <mailto:bbolker at gmail.com </mc/compose?to=bbolker at gmail.com>>>>
>>>> *Para:* r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>
>> <mailto:r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>>
>>> <mailto:r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>
>> <mailto:r-sig-mixed-models at r-project.org
> </mc/compose?to=r-sig-mixed-models at r-project.org>>>
>>>> *Enviado:* mar,14 junio, 2011 21:09
>>>> *Asunto:* Re: [R-sig-ME] Problem with computing gr and false
> convergence
>>>>
>>>> On 06/14/2011 02:57 PM, Iker Vaquero Alba wrote:
>>>>>
>>>>> Dear R-users:
>>>>>
>>>>> I am fitting a model with quite many terms, and I'm having a lot of
>>>>> problems. Either I get:
>>>>>
>>>>>
>>>>> "In mer_finalize(ans) : gr cannot be computed at initial par (65)"
>>>>>
>>>>> or false convergence problems, when the model is a little bit
>>>>> simpler. I attach the most complex one with "verbose=TRUE" to
> see if
>>>>> you can help me detectt where the problem is:
>>>>>
>>>>>
>>>>>
>>>>
>
>
> hatchcoltailmodel1<-lmer(nhatch~sex+mod+briventral+brithr+tlength+cond+briventral:brithr+briventral:tlength+briventral:cond+briventral:sex+briventral:mod+
>>>>>brithr:tlength+brithr:cond+brithr:sex+brithr:mod+tlength:cond+tlength:sex+tlength:mod+
>>>>>
>>>>
>
>
> cond:sex+cond:mod+(1|site/pair)+(1|year),family=poisson,na.action=na.omit,verbose=TRUE)
>>>>>
>>>>> 0: nan: 1.06217 0.584705 0.261488 -15.3440 -0.757161
>>>>> 16.9953 3.85760 1.53215 1.30924 5.14768 -0.0259057 0.0675839
>>>>> 0.172804 9.56398 0.000488889 -0.000119170 0.0156963 0.00445958
>>>>> -0.109915 -0.0119783 -0.0113610 -0.00132653 0.00594852 -2.29204e-05
>>>>> -0.0678690 0.00288760 0.321050 -0.0107392 0.0195566 -0.00242138
>>>>> -0.0538382 -0.0866086 -0.00632440 -0.168793 -0.0123981 0.00379237
>>>>> -0.00313671 -0.0228579 0.504940 nan -0.722186 -1.31712
> -0.618912
>>>>> -1.25492 Mensajes de aviso perdidos In mer_finalize(ans) : gr
> cannot
>>>>> be computed at initial par (65)
>>>>>
>>>>>
>>>>
>>>> Very hard to say without a reproducible example. Can you post the
>>>> data somewhere?
>>>>
>>>> How big is your data set?
>>>>
>>>> If any of your variables are continuous, consider centering and
>>>> scaling them (e.g. using scale())
>>>>
>>>>
>>>> Other than that, I only have a couple of coding style suggestions.
>>>>
>>>> 1. I *think* but am not sure that your very long model above is
>>>> equivalent to
>>>>
>>>> (briventral+brithr+tlength+cond+sex+mod)^2-mod:sex
>>>>
>>>> (all main effects + all pairwise interactions except mod:sex)
> but you
>>>> should of course check that (and the order might not be the same
> as the
>>>> model you have above).
>>>>
>>>> 2. It is best to use the 'data' argument to specify a data frame
>>>> (rather than attach()ing or having the variables floating around
> in the
>>>> workspace)
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>
>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>>
>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>
>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>>>
>>>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>
>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>>
>>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>
>> <mailto:R-sig-mixed-models at r-project.org
> </mc/compose?to=R-sig-mixed-models at r-project.org>>>> mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iEYEARECAAYFAk35HxIACgkQc5UpGjwzenMsgQCgllj8hR358PTj/2bRU0twJNoo
E1EAoJnvgEEkZgMrG9UPnY+elYKYteRS
=akZ5
-----END PGP SIGNATURE-----
More information about the R-sig-mixed-models
mailing list