Extract zero inflation and one inflation coefficients from brms ZOIB model

Interfaces brms
1

Is there a way to extract the zero-inflation and one-inflation parameters from a ZOIB model?

i have seen this: Interpretation of Zero-One Inflated Beta Regression Parameters

But am still unsure how to do it.

Here are a couple examples of the model specifications I am running.


f <- bf(
  neg ~ prpc*phase + 
    (1|pid) + 
    (1|vignette_num),
  phi ~ prpc*phase + 
    (1|pid) + 
    (1|vignette_num),
  zoi ~ prpc*phase + 
    (1|pid) + 
    (1|vignette_num),
  coi ~ prpc*phase + 
    (1|pid) + 
    (1|vignette_num)
)


f <- bf(
  rb ~ prpc + 
    (1|pid) + 
    (1|vignette_num),
  phi ~ prpc + 
    (1|pid) + 
    (1|vignette_num),
  zoi ~ prpc + 
    (1|pid) + 
    (1|vignette_num),
  coi ~ prpc + 
    (1|pid) + 
    (1|vignette_num)
)
  • Operating System: macos 14.6
  • brms Version:2.21.0
2

Hello! In the formulas you have there, the zero-one inflation (zoi) and conditional one inflation (coi) terms are both modeled quantities. So I think you first need to get the expected probabilities of zoi and coi:

# Let fit be a fitted {brms} model
# Let newdata be a data frame with predictor values of interest
zoi <- posterior_epred(fit, newdata, dpar = "zoi")
coi <- posterior_epred(fit, newdata, dpar = "coi")

I think you can then calculate the probabilities of one inflation (oi) and zero inflation (zi) as follows:

oi <- zoi * coi
zi <- zoi * (1 - coi)

My understanding is that zoi is the probability of there being an ‘inflated’ 0 or 1 and coi is the probability that this ‘inflated’ number is 1 (hence conditional). This in turn means that the conditional zero inflation probability should be 1 - coi. So the probability of one inflation overall is oi = zoi * coi (same idea for zi). Does this make sense and/or work for you?

2 Likes
3

Thank you kthayashi! I am specifically interested in understanding what is the effect of the predictor prpc on the odds of getting a 0 (versus >0) and the effect of prpc on the odds of getting a 1 (versus <1). Is there a way I can get that?

4

Sorry, I’m not totally sure how to formalize those estimands in the context of this model. Your quantities of interest (prob or odds of 0, prob or odds of 1) are not encoded directly as parameters in this model, I think. There might be a way to mathematically derive an expression for what you want, but this is beyond me.

If your inference does not need to rely on having exact coefficients for your effects, I would personally plot the effects of interest by expanding the above code like this:

# Predict zoi and coi at varying values of x1 and fixed x2
# Replace x1 and x2 with your own predictors with appropriate values
newdata <- data.frame(x1 = -3:3, x2 = 0)
zoi <- posterior_epred(fit, newdata, dpar = "zoi")
coi <- posterior_epred(fit, newdata, dpar = "coi")

# Compute oi and zi
oi <- zoi * coi # Pr(y = 1)
zi <- zoi * (1 - coi) # Pr(y = 0)
noz <- 1 - (oi + zi) # Pr(0 < y < 1)
oi_mean <- apply(oi, 2, mean)
oi_qi <- apply(oi, 2, quantile, c(0.025, 0.975))
zi_mean <- apply(zi, 2, mean)
zi_qi <- apply(zi, 2, quantile, c(0.025, 0.975))
noz_mean <- apply(noz, 2, mean)
noz_qi <- apply(noz, 2, quantile, c(0.025, 0.975))

# Plot probabilities (convert to odds of log-odds if desired)
par(mfrow = c(1, 3))
# oi
plot(x = NULL, y = NULL, xlim = c(-3, 3), ylim = c(0, 0.5), xlab = "x1", ylab = "Pr(y = 1)")
lines(newdata$x1, oi_mean, type = "b", pch = 16)
lines(newdata$x1, oi_qi[1, ], type = "b", lty = 2)
lines(newdata$x1, oi_qi[2, ], type = "b", lty = 2)
# zi
plot(x = NULL, y = NULL, xlim = c(-3, 3), ylim = c(0, 0.5), xlab = "x1", ylab = "Pr(y = 0)")
lines(newdata$x1, zi_mean, type = "b", pch = 16)
lines(newdata$x1, zi_qi[1, ], type = "b", lty = 2)
lines(newdata$x1, zi_qi[2, ], type = "b", lty = 2)
# noz
plot(x = NULL, y = NULL, xlim = c(-3, 3), ylim = c(0.5, 1), xlab = "x1", ylab = "Pr(0 < y < 1)")
lines(newdata$x1, noz_mean, type = "b", pch = 16)
lines(newdata$x1, noz_qi[1, ], type = "b", lty = 2)
lines(newdata$x1, noz_qi[2, ], type = "b", lty = 2)

Screenshot 2024-08-18 at 12.44.24 PM
Screenshot 2024-08-18 at 12.44.24 PM1688×550 81.5 KB

In your case, you could do something like this for prpc as I did here for x1. Below is code to simulate this example:

R code 
library(brms)
library(extraDistr)

# Simulate predictors (on unit scale)
n <- 1000
x1 <- rnorm(n, 0, 1)
x2 <- rnorm(n, 0, 1)

# Set 'true' parameters (a = intercept, b = slope)
a <- 0
b_x1 <- 0.2
b_x2 <- -0.3
b_x12 <- 0.1
phi_a <- 2
zoi_a <- -2
zoi_b_x1 <- 0.3
zoi_b_x2 <- 0.1
zoi_b_x12 <- 0
coi_a <- 1
coi_b_x1 <- 0.2
coi_b_x2 <- 0
coi_b_x12 <- 0

# Simulate zero-one inflation
mu_zoi <- zoi_a + zoi_b_x1*x1 + zoi_b_x2*x2 + zoi_b_x12*x1*x2
p_zoi <- inv_logit_scaled(mu_zoi)
is_zoi <- rbinom(n, 1, p_zoi)

# Simulate ones (conditional on zero-one inflation)
mu_coi <- coi_a + coi_b_x1*x1 + coi_b_x2*x2 + coi_b_x12*x1*x2
p_coi <- inv_logit_scaled(mu_coi)
is_coi <- rbinom(n, 1, p_coi)

# Compute phi
mu_phi <- phi_a
phi <- exp(phi_a)

# Simulate outcome (y) with zero-one inflation
mu <- a + b_x1*x1 + b_x2*x2 + b_x12*x1*x2
p <- inv_logit_scaled(mu)
y <- rprop(n, size = phi, mean = p)
y[is_zoi == 1 & is_coi == 1] <- 1
y[is_zoi == 1 & is_coi == 0] <- 0

# Prepare data
d <- data.frame(x1, x2, y)

# Prepare formula
f <- bf(
  y ~ x1*x2,
  phi ~ 1,
  zoi ~ x1*x2,
  coi ~ x1*x2,
  family = zero_one_inflated_beta()
)

# Fit model
fit <- brm(formula = f, data = d, chains = 4, cores = 4)

# Check model summary
fit

I hope this helps somewhat! Someone else might be able to provide you with a better solution.

1 Like