It’s quite straightforward. Instead of using the mean of each parameter, use the full posterior (a vector of values) for each parameter to compute the (full posterior distribution of the) conditional effects. Then take the mean of that posterior distribution of conditional effects. This is what I mean by “take the mean of a function of the posterior distribution”. You can extract the posterior samples of each parameter directly from your stanfit object. See also the rstan documention (assuming you’re using stan, like you indicate above): Extract samples from a fitted Stan model — extract • rstan