neg_binomial_2_log_lpmf: Precision parameter is 0, Even After Changing Prior

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](https://github.com/flatironinstitute/q2-matchmaker/issues/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

No repository file or test is named. Start by reproducing the custom TableModel Stan model and inspecting parameter construction around add_parameters(param_dict), then trace initialization and generated-quantities behavior for inv_disp. Done means identifying a reproducible cause and documenting or implementing a repository-appropriate fix.

Written by the indexing model from the issue text.