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

I get this error several times until the model fails:

Chain 4 Exception: Exception: neg_binomial_2_log_glm_lpmf: Failures variables[15] is -1.45964e+09, but must be nonnegative! (in ‘C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395423e33d2a.stan’, line 28, column 4 to column 88) (in ‘C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395423e33d2a.stan’, line 85, column 4 to column 120)

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!

Any help to solve this error would be much appreciated.

Lines 24-30:

    for (n in 1:N) {
      // add more terms to the linear predictor
      int nn = n + start - 1;
      mu[n] += r_1_1[J_1[nn]] * Z_1_1[nn];
    }
    ptarget += neg_binomial_2_log_glm_lpmf(Y[start:end] | XQ[start:end], mu, bQ, shape);
    return ptarget;
  }

Lines 83-87:

model {
  // likelihood including constants
  if (!prior_only) {
    target += reduce_sum(partial_log_lik_lpmf, seq, grainsize, Y, bQ, Intercept, XQ, offsets, shape, J_1, Z_1_1, r_1_1);
  }

Full model:

// full model
// generated with brms 2.16.4
functions {
  /* integer sequence of values
   * Args: 
   *   start: starting integer
   *   end: ending integer
   * Returns: 
   *   an integer sequence from start to end
   */ 
  int[] sequence(int start, int end) { 
    int seq[end - start + 1];
    for (n in 1:num_elements(seq)) {
      seq[n] = n + start - 1;
    }
    return seq; 
  } 
  // compute partial sums of the log-likelihood
  real partial_log_lik_lpmf(int[] seq, int start, int end, data int[] Y, vector bQ, real Intercept, data matrix XQ, data vector offsets, real shape, data int[] J_1, data vector Z_1_1, vector r_1_1) {
    real ptarget = 0;
    int N = end - start + 1;
    // initialize linear predictor term
    vector[N] mu = Intercept + rep_vector(0.0, N) + offsets[start:end];
    for (n in 1:N) {
      // add more terms to the linear predictor
      int nn = n + start - 1;
      mu[n] += r_1_1[J_1[nn]] * Z_1_1[nn];
    }
    ptarget += neg_binomial_2_log_glm_lpmf(Y[start:end] | XQ[start:end], mu, bQ, shape);
    return ptarget;
  }
}
data {
  int<lower=1> N;  // total number of observations
  int Y[N];  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  vector[N] offsets;
  int grainsize;  // grainsize for threading
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  int<lower=1> J_1[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_1;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int Kc = K - 1;
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  // matrices for QR decomposition
  matrix[N, Kc] XQ;
  matrix[Kc, Kc] XR;
  matrix[Kc, Kc] XR_inv;
  int seq[N] = sequence(1, N);
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
  // compute and scale QR decomposition
  XQ = qr_thin_Q(Xc) * sqrt(N - 1);
  XR = qr_thin_R(Xc) / sqrt(N - 1);
  XR_inv = inverse(XR);
}
parameters {
  vector[Kc] bQ;  // regression coefficients at QR scale
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> shape;  // shape parameter
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  vector[N_1] z_1[M_1];  // standardized group-level effects
}
transformed parameters {
  vector[N_1] r_1_1;  // actual group-level effects
  real lprior = 0;  // prior contributions to the log posterior
  r_1_1 = (sd_1[1] * (z_1[1]));
  lprior += student_t_lpdf(bQ | 3, 0, 1);
  lprior += student_t_lpdf(Intercept | 3, 0, 1);
  lprior += gamma_lpdf(shape | 0.01, 0.01);
  lprior += student_t_lpdf(sd_1 | 3, 0, 1)
    - 1 * student_t_lccdf(0 | 3, 0, 1);
}
model {
  // likelihood including constants
  if (!prior_only) {
    target += reduce_sum(partial_log_lik_lpmf, seq, grainsize, Y, bQ, Intercept, XQ, offsets, shape, J_1, Z_1_1, r_1_1);
  }
  // priors including constants
  target += lprior;
  target += std_normal_lpdf(z_1[1]);
}
generated quantities {
  // obtain the actual coefficients
  vector[Kc] b = XR_inv * bQ;
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
  // additionally sample draws from priors
  real prior_bQ = student_t_rng(3,0,1);
  real prior_Intercept = student_t_rng(3,0,1);
  real prior_shape = gamma_rng(0.01,0.01);
  real prior_sd_1 = student_t_rng(3,0,1);
  // use rejection sampling for truncated priors
  while (prior_shape < 0) {
    prior_shape = gamma_rng(0.01,0.01);
  }
  while (prior_sd_1 < 0) {
    prior_sd_1 = student_t_rng(3,0,1);
  }
}

The error:

Chain 4 Exception: Exception: neg_binomial_2_log_glm_lpmf: Failures variables[15] is -1.45964e+09, but must be nonnegative!

Is suggesting that the outcome count for your model contains unitialised values, which is pretty hard to do for arguments that are data. Can you post the brms code that you’re using here? Or have you modified this from what was generated by brms?

Additionally, do you get the same error if you run the brms model without threading?

1 Like

Thanks @andrjohns. Your input helped me find out that the errors I am getting are because of machine precision issues: my 15th y value exceeds the max value of an integer in R (2147483647).

How can I get around this effectively and efficiently given my model above and get the same inferences?

Does that mean that your y value is not currently representable in R before being passed to brms/Stan, or that it’s fine in R but breaks when you go to brms/Stan?

It is fine in R as a double, but not as an integer. It becomes NA when I try to coerce it to an integer. So, when I run the model, brms does not recognize it as an integer, and the errors appear. Using rstan instead of cmdstan as backend, the error is more informative, saying that some value in y is not int.

Hmm, this is a difficult one, because the maximum representable integer is the same for both R and c++.

A common alternative in c++ is to use ‘unsigned’ integers, essentially saying that we don’t need the negative numbers and can just increase the number of representable positive integers. However, our c++ generating infrastructure just generates integers.

@rok_cesnovar A fairly useful addition in the future could be having stanc3 generate c++ that uses size_t/uint_x if the stan code specifies int<lower=0>, we’d just need to make sure that the distributions don’t attempt to cast/coerce the inputs to integer.

For your model here, an immediately accessible workaround could be to pass your count outcome as a numeric variable (so that the integer limitations are avoided), and use the custom_family framework to re-define the neg_binomial_2_log family to take real outcomes.

Here’s an implementation of that which appears to work:

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.16.11). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar

set.seed(2022)

brm_data <- data.frame(y = sample(10:100,5,replace=F),
                       denom = sample(10:100,5,replace=F),
                       x1 = rnorm(5,0,1),
                       x2 = rnorm(5,3,2))

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]")


