Generate Posterior predictive distribution at every step in the MCMC chain for a hierarchical regression model

Question

I'm trying to fit a Bayesian Hierarchical regression model with a random correlated coefficients using R ,I'm using data having 160 groups (schools) to fit a model of math score as a function of one regressor .

The within group sampling model can be expressed as :

$$Y_j\sim MVN(X_j\beta_j,\sigma^2I)$$

where $Y_1...,Y_m$ the group specific data vectors being conditionally independent given $\beta_1...\beta_m$ . To describe the across- group heterogeneity I use the normal Hierarchical model , so that :

$$ \beta_1...\beta_m \sim i.i.d MVN(\theta,\Sigma)$$ where the values of $\theta$ and $\Sigma$ are fixed but unknown parameter to be estimated , so that we can write the between group sampling model as :

$$\beta_j=\theta+\gamma_j$$ where $$\gamma_1,...,\gamma_m \sim i.i.d MVN(0,\Sigma)$$ So more generally
$$ Y_i{_j}=\beta_j^{T}x_i{_j}+\epsilon_i{_j}$$ $$ =\theta^{T}x_i{_j} +\gamma_j^{T}x_i{_j} + \epsilon_i{_j}$$

Where $\theta$ is referred to as "fixed " effect and $\gamma_j$'s are called the random effect .Then the posterior distribution to be computed is : $$p(\beta_1,...,\beta_m,\theta,\Sigma,\sigma^2|X_1,...X_m,y_1,...,y_m)$$

I used semi-conjugate priors to approximate the posterior using MCMC (Gibbs sampler) so that $$\theta \sim MVN(\mu_0,\Lambda _0)$$ $$\Sigma\sim inverse.Wishart (\eta_0,S_0^{-1})$$ $$\sigma^2\sim inverse.Gamma(\nu_0/2,\nu_0\sigma_0/2)$$

Then according to Gelman and Hoff the full conditional distributions will be as follows :

The full conditionals of $\beta_j$'s $$(\beta_j,y_j,X_j,\theta,\Sigma)\sim MVN $$ with

$$Var(\beta_j,y_j,X_j,\theta,\Sigma)= (\Sigma^{-1}+{X_j}^{T}X_j/\sigma^2)^{-1}$$ and $$E(\beta_j,y_j,X_j,\theta,\Sigma)= (\Sigma^{-1}+{X_j}^{T}X_j/\sigma^2)^{-1}(\Sigma^{-1}\theta+ {X_j}^{T}y_j/\sigma^2)$$

The full conditional of $\theta$ is $$(\theta|\beta_1,...,\beta_m,\Sigma)\sim MVN(\mu_m,\Lambda _m)$$ with $$\Lambda_m=(\Lambda_0^{-1}+m\Sigma^{-1})^{-1}$$ and $$\mu_m=\Lambda_m(\Lambda_0^{-1}\mu_0+m\Sigma^{-1}\bar{\beta})$$ where $$\bar{\beta}=\frac{1}{m}\beta_j$$

The full conditional of $\Sigma$ is $$\Sigma|\theta,\beta_1,...,\beta_m \sim W^{-1}(\eta_0+m,[S_0+S_{\theta}]^{-1})$$ where $$S_{\theta}=\sum_{j=1}^{m}(\beta_j-\theta)(\beta_j-\theta)^{T}$$ That means $S_{\theta}$ depends on $\theta$ and must be recomputed each time $\theta$ is updated .

Finally the full conditional of $\sigma^2$ is :

$$\sigma^2\sim G^{-1}([\nu_0+\sum n_j]/2,[\nu_0{\sigma^2{_0}}+SSR]/2)$$ where $$ SSR =\sum_{j=1}^{m}\sum_{i=1}^{n_j}(y_i{_j}-\beta_j^{T}x_i{_j})^2$$ That also means that SSR must be recomputed in each scan of Gibbs Sampler before $\sigma^2$ is updated .

Now my first question is , how can I generate the posterior predictive distribution based on the data ? is that means that I have to generate posterior predictive distribution for each group based on the data within the group , such that : $$p(\tilde{y_j}|y_j)=\int p(\tilde{y_j}|\beta_j,\theta,\Sigma,\sigma^2)p(\beta_j,\theta,\Sigma,\sigma^2|y_j)d\beta_j,\theta,\Sigma,\sigma^2$$

If that right , that means we will generate a posterior predictive distribution for each group , and At every step in the MCMC chain, we will use the parameter values at that step to randomly generate a simulated value of $y_j$ from the model for each group .

My second question is , If we ignore the multilevel model and run a simple bayesian linear regression model and generate a posterior predictive distribution form this simple model , how can we compare the posterior predictive distributions of both models , since in the multilevel model we have generated 160(number of groups) posterior predictive distributions ?

Answer 1

Sorry for my sloppy answer in advance.

You would just predict the value und sample from it in your Gibbs sampler. If you write your sampler yourself, otherwise I would suggest to use JAGS.

You can predict for each of the 160 groups individually, otherwise you would use the prior to predict, that you estimated from this analysis. I thought posterior predictive is the probability of p(y+1|X,Y,w) you have not 160 predictive distributions.

I think you don't compare posterior predictive distributions, you often want to compare the prior of the hierarchical model, to the estimated weights of the non-hierarchical model. Because the assumption is that effects get estimated "better" when modeling heterogeneity. Of course, you can predict unseen data and compare MSE and other stuff, too.