Errors producing a prior predictive distribution using smoothers in stan_gamm4 with rstanarm

I am having difficulty generating a prior predictive distribution for my count model when I try to add smooth predictors to the model.

I have included a reproducible example below where I generate the input data (much more concise than just pasting in the output of the data generation).

library(tidyverse)
library(rstanarm)


set.seed(42)
stan_seed <- 42


sector_vals <- c("Technology", "Industrials", "Communications")
mktcap_vals <- c("QN1", "QN2", "QN3", "QN4", "QN5")
btm_vals    <- c("QN1", "QN2", "QN3", "QN4", "QN5")


create_company_details <- function(n) {
  details_tbl <- tibble(
    sector      = sample(sector_vals, n, replace = TRUE),
    mktcap_qn   = sample(mktcap_vals, n, replace = TRUE),
    btm_qn      = sample(btm_vals,    n, replace = TRUE),
    prev3_cases = sample(c("0", "1", "2", "3+"), n,
                         prob = c(0.90, 0.06, 0.03, 0.01),
                         replace = TRUE),
    obs_cases   = sample(0:2, size = n,
                         prob = c(0.90, 0.08, 0.02),
                         replace = TRUE)
  )

  return(details_tbl)
}

example_data_tbl <- tibble(year = 2000:2018) %>%
  mutate(n_company = rnbinom(n(), mu = 5000, size = 50),
         data      = map(n_company, create_company_details)
    ) %>%
  select(-n_company) %>%
  unnest(data)


temp_prior_stanreg <- stan_gamm4(
  obs_cases ~ s(year, k = 4) + sector + mktcap_qn + btm_qn + prev3_cases,
  data     = example_data_tbl,
  random   = NULL,
  family   = neg_binomial_2(),
  iter     = 1000,
  chains   = 8,
  QR       = TRUE,
  prior_intercept = normal(location = -4, scale = 1.0, autoscale = FALSE),
  prior_smooth    = normal(location =  0, scale = 1.0, autoscale = FALSE),
  prior_aux       = normal(location =  0, scale = 1.0, autoscale = FALSE),
  prior           = normal(location =  0, scale = 1.0, autoscale = FALSE),
  prior_covariance = decov(),
  prior_PD = TRUE,
  seed     = stan_seed
)

The “generated quantities” appears to be creating zero or negative values for the count RNG.

starting worker pid=6106 on localhost:11285 at 11:00:49.583
starting worker pid=6107 on localhost:11285 at 11:00:49.591
starting worker pid=6104 on localhost:11285 at 11:00:49.591
starting worker pid=6103 on localhost:11285 at 11:00:49.591
starting worker pid=6105 on localhost:11285 at 11:00:49.592
starting worker pid=6110 on localhost:11285 at 11:00:49.592
starting worker pid=6108 on localhost:11285 at 11:00:49.592
starting worker pid=6109 on localhost:11285 at 11:00:49.594

SAMPLING FOR MODEL 'count' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.2e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)


SAMPLING FOR MODEL 'count' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 2.3e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 2: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)

SAMPLING FOR MODEL 'count' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 2.5e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)

SAMPLING FOR MODEL 'count' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 6.3e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)


SAMPLING FOR MODEL 'count' NOW (CHAIN 5).
Chain 5: 
Chain 5: Gradient evaluation took 3.2e-05 seconds
Chain 5: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
Chain 5: Adjust your expectations accordingly!
Chain 5: 
Chain 5: 
Chain 5: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 5: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 5: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 5: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 5: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 5: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 5: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)

SAMPLING FOR MODEL 'count' NOW (CHAIN 6).
Chain 6: 
Chain 6: Gradient evaluation took 2.4e-05 seconds
Chain 6: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
Chain 6: Adjust your expectations accordingly!
Chain 6: 
Chain 6: 
Chain 6: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 6: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 6: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 6: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 6: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 6: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 6: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)

SAMPLING FOR MODEL 'count' NOW (CHAIN 7).
Chain 7: 
Chain 7: Gradient evaluation took 3e-05 seconds
Chain 7: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
Chain 7: Adjust your expectations accordingly!
Chain 7: 
Chain 7: 
Chain 7: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 7: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 7: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 7: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 7: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 7: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 7: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)

