# Setting prior on inverse square root of negative binomial's scale parameter

I’m working on a set of hurdle Poisson models in `brms`, but my posterior predictive checks are indicating that there’s some overdispersion that isn’t being captured properly. I’d like to switch to a hurdle negative binomial, but I’m having trouble specifying the `shape` prior.

I’d like to place a half-normal prior on a transformation of `shape` (i.e., `1/sqrt(shape) ~ normal(0, 1)` as described in the prior choice rec. wiki), but I’m not sure of the best way to do this. Defining a custom response distribution is one possibility, but that seems like a good opportunity to introduce unnecessary errors into the code, so I’m hoping there’s a simpler option.

When trying to implement this as a custom brms family, I encountered trouble getting the link functions and/or bounds of the hurdle parameter to work properly. I’ve defined the custom distribution as:

``````hurdle_negbinomial_scaled <- brms::custom_family(
"hurdle_negbinomial_scaled",
dpars = c("mu", "disp","hu"),
lb = c(0, 0, 0),
ub = c(NA, NA, 1),
type = "int")
``````

Support for the three dpars should be (0, Inf) for `mu`, (0, Inf) for `disp` (the excess dispersion) and (0, 1) for `hu` (the hurdle probability). My model has predictors for both `mu` and `hu`; however, the generated stancode doesn’t seem to apply the inverse-link function to hu.

``````// generated stancode, from the model block
// likelihood including all constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = Intercept + rep_vector(0.0, N);
// initialize linear predictor term
vector[N] hu = Intercept_hu + Xc_hu * b_hu;
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += r_1_1[J_1[n]] * Z_1_1[n];
}
for (n in 1:N) {
// add more terms to the linear predictor
hu[n] += r_2_hu_1[J_2[n]] * Z_2_hu_1[n] + r_3_hu_1[J_3[n]] * Z_3_hu_1[n];
}
for (n in 1:N) {
// apply the inverse link function
mu[n] = exp(mu[n]);
// I think there should be a similar section for hu[n]
}
for (n in 1:N) {
target += hurdle_negbinomial_scaled_lpmf(Y[n] | mu[n], disp, hu[n]);
}
}
``````

Unsurprisingly, this results in an error

``````4: Rejecting initial value:
Chain 4:   Error evaluating the log probability at the initial value.
Chain 4: Exception: Exception: bernoulli_lpmf: Probability parameter is 3.87143, but must be in the interval [0, 1]  (in 'model2631d11348284_a482d4a5e2a72a5f6c214b59315d36db' at line 15)
(in 'model2631d11348284_a482d4a5e2a72a5f6c214b59315d36db' at line 179)
``````

When I replace my custom family for the default `hurdle_negbinomial` family, The generated code snippet is pretty much the same, except that it calls a version of the lpmf that handles the scaling itself (e.g., `hurdle_neg_binomial_log_logit_lpmf`).

Can anyone help me figure out what I’m doing wrong?

Postscript: Here’s the custom stan code I used, which was directly adapted from hurdle negative binomial code included w/ brms.

``````/* hurdle negative binomial log-PDF of a single response
* Args:
*   y: the response value
*   mu: mean parameter of negative binomial distribution
*   disp: inverse scaled shape parameter of negative binomial distribution
*   hu: hurdle probability
* Returns:
*   a scalar to be added to the log posterior
*/
real hurdle_negbinomial_scaled_lpmf(int y, real mu, real disp, real hu) {
real phi = 1/ disp^2;
if (y == 0) {
return bernoulli_lpmf(1 | hu);
} else {
return bernoulli_lpmf(0 | hu) +
neg_binomial_2_lpmf(y | mu, phi) -
log1m((phi / (mu + phi))^phi);
}
}
/* hurdle negative binomial log-PDF of a single response
* logit parameterization for the hurdle part
* Args:
*   y: the response value
*   mu: mean parameter of negative binomial distribution
*   disp: inverse scaled phi  parameter of negative binomial distribution
*   hu: linear predictor of hurdle part
* Returns:
*   a scalar to be added to the log posterior
*/
real hurdle_negbinomial_scaled_logit_lpmf(int y, real mu, real disp, real hu) {
real phi = 1/ disp^2;
if (y == 0) {
return bernoulli_logit_lpmf(1 | hu);
} else {
return bernoulli_logit_lpmf(0 | hu) +
neg_binomial_2_lpmf(y | mu, phi) -
log1m((phi / (mu + phi))^phi);
}
}
/* hurdle negative binomial log-PDF of a single response
* log parameterization for the negative binomial part
* Args:
*   y: the response value
*   eta: linear predictor for negative binomial distribution
*   phi phi parameter of negative binomial distribution
*   hu: hurdle probability
* Returns:
*   a scalar to be added to the log posterior
*/
real hurdle_negbinomial_scaled_log_lpmf(int y, real eta, real disp, real hu) {
real phi = 1/ disp^2;
if (y == 0) {
return bernoulli_lpmf(1 | hu);
} else {
return bernoulli_lpmf(0 | hu) +
neg_binomial_2_log_lpmf(y | eta, phi) -
log1m((phi / (exp(eta) + phi))^phi);
}
}
/* hurdle negative binomial log-PDF of a single response
* log parameterization for the negative binomial part
* logit parameterization for the hurdle part
* Args:
*   y: the response value
*   eta: linear predictor for negative binomial distribution
*   disp: inverse scaled phi  parameter of negative binomial distribution
*   hu: linear predictor of hurdle part
* Returns:
*   a scalar to be added to the log posterior
*/
real hurdle_negbinomial_scaled_log_logit_lpmf(int y, real eta, real disp, real hu) {
real phi = 1/ disp^2;
if (y == 0) {
return bernoulli_logit_lpmf(1 | hu);
} else {
return bernoulli_logit_lpmf(0 | hu) +
neg_binomial_2_log_lpmf(y | eta, phi) -
log1m((phi / (exp(eta) + phi))^phi);
}
}
// hurdle negative binomial log-CCDF and log-CDF functions
real hurdle_negbinomial_scaled_lccdf(int y, real mu, real disp, real hu) {
real phi = 1/ disp^2;
return bernoulli_lpmf(0 | hu) + neg_binomial_2_lccdf(y | mu, phi) -
log1m((phi / (mu + phi))^phi);
}
real hurdle_negbinomial_scaled_lcdf(int y, real mu, real disp, real hu) {
real phi = 1/ disp^2;
return log1m_exp(hurdle_negbinomial_scaled_lccdf(y | mu, phi, hu));
}

``````

Looks like a bug. Let me check.

It was indeed a bug. Should now be fixed in the github version of brms.

2 Likes