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 ?


1 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.


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.