As someone who just implemented hierarchical probabilistic PCA in stan, I agree that it takes finesse, but it is no means impossible. Doing this sort of work efficiently in stan seems to require a some degree of understanding about how the sampler works. It also may require really thinking through your model. It saves you from deriving your own conditional distributions and writing a gibbs sampler, but you're going to have to do some analysis if you want to fit models of certain complexity.

KL divergence minimization (variational inference) is typically a weak approximation to the model you specified. I have seen it produce inferences on simulated which are just plain wrong. These "wrong" models are still often good predictors, so whether variational inference will work well for you depends on whether you care about making valid inferences or just doing prediction.