Mean of predictive distribution I observe independent, Poisson-distributed data $ D = \{x_1, ... x_n \} $ with mean parameter $ \mu $, i.e., $$x_i\stackrel{\text{iid}}{\sim}\mathcal{P}(\mu)$$ Over $ \mu $ I assume $ Gamma(\alpha_0, \beta_0) $ as a prior (where $ \alpha_0 $ and $ \beta_0 $ are some given hyperparameters). The posterior distribution is $$\mu|D \sim Gamma\left(a_0 + \sum\limits_{n=1}^N x_n, \beta_0 + N\right) $$ I am interested in the mean and variance of the predictive distribution $p(x_{N+1} | D)$. 
I know it's possible to show that the predictive distribution is a negative binomial, but, since I am only interested in mean and variance there should be a much quicker way by using identities, the law of total expectation and the law of total variance. However, I get confused when trying to apply them here.
The posterior predictive distribution for $ x_{N+1} $ is (leaving out conditioning on hyperparameters): $$ p(x_{N+1}|D) = \int_{\mu} p(x_{N+1}|\mu) \, p(\mu|D) \operatorname{d}\!\mu $$
That is equivalent to the expectation:
$$ p(x_{N+1}|D) = \mathbb{E}_{\mu|D}\Big[p(x_{N+1}|\mu)\Big] $$
The mean of the posterior distribution then is:
$$ \mathbb{E}_{x_{N+1}|D} \Big[\mathbb{E}_{\mu|D}\Big[p(x_{N+1}|\mu)\Big]\Big]$$
which is an iterated expectation to which I would like to apply the law of total expectation.
The law states: $$ \mathbb{E}_X[X] = \mathbb{E}_Y[\mathbb{E}_{X|Y}[X|Y]] $$
Looking at the indices however, I don't know how to proceed from here: That is I cannot map my indices properly to the ones of the law of total expectation.
 A: While writing$$p(x_{N+1}|D)=\mathbb{E}_{\mu|D}[p(x_{N+1}|\mu)]$$is formally correct, this leads to confusion when applying the Law of Total Expectation, in that the formula $$\mathbb{E}_{x_{N+1}|D} \Big[\mathbb{E}_{\mu|D}\Big[p(x_{N+1}|\mu)\Big]\Big]$$ does not make sense.
What you want to find is $\mathbb{E}_{x_{N+1}|D}[X_{N+1}]$ and $\text{var}_{x_{N+1}|D}[X_{N+1}]$. So, if you apply the Law of Total Expectation, you get
$$\mathbb{E}_{x_{N+1}|D}[X_{N+1}]=\mathbb{E}_{\mu|D}\left\{\mathbb{E}_{x_{N+1}|\mu,D}[X_{N+1}]\right\}=\mathbb{E}_{\mu|D}\left\{\mu\right\}=\frac{a_0 + \sum_{n=1}^N x_n}{\beta_0 + N}\,.$$
Similarly, for the variance,
\begin{align*}
\text{var}_{x_{N+1}|D}[X_{N+1}] &= \mathbb{E}_{\mu|D}\left\{\text{var}_{x_{N+1}|\mu,D}[X_{N+1}]\right\}+\text{var}_{\mu|D}[\mathbb{E}_{x_{N+1}|\mu,D}\left\{X_{N+1}\right\}]\\
&= \mathbb{E}_{\mu|D}\left\{\mu\right\}+\text{var}_{\mu|D}[\mu]\\
&= \frac{a_0 + \sum_{n=1}^N x_n}{\beta_0 + N}+\frac{a_0 + \sum_{n=1}^N x_n}{(\beta_0 + N)^2}\\
&= \frac{a_0 + \sum_{n=1}^N x_n}{\beta_0 + N}\,\frac{\beta_0 + N+1}{\beta_0 + N}
\end{align*}
