Right, I poked around a bit but I *think* closed-form marginalisation is out of reach. Nevertheless, we can do a few calculations to help us understand the situation a little bit better.

Consider the following hierarchical model:

\begin{align*}
W &\sim \mathcal{P}(\lambda), \\
\lambda &= \exp(\mu + Z), \\
Z &\sim \mathcal{N}(0, \sigma^2),
\end{align*}

with \mu \in \mathbb{R} and \sigma \in \mathbb{R}^+.

We would like to compute E[W] and \text{Var}(W).

We will use the laws of total expectation and total variance

\begin{align*}
E[X] &= E_Y[E_{X|Y}[X|Y]], \\
\text{Var}(X) &= E_Y[\text{Var}_{X|Y}(X|Y)] + \text{Var}_Y(E_{X|Y}[X|Y]).
\end{align*}

So let us first compute E[W|Z] and \text{Var}(W|Z).

From the Poisson distribution we know that if U \sim \mathcal{P}(\nu) then E[U] = \text{Var}(U) = \nu, hence E[W|Z = z] = \text{Var}(W|Z = z) = \exp(\mu + z).

Next we need

\begin{align}
E[W] &= E_Z[\exp(\mu + z)] = \int_{-\infty}^\infty \exp(\mu + z) \frac{\exp(-z^2/2\sigma^2)}{\sqrt{2\pi\sigma^2}}dz, \\
&= \exp(\mu + \sigma^2/2).
\end{align}

Also,

\begin{align}
\text{Var}(W) &= E_Z[\exp(\mu + z)] + \text{Var}_Z(\exp(\mu + z)),\\
&= \exp(\mu + \sigma^2/2) + E_Z[\left(\exp(\mu + z)\right)^2] - \left(E_Z[\exp(\mu + z)]\right)^2, \\
&= \exp(\mu + \sigma^2/2) + \left(\exp(\sigma^2)-1\right)\exp(2\mu + \sigma^2).
\end{align}

## Sanity check

```
set.seed(666)
N <- 1E6
mu <- log(5)
sigma2 <- .5
Z <- rnorm(N, mean = 0, sd = sqrt(sigma2))
W <- rpois(N, lambda = exp(mu + Z))
mean(W); exp(mu + sigma2/2)
[1] 6.416886
[1] 6.420127
var(W); (exp(sigma2)-1) * exp(2*mu + sigma2) + exp(mu + sigma2/2)
[1] 33.14435
[1] 33.15914
```

## Conclusion

No wonder things go awry when the normal random effectâs variance grows, look at what happens to the variance of W when \sigma^2 grows. I donât quite know what do to from here other than maybe fix the linear predictor (\mu in the calculations above) and see what the induced distribution on the first two moments looks like under various priors for the random effect variance, \sigma^2.