[R] Poisson Lognormal

[Robert Ruser]

Hi R Users,
I'm going to estimate via. ML the parameters in Poisson Lognormal
model. The model is:

x | lambda ~ Poisson(lambda)
lambda ~ Lognormal(a,b)

Unfortunately, I haven't found a useful package allowing for such
estimation. I tried to use "poilog" package, but there is no equations
and it's hard to understand what exactly this package really does.
Using it I get the incorrect estimators.
I thing that I could use any package that allows to estimate
generalized linear mixed model, because above models is equivalent to:

x | exp(lambda) ~ Poisson (exp(lambda))
exp(lambda) ~ Normal(c,d)

Does it exist any package that can estimate it? May be you know a
package that do Gauss-Hermite quadrature for estimation or simply do
estimation for the first model?

Robert