Exception: neg_binomial_2_log_glm_lpmf: Failures variables[15] is -1.45964e+09

Modeling
, ,
22

Unfortunately, at those extreme values, ensuring that the neg. binomial density is numerically well behaved becomes challenging. The ordering of operations seen at math/neg_binomial_2_lpmf.hpp at develop · stan-dev/math · GitHub should probably be an improvement (but it is a bit cryptic code and was never tested outside the range of integeres)

I propose a different workaround: If all your counts are very large, then there should be very little difference between using Gamma-distributed response and a neg.binomial response (Neg. binomial is a Poisson noise over a gamma-distributed variable, with large means, the Poisson noise is negligible). Gamma works on doubles so should not have this type of problems (and also has simple scaling properties). The shape parameter used by brms should correspond to the phi parameter of the neg. binomial.

A quick example

large_data <- data.frame(y = rnbinom(n = 100, mu = 50000, size = 2))
brm(y ~ 1, family = "negbinomial", data = large_data)

gives:

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    10.79      0.07    10.65    10.94 1.00     2735     2316

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     1.90      0.25     1.46     2.40 1.00     3448     2681

while

brm(y ~ 1, family = "gamma", data = large_data)

gives

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    10.80      0.07    10.66    10.94 1.00     2933     2559

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     1.90      0.25     1.43     2.43 1.00     2829     2084

Best of luck with your model!

4 Likes
23

Thanks. There are other considerations that make nb2 desirable, including the need to account for the exposure variable, etc.

24

I don’t know why that doesn’t work. It works fine for me. I even dialed up the values to

brm_data <- data.frame(y = 2.83533e+30,
                       denom = 2.83533e+29)

mod_cont <- brm(y | vreal(denom) ~ 1,
                data = brm_data,
                family = nb_continuous,
                stanvars = stan_lpmf,
                backend = "cmdstanr",
                control = list(adapt_delta = 0.99),
                refresh = 0, silent = 2)

edit: adding brms and cmdstan versions

# Loading 'brms' package (version 2.16.3). 
cmdstanr::cmdstan_version()
# [1] "2.29.1"
5 Likes
25

Which ones? Things like offset should work exactly the same for gamma. At those values, they are basically indistinguishable and you can use gamma_lpdf(y | shape, shape / mu) as a quite close approximation to neg_binomial2_lpmf(y | mu, shape):

mu = 1e5
shape = 3
y = round(seq(0.5e5, 2e5, length.out = 100))
plot(y, dgamma(y, shape = shape, rate = shape/mu, log = TRUE))
lines(y, dnbinom(y, mu = mu, size = shape, log = TRUE))

image

Or we can look at the difference:

plot(y, dgamma(y, shape = shape, rate = shape/mu, log = TRUE) - 
  dnbinom(y, mu = mu, size = shape, log = TRUE))

(note that the scale is know extremely small):

image

4 Likes
26

I just tried again:

library(brms)

set.seed(2022)

neg_binomial_2_log_continuous <- "
  real neg_binomial_2_log_continuous_lpdf(real y, real mu_input,
                                          real phi_input, real offset_arg) {

    real phi = phi_input * offset_arg;
    real log_phi = log(phi);

    real log_offset = log(offset_arg);
    real eta = log(mu_input) + log_offset;
    real exp_eta = exp(eta);
    real eta_over_exp_eta_phi = -log_phi - log1p(exp_eta);
    real exp_eta_over_exp_eta_phi = exp(eta_over_exp_eta_phi);
    real log1p_exp_eta_m_logphi = log1p_exp(eta - log_phi);
    real y_plus_phi = y + phi;

    real logp = 0;
    logp += lchoose(y_plus_phi - 1, y);
    logp += y * eta;
    logp += -phi * log1p_exp_eta_m_logphi - y * (log_phi + log1p_exp_eta_m_logphi);

    return logp;
  }
"
stan_lpmf <- stanvar(scode = neg_binomial_2_log_continuous,
                     block = "functions")

nb_continuous <- custom_family(
  name = "neg_binomial_2_log_continuous",
  dpars = c("mu", "shape"),
  links = "log",
  lb = c(NA, 0),
  type = "real",
  vars = "vreal1[n]")

