I have two models. One is a hierarchical Bayesian model that estimates parameters $p_{i, j}$ for group $i$ and sex $j$. This model is set up hierarchically because of the inherent structure of the groups and sexes, and does not work well when the groups are estimated separately. Practically, it's set up in JAGS in R.

The other model is extremely complicated and is currently optimized via maximum likelihood estimation. This model takes in $p_j$ as known parameters (as well as many other knowns) and spits out the real information that I'm interested in for a single group $i$. This model is written in R and uses optim.

Generally it is not good practice to stick an estimated parameter into a model. I'd like to reparameterize the models so they're compatible with each other and merge them together so $p_j$ is estimated together with everything in the second model. My question is, practically how can I take an extremely complicated model and "add in" my Bayesian model? Could I somehow extract the likelihood value in the Bayesian model and multiply it to the overall likelihood in the larger second model?

  • $\begingroup$ Why do you want to "add" them..? $\endgroup$
    – Tim
    Commented Sep 22, 2017 at 14:21
  • $\begingroup$ The information of interest comes out of the second, complicated model. The parameter coming out of the first model feeds into the second one. But it would be preferable if all parameters were estimated together, not in this two-step model process. $\endgroup$
    – Peter
    Commented Sep 22, 2017 at 14:24
  • $\begingroup$ So why not building the "two-step" model all within the Bayesian framework..? $\endgroup$
    – Tim
    Commented Sep 22, 2017 at 14:27
  • $\begingroup$ Because the second model is incredibly long and complicated and it would take me probably about 3 weeks to properly implement it. I also didn't write the long and complicated one, so while I have the code and a short 1-page document describing the model, using and understanding someone else's code has a steep learning curve. Basically, I'm lazy, and I was hoping for an easier solution. $\endgroup$
    – Peter
    Commented Sep 22, 2017 at 14:57

1 Answer 1


Perhaps the following framework will be useful. In the first stage you have the posterior distribution for some unobserved parameters: \begin{equation} p(\theta|y) \end{equation} where $y$ denotes the data and $\theta$ denotes the parameters. (These would be the $p_{ij}$ according to the notation in the question, but I want to use $p$ to stand for probability densities.) In the second stage you have a deterministic function that takes $\theta$ as an input and returns the variable of interest $x$ as an output, \begin{equation} x = f(\theta) . \end{equation} Note that the uncertainty in $\theta$ introduces uncertainty into $x$. Let $\{\theta^{(r)}\}_{r=1}^R$ denote draws from $p(\theta|y)$. The corresponding draws of $x$ from $p(x|y)$ are $\{x^{(r)}\}_{r=1}^R$, where $x^{(r)} = f(\theta^{(r)})$.

The follow-up question might be: How to use the draws from $p(x|y)$? The answer depends on the "decision problem" you are trying to solve. For example, if your loss function were quadratic, then you would compute the posterior mean: \begin{equation} \widehat x = \int x\,p(x|y)\,dx \approx \frac{1}{R} \sum_{r=1}^R x^{(r)} . \end{equation} There is a shortcut which may or may not be very good (depending on how nonlinear $f(\theta)$ is). Let \begin{equation} \widetilde x = f(\widehat \theta) , \end{equation} where $\widehat \theta = \frac{1}{R} \sum_{r=1}^R \theta^{(r)}$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.