I thought that since the conjugate prior of Poisson distribution is gamma, we needed to use that when assigning prior distributions to the beta coefficient. MCMCpack and rstanarm both specify a multivariate normal distribution. Can someone explain the reasoning behind that?
1 Answer
If $Y_i \stackrel{ind}{\sim} Poisson(\lambda)$ and $\lambda \sim Gamma(a, b)$, then $$\lambda | Y_1, \ldots, Y_n \sim Gamma\left(\sum_{i = 1}^n y_i + a, n + b\right)$$ as you pointed out.
However, a Poisson regression is set up in a different way:
If $Y_i \stackrel{ind}{\sim} Poisson(\lambda)$, a regression model for $\lambda$ can defined as $$ log(\lambda) = \beta_0 + \beta_1 X_{i,1} + \ldots + \beta_p X_{i,p} $$ where $\mathbf{\beta} = [\beta_0, \beta_1, \ldots, \beta_p]$ is a vector of coefficients. There is no conjugate prior distribution for the $\beta$ in the Poisson regression.
You could define the gamma prior distribution for each $\beta_i$ as you proposed, but then we would state that each $\beta_i$ cannot be negative because the gamma distribution is only defined for positive values.
This could be a problem because you are not allowing that a covariable has a negative relationship with the average number of counts $\lambda$. For example, if $Y_i$ is the number of children in country $i$ and $X_{i,1}$ is the wealth of country $i$, then $\beta_1$ could be negative indicating that people have fewer children as their countries get richer.
The only case that we have conjugate prior for $\beta$ is the normal linear regression, where the normal prior distribution for $\beta$ is conjugate with the normal distribution for the response variable.
I guess the normal prior distribution for $\beta$ in the Poisson regression just follows the normal linear regression. In addition, the hyperparameters of Normal prior distribution are easy to interpret, a less informative prior (i.e., high variance) is easily defined and the normal distribution is defined on the real line, which means that the coefficients $\beta_i$ can assume any real value.