Hi, I work in healthcare and the data is the number of scripts written by each doctor. Naturally, Poisson makes sense. But during the data collection process, the number of scripts is not integer but just non-negative numbers! I wrote a custom likelihood function using optim in R and the results make sense. The key is to use the gamma function instead of the factorial, which does not even affect the coefficients optimization.
I am trying to do the same thing in Stan. But can someone suggest any further optimization? Thanks.
// Xb = X*beta; X is a [NxK] matrix and beta is [Kx1] vector. This term is created beforehand
real cust_PoissonReg_lpdf(real[] Y, vector Xb) {
return sum( to_vector(Y) .* Xb - exp(Xb) );
}
If the data aren’t integer counts, it sounds like the data were somehow summarized. Can you elaborate? Are these perhaps counts per person or per time unit? Is there anyway of converting back to counts and the associated number of person or amount of time? If not, it doesn’t make sense to model the data as Poisson-distributed. Instead use a strictly positive distribution like log-normal or gamma distribution. And if you want, impose additional structure on the variance as function of the mean (as you are now doing in a kind of off way with this faux-Poisson).
These are counts per person. They are just not integers. I am less worried about using Poisson (or Negative binomial) to model non-integers, as I have seen others (for example, Dr. Andrew Gelman in one of his old blog posts) use binomial distribution with non-integers.