SAMPLING FOR MODEL 'count' NOW (CHAIN 8).
Chain 8: 
Chain 8: Gradient evaluation took 3.6e-05 seconds
Chain 8: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
Chain 8: Adjust your expectations accordingly!
Chain 8: 
Chain 8: 
Chain 8: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 8: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 8: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 8: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 8: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 8: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 8: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 5: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 7: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 1: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 6: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 6: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 7: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 7: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 1: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 1: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 7: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 2: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 2: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 8: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 5: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 6: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.044179 seconds (Warm-up)
Chain 1:                17.2622 seconds (Sampling)
Chain 1:                17.3064 seconds (Total)
Chain 1: 
Chain 3: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 7: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 6: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.042655 seconds (Warm-up)
Chain 2:                17.6441 seconds (Sampling)
Chain 2:                17.6867 seconds (Total)
Chain 2: 
Chain 5: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 8: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 6: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 6: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 6: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 7: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 8: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 7: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.043791 seconds (Warm-up)
Chain 3:                19.1233 seconds (Sampling)
Chain 3:                19.1671 seconds (Total)
Chain 3: 
Chain 6: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 5: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.098138 seconds (Warm-up)
Chain 4:                19.4556 seconds (Sampling)
Chain 4:                19.5537 seconds (Total)
Chain 4: 
Chain 7: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 8: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 6: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 6: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 6: 
Chain 6:  Elapsed Time: 0.081746 seconds (Warm-up)
Chain 6:                17.0838 seconds (Sampling)
Chain 6:                17.1656 seconds (Total)
Chain 6: 
Chain 7: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 5: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 5: 
Chain 5:  Elapsed Time: 0.085009 seconds (Warm-up)
Chain 5:                19.2444 seconds (Sampling)
Chain 5:                19.3294 seconds (Total)
Chain 5: 
Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Exception: poisson_rng: Rate parameter is 0, but must be > 0!  (in 'model_count' at line 205)

Chain 8: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 8: 
Chain 8:  Elapsed Time: 0.07378 seconds (Warm-up)
Chain 8:                16.2533 seconds (Sampling)
Chain 8:                16.3271 seconds (Total)
Chain 8: 
Chain 7: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 7: 
Chain 7:  Elapsed Time: 0.074383 seconds (Warm-up)
Chain 7:                17.9131 seconds (Sampling)
Chain 7:                17.9874 seconds (Total)
Chain 7: 
Error in checkForRemoteErrors(val) : 
  8 nodes produced errors; first error: invalid class “stanfit” object: The following variables have undefined values:  alpha[1],The following variables have undefined values:  mean_PPD

I imagine I need to adjust the various priors to prevent this from happening, but I am not sure how to proceed?

My old version of the model had the year as a categorical and I used stan_glm without the smoothing on year, and that worked fine, but I am not sure what to do for this case where I want to use smoothing.

The code is running on a Docker instance based on rocker/verse:4.0.0 with rstanarm 2.19.3

Thanks in advance for all your help!

Happy to give more details if needed.

You can try changing the seed as in this thread here:

I know a lot folks don’t run with seed values initially.

You might want to also set the init value to 0 or 1 I can’t remember with the poosson_rng what the init should be at.

There were a few typos in my reprex, so I fixed those.

No luck so far changing seeds, but I thought with the init values you just specify an interval for the unconstrained space - either way, this all sounds fairly brittle so I suspect I will need to change my priors to sort the problem properly.

1 Like

Just bumping this thread in case someone has some pointers for this - it may not be doable as it stands.

@kaybenleroll Thanks for bumping the thread. When prior_PD = TRUE I think we can turn off the part of generated quantities where the error is occurring without affecting anything. I’ll fix this on GitHub, but there’s a hack you can use in the meantime (I think, can you let me know if this works?). Right before running your code to create temp_prior_stanreg can you try running this?

stan_glm.fit2 <- rstanarm::stan_glm.fit
args <- formals(stan_glm.fit2)
args$mean_PPD <- FALSE
formals(stan_glm.fit2) <- args
assignInNamespace("stan_glm.fit", stan_glm.fit2, ns = "rstanarm")
# now run your code for temp_prior_stanreg

Thanks @jonah - I’ll try that in the next hour or so and let you know how I get on.

@jonah - I think the second line didn’t need the rstanarm:: part due to the assignment in the first line, but that did seem to do the trick - I can now run it to get a prior predictive distribution.

1 Like

Oh, you’re right, thanks. Just edited my post to fix that.

Great, glad that works! I fixed this on GitHub so it shouldn’t be necessary to use this workaround when the next release comes out.

Brilliant - I’ll use the workaround for now and keep an eye out for the next release.