Estimating phi, zoi, and coi in a multivariate zoib

Hello! I am familiar with rstan and brms, but less so with zero-one inflated beta (zoib) models so the misunderstanding of this question might come from there.

I have found that I can run a univariate and multivariate zoib, great! However, I am running into an error when i try to estimate the phi, zoi, and coi parameters in the multivariate model. Any thoughts or tips would be appreciated!

This first block of code creates a toy data set, declares the univariate zoib formula and runs the model succesfully in brms

library(brms)
library(cmdstanr)

df = data.frame(y = c(0, 1,1,1, 0, 1),
                y2 = c(1, 0, 0, 0, 1, 0),
                x = c(1.9, 3.8, 3.0, 4.0, 2.0, 3.7),
                group = c(0, 0, 1, 1, 1, 0),
                person = c(1,2,3,4,5, 6))

#univariate zoib
zoib_uni_model <- bf(
  y ~ x + group + (1|person),
  phi ~ x + group + (1|person),
  zoi ~ x + group + (1|person),
  coi ~ x + group + (1|person)
)


fit_uni = brm(
  zoib_uni_model,
  data = df,
  cores = 4, 
  iter = 100,
  chains = 1,
  seed = 1,
  family = zero_one_inflated_beta(),
  backend = "cmdstanr")
)

To run the multivariate zoib succesfully I can use the following formula and call to brms

zoib_mv_model1 <- mvbf(
  y ~ x + group + (1|person),
  y2 ~ x + group + (1|person)

fit_mv = brm(
  zoib_mv_model1,
  data = df,
  cores = 4, 
  iter = 100,
  chains = 1,
  seed = 1,
  family = zero_one_inflated_beta(),
  backend = "cmdstanr")

**However, if i add back in the model specific parameters and run the following model I get the following error message: **

Error: The following variables can neither be found in 'data' nor in 'data2':
'phi', 'zoi', 'coi'

MVBF model with model specific parameters

#multivariate zoib
zoib_mv_model 2<- mvbf(
  y ~ x + group + (1|person),
  y2 ~ x + group + (1|person),
  phi ~ x + group + (1|person),
  zoi ~ x + group + (1|person),
  coi ~ x + group + (1|person))



fit_mv =  brm(
  zoib_mv_model2,
  data = df,
  cores = 4, 
  iter = 100,
  chains = 1,
  seed = 1,
  family = zero_one_inflated_beta(),
  backend = "cmdstanr")

Try specifying it like this:

zoib_m1 <- bf(y ~ x + group + (1|person), phi ~ x + group + (1|person), zoi ~ x + group + (1|person), coi ~ x + group + (1|person)) + zero_one_inflated_beta()
zoib_m2 <- bf(y2 ~ x + group + (1|person), phi ~ x + group + (1|person), zoi ~ x + group + (1|person), coi ~ x + group + (1|person)) + zero_one_inflated_beta()

fit_mv =  brm(
  zoib_m1 + zoib_m2,
  data = df,
  iter = 100,
  chains = 1)

Note that you can also get correlations between varying intercepts by doing (1|p|person). Also, I know this was just a toy data, but it seems odd to fit a zoib model to data that is just binary 0 and 1.

Thanks so much JD, will try this now.

Yeah, the data i mustered up was pre coffee , agree that zoib would not be a great choice for binary output data.

Re: ```(1|p|person), can you elaborate on what the p represents in the random effect structure?

See ?brmsformula, specifically:

Group-level terms

Multiple grouping factors each with multiple group-level effects are possible. (Of course we can also run models without any group-level effects.) Instead of | you may use || in grouping terms to prevent correlations from being modeled. Equivalently, the cor argument of the gr function can be used for this purpose, for example, (1 + x || g) is equivalent to (1 + x | gr(g, cor = FALSE)) .

It is also possible to model different group-level terms of the same grouping factor as correlated (even across different formulas, e.g., in non-linear models) by using |<ID>| instead of | . All group-level terms sharing the same ID will be modeled as correlated. If, for instance, one specifies the terms (1+x|i|g) and (1+z|i|g) somewhere in the formulas passed to brmsformula , correlations between the corresponding group-level effects will be estimated. In the above example, i is not a variable in the data but just a symbol to indicate correlations between multiple group-level terms. Equivalently, the id argument of the gr function can be used as well, for example, (1 + x | gr(g, id = "i")) .