Conditional expectation of hierarchical model parameters via marginalisation 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?
 A: I received some guidance on this topic so I'd like to return and add an answer. Bear in mind I don't have a formal math/stats background so my notation is  likely to be a bit off.
I think that my intuition of the nested iteration approach was mostly correct (I thought this was too simple!). For this Dirichlet-Multinomial model, the nested iterated expection is as follows:
$$
\mathbb{E}[X_{ij}|\mu_i, \sigma_i] =  \mathbb{E}\left[\mathbb{E}[X_{ij}|p_{ij}]\right] \\
\text{where} \\
\mathbb{E}[X|p_{ij}] = \mathbb{E}\left[\mathbb{E}[p_{ij}|\alpha_{ij}]\right] \\
\text{where} \\
\mathbb{E}[p_{ij}|\alpha_{ij}] = \mathbb{E}\left[\mathbb{E}[\alpha_{ij}|\mu_i, \sigma_i]\right]
$$
The expectation of these distributions are closed form and widely known (eg Wikipedia).
$$
\mathbb{E}[X_{ij}|p_{ij}] = np_{ij} \ \text{where n is the sum of} \ X_j \\
\mathbb{E}[X_{ij}|\alpha_{ij}] = n\mathbb{E}[\alpha_{ij}] = n\frac{\alpha_{ij}}{\sum_{i=1}^q \alpha_{ij}}
$$
From here I would like to continue and plug in the expectation of the log-normal distribution, but it is not that easy in this case because we of course have dependence between the $\alpha_{ij}$s.
As noted in two recent papers (1 and 2), we can instead approximate this model. Since the marginal distribution of multinomial is binomial, we can approximate the model using a Poisson or negative binomial distribution for each $i$ as long as $q$ is large. When $q$ is small, I may have to return and attempt to solve this more formally.
