Intensity function in Poisson random effect model I have a somewhat general question about intensity functions in Poisson random effect models. 
Consider the Poisson random effects model in which conditional on a random effect $u$, an individual experiences events according to a Poisson process with intensity function $u\rho(t)$. Furthermore, suppose $u$ has a gamma density $g(u)$, with mean 1 and variance $\phi$. Denote $N(t)$ the number of events, $H(t)$ the history, and $\rho(t)= \mu^{\prime}(t)$. 
My goal is to show 
$$ \lambda(t|H(t)) = \left(\frac{1+\phi N(t^{-})}{1+\phi \mu(t)}\right) \rho(t). $$
Here, we have the intensity $\lambda(t|H(t))$ is given by 
$$ \lambda(t|H(t)) = \lim\limits_{\Delta t \to 0}\frac{P(\Delta N(t)=1|H(t))}{\Delta t}. $$
My first thought was the following:
$$
 \lambda(t|H(t)) = \int \lambda(t|H(t),u) g(u)du
= \int u\rho(t) g(u)du = \rho(t).
$$
This is obviously wrong, but I am not sure why. I figured that $\lambda(t|H(t),u)=u\rho(t)$ because conditional on the random effect $u$, we have a Poisson process with intensity $u\rho(t)$, but I could be mistaken.
Edit: Upon further thinking about it, I assume the following statement is incorrect: 
$$
 \lambda(t|H(t)) = \int \lambda(t|H(t),u) g(u)du.
$$
The following should hold instead:
$$
 \lambda(t|H(t)) = \int \lambda(t|H(t),u) g(u|H(t))du.
$$
Any thoughts?
 A: What follows is my solution. 
For small $\Delta t$ with $\Delta N(t)=N(t, t+\Delta t) = N(t+\Delta t)-N(t)$ holds:
$$ P(N(t,t+\Delta t)=1|H(t))
\\=\frac{P(\Delta N(t)=1, H(t))}{P(H(t))}
\\=\frac{\int P(\Delta N(t)=1|H(t),u)P(H(t)|u)g(u)du}{\int  P(H(t)|u)g(u) du}
\\=\frac{\int u\cdot  \rho(t) \cdot \Delta t \cdot P(H(t)|u)g(u) du}{\int  P(H(t)|u)g(u) du} .$$
The history $H(t)$ is defined by $H(t)=\{N(s):0\leq s <t\}$. Conditioned on $u$, the history $H(t)$ follows a Poisson distribution, i.e. 
$$ P(H(t)|u)=P(N(t^{-})|u) =\exp(-\mu(t) u) \frac{(\mu(t) u)^{N(t^{-})}}{N(t^{-})!}. $$
Since $g(\cdot)$ is the density of Gamma distribution, $\int  P(H(t)|u)g(u) du$ is the probability that a negative binomial distribution has the value $N(t^{-})$.
Moreover, is is straight forward to show 
$$\int u P(H(t)|u)g(u) du = \frac{N(t^{-})+1}{\mu(t)}P_{NB}(N(t^{-})+1),$$ 
where $P_{NB}(\cdot)$ is the probability mass function of a negative binomial distribution. 
Finally, it follows that 
$$P(N(t,t+\Delta t)=1|H(t))\\
= \frac{N(t^{-})+1}{\mu(t)} \frac{P_{NB}(N(t^{-})+1)}{P_{NB}(N(t^{-}))}\rho(t)\Delta t\\
= \frac{1+\phi N(t^{-})}{1+\phi \mu(t)} \rho(t)\Delta t$$
and therefore the we obtain the intensity 
$$\lambda(t|H(t))= \frac{1+\phi N(t^{-})}{1+\phi \mu(t)} \rho(t).$$ 