mod_cont <- brm(y | vreal(denom) ~ x1 + x2,
                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 0.4 seconds.
#> Chain 2 finished in 0.5 seconds.
#> Chain 3 finished in 0.4 seconds.
#> Chain 4 finished in 0.4 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 0.5 seconds.
#> Total execution time: 2.4 seconds.

mod_disc <- brm(y | rate(denom) ~ x1 + x2,
                data = brm_data,
                family = negbinomial(),
                backend = "cmdstanr",
                control = list(adapt_delta = 0.99),
                refresh = 0, silent = 2)
#> Running MCMC with 4 sequential chains...
#> 
#> Chain 1 finished in 0.5 seconds.
#> Chain 2 finished in 0.4 seconds.
#> Chain 3 finished in 0.5 seconds.
#> Chain 4 finished in 0.5 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 0.5 seconds.
#> Total execution time: 2.6 seconds.

summary(mod_cont)
#>  Family: neg_binomial_2_log_continuous 
#>   Links: mu = log; shape = identity 
#> Formula: y | vreal(denom) ~ x1 + x2 
#>    Data: brm_data (Number of observations: 5) 
#>   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     1.70      1.36    -0.13     5.28 1.01      537      303
#> x1           -1.05      0.64    -2.68    -0.18 1.01      485      299
#> x2           -0.14      0.24    -0.58     0.24 1.01      714      382
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> shape     0.18      0.22     0.01     0.65 1.01      328      364
#> 
#> 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).
summary(mod_disc)
#>  Family: negbinomial 
#>   Links: mu = log; shape = identity 
#> Formula: y | rate(denom) ~ x1 + x2 
#>    Data: brm_data (Number of observations: 5) 
#>   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     1.70      1.12     0.13     4.75 1.01      760      442
#> x1           -1.02      0.57    -2.60    -0.20 1.01      738      351
#> x2           -0.16      0.16    -0.54     0.11 1.01     1351      561
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> shape     0.19      0.21     0.01     0.78 1.01      417      792
#> 
#> 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).

Created on 2022-03-24 by the reprex package (v2.0.1)

4 Likes

Thanks for this solution. I think we are making progress. Now I get this error:

Semantic error in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-39542fd962ab.stan', line 57, column 74 to column 80:
   -------------------------------------------------
    55:      for (n in 1:N) {
    56:        int nn = n + start - 1;
    57:        ptarget += neg_binomial_2_log_continuous_lpdf(Y[nn] | mu[n], shape, vreal1[nn]);
                                                                                   ^
    58:      }
    59:      return ptarget;
   -------------------------------------------------

Identifier 'vreal1' not in scope.mingw32-make.exe: *** [make/program:47: C:\Users\x\AppData\Local\Temp\Rtmp2D4vAA\model-39542fd962ab.hpp] Error 1

Hm, it looks like brms isn’t exporting custom family family parameters to the reduce_sum function.

Does it run if you don’t use the threading argument?

Using rstan as backend, the error message is:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "vreal1" does not exist.
 error in 'model395419756840_df84cf1cd3b42910b45eaa2fc51dbe07' at line 89, column 78
  -------------------------------------------------
    87:     }
    88:     for (n in 1:N) {
    89:       target += neg_binomial_2_log_continuous_lpdf(Y[n] | mu[n], shape, vreal1[n]);
                                                                                     ^
    90:     }
  -------------------------------------------------

 Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'df84cf1cd3b42910b45eaa2fc51dbe07' due to the above error.

