atahk/pscl

Zeroinfl may crash while computing SE.logtheta

Closed this issue · 4 comments

dstoeckel commented

Hi, first of all thank you for the pscl package. I am using it directly and via mpath.

However, I have a numerical problem that pops up from time to time in the current CRAN version (1.5.2).

zeroinfl calls optim to minimize the log-likelihood of the zero-inflated model. It then uses the inverse of the inferred Hessian (Fisher Information matrix) to compute the standard error of theta. If the Hessian is singular, this crashes. However, the inferred Hessian can simply become singular due to numerical issues in the optimization procedure (method="BFGS"), which crashes the solve call:

Error in solve.default(as.matrix(fit$hessian)) :
  system is computationally singular: reciprocal condition number = 3.08817e-32

The problem disappears if

  • The data is modified enough
  • The method is switched to
    • L-BFGS-B: Fast, but leads to warnings due to reltol being defined
    • CG: Significantly slower
  • EM=TRUE is used (Much slower)

The following code reproduces the problem:

y <- c(8, 73, 14, 0, 0, 40, 2, 0, 0, 0, 11, 0, 0, 0, 0, 12, 2, 0, 0, 7, 0, 78, 0, 0, 20, 0, 11, 0, 0, 0, 0, 8, 0, 0, 0, 0, 4, 0,
0, 0, 0, 3, 0, 72, 50, 3, 0, 0, 0, 0, 9, 0, 0, 0, 0, 13, 0, 1, 2, 0, 0, 0, 4, 0, 0, 10, 30, 0, 0, 5, 0, 9, 0, 8, 0, 0, 0, 1,
0, 18, 0, 7, 40, 0, 0, 0, 11, 0, 7, 9, 3, 40, 10, 0, 4, 35, 11, 4, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 33, 10, 0, 3, 0, 0, 35, 0, 0,
0, 0, 5, 0, 0, 0, 0, 3, 0, 0, 16, 0, 0, 9, 15, 6, 7, 0, 0, 30, 0, 9, 0, 0, 0, 30, 10, 0, 0, 0, 0, 0, 0, 9, 12, 5, 31, 0, 0,
0, 2, 0, 0, 7, 20, 28, 4, 0, 0, 0, 0, 0, 0, 0, 0, 64, 0, 7, 0, 0, 0, 5, 0, 0, 0, 30, 0, 0, 0, 5, 0, 0, 0, 4, 0, 0, 0, 10, 0,
0, 0, 31, 0, 7, 0, 0, 1, 0, 7, 0, 28, 0, 6, 0, 0, 0, 0, 0, 2, 0, 15, 30, 0, 0, 0, 0, 0, 0, 0, 23, 5, 0, 0, 0, 5, 0, 0, 0, 0,
6, 15, 0, 3, 19, 0, 3, 0, 0, 50, 7, 0, 0, 4, 0, 0, 30, 0, 5, 0, 0, 0, 0, 0, 0, 0, 9, 56, 0, 0, 0, 2, 0, 0, 0, 24, 5, 0, 0,
0, 0, 0, 0, 0, 0, 0, 11, 0, 0, 5, 29, 0, 0, 0, 3, 0, 1, 28, 0, 0, 0, 0, 0, 0, 0, 0, 21, 0, 0, 8, 0, 26, 0, 0, 8, 0, 0, 25, 7,
0, 0, 0, 0, 1, 0, 0, 0, 0, 26, 8, 0, 0, 0, 1, 0, 0, 0, 36, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 36, 1, 0, 49, 0, 2, 0, 7,
0, 0, 3, 0, 0, 31, 9, 0, 0, 0, 0, 0, 0, 9, 0, 41, 0, 2, 6, 0, 0, 0, 4, 0, 0, 31, 4, 0, 0, 0, 0, 2, 0, 1, 55, 0, 20, 3, 0, 0,
0, 5, 0, 0, 31, 0, 2, 0, 0, 0, 2, 12, 0, 0, 0, 6, 15, 0, 0, 23, 0, 0, 8, 22, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 6, 0, 22,
0, 15, 0, 0, 0, 22, 6, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 22, 13, 0, 0, 0, 0, 5, 20, 22, 0, 0, 3, 10, 0, 0, 0, 0, 0, 2, 0,
7, 0, 0, 36, 3, 0, 0, 0, 0, 2, 36, 0, 0, 4, 0, 0, 8, 12, 0, 0, 0, 0, 6, 0, 0, 36, 0, 16, 0, 0, 0, 5, 36, 0, 0, 0, 0, 10, 0,
0, 0, 0, 0, 0, 0, 0, 0, 41, 35, 0, 0, 0, 0, 0, 35, 0, 1, 0, 0, 0, 5, 0, 0, 0, 4, 0, 0, 0, 26, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0,
5, 25, 0, 5, 0, 0, 0, 5, 0, 0, 0, 7, 25, 0, 0, 5, 0, 0, 0, 0, 4, 25, 0, 20, 8, 0, 0, 6, 0, 0, 0, 0, 0, 46, 9, 0, 0, 0, 0, 0,
0, 7, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 50, 0, 13, 0, 0, 0, 0, 6, 30, 0, 0, 0, 0, 4, 0, 0, 0, 0, 6, 0, 30, 0, 30, 4,
0, 4, 0, 0, 30, 0, 4, 0, 0, 0, 0, 0, 0, 2, 0, 0, 30, 9, 0, 0, 0, 11, 0, 0, 0, 0, 34, 0, 0, 0, 0, 0, 0, 0, 25, 0, 0, 0, 0, 10,
0, 32, 0, 0, 3, 0, 0, 0, 0, 0, 5, 30, 0, 0, 0, 0, 1, 0)
fit <- pscl::zeroinfl(y~1|1, dist="negbin")

