Maybe I am missing something (and maybe you resolved your issue in the meantime), but what would be wrong about fitting only the data without missingness and using posterior_predict(fit, data = orig_data %>% filter(is.na(Y), allow_new_levels = TRUE) to estimate the uncertainty you have about the missing Y values? This just takes the fitted uncertainty in your factor levels and draws the coefficient for the combinations of predictors not seen in the non-missing data using this uncertainty. This needs to assume that the unobserved combinations are in some sense “from the same population” as the observed ones, so it won’t help you if there is systematic bias in the unobserved. But it might be a good start…
Hope that helps!