That is IMHO reasonable. It is definitely almost exactly the thing I did when evaluating the effects in Markov model in my recent work and I couldn’t figure out anything better. Few changes I would suggest:
- I think it is more useful to not actually compute Y, but only the probabilities P(Y) of observing each Y value at each time - any quantity of interest can be computed easily from this distribution and you avoid some issues (like impossibility to convert to log-odds if some event has low probability and thus is never observed in your particular Y). Maybe this is what you discuss in the adendum? (I don’t think I undestood the addendum completely)
- I don’t see how step 5 is useful - if you have enough posterior draws than IMHO computing just one simulation per draw is enough. (the consideration would be IMHO different if draws were particularly expensive and you would have much smaller number of posterior draws).
The obvious problem is that with this algorithm you get the effect estimate for exactly the distribution of patient covariates you have observed - which might (and might not) be desirable. If I remember correctly, I think @stablemarkets advised in his StanCon talk to also do a “Bayesian bootstrap”, i.e. reweight the patient-level differences in each sample by a sample drawn from Dirichlet(1,1,1,...,1). This way you generalize beyond the set of patients you have actually observed to the hypothethical population your patients were sampled from (hope I am not misrepresenting the argument). This should IMHO be quite similar to resampling the patients with replacement in each sample, but I am not 100% sure about this. I would however be surprised if this bootstrap made a notable difference when you have 1000 patients.
The problem actually seems similar to the model I used for some Covid data (described at: Fitting HMMs with time-varying transition matrices using brms: A prototype)
Best of luck with your model!