What brms syntax are you using?

brms sintax:

brm(formula,family=nb_continuous,stanvars = stan_lpmf,
     			data=dat,
                iter=10000,
				chains=4, seed=2020, control=list(
				adapt_delta=
				0.999, 
                max_treedepth=
                15), 
				backend = "rstan",
				cores = cores,
				inits = "random",				
				sample_prior = "yes",
				save_pars = save_pars(latent = TRUE, all = TRUE),
				prior = prior)

And forformula?

formula = bf(y ~ 1 + offset(x1) + (1|x2) + x3, decomp = "QR")

Ah there it is. Instead of offset you need to use vreal, have a look at the example code I used above

1 Like

I see, thanks. Now I get this error:

Chain 4 Exception: Exception: Exception: binomial_coefficient_log: first argument is -1.35428e+09, but must be greater than or equal to -1 (in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395410233d3.stan', line 33, column 4 to column 39) (in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395410233d3.stan', line 57, column 6 to column 86) (in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395410233d3.stan', line 116, column 4 to column 119)
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!

Hmm, that’s also indicative of some overflow issues. What is the largest value of your outcome? The largest number c++ can represent as a double (i.e., real) is:

#include <iostream>
#include <limits>

int main(){
    std::cout << std::numeric_limits<double>::max() << std::endl;
    return 0;
}

> 1.79769e+308
1 Like

max(y)=2.83533e+09

Can you try running the model with subsets of data and with fewer/no predictors to identify what inputs are causing the issue? I can run the model with that sized outcome with no issues:

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)

Running MCMC with 4 chains, at most 8 in parallel...

Chain 2 finished in 0.4 seconds.
Chain 3 finished in 0.4 seconds.
Chain 1 finished in 0.8 seconds.
Chain 4 finished in 2.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.9 seconds.
Total execution time: 2.3 seconds.

On my machine:

##y = 2.83533e+08 works:

brm_data <- data.frame(y = 2.83533e+08,
                       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)

Running MCMC with 4 chains, at most 8 in parallel...

Chain 2 finished in 0.6 seconds.
Chain 4 finished in 0.6 seconds.
Chain 1 finished in 0.7 seconds.
Chain 3 finished in 4.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 1.4 seconds.
Total execution time: 4.3 seconds.

##But y = 2.83533e+09 does not work:

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)

Chain 4 Exception: Exception: Exception: binomial_coefficient_log: first argument is -1.4211e+09, but must be greater than or equal to -1 (in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395410233d3.stan', line 19, column 4 to column 39) (in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395410233d3.stan', line 55, column 6 to column 83) (in 'C:/Users/x/AppData/Local/Temp/Rtmp2D4vAA/model-395410233d3.stan', line 116, column 4 to column 119)
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!

Can I ask for help on this issue.

Pinging: @andrjohns, @avehtari, @betanalpha, @bgoodri, @Bob_Carpenter, @martinmodrak, @maxbiostat, @paul.buerkner, @spinkney, et al.