I have a multivariate model with two outcomes, Y_1 and Y_2:

- E[Y_1] is an average rate (proportion)
- E[Y_2] is an average amount

The goal requires E[Y_1] and E[Y_2] conditional on certain combinations of covariates, and averaged across others, but also E[Y_1 Y_2] conditional on those same combinations.

My understanding is that `marginaleffects::avg_predictions()`

can compute E[Y_1|\dots] and E[Y_2|\dots], but I don’t see how it can deliver E[Y_1Y_2|\dots]. (If it can, great!) Y_1 and Y_2 are both positive, and positively correlated, so the product E[Y_1|\dots] E[Y_2|\dots] will be biased low compared to E[Y_1 Y_2|\dots].

Wondering which of these approaches might seem best, or if there is a fourth option:

- Get
`avg_predictions()`

to yield the desired product - Modify the
`generated quantities`

block using`brm(..., stanvars = ...)`

(see below) - Sample from the paired
`draws()`

and calculate, then summarize, the products from the paired draws

Here is an example using a toy dataset and a bit of pseudo-ish code inside `stanvars()`

that doesn’t work, but hopefully gives the idea for approach #2:

```
require(brms)
require(cmdstanr)
bf_mpg <- bf(mpg ~ cyl)
bf_am <- bf(am ~ disp) + bernoulli()
mod <- brm(
bf_mpg + bf_am,
stanvars = stanvar(
scode = "real foo = mpg .* am;", # What should this be?
block = "genquant")
backend = "cmdstanr",
data = mtcars)
```

Solutions or ideas? Thank you in advance for your time and attention!