Bayesian (Stan) advantage for estimating small variance components in a multilevel model Imagine you have some data with one continuous predictor x, a categorical predictor A, and some dependent variable y. You have some amount (enough) data for several subjects S across all conditions. You fit a model:
y ~ x * A + (x*A|S)
using lme4 and brms (using Stan) in R. The lme4 model returns the following error "boundary (singular) fit: see help('isSingular')" meaning the model had problems estimating very small variance components. The Bayesian model will not usually have problems fitting the same model and Gelman suggests keeping 'random effects' with small variances in our models so that we can say that they have a small variance.
My question is: What is it specifically about Stan that allows it to fit models with small variance components with fewer problems than lme4? I'm not saying estimating these components is not more difficult or does not cause problems, I just mean that small variance components are not generally considered problematic for a Bayesian model fit using Stan. Is it:

*

*The Hamiltonian Monte Carlo moves around in such a way that boundaries (e.g. 0) don't cause as much of a problem?


*Priors somehow help the posterior estimate these small values?


*Random sampling itself offers an advantage over trying to find the single 'best' value?


*Some combination of the three? Something else?
 A: When estimates on the boundary of the parameter space (here zero for a variance) cause problems for frequentist methods that rely on asymptotic methods, this is not an issue with the model (likelihood) itself, but rather with how you are estimating parameters/fitting the model. In some circumstances (e.g. binomial data with 10 cases out of 100 for one factor level vs. 0 out of 100 for another factor level), there's even purely frequentist solutions (e.g. in the example I gave in brackets things like exact logistic regression). I don't think there's anything like that for random effects models, but perhaps there is in some specific circumstances.
In all of these cases Bayesian methods with proper prior distributions can help, too. These priors could either be very vague (but still a proper distribution) to not influence inference too much, could reflect what you know about the problem, and, yes, they could allow for the possibility of the parameter value being exactly 0. The effect of the prior distribution, will of course be to pull values away from zero.
The only thing that tends to cause issues for Hamiltonian Monte Carlo (and thus, Stan) would be having a discrete parameter such as a latent parameter for "Is variance 0 or not?". The issue is that the algorithm (in its standard version) needs derivatives with respect to the parameter, which would obviously be an issue. One can avoid such issues though, if one can integrate out the parameter or find an approximation (e.g. instead of a prior with a point mass at zero, just concentrate a lot of prior weight near zero at values that are for all practical intents and purposes the same thing as zero).
Without a proper prior distribution, you don't know that you would have a proper posterior distribution (and there's examples in the literature where people didn't notice that their MCMC sampler happily sampled away at a problem, for which it could not possibly provide proper MCMC samples).
(Pseudo-)random sampling from a posterior distribution is not a thing that is about a single point estimate vs. characterizing the distribution. It's mostly used, because Bayesian models of realistic complexity usually don't allow us to analytically calculate or with minimization capture the posterior distribution or even to approximate its shape well around its mode. You can still afterwards give some suitable point estimate (e.g. posterior median), it's just that the samples from the posterior distribution also let you easily describe the posterior distribution very nicely (and it does not matter whether it's some nice normal shape or anything like that, which the standard frequentist random effects model fit with lme4 would assume for the sampling distribution). Additionally, you gain a few bonuses like being able to easily describe the posterior for transformations of the model parameters (which is often tricky for maximum likelihood methods).
