Sure. I am thinking about implementing it myself. I saw from this post that for making predictions on new data, I can sample values for each parameter (in my case they are alpha, alpha_i, theta, and pie) from their posterior distribution respectively and then sample yrep from the ZINB distribution with sampled alpha, alpha_i, theta, and pie. Is this how the posterior_predict function works? If yes, then I can add the constraints to the first step to customize the sampling process. I tried to look for the source code of posterior_predict but the posterior_predict.R file on Github doesn’t include all the details…