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"),
links = c("log", "log", "logit"),
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));
}