# Metropolis-Hastings in a Bayesian Hierarchical model

I am trying to estimate a Bayesian Hierarchical model using the random-walk Metropolis-Hastings algorithm. While in a non-Hierarchical model, the algorithm is staight-forward, I am not sure I am getting it right for the Hierarchical. Below I describe the algorithm for the Hierarchical model, as I have understood it.

Consider the following: Assume we have data from N subjects who complete a task with M trials. Based on this dataset, we want to estimate a very simple Bayesian Hierarchical model with two parameters, $r$ and $s$. We make the assumption that $r_i\sim N(\mu_r,\sigma_r)$ and $s_i\sim N(\mu_s,\sigma_s)$, for each subject $i$, with population parameters $\mu_r,\sigma_r, \mu_s, \sigma_s$. We also assume that $\mu_r\sim N(0,1),\mu_s\sim N(0,1)$ and $\sigma_r\sim Unif(0,3), \sigma_s\sim Unif(0,3)$.

For the Metropolis-Hastings algorithm, I start with some initial values for $\mu_r,\sigma_r, \mu_s, \sigma_s$. Based on these values, I generate random samples for $r$ and $s$ for all the subjects. Using these values, I calculate the total posterior, as the sum of the log-likelihood of the data for every subject $i$, given $r_i$ and $\sigma_i$ plus the log-likelihood of the individual posteriors, plus the log-likelihood of the posterior for the population parameters.

Then I start the iterations of the chain. At each iteration, I generate a proposal for all $\mu_r,\sigma_r, \mu_s, \sigma_s$ by adding a random normal variable $N(0,k)$ to the values of the previous iteration (k is the tuning parameter). I then generate random samples for $r$ and $s$ for all the subjects, based on the proposal values. Then, I loop over all subjects and calculate the individual log-likelihood +log posterior for both the new and old values. Then with probability equal to the ratio of the likelihoods, I keep the proposal values, otherwise I keep the old values.

As a final step, I compare the total log-posterior (all subjects+log of posteriors for the population parameters) for both the old and the new values. Again, with probability equal to the ratio of the two, I keep the proposal, otherwise I keep the old values.

Then repeat the process for all the iterations.

I have mainly two questions, with the obvious one: Does it make any sense or I got it totally wrong?

1.For the subject level parameters, do I need to suggest a proposal by adding a random variable as in the population parameters or is it fine as is described above?

2.Is it correct to check the ratio for both the individual posteriors and the population ones?

• Without the complete maths formulas that describe what the code is doing, rather than doing it by words only, it is difficult to reply. However, adding a final accept/reject step sounds incorrect. To have two acceptance steps you need to make two proposal steps. Aug 12, 2018 at 19:48

If $x_i$ denotes the observation for subject $i$, and if $f(x_i|r_i,s_i)$ denotes the density associated with it, assuming independence between the subjects, the likelihood is$$\prod_i f(x_i|r_i,s_i)$$and the posterior is$$p(\mu_r,\sigma_r, \mu_s, \sigma_s) \prod_i f(x_i|r_i,s_i)\times p(r_i|\mu_r,\sigma_r)p(s_i|\mu_s, \sigma_s)$$Making a single proposal on the whole vector$$(r_1,s_1,...,\mu_r,\sigma_r, \mu_s, \sigma_s)$$with density$$q(\mu_r,\sigma_r, \mu_s, \sigma_s|\mu_r^-,\sigma_r^-, \mu_s^-, \sigma_s^-)\prod_i p(r_i|\mu_r,\sigma_r)p(s_i|\mu_s, \sigma_s)$$where $(\mu_r^-,\sigma_r^-, \mu_s^-, \sigma_s^-)$ denotes the previous values leads to a Metropolis-Hastings acceptance ratio of$$\frac{p(\mu_r,\sigma_r, \mu_s, \sigma_s) \prod_i f(x_i|r_i,s_i)\times p(r_i|\mu_r,\sigma_r)p(s_i|\mu_s, \sigma_s)}{p(\mu_r^-,\sigma_r^-, \mu_s^-, \sigma_s^-) \prod_i f(x_i|r_i^-,s_i^-)\times p(r_i^-|\mu_r^-,\sigma_r^-)p(s_i^-|\mu_s^-, \sigma_s^-)}$$ $$\times \frac{ q(\mu_r^-,\sigma_r^-, \mu_s^-, \sigma_s^-|\mu_r,\sigma_r, \mu_s, \sigma_s)\prod_i p(r_i^-|\mu_r^-,\sigma_r^-)p(s_i^-|\mu_s^-, \sigma_s^-) }{q(\mu_r,\sigma_r, \mu_s, \sigma_s|\mu_r^-,\sigma_r^-, \mu_s^-, \sigma_s^-)\prod_i p(r_i|\mu_r,\sigma_r)p(s_i|\mu_s, \sigma_s)}$$or$$\frac{p(\mu_r,\sigma_r, \mu_s, \sigma_s) \prod_i f(x_i|r_i,s_i)}{p(\mu_r^-,\sigma_r^-, \mu_s^-, \sigma_s^-) \prod_i f(x_i|r_i^-,s_i^-)}$$accounting for the symmetry in the proposal $q(\cdot|\cdot)$. This is the only ratio involved in the algorithm.