negative binomial models using glmmTMB & glmer.nb

[Jesse Delia]

Dear list,

Hope all is well. Im running into issues trying to fit negative binomial mixed models.  I used two packages, neither are converging. I've googled a bunch and checked the list, cant seem to find a solution... I'm hoping someone might have ideas? Below is a reproducible example using my data. 

The model is relatively straight forward: comparing days until feeding for tadpoles between two treatments, while controlling for clutch of origin/sibling random effects. One thing -- there are decimals in the dataset because (0.5 days), however, i converted them to integers in the example below, but still ran into issues. 

Thanks for your time,

-Jesse

#####

library(glmmTMB); library(MASS); library(lme4)

data<-structure(list(clutch.ID = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 

3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 

8L, 8L, 9L, 9L, 9L, 10L, 10L, 10L, 11L, 11L, 11L, 12L, 12L, 12L, 

13L, 13L, 13L, 14L, 14L, 14L, 15L, 15L, 15L, 16L, 16L, 16L, 17L, 

17L, 17L, 18L, 18L, 18L, 19L, 19L, 19L, 20L, 20L, 20L, 21L, 21L, 

21L), .Label = c("CGH10", "CGH11", "CGH12", "CGH13", "CGH15", 

"CGH16", "CGH18", "CGH19", "CGH2", "CGH21", "CGH22", "CGH23", 

"CGH24", "CGH4", "CGH5", "CGH5.1", "CGH6", "CGH6.1", "CGH7", 

"CGH7.1", "CGH9"), class = "factor"), treatment = structure(c(1L, 

1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 

2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 

2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 

1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L), .Label = c("Early", 

"Late"), class = "factor"), feedinghatch = c(5.5, 5.5, 5.5, 6, 

6, 6, 6, 6, 6.5, 7, 7, 7, 0.5, 1, 1, 0.5, 1, 1, 2.5, 2.5, 2.5, 

2, 2, 2, 5, 5, 5, 1.5, 1.5, 1.5, 0.5, 0.5, 0.5, 1, 1, 1, 2, 2, 

2, 4.5, 4.5, 4.5, 6.5, 6.5, 6.5, 0.5, 0.5, 0.5, 6.5, 6.5, 6.5, 

0.5, 0.5, 0.5, 5, 5, 5, 0.5, 0.5, 0.5, 5, 5, 5), feedingdpo = c(13, 

13, 13, 13.5, 13.5, 13.5, 13.5, 13.5, 14, 14.5, 14.5, 14.5, 15, 

15.5, 15.5, 15, 15.5, 15.5, 17, 17, 17, 16.5, 16.5, 16.5, 12.5, 

12.5, 12.5, 16, 16, 16, 15, 15, 15, 15.5, 15.5, 15.5, 16.5, 16.5, 

16.5, 12, 12, 12, 14, 14, 14, 15, 15, 15, 14, 14, 14, 15, 15, 

15, 12.5, 12.5, 12.5, 15, 15, 15, 12.5, 12.5, 12.5)), .Names = c("clutch.ID", 

"treatment", "feedinghatch", "feedingdpo"), class = "data.frame", row.names = c(NA, -63L))

NB <-glmer.nb(feedingdpo ~ treatment + (1 | clutch.ID), data = data)

#Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GHrule(0L), compDev = compDev,  : (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

#There were 50 or more warnings (use warnings() to see the first 50)

warnings(NB)

#Warning messages:

#1: In (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L,  ... : non-integer x = 13.500000 Error in cat(list(...), file, sep, fill, labels, append) : argument 2 (type 'S4') cannot be handled by 'cat'

TMB<-glmmTMB(feedingdpo ~ treatment+ (1 | clutch.ID), family="nbinom2", data=data)

#Warning messages:

#1: In fitTMB(TMBStruc) : Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

#2: In fitTMB(TMBStruc) : Model convergence problem; false convergence (8). See vignette('troubleshooting')

summary(NB)

#Error in summary(NB) : object 'NB' not found

summary(TMB)

#Family: nbinom2  ( log )

#Formula:          feedingdpo ~ treatment + (1 | clutch.ID)

#Data: data

 #    AIC      BIC   logLik deviance df.resid 

  #    NA       NA       NA       NA       59 

#Random effects:

