Hello, I have individual-level predictions produced by a non-parametric Bayesian model. Specifically I have, for each observation, a posterior distribution of the individual Relative Risks associated with a simulated change in a predictor (relative individual treatment effect, rITE).
My goal is to find subgroups in the rITE (I use a tree on the median posterior values), and then to evaluate the posterior distribution of the rITE in each subgroup. My problem is that these groups can be small and I’d like to regularize the estimates modeling them as a random effects. That is, in a brms like formula: rITE ~ (1|group) with a lognormal model.
My question is: how do model this considering that I have a distribution for each observation? I’d like to use the information in the whole distribution without simplifying it taking for example mean and std dev.
You have at least three options:
- Re-fit the original model and the new model jointly. If this is feasible, it is the cleanest solution.
- Obtain a parametric approximation to the posterior from the first model, and fit the second model as a hierarchical model wherein the true values are sampled from this parametric approximation. Note that it is important that the parametric approximation adequately capture the joint posterior, and independent approximations of the marginal posteriors for each quantity of interest may or may not be an adequate joint approximation.
- Iteratively fit the second model using as data many independent draws from the posterior of the first model. Combine the posterior chains from these many fits for inference. Unlike the previous two options, this option does not allow the second model to inform the posterior expectations from the first model, which can be a feature or a bug depending on your exact use case. One challenge of this method is that you might find that different parameterizations of the second model might be necessary to yield accurate computation over different posterior iterations of the first model. Thus, if you’re unlucky you might need to use multiple different parameterizations of the second model (ensuring that they encode exactly equivalent models). And if you’re really unlucky, you might find that some realizations from the first model’s posterior are difficult or impossible to reparametrize in such a way as to yield correct computation. If that happens, you can’t safely proceed.