[R-sig-ME] glmmADMB fails to fit poisson data

[Vladimir Trifonov]

Hi,

Why is the following  “trivial” model fit failing? After all, we are fitting a Poisson GLM on a Poisson data.

--------
> glmmadmb(rpois(100, lambda=8000)~1, family="poisson")
 Hessian is 0 in row 1
 This means that the derivative if probably identically 0  for this parameter
Error in matrix inverse -- matrix singular in inv(dmatrix)
Estimated covariance matrix may not be positive definite
 1e+20
Estimated covariance matrix may not be positive definite
 1e+20

GLMM's in R powered by AD Model Builder:

  Family: poisson
  link = log

Fixed effects:
  Log-likelihood: -2302.59
  AIC: 4607.18
  Formula: rpois(100, lambda = 8000) ~ 1
(Intercept)
          0

Number of observations: total=100
--------

I tried varying the number of samples (100 above), the mean (8000 above) and this fails. Here is another attempt (with a different kind of fail)

--------
> glmmadmb(rpois(100, lambda=50)~1, family="poisson")
Parameters were estimated, but standard errors were not: the most likely problem is that the curvature at MLE was zero or negative
Error in glmmadmb(rpois(100, lambda = 50) ~ 1, family = "poisson") :
  The function maximizer failed (couldn't find parameter file) Troubleshooting steps include (1) run with 'save.dir' set and inspect output files; (2) change run parameters: see '?admbControl';(3) re-run with debug=TRUE for more information on failure mode
In addition: Warning message:
running command './glmmadmb -maxfn 500 -maxph 5 -noinit -shess' had status 1
--------

After some trial and error, the only thing I found that “works” (with a warning about the covariance matrix) consistently is to set a “mysterious” poiss_prob_bound = FALSE

--------
> glmmadmb(rpois(100, lambda=50)~1, family="poisson", admb.opts=admbControl(poiss_prob_bound=FALSE))
Estimated covariance matrix may not be positive definite
 0.020141
Estimated covariance matrix may not be positive definite
 0.020141

GLMM's in R powered by AD Model Builder:

  Family: poisson
  link = log

Fixed effects:
  Log-likelihood: -327.749
  AIC: 657.498
  Formula: rpois(100, lambda = 50) ~ 1
(Intercept)
   3.904998

Number of observations: total=100
--------

While this is a “solution”, I would expect this kind of fit to work “off-the-box” (it certainly does with the standard glm function). Can someone comment on this.  I am not an expert on GLMs and needed something that would fit ZINB + random effects (admittedly the deep end of GLMs) and was excited to find glmmADMB, but using it has proven rather frustrating because it often fails (with a dramatic stop that does not reach the top level and so cannot be caught with tryCatch, nonetheless) and after I found the above example, I am not so sure it is always for a good reason.

Thanks,
Vladimir



Session:

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] C

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

other attached packages:
[1] glmmADMB_0.8.3.3 MASS_7.3-31

loaded via a namespace (and not attached):
 [1] Matrix_1.1-3    R2admb_0.7.13   Rcpp_0.12.3     coda_0.17-1
 [5] grid_3.1.0      lattice_0.20-29 magrittr_1.5    nlme_3.1-117
 [9] plyr_1.8.1      stringi_0.5-5   stringr_1.0.0   tools_3.1.0

> glmmADMB:::get_bin_version()
ADMB Program: ./glmmadmbADMB-11.5 safe libraries compiled with GNU C++ 4.9.0 (64bit)Copyright (c) 2008-2015 ADMB Foundation and Regents of the University of CaliforniaBuild date: Jun  9 2016