I'm fitting a Bayesian Hierarchical logistic model with a random intercept, having data at two levels: subjects and hospital,




For the "subjects" equation I have 3-5 variables, and for the "group" ones, 5-8 variables. I've fitted 6 models including or not some variables, adding or not random effect, and adding or not a model for the intercept.

My question is about the posterior predictive checking:

For each model (in the picture, only the 6th),in black there is the proportion of 1's in the response variable y (is a status variable with only 2 possible values 0 or 1) for each centre; the red points come from:

  1. Sample from the posterior $p(\beta's,a's,\sigma_{a}|y)$
  2. Simulate from the posterior predictive $p(y^{rep}|\beta's,a's,\sigma_{a},y)$.
  3. From the posterior predictive I get $y^{rep}_{ij}\in[0,1]$, but my data $y$ consists only of 0's & 1`s.
  4. If $y^{rep}_{ij}\geq 0.5 \rightarrow y^{rep}_{ij}=1$ else $y^{rep}_{ij}=0$. Then for each centre $j$ I have a vector of 0's and 1's (as in original data $y$) and I compute the proportion of 1's and I get 1 red point, then I repeat 30 times to get 30 red points for each centre.

Is the 4th step statistically correct? If not, there is a way to go from the $p_{ij}$ to the data $y$ in order to compare?

Thank you in advance.

  • $\begingroup$ Can you please write down the posterior predictive that you are using? I presume the likelihood $p(y | .) \sim Bern(p_{ij})$? If so then the posterior predictive should also have a binary domain. I suspect your posterior predictive is about $p_{ij}$ not $Y^*$. $\endgroup$ – user130957 Sep 13 '16 at 18:56

Your Answer


Browse other questions tagged or ask your own question.