neg_binomial_2_log_lpmf: Precision parameter is 0, Even After Changing Prior

[stefanom64]

Hi BIRDMAn devs 👋

We are analyzing microbiome count data using a custom negative binomial TableModel in BIRDMAn, which we find to be a very interesting and powerful tool for this kind of analysis — thank you for developing and maintaining it!

However, we’re encountering a persistent issue during sampling and would appreciate your insights or suggestions.

We are fitting a custom Stan model using neg_binomial_2_log_lpmf() for microbiome count data (samples × taxa), via the TableModel interface. The design includes covariates like Diet * Day + Plasma_Glucose + HSI, and we also account for sequencing depth as an offset.

Following q2-matchmaker issue #24, we updated the prior on the intercept beta_0 to ensure proper model compilation.

We're passing these priors when building the model:

param_dict = { "depth": np.log(filt_tbl.sum(axis="sample")), "B_p": 2.5, "inv_disp_sd": 3.0, "beta_0_mu": -4.54, "beta_0_sd": 1.47 } nb_model.add_parameters(param_dict)

Despite that fix, we consistently get this Stan error across chains:
`Exception: neg_binomial_2_log_lpmf: Precision parameter is 0, but must be positive finite! (in 'negative_binomial_model.stan', line 44)

During posterior predictive sampling in the generated quantities block, we also observed overflows when calling neg_binomial_2_log_rng(). To mitigate this, we manually clamp the lam_clr value in Stan code using:

generated quantities {
  array[N, D] int y_predict;
  array[N, D] real log_lhood;

  for (n in 1:N) {
    for (i in 1:D) {
      // Clamp lam_clr to prevent overflow in RNG
      real safe_lam = fmin(lam_clr[n, i], 13.8);  // exp(13.8) ≈ 1e6
      y_predict[n, i] = neg_binomial_2_log_rng(safe_lam, inv_disp[i]);
      log_lhood[n, i] = neg_binomial_2_log_lpmf(y[n, i] | lam_clr[n, i], inv_disp[i]);
    }
  }

I share with you the Stan code used for analysis:

data {
  int N;               // number of samples
  int D;               // number of taxa/features
  int p;               // number of covariates
  matrix[N, p] x;               // design matrix (Diet * Time + Plasma_Glucose + HSI)
  vector[N] depth;              // log sequencing depth
  array[N, D] int y;            // microbial counts

  // Priors
  real B_p;            // SD for beta coefficients
  real inv_disp_sd;    // SD for lognormal prior on inverse dispersion
  real beta_0_mu;               // Prior mean for intercepts
  real beta_0_sd;      // Prior SD for intercepts
}

parameters {
  row_vector[D-1] beta_0;                 // intercepts (excluding baseline taxon)
  matrix[p, D-1] beta_x;                  // covariate coefficients
  vector[D] inv_disp;            // inverse dispersion
}

transformed parameters {
  matrix[N, D-1] lam;                     // log-rates (excluding baseline)
  matrix[N, D] lam_clr;                   // full CLR log-rate including baseline

  lam = x * beta_x;                       // fixed effects
  for (n in 1:N) {
    lam[n] += beta_0 + depth[n];          // add intercepts + log-depth offset
  }

  lam_clr = append_col(to_vector(rep_array(0, N)), lam);  // first component fixed at 0
}

model {
  // Priors
  beta_0 ~ normal(beta_0_mu, beta_0_sd);           // well-informed prior on intercepts
  to_vector(beta_x) ~ normal(0, B_p);              // covariate priors
  inv_disp ~ lognormal(0, inv_disp_sd);            // dispersion

  // Likelihood
  for (n in 1:N) {
    for (i in 1:D) {
      target += neg_binomial_2_log_lpmf(y[n, i] | lam_clr[n, i], inv_disp[i]);
    }
  }
}

generated quantities {
  array[N, D] int y_predict;
  array[N, D] real log_lhood;

  for (n in 1:N) {
    for (i in 1:D) {
      // Clamp lam_clr to prevent overflow in RNG
      real safe_lam = fmin(lam_clr[n, i], 13.8);  // exp(13.8) ≈ 1e6
      y_predict[n, i] = neg_binomial_2_log_rng(safe_lam, inv_disp[i]);
      log_lhood[n, i] = neg_binomial_2_log_lpmf(y[n, i] | lam_clr[n, i], inv_disp[i]);
    }
  }
}
 Help requested

-Is there any way inv_disp could still be initialized to zero or become numerically unstable, despite using  and a lognormal prior?
-Should we further constrain it (e.g., )?
-Are there BIRDMAn-specific initialization tricks for custom TableModels that could help?

Thank you again for your support and for creating such a great tool!

Sincerely,
Stefano

[mortonjt]

It maybe worthwhile to initialize the dispersion factors by exploiting the mean vs variance relationship. This is described in the DEseq2 paper. It would involve computing means / variances of the counts across all samples, computing linear regression per feature, and initializing the dispersion factors from the regression coefficients. We found it made a difference when initializing the models

[stefanom64]

Hi @mortonjt,
thank you for your feedback. I will try that and I will update you!
cheers

[stefanom64]

Gi @mortonjt,
I tried implementing your suggestions, but I still encountered a few warning messages in the console during the fit_model() step.

At the very beginning of sampling, I get:
Exception: neg_binomial_2_log_lpmf: Precision parameter is 0, but must be positive finite! (in 'negative_binomial_model.stan', line 45, column 6 to column 78)
And later during sampling:
Exception: neg_binomial_2_log_rng: Random number that came from gamma distribution is 1.39434e+09, but must be less than 1073741824.000000 (in 'negative_binomial_model.stan', line 58, column 6 to column 75)
Here is how I initialized the dispersion values:

smoothed_f = np.clip(a0_f + a1_f / np.clip(means_f.values, 1e-12, None), 1e-8, None)
inv_disp_init = (1.0 / smoothed_f).astype(float)

Then I passed them to the model as follows:

inv_series = pd.Series(inv_disp_init, index=means_f.index)
inv_disp_init = inv_series.reindex(features).fillna(inv_series.median()).to_numpy()

inv_disp_init = np.asarray(inv_disp_init, dtype=float)
inv_disp_init = np.clip(inv_disp_init, 1e-8, None)

D = len(nb_model.feature_names)
p = len(nb_model.colnames)

init_dict = {
    "inv_disp": inv_disp_init,                   # shape: (D,)
    "beta_0": np.zeros((D - 1,), dtype=float),   # shape: (D-1,)
    "beta_x": np.zeros((p - 1, D - 1), dtype=float),
}

inits_per_chain = [init_dict] * 4

fit = nb_model.fit_model(
    method="mcmc",
    num_draws=500,
    mcmc_warmup=500,
    mcmc_chains=4,
    mcmc_kwargs={
        "inits": inits_per_chain,
        "show_console": True
    },
    seed=42,
)

Despite the warning messages, the rhat values are all around 1 and the effective sample sizes (ESS) are large (over 3000). So I’m wondering:

Are these warnings during initialization and sampling harmful?

Could it be that they’re just artifacts of some extreme values, but the model still fits reasonably?

Would it help if I shared the full Jupyter notebook in case you'd like to review the entire pipeline?

Thanks in advance!

Stefano