Hi,
sorry for note getting to you earlier. Did you manage to move forward in the meantime?
A set of hints to debug models is at Divergent transitions - a primer. Particularly, I would try, if you can fit the Poisson part as a standalone model.
One thing that looks suspicious to me is the b0[i]/(1-b2[i]) term as this is not continuous in b2[i] - i.e. the value flips from “very large positive” to “very large negative” as b2 crosses 1. This is likely to make Stan (and any other sampler really) unable to traverse this region in the parameter space, so provided that there is at least one i for which there is non-negligible posterior mass of b2[i] close to 1 from either side, this IMHO should cause problems. Alternatively, I can see how that would lead to bimodal posterior for b2[i] (and other coefficients) as values close to 1 could lead to very low posterior density while some values a bit further on either side could be compatible with the data. Once again, bimodal posteriors will frustrate almost any sampler.
I do not understand the math of the model enough to judge to what extent can we say a priori that some of the values are non-sensical and exclude them. Alternatively, one could try to reparametrize the model so that linear predictors for both outcomes are smooth with respect to all the parameters, but I don’t immediately see such a parametrization. Maybe treating the cases b2[i] > 1 and b2[i] < 1 separately, making the model a mixture could also help.
It is unfortunately not uncommon that when people see problems when reimplementing a model in Stan that it turns out that the original sampler also didn’t work for the problem, but it failed silently while Stan fails loudly.
Best of luck with the model.