I am using cmdstanr and the posterior package. I would like to compare the distribution of parameters \theta_1, \dots, \theta_7 conditional on different values for a set of other parameters \alpha_1, \dots, \alpha_{100}. What is the quickest way to do this in Stan?

Assuming that \alpha_1,... are also parameters that are estimates jointly with \theta_1,... I think this is something you will need to do outside of the Stan run in post-processing.

But I think a some additional background is needed to give a clearer and more certain response.

The parameters are estimated jointly. In the simplest case, say we only have one \alpha and one \beta and the model is y \sim \alpha \cdot x + \theta \cdot x. So here I would want to have p(\theta | \alpha = a_1) or p(\theta | \alpha = a_2) etc, i.e., the effect of x on y conditional on \alpha \in (a_1, a_2).

My actual model is more complex as it has more parameters and the functional relationship is not just linear. I would basically be looking for a command gather_draws(model_fit, parameter = theta, alpha = a_1). That is, get posterior draws of \theta given \alpha is a fixed and and specified value.

Does that make it clearer regarding what I would need?

I am unfortunately a data.table person, which means that I use subset_draws to extract parameters of interest, which I then put into a data.table, before I use melt, cast, and other functions to calculate what I need. So I can’t help much with gather_draws (which is from tidybayes) and the like.

Here is how you could do this with posterior and data.table:

The reason I asked about the joint estimation is that I was wondering how one would deal with the fact that one can’t calculate basic indices like Rhat for these conditional estimates. I also don’t see how one could use posterior or tidybayes functions to get statistics for conditional thetas. Lastly, I am uncertain how to deal with the fact that all alpha value will not have the same frequency in the posterior.

Ah yeah that is about the way I would have done it in dplyr. But then I would have the same Issue as you mention with the unequal frequency in the posterior.

One question regarding your code:

subset_draws("theta | a_1", regex = T)

So these are draws of theta conditional on a_1 or just joint draws of theta and a_1?

theta_a =
draws %>%
subset_draws("theta|a", regex = T) %>%
as_draws_df()

this gives you a data.frame with 2 columns, one for theta and one for a. And yes you will have unequal frequencies in the posterior, which is why I would be hesitant to do this.

One way around this could be to model theta as a function of a.