## 2 Method
### Model Distributions  
The negative binomial, gamma-Poisson model, is equivalent to the Poisson distributions with mean parameter \(\mu\), where the \(\mu\) parameter follows a gamma distribution with shape parameter \(\alpha\), and rate parameter of \(\alpha/\mu\). Although maximum likelihood has been used extensively in the literature by the generalized linear model (GLM), the Bayesian model especially favors the data, if there is prior expert knowledge about the parameters of interest, or when especially there is a small sample size.  
Although we did not have any prior knowledge about the possible distributions' parameters of interest, the hierarchical technique was used due to a relatively low number of observations, and ease of implementing the hierarchical technique in the Bayesian context. In this paper, a hierarchical approach has been developed where there were 3 levels of hierarchy due to the factor of cluster.  
Often, the response variable is a count response, which makes the algorithm for estimating the parameters more challenging, compared with a binary response like Probit or logit models. For this model we have a set of covariates as \(x=\big{[}1,x_{1},\ldots,x_{p}\big{]}^{T}\), and coefficients as \(\beta=[\beta_{0},\ldots,\beta_{p}]^{T}\). Compared with the Poisson model, where the mean and variance are both equal to \(\lambda\), in this study due to reasons such as differences across the individual observations, crash count, the model would be characterized as overdispersed, having a variance greater than the mean. To address the shortcoming of the Poisson model, the NB model could be used by placing a gamma distribution prior with the shape of \(\alpha\) and scale of \(\frac{p}{1-p}\) so we would have:  
\[f(r,p)=\int_{0}^{\infty}f_{Poisson(\lambda)}f_{Gamma}\bigg{(}\alpha,\frac{1-P} {P}\bigg{)}d\lambda, \tag{1}\]  
where \(\lambda\) is the only parameter of Poisson, while there are two parameters of shape and scale for gamma.  
The probability density function (PDF) of the Poisson-gamma mixture, from Equation (1), would be written as:  
\[P(Y=y_{i}|\lambda_{i},\alpha)=\frac{\Gamma(y_{i}+\alpha)}{\Gamma(y_{i}+1) \Gamma(\alpha)}\bigg{(}\frac{\alpha}{\alpha+\lambda_{i}}\bigg{)}^{\alpha} \bigg{(}\frac{\lambda_{i}}{\alpha+\lambda_{i}}\bigg{)}^{y_{i}}, \tag{2}\]  
where \(\Gamma(.)\) is the gamma function, \(P\) is a probability parameter, and \(Y\) is assumed to follow the negative binomial (NB) distribution through the log link function.  
Now based on the above, the posterior probability of beta conditioned on data could be obtained as:  
\[P(\beta|data)=\frac{P(\beta)P(data|\beta)}{\int P(\beta)d\beta}, \tag{3}\]  
As the equation \(P(data|\beta)\) in the numerator is intractable, after denoting distributions to the parameter's prior, a simulation-based approach would be used for estimating the posterior distributions of \(\beta\).  
The log likelihood function for finding the model parameters of the above equation could be written as follows [3]:  
\[\mathcal{L}=\sum_{i=1}^{n}\{\ln\Big{[}\Gamma\Big{(}y_{i}+\alpha^{-1}\Big{)} \Big{]}-\ln[\Big{(}\Gamma\Big{(}\alpha^{-1}\Big{)}\Big{]}-\ln[\Gamma(y_{i}+1) ]-\alpha^{-1}\ln(1+\alpha\mu_{i})-y_{i}\ln(1+(\alpha\mu_{i}))+y_{i}\ln(\alpha )+y_{i}\ln(\mu_{i}) \tag{4}\]  
The above is the key for the parameters' estimates in the context of maximum likelihood estimation, which would be solved with the help of gradient and Hessian. In the Bayesian context the process is based on sampling, which an upcoming subsection would talk about.