brm_data <- data.frame(y = 2.83533e+09,
                       denom = 2.83533e+07)

mod_cont <- brm(y | vreal(denom) ~ 1,
                data = brm_data,
                family = nb_continuous,
                stanvars = stan_lpmf,
                backend = "cmdstanr",
                control = list(adapt_delta = 0.99),
                refresh = 0, silent = 2)

Here is the full error message from chain 4 up to the end:

Chain 4 Rejecting initial value:
Chain 4 Error evaluating the log probability at the initial value.
Chain 4 Exception: Exception: binomial_coefficient_log: first argument is -1.41176e+09, but must be greater than or equal to -1 (in ‘C:/Users/x/AppData/Local/Temp/RtmpsXf9uZ/model-4b3c72c2e9e.stan’, line 19, column 4 to column 39) (in ‘C:/Users/x/AppData/Local/Temp/RtmpsXf9uZ/model-4b3c72c2e9e.stan’, line 55, column 6 to column 83)
Chain 4 Exception: Exception: binomial_coefficient_log: first argument is -1.41176e+09, but must be greater than or equal to -1 (in ‘C:/Users/x/AppData/Local/Temp/RtmpsXf9uZ/model-4b3c72c2e9e.stan’, line 19, column 4 to column 39) (in ‘C:/Users/x/AppData/Local/Temp/RtmpsXf9uZ/model-4b3c72c2e9e.stan’, line 55, column 6 to column 83)
Chain 4 Initialization between (-2, 2) failed after 100 attempts.
Chain 4 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
Warning: Chain 4 finished unexpectedly!

Warning: Use read_cmdstan_csv() to read the results of the failed chains.
Error in cmdstanr::read_cmdstan_csv(out$output_files(), variables = “”, :
Assertion on ‘files’ failed: No file provided.
In addition: Warning messages:
1: All chains finished unexpectedly! Use the $output(chain_id) method for more information.

2: No chains finished successfully. Unable to retrieve the fit.

sessionInfo():

R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22000)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

other attached packages:
[1] brms_2.16.4 Rcpp_1.0.6

loaded via a namespace (and not attached):
[1] nlme_3.1-152 matrixStats_0.59.0 xts_0.12.1
[4] threejs_0.3.3 rstan_2.21.2 tensorA_0.36.2
[7] tools_4.1.0 backports_1.2.1 utf8_1.2.1
[10] R6_2.5.0 DT_0.18 mgcv_1.8-35
[13] DBI_1.1.1 projpred_2.0.2 colorspace_2.0-2
[16] withr_2.4.2 tidyselect_1.1.1 gridExtra_2.3
[19] prettyunits_1.1.1 processx_3.5.2 Brobdingnag_1.2-6
[22] curl_4.3.2 compiler_4.1.0 cli_3.0.0
[25] shinyjs_2.0.0 colourpicker_1.1.0 posterior_1.2.0
[28] scales_1.1.1 dygraphs_1.1.1.6 checkmate_2.0.0
[31] mvtnorm_1.1-2 ggridges_0.5.3 callr_3.7.0
[34] stringr_1.4.0 digest_0.6.27 StanHeaders_2.21.0-7
[37] minqa_1.2.4 base64enc_0.1-3 pkgconfig_2.0.3
[40] htmltools_0.5.1.1 lme4_1.1-27.1 fastmap_1.1.0
[43] htmlwidgets_1.5.3 rlang_0.4.11 shiny_1.6.0
[46] farver_2.1.0 generics_0.1.0 zoo_1.8-9
[49] jsonlite_1.7.2 crosstalk_1.1.1 gtools_3.9.2
[52] dplyr_1.0.7 distributional_0.2.2 inline_0.3.19
[55] magrittr_2.0.1 loo_2.4.1 bayesplot_1.8.1
[58] Matrix_1.3-3 munsell_0.5.0 fansi_0.5.0
[61] abind_1.4-5 lifecycle_1.0.0 stringi_1.6.2
[64] MASS_7.3-54 pkgbuild_1.2.0 plyr_1.8.6
[67] grid_4.1.0 parallel_4.1.0 promises_1.2.0.1
[70] crayon_1.4.1 miniUI_0.1.1.1 lattice_0.20-44
[73] splines_4.1.0 knitr_1.33 ps_1.6.0
[76] pillar_1.6.1 igraph_1.2.6 boot_1.3-28
[79] markdown_1.1 shinystan_2.5.0 reshape2_1.4.4
[82] codetools_0.2-18 stats4_4.1.0 rstantools_2.1.1
[85] glue_1.4.2 V8_3.4.2 RcppParallel_5.1.4
[88] nloptr_1.2.2.2 vctrs_0.3.8 httpuv_1.6.1
[91] gtable_0.3.0 purrr_0.3.4 assertthat_0.2.1
[94] ggplot2_3.3.5 xfun_0.24 mime_0.11
[97] xtable_1.8-4 coda_0.19-4 later_1.2.0
[100] rsconnect_0.8.18 tibble_3.1.2 shinythemes_1.2.0
[103] tinytex_0.32 gamm4_0.2-6 cmdstanr_0.4.0
[106] ellipsis_0.3.2 bridgesampling_1.1-2

