Bayesian hierarchical modeling of a two-stage process model (Planning stages)

Hey everyone,

I’m a new Stan user, and although I’ve read through the user manual, done a tutorial, and looked at a few sample programs, I want to make sure I understand how to properly specify a model in Stan for Bayesian hierarchical modeling.

There is no code to share yet because the questions I have are of a very fundamental nature, and I want to ensure that I’m not wasting my lab’s time jumping into things without knowing what I’m doing.

I’m working with a two-stage process model, which has already been specified before in MATLAB and python for parameter estimation by MLE. Now we want to move on to Bayesian parameter estimation, specifically using hierarchical modeling because the experiments we perform are, by their nature, data poor.

Is the correct way to implement our model as a Bayesian hierarchical model to write its likelihood function (in the functions section), specify that the individual-level data is assumed to be drawn from the distribution defined by that function given some parameters (in the model section), and to specify that those parameters are drawn from a group-level distribution that has its own hyperparameters that also need to be fit (in the model section)?

I appreciate any help or feedback! Thanks!


If your model requires a custom likelihood function, then yes the functions block is the right place to put it. A walk-through of how to do this that I like can be found here. It’s a re-implementation of the exponential distribution, so fairly familiar/simple to keep the focus on Stan implementation rather than other details.

Notice that I said *re-*implementation there. Just in case you’ve missed it, Stan has a huge number of distributions already implemented, tested/validated, and optimized for stability and efficiency. These are documented in the Stan Functions Reference. Depending on what a “two-stage process model” means in your discipline, a lot of the heavy lifting may have already been done for you.

Yes and yes. You’ll use a sampling function here that increments the target log density. This block looks the most like the math-stats way of writing a BHM.

You can use the data block to pass in values for your (hyper)priors. The transformed data and transformed parameters can also be helpful if you are expecting to ingest data or priors in a format that needs to be transformed by a function to interact with the core of your model.

Best of luck getting started. The folks on this forum have been super helpful as I’ve learned to use Stan. Please don’t hesitate to ask if you get stuck.


Thank you for your reply! Knowing I’m on the right track is very helpful!