Initialization failure in negative-binomial model with exact and right-censored likelihood contributions

Modeling
,
1

I’m trying to fit a large spatial-temporal Bayesian model with close to 1,700 areas across 10 years in Stan using a negative-binomial model with exact and right-censored likelihood contributions. Predictors include meterological variables built using DLNM statistical framework.

Here is the Stan code of the model ```


data {
  // ---------------------------------------------------------------------------
  // Exact observed data
  // ---------------------------------------------------------------------------

  int<lower=0> N;
  // Number of exactly observed observations

  array[N] int<lower=0> Y;
  // Exactly observed case counts

  int<lower=2> K;
  // Intercept and at least one basis column

  matrix[N, K] X;
  // Design matrix including the intercept column

  vector[N] log_offset;
  // Log population offset

  int<lower=1> N_month_lfo;
  // Number of calendar months

  int<lower=1> S;
  // Number of subcounties

  array[N] int<lower=1, upper=S> subcounty_id;
  // Subcounty index for exact observations

  // ---------------------------------------------------------------------------
  // Right-censored observed data
  // ---------------------------------------------------------------------------

  int<lower=0> N_cens;
  // Number of right-censored observations

  array[N_cens] int<lower=1> Y_cens;
  // Censored lower-bound counts

  matrix[N_cens, K] X_cens;
  // Design matrix for censored observations

  vector[N_cens] log_offset_cens;
  // Log population offset for censored observations

  array[N_cens] int<lower=1, upper=S> subcounty_id_cens;
  // Subcounty index for censored observations

  // ---------------------------------------------------------------------------
  // Held-out exact data for M-step-ahead LFO-CV
  // ---------------------------------------------------------------------------

  int<lower=0> N_test;
  // Number of held-out exact observations

  array[N_test] int<lower=0> Y_test;
  // Held-out exact case counts

  matrix[N_test, K] X_test;
  // Held-out exact design matrix

  vector[N_test] log_offset_test;
  // Held-out exact log population offset

  array[N_test] int<lower=1, upper=N_month_lfo> month_id_test;
  // Held-out exact LFO-CV month index

  array[N_test] int<lower=1, upper=S> subcounty_id_test;
  // Held-out exact subcounty index

  // ---------------------------------------------------------------------------
  // Held-out right-censored data for M-step-ahead LFO-CV
  // ---------------------------------------------------------------------------

  int<lower=0> N_test_cens;
  // Number of held-out censored observations

  array[N_test_cens] int<lower=1> Y_test_cens;
  // Held-out censored lower-bound counts

  matrix[N_test_cens, K] X_test_cens;
  // Held-out censored design matrix

  vector[N_test_cens] log_offset_test_cens;
  // Held-out censored log population offset

  array[N_test_cens] int<lower=1, upper=N_month_lfo>
    month_id_test_cens;
  // Held-out censored LFO-CV month index

  array[N_test_cens] int<lower=1, upper=S>
    subcounty_id_test_cens;
  // Held-out censored subcounty index

  // ---------------------------------------------------------------------------
  // LFO-CV indices
  // ---------------------------------------------------------------------------

  array[N] int<lower=1, upper=N_month_lfo> month_id;
  // LFO-CV month index for exact observations

  array[N_cens] int<lower=1, upper=N_month_lfo> month_id_cens;
  // LFO-CV month index for censored observations

  int<lower=0, upper=1> prior_only;
  // 1 = prior predictive sampling; 0 = fit to observations

  // ---------------------------------------------------------------------------
  // Priors
  // ---------------------------------------------------------------------------

  vector[K - 1] prior_b_mean;
  // Prior means for fixed-effect coefficients

  vector[K - 1] prior_b_sd;
  // Prior standard deviations for fixed-effect coefficients

  real prior_intercept_mean;
  // Prior mean for the intercept

  real prior_intercept_sd;
  // Prior standard deviation for the intercept

  real prior_shape;
  // Rate parameter for the exponential prior on NB shape

  real<lower=0> prior_sigma_subcounty;
  // Prior standard deviation for between-subcounty variation
}

parameters {
  // ---------------------------------------------------------------------------
  // Fixed effects
  // ---------------------------------------------------------------------------

  vector[K - 1] b;
  // Fixed-effect coefficients

  real Intercept;
  // Overall model intercept

  real<lower=0> shape;
  // Negative-binomial overdispersion/precision parameter

  // ---------------------------------------------------------------------------
  // Subcounty effects
  // ---------------------------------------------------------------------------

  vector[S] z_subcounty;
  // Standardized subcounty effects

  real<lower=0> sigma_subcounty;
  // Between-subcounty standard deviation
}

transformed parameters {
  // ---------------------------------------------------------------------------
  // Non-centred subcounty random intercepts
  // ---------------------------------------------------------------------------

  vector[S] subcounty_effect =
    sigma_subcounty * z_subcounty;
}

model {
  // ---------------------------------------------------------------------------
  // Priors
  // ---------------------------------------------------------------------------

  b ~ normal(prior_b_mean, prior_b_sd);
  Intercept ~ normal(prior_intercept_mean, prior_intercept_sd);
  shape ~ exponential(prior_shape);
  z_subcounty ~ std_normal();
  sigma_subcounty ~ normal(0, prior_sigma_subcounty);

  // ---------------------------------------------------------------------------
  // Likelihood
  // ---------------------------------------------------------------------------

  if (!prior_only) {
    vector[N] eta =
      Intercept +
      subcounty_effect[subcounty_id] +
      log_offset +
      X[, 2:K] * b;

    if (N > 0) {
      Y ~ neg_binomial_2_log(eta, shape);
    }

    if (N_cens > 0) {
      vector[N_cens] eta_cens =
        Intercept +
        subcounty_effect[subcounty_id_cens] +
        log_offset_cens +
        X_cens[, 2:K] * b;

      for (n in 1:N_cens) {
        target += neg_binomial_2_lccdf(
          Y_cens[n] - 1 |
          exp(eta_cens[n]),
          shape
        );
      }
    }
  }
}

generated quantities {
  // ---------------------------------------------------------------------------
  // Training-data pointwise log likelihoods
  // ---------------------------------------------------------------------------

  vector[N] log_lik_obs;
  // Exact-observation log likelihoods

  vector[N_cens] log_lik_cens;
  // Censored-observation log likelihoods

  vector[N_month_lfo] log_lik_month;
  // Training log likelihood summed by calendar month

  // ---------------------------------------------------------------------------
  // Held-out pointwise log likelihoods
  // ---------------------------------------------------------------------------

  vector[N_test] log_lik_test;
  // Held-out exact log likelihoods

  vector[N_test_cens] log_lik_test_cens;
  // Held-out censored log likelihoods

  vector[N_month_lfo] log_lik_test_month;
  // Held-out log likelihood summed by calendar month

  // ---------------------------------------------------------------------------
  // Posterior predictive counts
  // ---------------------------------------------------------------------------

  array[N] int Y_rep;
  // Posterior predictive exact training counts

  array[N_cens] int Y_rep_cens;
  // Posterior predictive latent counts for censored training rows

  array[N_test] int Y_rep_test;
  // Posterior predictive exact held-out counts

  array[N_test_cens] int Y_rep_test_cens;
  // Posterior predictive latent counts for censored held-out rows

  // ---------------------------------------------------------------------------
  // Posterior expected counts
  // ---------------------------------------------------------------------------

  vector[N] mu_hat;
  // Posterior expected counts for exact training observations

  vector[N_cens] mu_hat_cens;
  // Posterior expected latent counts for censored training observations

  vector[N_test] mu_hat_test;
  // Posterior expected counts for exact held-out observations

  vector[N_test_cens] mu_hat_test_cens;
  // Posterior expected latent counts for censored held-out observations

  // Initialize monthly log-likelihood totals
  log_lik_month = rep_vector(0, N_month_lfo);
  log_lik_test_month = rep_vector(0, N_month_lfo);

  {
    // -------------------------------------------------------------------------
    // Exact training observations
    // -------------------------------------------------------------------------

    vector[N] eta =
      Intercept +
      subcounty_effect[subcounty_id] +
      log_offset +
      X[, 2:K] * b;

    if (N > 0) {
      for (n in 1:N) {
        mu_hat[n] = exp(eta[n]);

        Y_rep[n] =
          neg_binomial_2_log_rng(eta[n], shape);

        log_lik_obs[n] =
          neg_binomial_2_log_lpmf(
            Y[n] | eta[n], shape
          );

        log_lik_month[month_id[n]] +=
          log_lik_obs[n];
      }
    }

    // -------------------------------------------------------------------------
    // Right-censored training observations
    // -------------------------------------------------------------------------

    if (N_cens > 0) {
      vector[N_cens] eta_cens =
        Intercept +
        subcounty_effect[subcounty_id_cens] +
        log_offset_cens +
        X_cens[, 2:K] * b;

      for (n in 1:N_cens) {
        mu_hat_cens[n] = exp(eta_cens[n]);

        Y_rep_cens[n] =
          neg_binomial_2_log_rng(
            eta_cens[n], shape
          );

        log_lik_cens[n] =
          neg_binomial_2_lccdf(
            Y_cens[n] - 1 |
            exp(eta_cens[n]),
            shape
          );

        log_lik_month[month_id_cens[n]] +=
          log_lik_cens[n];
      }
    }

    // -------------------------------------------------------------------------
    // Exact held-out observations
    // -------------------------------------------------------------------------

    if (N_test > 0) {
      vector[N_test] eta_test =
        Intercept +
        subcounty_effect[subcounty_id_test] +
        log_offset_test +
        X_test[, 2:K] * b;

      for (n in 1:N_test) {
        mu_hat_test[n] = exp(eta_test[n]);

        Y_rep_test[n] =
          neg_binomial_2_log_rng(
            eta_test[n], shape
          );

        log_lik_test[n] =
          neg_binomial_2_log_lpmf(
            Y_test[n] | eta_test[n], shape
          );

        log_lik_test_month[month_id_test[n]] +=
          log_lik_test[n];
      }
    }

    // -------------------------------------------------------------------------
    // Right-censored held-out observations
    // -------------------------------------------------------------------------

    if (N_test_cens > 0) {
      vector[N_test_cens] eta_test_cens =
        Intercept +
        subcounty_effect[subcounty_id_test_cens] +
        log_offset_test_cens +
        X_test_cens[, 2:K] * b;

      for (n in 1:N_test_cens) {
        mu_hat_test_cens[n] =
          exp(eta_test_cens[n]);

        Y_rep_test_cens[n] =
          neg_binomial_2_log_rng(
            eta_test_cens[n], shape
          );

        log_lik_test_cens[n] =
          neg_binomial_2_lccdf(
            Y_test_cens[n] - 1 |
            exp(eta_test_cens[n]),
            shape
          );

        log_lik_test_month[
          month_id_test_cens[n]
        ] += log_lik_test_cens[n];
      }
    }
  }
}

I’m using Pathfinder for initialization as

init <- if (prior_only) NULL else model$pathfinder(data = sd, seed = 124, num_paths = 4, output_dir = "cmdstan_csv")

The error arises from the right-censored likelihood. The call to neg_binomial_2_lccdf(Y_cens[n] - 1 | exp(eta_cens[n]), shape) either evaluates the probability as numerically zero during initialization or causes an overflow in the internal hypergeometric_pFq calculation during sampling. Pathfinder finds a finite initial point, but Chain 1 subsequently fails and CmdStanR cannot retrieve the fit metadata.

Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can't start sampling from this initial value.
Path [3] :Initial log joint density = -571168.607012
Chain 1 Exception: Error in function boost::math::hypergeometric_pFq(double,double,double): numeric overflow (in '/var/folders/qf/bzkp2r5d0zz_5lxfrbnc1wdh0000gn/T/Rtmplay7A8/model-7de54f298a8e.stan', line 124, column 8 to line 126, column 10)

Is there a way to resolve this? Thanks in advance!

2

A similar underflow seems to have been reported and solved (by @sakrejda) for left-censored data and the neg_binomial_2_lcdf function:

Not sure if it matters or whether the fix was also applied to the lccdf function at the time as well?

3

Sorry about that. Our cdfs and ccdfs are generally not very stable as most of them work on the linear scale and take a log because nobody’s worked out robust and efficient expansions of them. So this is mostly doing what @LuuC suggests. If there’s something that should be patched in Stan, you can open an issue in stan-dev/math and we’ll try to fix it.

How did you conclude that Pathfinder found a finite initial point? I ask because it looks like the sampler is rejecting the initial value.

This isn’t relevant for fixing your initialization problems, but you can make a lot of this more efficient by vectorizing operations. For example, this

    if (N_cens > 0) {
      vector[N_cens] eta_cens =
        Intercept +
        subcounty_effect[subcounty_id_cens] +
        log_offset_cens +
        X_cens[, 2:K] * b;

      for (n in 1:N_cens) {
        target += neg_binomial_2_lccdf(
          Y_cens[n] - 1 |
          exp(eta_cens[n]),
          shape
        );
      }
    }

can just be

transformed data {
  matrix[N_cens, K - 1] X_cens_2 = X_cens[ , 2:K];
  array[N_cens] Y_cens_minus_1 = Y_cens;
  for (n in 1:N_cens) {
    y_cens_minus_1[n] -= 1;
  }
}

vector[N_cens] eta_cens
   = Intercept
     + subcounty_effect[subcounty_id_cens]
     + log_offset_cens
     + X_cens_2 * b;
target += neg_binomial_2_lccdf(
    Y_cens_minus_1 |
    Intercept +  subcounty_effect[subcounty_id_cens] + log_offset_cens  + ba *  X_cens_2K,
    shape);
4

This PR Evaluate the incomplete beta complement without cancellation in the lccdf and Student t cdf functions by avehtari · Pull Request #3433 · stan-dev/math · GitHub might fix this, at least it is improving neg_binomial_2_lccdf() a lot. Hopefully merged before the next release.

There are now 10 stan-dev/math PRs greatly improving our cdfs, lcdfs and lccdfs

5

While waiting for these updates to be merged and released, is there anything that could be manually coded as a custom function. I was looking at a PR about the NB lccdf function, but that is only for the special case when the data (count) equals zero.

6

Probably yes, but it is possible that even if the value would be finite, the autodiffed gradient might be non-finite. Many of the numerical stability fixes have been specifically for gradients and those can be fixed only in C++. If you need these before next release, it is possible to build CmdStan using specific github branches and get these updates (getting updates from different PRs before they are merged is a bit more complicated)