1 Like
27

Thanks. I will continue exploring gamma, but nb2 is the desirable model for my case. I tested gamma early in the project and found that with almost the same specification as nb2, nb2 performed better. Also, it looks like gamma could be challenging to interpret in my case, etc.

1 Like
28

and what’s your version of cmdstan? Running the above

29 
cmdstanr::cmdstan_version()
"2.27.0"
30

can you update cmdstan by running

cmdstanr::install_cmdstan(cores = 4)
1 Like
31

Thanks. Updated cmdstan to 2.29.2. Now it seems to work BUT:

mod_cont <- brm(y | vreal(denom) ~ 1,
                 data = brm_data,
                 family = nb_continuous,
                 stanvars = stan_lpmf,
                 backend = "cmdstanr",
                 control = list(adapt_delta = 0.99),
                 refresh = 0, silent = 2)

Running MCMC with 4 sequential chains…

Chain 1 finished in 5.5 seconds.
Chain 2 finished in 0.9 seconds.
Chain 3 finished in 0.6 seconds.
Chain 4 finished in 4.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.8 seconds.
Total execution time: 12.2 seconds.

Warning: 8 of 4000 (0.0%) transitions ended with a divergence.
This may indicate insufficient exploration of the posterior distribution.
Possible remedies include:

  • Increasing adapt_delta closer to 1 (default is 0.8)
  • Reparameterizing the model (e.g. using a non-centered parameterization)
  • Using informative or weakly informative prior distributions

1079 of 4000 (27.0%) transitions hit the maximum treedepth limit of 10 or 2^10-1 leapfrog steps.
Trajectories that are prematurely terminated due to this limit will result in slow exploration.
Increasing the max_treedepth limit can avoid this at the expense of more computation.
If increasing max_treedepth does not remove warnings, try to reparameterize the model.

mod_cont

Family: neg_binomial_2_log_continuous
Links: mu = log; shape = identity
Formula: y | vreal(denom) ~ 1
Data: brm_data (Number of observations: 1)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 12.20 8.57 4.60 26.59 1.73 7 144

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.40 7.22 0.00 21.31 1.95 6 22

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Warning messages:
1: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results! We recommend running more iterations and/or setting stronger priors.
2: There were 8 divergent transitions after warmup. Increasing adapt_delta above may help. See Runtime warnings and convergence problems

Now, testing the actual data…

2 Likes
32

Thanks, @spinkney, @andrjohns, @martinmodrak, for your help. I think my problem has been solved now (using nb_continuous and after updating cmdstan to 2.27.0).

I can fit a nb2 model without divergence and the estimates are technically sound. Perhaps it is difficult to imagine the relief that this brings.

Now, which post should I select as the answer?

Additionally, can I ask additional help to solve this issue when using brms::fitted (or should I open another post):

Error in get(out, family$env) : 
object 'posterior_epred_neg_binomial_2_log_continuous' not found
2 Likes
33

When using brms's post-processing facilities with custom families you also need to define the necessary functions for brms to use, have a look at the vignette for more background on how to do this: Define Custom Response Distributions with brms • brms

3 Likes
34

Done. This link you shared is all I needed for the last part. Again, thanks, @andrjohns.

2 Likes