Problem with computing gr and false convergence

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)

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?

   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)

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