Firstly apologies, I have fairly limited mathematical skills so there is a good chance that my question is simple or obvious.
I have a model in which I want to calculate the conditional expectation of the data given some higher-level model parameters, in order to gain an intuition on how the higher level parameters relate to the data. This would be done by marginalising out the intermediate parameters.
It was suggested to me to use nested iterated expectations to get to this point but I am having difficulty understanding exactly how to do this. Here is the model for clarity, though I am happy to deal with arbitrary distributions and (attempt to) work through the math myself. Here $X$ is a matrix of integers, $i\in\{1,\dots, q\}, j\in\{1,\dots,r\}$. $$ X_{ij} \sim \text{Multinomial}(p_{\cdot j}) \\ p_{\cdot j} = (p_{ij}, \dots, p_{gj}) \sim \text{Dirichlet}(\alpha_{\cdot j}) \\ \alpha_{\cdot j} = (\alpha_i, \dots, \alpha_g) \\ \alpha_{ij} \sim \text{log-Normal}(\mu_i, \sigma^2_i) $$
In this case I would like to know the form of $\mathbb{E}\left[X_{ij} | \mu_i, \sigma^2_i\right]$. The "obvious" or brute-force way to do this would (I think) be to write the joint distribution $\Pr(X_{ij}, p_{ij}, \alpha_{ij}, \mu_i, \sigma^2_i)$ and integrate $\mathrm{d}p$ $\mathrm{d}\alpha$ but this is likely intractable.
My understanding is that the nested iterated expectation approach would amount to evaluating $$ \mathbb{E}\left[X_{ij}\right] = \mathbb{E}\left[\mathbb{E}\left[X_{ij}|p_{ij}\right]\right] \text{where} \\ \mathbb{E}\left[X_{ij}|p_{ij}\right] = \mathbb{E}\left[\mathbb{E}\left[p_{ij}|\alpha_{ij}\right]\right] $$ ...and so on. Is this a correct intuition for this approach? Alternatively, is there another way to approach the problem?