Proposed solutions:

  • tryCatch the error from the solve call and attempt to re-solve using a different method if the Hessian is singular.
  • Allow unsetting reltol such that L-BFGS-B can be used.
atahk commented

Thanks for this!

The issue in your example seems to be that the Hessian actually is singular because optim is finding the "wrong" solution (a boundary solution that fits the large number of zeros with the count model). As a result, you can get it to perform properly by giving it better start values. For example, you can use the estimates of a hurdle model, which are much closer to the solution, for the start values. Thus, the following works without error for me:
fit.hurdle <- pscl::hurdle(y~1|1, dist="negbin")
fit <- pscl::zeroinfl(y~1|1, dist="negbin", start=fit.hurdle$coefficients)

Obviously, the zeroinfl function should handle this situation better. I don't think using L-BFGS-B or computing the Hessian differently are the best solutions here because the Hessian actually should be singular (although I think we might want to compute the Hessian differently anyway). But we should at least catch this error. Maybe we should issue a warning with that and suggest different start values? Or maybe the default should be to get better start values rather than using a vector of zeros (e.g., to get start values for the count side from a truncated count model and observations with non-zero responses, as happens with hurdle.

dstoeckel commented

Thank you for your answer! I agree that switching to L-BFGS-B or anything else as a default would only be a band aid for this particular case and not a solution in general. The proper solution is probably to get better starting values and (as this will make the problem only less likely) to catch a singular Hessian + issue a warning.

For my particular case I will have to see whether I can convince the mpath package to take different starting values (It is computing a whole LASSO path, so the probabilities may change a bit ...). But this is my problem 🙁

simonjackman commented
atahk commented

I added error handling with the Hessian and improved the starting values a bit. This at least fixes the example above (both in that it finds the correct solution with the new default starting values and it handles the singular Hessian properly if given the old starting values). It's also worked fine on all the examples I've tried. It's not on CRAN yet and should probably be tested a bit more first. But, until then, it can be installed from the github version using devtools.