Hi, my model is as follows :

model =

bf(Y1 ~ IV1+…IVn+(1|rand|individual)+factor(time), family=bernoulli(logit)) +

bf(Y2 ~ IV1+…IVn+(1|rand|individual)+factor(time), family=bernoulli(logit)) +

bf(DV1 | subset(Y1) ~ IV1+…IVn+(1|rand|individual)+factor(time), family=gaussian()) +

bf(DV2 | subset(Y1) ~ IV1+…IVn+(1|rand|individual)+factor(time), family=gaussian())

Prior_int = c(

set_prior(“normal(0,5)”, resp = "Y1, class = “Intercept”),

set_prior(“normal(0,5)”, resp = “Y2”, class = “Intercept”),

set_prior(“normal(0,5)”, resp = “DV1”, class = “Intercept”),

set_prior(“normal(0,5)”, resp = “DV2”, class = “Intercept”) )

fit = brm (model + set_rescor(FALSE), data=DATASET,

warmup = 2000, iter = 4000, seed = 2021,

chains = 4, control = list(adapt_delta =.80),

cores = 4, prior = c(Prior_int), inits = 0)

)

What I want to do is to obtain posterior samples from different values of `IV1`

(e.g., IV1=0, IV1=50, IV1=100) with each of the outcomes (e.g., Y1,Y2,…) and put them in a plot to compare the fitted values.

I found that I could get population-level (fixed effects) posterior draws using “fixef(fit, summary=F)”.

And I want to know

- How could I get 95% credible interval when calculating the fitted values using fixef for all of the draws (as I know, it is not possible to use the ‘prob’ option when ‘summary=F’ is designated).
- I saw a posting that recommends to use the function ‘fitted’ , but when I use the function error comes out (see below)

**> fitted(fit_main, newdata=data.frame(IV1=xx), scale = “response”, summary=FALSE)**

*Error: Argument ‘resp’ must be a single variable name for models using addition argument ‘draw_ids’.* - how to draw a plot for comparing the fitted value of two outcomes (e.g., comparing the fitted value of 8,000 draws for Y1 and Y2 when IV1=0 in the same plot)

Thank you in advance for your help.