Why does brms parameterize ZOIB model with zoi and coi?

Interfaces brms
1

In statistical papers I read about ZOIB models (e.g. liu and eugenio, 2016), if I am understanding correctly, it is parameterized with alpha = p(y=0), gamma = p(y = 1 | y != 0) and a beta distribution conditional on y !=0 and y != 1. In contrast, brms seems to parameterize the model with zoi = p(y = 0 or y = 1), coi = p(y = 1 | y = 0 or y = 1), and beta distribution.

I can’t imagine when that latter parameterization with the p(0 or 1) would be conceptually preferable. Is there some reason, perhaps related to model fitting, that I am missing?

2

Sorry this hasn’t been answered sooner. I don’t know brms, but these two alternatives are going to produce the same answers with different parameterizations. I think if you actually write out the probabilities for y = 0, y = 1, and y >= 2 it will make more sense.

Your approach takes ZI to be probability of zero inflation and COI to be probability of one inflation given there’s not zero inflation. That works out to something like this:

Pr[y = 0] = ZI + (1 - ZI - COI) * B(0 | ...)
Pr[y = 1] = COI + (1 - ZI - COI) * B(1 | ...)
Pr[y = n] = (1 - ZI - COI) * B(n | ...)

and the alternative with zero-or-one inflation ZOI and COI to be probability of 1 inflation if it’s 0 or 1 inflated, would be:

Pr[y = 0] = ZOI * (1 - COI) + (1 - ZOI) * B(0 | ...)
Pr[y = 1] = ZOI * COI + (1 - ZOI) * B(0 | ...)
Pr[y = n] = (1 - ZOI) * B(n | ...)

Seems like six of one, half dozen of the other to me. The brms approach does a bit less work when the value isn’t 0 or 1 and a bit more work when it is 0 or 1.

5 Likes
3

I don’t think I had a good reason back in the days to implement it with ZOI COI parameterization. I wanted to ensure that the sum of the probabilities does never leave the 0-1 interval as could be the case with ZI OI parameterization. But ZI - COI seems to work just fine too. If you want the latter in brms, you can make a custom family implementing it.

2 Likes
4

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)
}
8 Likes
5

Thank you for sharing your custom family code! I am sure it will be helpful also to others!