Hi all,
I’ve fitted a multi-level model in brms and I want to know the MAP parameter values. My purpose is sampling from the MAP random effects structure to get the best guess of the predictive distribution without any estimation uncertainty. There is already great functionality to simulated new datasets from the random effects from randomly drawn posterior samples using brms::predict and sample_new_levels = "gaussian". But, as I understand it, this incorporates the model uncertainty baked into the posterior samples. I’d like to do this for the MAP parameters.
I’ve already tried a couple of ways.
-
taking the MAP estimates of posterior samples using bayestestR::map_estimate. But when I simulate the predictive distribution by constructing the covariance matrix from the MAP sds and correlations, the resulting covariance matrix is not positive definite. I’m wonder whether this is because taking the MAP of each parameter value means that the combined parameters are no longer jointly-credible.
-
extracting the stanmodel from brms using mod[["fit"]]@stanmodel and using rstan::optimizing. But this throws an error about variables not existing.
I’ve also wondered whether I’m over-complicating things and it could be as simple as taking the posterior sample with maximum lp_ (if I interpret correctly that this would be the sample where the parameters match the model the most)?
Any direction appreciated.
Thanks for your time.
Callum
Why would you like to get MAP estimates? Using posterior means or medians are usually much more stable and overall favorable point estimates. The latter two are easiy to extract but just computing means or medians over the posterior draws.
Hi Paul, thanks for replying.
Maybe I’m misunderstanding and taking the means/medians is the best options. I’ll explain more specifically my issue (apologies if some of this is already clear from the first post).
My model is multi-level and predicts a normally distributed response. I’m specifically interested in the variability of the model response predictions. As far as I understand it, the fully predictive intervals have three influences: the response spread as given by sigma (within participant variability), the random effects (between participant variability), and the model/estimation uncertainty (given by the variation of parameters values across the posterior samples).
I’d like to assess how much these sources of variability contribute to the fully predictive intervals. So would like to get predictive intervals for the model’s best guess at the fixed effects (as the lowest variability), then add the best guess for the random effects, then add the model uncertainty.
The (fixed or random) effects + model uncertainty is easy because you just integrate across the posterior samples. But I’m puzzled as to how you get the best guess of the parameterisation, ignoring the model uncertainty.
I understand (perhaps wrongly) that each posterior sample’s parameter values are jointly-credible, in the sense that they have the parameter covariance structure baked into the sample (and probably for other reasons that I don’t appreciate). So you can simulate a new hypothetical dataset from each sample (i.e. the discussion here: https://github.com/paul-buerkner/brms/issues/191). Wouldn’t taking means or medians over the posterior draws mean that the resulting combination of values was no longer jointly-credible? So simulating a new dataset from the mean/median parameter values would be biased in some way?
This reasoning led me to trying the get the sample that represents the model’s best guess, which in turn led me on to trying to get the MAP values. But perhaps this is the wrong way to go?
Using MAP is the wrong conclusion, but your idea is reasonable in principle, if I understand it correctly. What is important (and what I understand you mean) is that all transformations need to be done for each individual draw and then, only after draws for the quantity of final interest are computed, you can do summaries, for example, means, median, or quantiles.