#Conditional model:

 #Groups    Name        Variance  Std.Dev.

 #clutch.ID (Intercept) 2.756e-11 5.25e-06

#Number of obs: 63, groups:  clutch.ID, 21

#Overdispersion parameter for nbinom2 family (): 4.56e+08 

#Conditional model:

  #            Estimate Std. Error z value Pr(>|z|)    

#(Intercept)    2.58148    0.05022   51.40   <2e-16 ***

#treatmentLate  0.16909    0.06677    2.53   0.0113 *  

#i checked the vignette, but the only solution i found was rescaling the predictor. however, im not sure if that makes sense considering i only have two treatments-- i tried it but got one of the same errors (but not both):

data$treat<-as.numeric(data$treatment)

data$scaletreatment<-scale(data$treat, center = TRUE, scale = TRUE)

TMBS<-glmmTMB(feedingdpo ~ scaletreatment+ (1 | clutch.ID), family="nbinom2", data=data)

#i tried an optimizer for the glmer.nb, but got the same error:

NB <-glmer.nb(feedingdpo ~ treatment + (1 | clutch.ID), data = data,  control=glmerControl( optimizer ="bobyqa",optCtrl=list(maxfun=1e5)) )

 

#i thought maybe the issue might be the decimal entries for 0.5 days, so i tired converting them to an integer as number of checks (2) per day. But that leads to other issues..

data$checksdpo<-data$feedingdpo*2 

NB2 <-glmer.nb(checksdpo ~ treatment + (1 | clutch.ID), data = data)

#Warning messages:

#1: In theta.ml(Y, mu, weights = object@resp$weights, limit = limit,  : iteration limit reached

#2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  : Model failed to converge with max|grad| = 0.00133527 (tol = 0.001, component 1)

#3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  : Model failed to converge with max|grad| = 0.00142708 (tol = 0.001, component 1)

#4: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  : Model failed to converge with max|grad| = 0.00217178 (tol = 0.001, component 1)

#5: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  : Model failed to converge with max|grad| = 0.00158022 (tol = 0.001, component 1)

TMB2<-glmmTMB(checksdpo ~ treatment+ (1 | clutch.ID), family="nbinom2", data=data)

#same warnings as other TMB model

[Justin Touchon]

Hi Jesse,

First off, thanks for using the STRI R Google group!  It only happens once a year or so, but it is very satisfying when it does happen.

I can get your example data to run in yet a different mixed model package, glmmADMB.  glmmADMB uses a totally different modeling algorithm, which is why it often works when lme4 does not.  It’s usually slower, but it is good to know about.  Note that it is more strict as far as needing integers for a negative binomial model, so this needs to be done with your “number of checks” variable.  Also note that last I checked, glmmADMB still has to be installed from R-Forge, so see here (http://glmmadmb.r-forge.r-project.org) for more info on installation.

Saludos,

Justin

> NB<-glmmadmb(checksdpo~treatment+(1|clutch.ID), family="nbinom", data=data)

> summary(NB)

Call:

glmmadmb(formula = checksdpo ~ treatment + (1 | clutch.ID), data = data, 

    family = "nbinom2")

AIC: 344.9 

Coefficients:

              Estimate Std. Error z value Pr(>|z|)    

(Intercept)     3.2746     0.0367   89.33  < 2e-16 ***

treatmentLate   0.1691     0.0489    3.46  0.00054 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of observations: total=63, clutch.ID=21 

Random effect variance(s):

Group=clutch.ID

             Variance    StdDev

(Intercept) 3.529e-07 0.0005941

Negative binomial dispersion parameter: 403.43 (std. err.: 0.0019028)

Log-likelihood: -168.451 

> NB.1<-glmmadmb(checksdpo~1+(1|clutch.ID), family="nbinom", data=data)

> anova(NB,NB.1)

Analysis of Deviance Table

Model 1: checksdpo ~ 1

Model 2: checksdpo ~ treatment

  NoPar  LogLik Df Deviance  Pr(>Chi)    

1     3 -174.47                          

2     4 -168.45  1   12.044 0.0005196 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Warning message:

In anova.glmmadmb(NB, NB.1) :

  rearranging models in order of increasing complexity

--
You received this message because you are subscribed to the Google Groups "STRI R users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stri-R-group...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.