Bayesian Poisson Regression with Gamma Prior Formulas

Question

Are there closed form formulas for the posterior and evidence of a Poisson-Gamma Bayesian regression model? I was not able to find anything that is accessible online. I am not sure for which model can we find such an expression. I was thinking of the following formulation (notation follow):

Likelihood is Poisson: $P(y|X, \theta) = \prod_{i=1}^m \frac{1}{y_i!}(x_i^t \theta)^{y_i} \exp(x_i^t \theta)$ and prior is Gamma: $P(\theta) = \prod_{j=1}^n \frac{\beta^{\alpha}}{\Gamma (\alpha)} \theta_j^{\alpha-1}\exp(-\beta\theta_j)$.

Notation: $m$ number of observations, $n$ number of parameters, $X \in \mathbb{R}^{m\times n}$ a design matrix with rows $x_i \in \mathbb{R}^n$, $\theta \in \mathbb{R}^{n}$ the parameter I would like to infer and $y \in \mathbb{R}^m$ data.

My questions:

1. Is there a closed form expression for the posterior / evidence for the above model?
2. Is there a closed form expression if we formulate the model as a GLM with (e.g.) with a $\log$ link?
3. If there are no closed form formulas - is there a regression model with a conjugate prior, for positive data and parameters?

Freely available references are great, a complete derivation is even better.

Answer 1

It depends on what you would call a "closed form expression". If you expand the powers of $y_i$ into of sum of monomials (in $\theta_j$) then you can write the likelihood as a sum of (potentially very large number of) components, where each component has the form of a product of Gamma distributions. So the posterior is also a mixture distribution of independent Gamma components.

Your prior is technically a conjugate prior in this case because it belongs to the same family (a mixture with one component).

Practically however this will not be very useful unless the numbers are really small - I believe the number of components will be roughly exponential in $\sum y_i$.

On the other hand it is rather straightforward to calculate the gradient of the likelihood/posterior w.r.t $\theta_j$, so you can use it to find the MLE/MAP using gradient-descent optimization methods.

Answer 2

There's a conjugate updating solution for Poisson regression, if your model terms are for categorical variables. I.e. whenever you can summarize the results with sufficient statistics (which in this case are number of events for each combination of categories and, if time periods differ between observations, also the total follow-up time). This is basically straight from the theorems on when conjugate priors exist.

The difficult thing is of course, whether you actually get a conjugate prior that is analytically tractable.

Chen and Ibrahim have some results on how to construct priors for generalized linear models (Chen, M.-H. and Ibrahim, J. G. (2003). Conjugate priors for generalized linear models. Statistica Sinica, 13(2):461–476.). I used those in Section 5 of my PhD thesis to provide the analytical solution for conjugate updating for exponential time-to-event with right-censoring (which is also the solution for Poisson regression, because up to constants the likelihoods are the same) comparing two groups (i.e. there is just a single categorical variable in the regression model) in randomized trial (fully analytically tractable, just a bit hard to get specific priors your might want for the regression coefficient - I proposed one particular solution for that in Section 5.3). This solution is equivalent to assuming independent Gamma priors for the rates in the treatment group and for the control group. I gave formulations in terms of hazard (or rate) ratio + control group (hazard) rate, as well as log-ratio + control rate, and while I didn't give it you should be able to get the log-ratio + log-control-rate parameterization via change of variables (no promises that that's nicely tractable). However, once you do a meta-analysis - i.e. you have two categorical variables, one with 2 levels (=treatment vs. control group) and one with many levels (different studies) - you still have conjugate solution, but you are left with an integral that, as far as I am aware, does not have a nice analytical solution and needs to be evaluated numerically.

I also once tried to do this for logistic regression and I believe that's actually less tractable than the Poisson situation (I'd love to hear otherwise, if someone knows a nice treatment of the topic). As far as I know, you are also in trouble, if you try to add random effects to commonly used GLMs. To be honest, I have not encountered many real-life problems, where conjugate solutions for regression models were really an option. Luckily we have really good MCMC samplers nowadays.