1

I am running what I suppose is the same mixed-effect model with a negative binomial distribution (log link) in both lme4 and the glmmTMB package in R. Code shown below:

mNB<-glmer.nb(size~ scale(Group1) + (1|PairID), weights=w, verbose=T, data=imp.NB)

which gives the warning

boundary (singular) fit: see ?isSingular

and summary

    summary(mNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(0.0515)  ( log )
Formula: size ~ scale(Group1) + (1 | PairID)
   Data: imp.NB
Weights: w

      AIC       BIC    logLik  deviance  df.resid 
-451160.3 -451137.1  225584.1 -451168.3      2458 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
  17.57   21.73   79.29  218.52 2655.68 

Random effects:
 Groups Name        Variance Std.Dev.
 PairID (Intercept) 0        0       
Number of obs: 2462, groups:  PairID, 183

Fixed effects:
               Estimate Std. Error  z value Pr(>|z|)    
(Intercept)   -4.491636   0.007533 -596.276  < 2e-16 ***
scale(Group1) -0.031555   0.007361   -4.287 1.81e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
scale(Grp1) -0.023
convergence code: 0
boundary (singular) fit: see ?isSingular

Meanwhile this gives no complaints and summary

t2NB<-glmmTMB(size~scale(Group1) + (1|PairID), weights=w, family=nbinom2, verbose=T, data=imp.NB, ziformula=~0)         

summary(t2NB)
     Family: nbinom2  ( log )
    Formula:          size ~ scale(Group1) + (1 | PairID)
    Data: imp.NB
    Weights: w

         AIC      BIC   logLik deviance df.resid 
     86472.4  86495.6 -43232.2  86464.4     2458 

    Random effects:

    Conditional model:
     Groups Name        Variance Std.Dev.
     PairID (Intercept) 0.3382   0.5816  
    Number of obs: 2462, groups:  PairID, 183

    Overdispersion parameter for nbinom2 family (): 1.02 

    Conditional model:
                  Estimate Std. Error z value Pr(>|z|)    
    (Intercept)   1.983829   0.046852   42.34   <2e-16 ***
    scale(Group1) 0.007061   0.034005    0.21    0.836

I have reason to believe the weights argument causes the singularity as discussed in this post. But I want to know why one converges and not the other? Can I trust the results from glmmTMB?

EDIT

Note: running both models without the weights gives nearly identical results. While running the lme4 model with weights gives 0 variance for the random effect.

Nebulloyd
  • 263
  • 3
  • 10

1 Answers1

1

These are fairly surprising results, and I would want to do more investigation before accepting either set of results ("the man with two watches never knows what time it is"). Some general thoughts:

  • weights are often problematic because (1) less well-exercised and (2) sometimes interpreted differently by different packages. Have you fitted the model with both packages without weights included to see if you get similar answers? (What is the meaning of the weights in your model, i.e. what are you trying to account for by using weights?)
  • you could try with a third package, e.g. GLMMadaptive (i.e. mixed_model(): see ?negative.binomial, and be aware that loading GLMMadaptive after lme4 will break the CRAN version of lme4::glmer.nb [development version works] due to the masking of MASS::negative.binomial.
  • almost everything about the models is very different (e.g. negative binomial parameter of 0.05 (glmer.nb) vs 1.02 (glmmTMB), very different log-likelihoods ...
  • if taken at face value, the log-likelihood for glmer.nb (225584.1) is much higher than that for glmmTMB (-43232.2), which would say that the glmer.nb fit is much better, but I'm not 100% sure I trust that the weights have been appropriately accounted for in the log-likelihoods.

There's a small chance that there are multiple optima, as in this example.

  • you could try different starting values (especially, reciprocally starting the glmer.nb fit at the glmmTMB values and vice versa).
  • with a bit more work, you could examine the log-likelihood of the model across a range of all four parameters (two fixed effects, one random effect, and one overdispersion parameters) and plot the results
Ben Bolker
  • 43,543
  • 2
    Running each model without weights gives almost identical summary outputs. If there's a specific parameter you're interested in I can dig it up? This brings your comment about the treatment of weights to the forefront. The weights should compensate for imperfect detection during the survey, modeled as probabilities in Distance. They have been rescaled (1-10). 1 represents observations in habitat with the highest chance of detection, and 10 should give observations in habitat with the lowest chance of detection greater weight. In the case of this species, I think the values go from 3-7. – Nebulloyd Oct 27 '19 at 22:02