# Probability of event at the individual level from Bayesian posterior mean distribution

This is the first time I've come across this type of problem. I really appreciate any comments or suggestions.

In a Bayesian model implemented through Stan, with populations as the units of analysis (e.g., hospitals or institutions) , I have estimated a continuous parameter denoted as $$\mu$$. In this context:

$$\mu_i \sim N(\theta, \tau^2)$$

where $$i$$ denotes the mean for the $$i^{th}$$ population (for $$i$$ = 1,...,$$k$$) and $$\tau^2$$ represents the between-population variance. Thus, $$\theta$$ represents the average of true means across populations.

Based on $$n$$ samples from the posterior mean distribution, I want to calculate the probability that a person/individual (across all populations) has a value greather than 50.

In other words, assuming that the value for the $$j^{th}$$ individual in the $$i^{th}$$ population is represented by:

$$\mu_{ji} \sim N(\mu_i, \sigma^2)$$

I want to calculate Pr($$\mu_{ji}>50$$)

where $$\sigma^2$$ represents the population-specific variance, assumed to be the same for all populations for simplicity.

To compute Pr($$\mu_{ji}>50$$), I calculated for each sample from the posterior distribution Pr($$\mu_{ji}>50$$) = 1-$$\Phi\left(\frac{50-\mu_{i}}{\sigma} \right)$$

and obtained the median across all $$k$$ samples

Is this a reasonable approach?

Any suggestions are welcome.

This is a hierarchical Bayesian model, hence $$\mu_{ij}$$ only depends on $$\mu_i$$, and $$\mathbb P(\mu_{ji}>50|D)=\mathbb E\left[ \mathbb P(\mu_{ji}>50\big|\mu_i)|D\right]$$ where $$D$$ denotes the actual sample (to be distinguished from $$\mu_{ij}$$, which is distributed from the predictive distribution). Now, $$\mathbb P(\mu_{ji}>50\big|\mu_i)=1-\Phi(\sigma^{-1}\{50-\mu_i\})$$ Therefore $$\mathbb P(\mu_{ji}>50|D)=1-\mathbb E\left[ \Phi(\sigma^{-1}\{50-\mu_i\}) | D \right]$$ which can be approximated by averaging over the simulated posterior sample.