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)
```

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.