[R-sig-ME] CHOLMOD error in lmer-specifying nested random effects

Fri Aug 7 20:31:11 CEST 2009

Hello,

a quick reply to get this started (you may get more expert feedback by other members of the list).

1. I'd expect sex to have an effect on nursing time (longer time for males?). Sex is definitively a fixed effect (see also the mailing list archives).

2. Specify a unique ID code for each cub and for each mother, then lmer will work out the nesting without the need to specify it explicitly (see posts by Doug Bates on this issue).

Thus I'd try something like:

ger<-lmer(Nurse.min~Mom.min+Mom.rank+Sex+(1|Cub)+(1|Mom),data=data,family=poisson,control=list(msVerbose=T))

If you don't have multiple cubs for each mother it might be that a model with only cubID might be more appropriate (but am guessing here). If the warning persists try also centrering/standardizing Mom.min.

Re the overdispersion issue have a look at the mailing list archive (e.g. postings by Ben Bolker).

HTH

Cheers,

Luca

----- Messaggio originale -----
Da: "Eli Swanson" <eliswanson at gmail.com>
A: r-sig-mixed-models at r-project.org
Inviato: Venerdì, 7 agosto 2009 18:17:06 GMT +01:00 Amsterdam/Berlino/Berna/Roma/Stoccolma/Vienna
Oggetto: [R-sig-ME] CHOLMOD error in lmer-specifying nested random effects

Hi all,

I'm working on an analyzing some data right now that's causing me some
difficulties.  I have 2 questions, and would really appreciate any
advice that people can offer.  I don't post data because there are
over 2000 cases.  sessionInfo() is posted at the end of my post.

The data was collected during focal animal surveys (FAS).  The
response variable is the number of minutes a cub was observed nursing
(Nurse.min).  The fixed predictors are 1) the number of minutes the
mother was present during the FAS (Mom.min), 2) the mother's social
rank (Mom.rank), and 3)sex of cub(Sex).  The random effects are 1)
Identity of mom, and 2) identity of cub, as there are a large number
of observations for every mother, and most of the cubs.  Cub is nested
within mom, and when I include sex in the model, my understanding is
that I have to nest it within cub.

Question the first:
I try to create the following model, leaving out sex for the moment,
because my understanding is that with only two cases, sex complicates
matters:

ger<-lmer(Nurse.min~Mom.min+Mom.rank+(1|Mom:Cub)+(1|Mom),data=data,family=poisson,control=list(msVerbose=T))

Is this model correctly specified?  I receive a CHOLMOD warning, as follows:

   11918.870: 0.289063 0.172748 -0.0419576 0.0575048 0.00979307
CHOLMOD warning: 7u‡e_
Error in mer_finalize(ans) :
  Cholmod error `not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 432

When I try to include Sex as a fixed effect, I'm even less sure I'm
specifying the model correctly but my model looks like this:

ger<-lmer(Nurse.min~Mom.min+Mom.rank+Sex+(1|Mom:Cub:Sex)+(1|Mom),data=data,family=poisson,control=list(msVerbose=T))

And my output:
  0:     11926.160: 0.289063 0.172748 0.172336 0.0572882 0.00311534
-0.244655 -0.683584
CHOLMOD warning: 7u‡e_
Error in mer_finalize(ans) :
  Cholmod error `not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 432

The model works if I specify it as follows, but I can't get it to work
with sex, and I think I'm specifying the random effects incorrectly
here:

ger<-lmer(Nurse.min~Mom.min+Mom.rank+(1|Mom:Cub),data=data,family=poisson)

summary(ger)
Generalized linear mixed model fit by the Laplace approximation
Formula: Nurse.min ~ Mom.min + Mom.rank + (1 | Mom:Cub)
   Data: data
   AIC   BIC logLik deviance
 11381 11404  -5687    11373
Random effects:
 Groups  Name        Variance Std.Dev.
 Mom:Cub (Intercept) 2.1313   1.4599
Number of obs: 2234, groups: Mom:Cub, 70

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.4669098  0.1962970   -7.47 7.84e-14 ***
Mom.min      0.0688668  0.0007541   91.32  < 2e-16 ***
Mom.rank     0.0676121  0.0069510    9.73  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr) Mom.mn
Mom.min  -0.134
Mom.rank -0.332  0.115

So, am I specifying these models correctly?  What could be the
problem?  I only have 1 observation for some of the cubs, and I
thought this might be the problem, but when I remove these cubs, the
error remains.

  I'm not sure if it makes a difference, but I do have slightly
overdispersed data (its 0-inflated i think, but not hugely so), hence
my second question:

When I use the quasipoisson family to estimate the model like so
(specifying the model the only way i could get it to work, even though
this may be incorrect)
:
 ger<-lmer(Nurse.min~Mom.min+Mom.rank+(1|Mom:Cub),data=data,family=quasipoisson,control=list(msVerbose=T))

Then, by my understanding, estimate the degree of overdispersion:

lme4:::sigma(ger)
 sigmaML
1.233926

How badly overdispersed is this?  I tried using MCMCglmm with a
"zipoisson", but its giving me strange results (in fact, estimates for
the fixed effects that appear to have an opposite sign of those I find
with lmer).  As I'm not an expert in Bayesian methods, I assume I'm
specifying something wrong there and would prefer to stick with lmer
for now if possible.  Can I just continue to use quasipoisson?  I know
that the standard errors are inflated, but if im not worried about
that does this mean that parameters are correctly estimated?

sessionInfo()
R version 2.9.0 (2009-04-17)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lme4_0.999375-28   Matrix_0.999375-24 lattice_0.17-22

loaded via a namespace (and not attached):
[1] grid_2.9.0

i would appreciate any help that people can offer,
Thank you very much,
Eli

-- 
Eli Swanson
Department of Zoology
Ecology, Evolutionary Biology, and Behavior Program
Michigan State University

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models