0

I am trying to fit zero inflated poisson random effects model with count data using glmmTMB. However, I received some errors and warnings.

head(data)

  count time study    
1     0  259     1 
2     0  199     1 
3     0  571     1 
4     0  927     1 
5     7  254     1 
6     0  877     1 

str(data)

 'data.frame':  959 obs. of  3 variables:
 $ count : int  0 0 0 0 7 0 0 0 0 0 ...
 $ time  : int  259 199 571 927 254 877 555 158 1014 705 ...
 $ study : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...

time is an offset term, and study is a random effect I have to add.

Error 1

If I add a random effect, then the following random-effect models report the same error

f1<-glmmTMB(count~1+offset(time)+1|study, data = data, family = poisson, ziformula = ~1|study)

f2<-glmmTMB(count~1+offset(time)+1|study, data = subset(data, study=="1"), family = poisson, ziformula = ~1|study)

NA/NaN function evaluationError in nlminb(start = par, objective = fn, gradient = gr, control = control$optCtrl) : 
 gradient function must return a numeric vector of length 4

Error 2

I thought I need to try a simpler model, and thus I dropped the random-effect. Unfortunately, I received a different error.

f3<-glmmTMB(count~1+offset(time), data = data, family = poisson, ziformula = ~1)

f4<-glmmTMB(count~1+offset(time), data = subset(data, study=="1"), family = poisson, ziformula = ~1)

NA/NaN function evaluationError in nlminb(start = par, objective = fn, gradient = gr, control = control$optCtrl) : 
  NA/NaN gradient evaluation

Warning 1

OK. I wonder how I can get an output without any error. So, I dropped the offset term instead of the random-effect, and thus obtained a warning instead of an error.

f5<-glmmTMB(count~1+1|study, data = data, family = poisson, ziformula = ~1|study)

f6<-glmmTMB(count~1+1|study, data = subset(data, study=="1"), family = poisson, ziformula = ~1|study)

Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729

"Good" fit with neither error nor warning

After a bunch of refits, I find that if I drop at least one of two random-effect terms and the offset term, models will report neither an error or a warning.

f7<-glmmTMB(count~1+1|study, data = data, family = poisson, ziformula = ~1)

f8<-glmmTMB(count~1, data = data, family = poisson, ziformula = ~1|study)

Question

What issues caused these errors and warining? If possible, how could I fit a ZIP mixed effect model and receive neither error nor warning?

Example

I have a reproducible example which produced the warning 1

count.ex<-rpois(500, 0.2)
study.ex<-as.factor(sample(1:5, 500, replace = TRUE))
time.ex<-rexp(500, 150)

fit.ex<-glmmTMB(count.ex~1+offset(time.ex)+1|study.ex, family = poisson, ziformula = ~1|study.ex)

Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
C.C.
  • 55
  • 9

1 Answers1

2

This is rather difficult without a reproducible example!

  • There is one obvious problem with your model: you've specified offset(time) in your model, you almost certainly want to use offset(log(time)). Your time values are large (the values I can see in your str() range from 500-1000; the offset enters the predicted value in the model as exp(offset), which is either ridiculously large (if offset < log(.Machine$double.xmax) == 709.78) or infinite.

  • You should be careful to protect random-effect terms with parentheses (e.g. you use ~1 + 1|study rather than ~1 + (1|study); I don't think it causes problems in the examples you show here, but in general the pipe (|) operator has low precedence, which means that e.g. ~1+a+b|study is equivalent to ~(1+a+b)|study ...

If I use

 ~count.ex~1+offset(log(time.ex))+(1|study.ex)

in your example, I don't get any warnings.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453