Thank you both very much! I really appreciate your responses. Dr. Buerkner, great to hear from you.
For posterity, here is the code I used to implement a custom family.
zi_coi_zoib <-
custom_family("zi_coi_zoib",
dpars = c("mu", "zi", "coi", "phi"),
links = c("logit", "logit", "logit", "log"),
lb = c(0, 0, 0, 0),
ub = c(1, 1, 1, NA),
type = "real")
stan_funs <- "
real zi_coi_zoib_lpdf(real y, real mu, real zi, real coi, real phi) {
row_vector[2] shape = [mu * phi, (1 - mu) * phi];
if (y == 0) {
return log(zi);
} else if (y == 1) {
return log1m(zi) + log(coi);
} else {
return beta_lpdf(y | shape[1], shape[2]) + log1m(zi) + log1m(coi);
}
}
"
stanvars <- stanvar(scode = stan_funs, block = "functions")
posterior_predict_zi_coi_zoib <- function(i, prep, ...) {
mu <- brms::get_dpar(prep, "mu", i = i)
zi <- brms::get_dpar(prep, "zi", i = i)
coi <- brms::get_dpar(prep, "coi", i = i)
phi <- brms::get_dpar(prep, "phi", i = i)
shape1 <- mu * phi
shape2 <- (1 - mu) * phi
# Generate uniform random variables
u <- runif(prep$ndraws)
# Compute cumulative probabilities
cdf_zero <- zi
cdf_one <- zi + (1 - zi) * coi
# Initialize output vector
y <- numeric(prep$ndraws)
# Assign outcomes based on cumulative probabilities
y[u < cdf_zero] <- 0
y[u >= cdf_zero & u < cdf_one] <- 1
# For the remaining cases, draw from the beta distribution
beta_indices <- which(u >= cdf_one)
y[beta_indices] <- rbeta(length(beta_indices), shape1, shape2)
return(y)
}