I often come across problems with overdispersed poisson-distributed data. For example time series of counts of animals in a given area, e.g. one count per year over 30 years. In these situations it makes sense to me to model these using a gaussian observation level error, a logarithmical link function and a poisson likelihood. The gaussian error would for example represent the effect of weather.
In practice I have the problem that the hyperparameter that controls the size of the observation level errors can vary over several orders of magnitude, which leads to sampling problems. (also with non-centered parameterisations)
I’ve thus recently started trying to approximately integrate out the observation level error, using two different approaches. Basically the idea is to create a custom _lpmf function that includes both the poisson and the gaussian component. But developing them is quite time intensive so I wonder if a solution might already be out there.
Do you know of any existing solutions?
For example I’ve heard of a laplace integrator that has or is beeing developed in stan, but couldn’t find much info about it. Also I wonder whether someone might already have coded a custom lmpf-function that uses an analytical approximation.
Just for completeness:
My first approach is to just use a, not very smart, grid based integrator.
The second approach is to approximate the distribution by a mixture of a negative binomial distribution and a lognormal distribution. This probably would work quite well for high numbers of animals, but for low numbers it fails because the lognormal distribution is continuous while my data is discrete.