0

I am running a two-part simulation based on a hurdle model. However, when bootstrapping and resampling, I am running into an error

Error in solve.default(as.matrix(fit_zero$hessian)) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0

After running traceback(), I am greeted with with this :

19: solve(as.matrix(fit_zero$hessian))
18: hurdle(formula = y ~ m + x, data = boot.data0, dist = "poisson", 
        zero.dist = "binomial") at #21
17: bootstrap(v)
16: FUN(X[[i]], ...)
15: lapply(X = X, FUN = FUN, ...)
14: sapply(integer(n), eval.parent(substitute(function(...) expr)), 
        simplify = simplify)
13: replicate(r, bootstrap(v), simplify = FALSE) at #32
12: unlist(replicate(r, bootstrap(v), simplify = FALSE)) at #32
11: matrix(unlist(replicate(r, bootstrap(v), simplify = FALSE)), 
        ncol = 8, byrow = TRUE) at #32
10: qnorm((sum(boot[, 3] > boot[, 1]) + sum(boot[, 3] == boot[, 1])/2)/r) at #3
9: bias.correct(matrix(unlist(replicate(r, bootstrap(v), simplify = FALSE)), 
       ncol = 8, byrow = TRUE), r) at #32
8: FUN(X[[i]], ...)
7: lapply(X = X, FUN = FUN, ...)
6: sapply(model_values, bootstrapper) at #35
5: FUN(newX[, i], ...)
4: apply(getParameters()[, 2:7], 1, function(parameters) {
       print(parameters)
       PROGRESS <- 0
       seed <- parameters["seed"]
       n <- parameters["n"]
       a <- parameters["a"]
       b <- parameters["b"]
       c <- parameters["c"]
       i <- parameters["i"]
       set.seed(seed)
       model_values <- replicate(iterations, model(n, a, b, c, i), 
           simplify = FALSE)
       bootstrapper <- function(v) {
           PROGRESS <<- PROGRESS + 1
           print(PROGRESS)
           bias.correct(matrix(unlist(replicate(r, bootstrap(v), 
               simplify = FALSE)), ncol = 8, byrow = TRUE), r)
       }
       boot.fit <- sapply(model_values, bootstrapper)
       boot.fit.matrix <- matrix(unlist(boot.fit), ncol = 4, byrow = TRUE)
       averaged <- apply(boot.fit.matrix, 2, mean)
       return(averaged)
   }) at #9
3: main()
2: as_mapper(.f)
1: possibly(main(), otherwise = "error") 

After a previous, vague stackoverflow question and a lot of research, I realized that it is entirely possible for the distribution to be all 0's or all 1's in that particular sample. While I cannot use that sample, I do not want the loop (or replicate in this case) to stop.

What I would like to know is if there is a way to create a condition such as "if the matrix is singular, sample the next iteration".

My randomly generated data comes from this part of the code:


gen.hurdle = function(n, a, b1, b2, c1, c2, i0, i1, i2){

  x = round(rnorm(n),3)
  e = rnorm(n)
  m = round(i0 + a*x + e, 3)

  lambda = exp(i1 + b1*m + c1*x)                       # PUT REGRESSION TERMS FOR THE CONTINUUM PART HERE; KEEP exp()
  ystar = qpois(runif(n, dpois(0, lambda), 1), lambda) # Zero-TRUNCATED POISSON DIST.; THE CONTINUUM PART

  z = i2 + b2*m  + c2*x                                # PUT REGRESSION TERMS FOR THE BINARY    PART HERE
  z_p = exp(z) / (1+exp(z))                            # p(1) = 1-p(0)
  tstar = rbinom(n, 1, z_p)                            # BINOMIAL DIST.         ; THE BINARY    PART

  y= ystar*tstar                                       # TWO-PART COUNT OUTCOME

  return(cbind(x,m,y,z,z_p,tstar))
}

# Returns the base model's powers, type 1 errors, and data
model <- function(n, a, b, c, i){
  #generate random data
  data  = data.frame(gen.hurdle(n, a, b, b, c, c, i, i, i))
  data0 = data.frame(gen.hurdle(n, a, 0, 0, c, c, i, i, i))

I have tried try() and trycatch() to no avail because I have many nested functions. However, I read a similar SO post that stated "Instead of using tryCatch you could simply calculate the determinant of the matrix with the function det. A matrix is singular if and only if the determinant is zero."

I need some help with setting this condition.

Nazim Kerimbekov
  • 4,712
  • 8
  • 34
  • 58
J Bacon
  • 27
  • 6
  • Sounds like you want TryCatch: http://mazamascience.com/WorkingWithData/?p=912 – jeremycg May 29 '19 at 18:58
  • @jeremycg I can give it another go for sure. One big question I have with TryCatch is does it replace the error'd sample? For instance, say I am trying to resample 50 times. Time 15 throws an error. Does it proceed and have 49 samples or skip the error message, sample the next iteration, and keep the n=50? I would prefer the latter to the former. – J Bacon May 29 '19 at 19:02
  • Please show some more of the relevant code; I'm sure that `gen.hurdle` alone works fine. Ideally we should be able to reproduce the error. – Julius Vainora May 29 '19 at 23:12

0 Answers0