# Hamiltonian Monte Carlo with Random Effect for Analyzing Cyclist Crash Severity
## 2 Method
### Hamiltonian Monte Carlo
#### 2.2.1 Model Formulation  
For estimation of \(\alpha\) in Equation (9), the log likelihood or \(log(f(\theta))\) is needed. The probability of the binary mixed effect model could be written as:  
\[Pr(y=1|X,Z)=\frac{exp(X\beta+Zu)}{1+exp(X\beta+Zu)},\ where\ u\sim N(0,G) \tag{11}\]  
where \(u\) follows multivariate normal distribution with mean 0 and variance covariance of \(G\), \(u\sim N(0,G)\), \(Z\) is a matrix where a column \(j\) in \(z_{j}\) would be filled by values of 1 if the value of random effect exists in \(j\) column, and 0 otherwise. \(u\) could be written as  
\[u=G^{0.5}\tau \tag{12}\]  
where \(\tau\) is the initial values of random effects, to the number of categories.  
As the square root of \(G\) is required for Equation (12), LDL matrix decomposition would be used. Therefore, we have:  
\[G=LDL^{T}=LD^{1/2}D^{1/2}L^{T} \tag{13}\]  
Thus,  
\[\sqrt{G}=LD^{1/2} \tag{14}\]  
where \(L\) is an identity matrix of \(I_{n}\), \(D^{1/2}\) is a matrix with diagonal values of \(\lambda\), where \(\lambda\) follows half-t student distribution [20]. Therefore, matrix of \(\sqrt{G}\) would be built based on Equation (14). Now, to finalize the \(u\) matrix, we need \(G^{1/2}\times\tau\), where \(\tau\) is the parameters to be estimated for random effect, and already the initial values being updated through the leapfrog algorithm.  
As the equation of half-t has a form of the power of \(-(\
u+1)/2\), transformation would be used to make it in a linear form by the help of Jacobian transformation. Due to transformation variables, and from Jacobian transformation, the processes, \(log\ P(\beta|\sigma_{\beta}^{2})\), \(log\ P(\tau)\) and \(log\ P(\xi)\), would be added to the log posterior. The likelihood of \(log\ \pi(\beta|\sigma_{\beta}^{2})\)follows the multinomial probability density function with variance of \(\sigma_{\beta}^{2}\) so it would be written as:  
\[log\ P(\beta|\sigma_{\beta}^{2})\rightarrow-\frac{\beta^{T}\beta}{2\sigma_{\beta }^{2}} \tag{15}\]  
\(\tau\) follows normal distribution probability density function (PDF) with \(N(0,1)\), so it would be written as  
\[og\ P(\tau)\rightarrow-0.5\tau^{T}\tau \tag{16}\]  
After transforming the variance covariance of \(\lambda\) to make it in a form of linear model through Jacobian transformation, the probability of transformed \(\lambda\), \(\zeta=log\ (\lambda)\) and based on PDF of half-t distribution would be written as:  
\[log\ P(\xi)\rightarrow-0.5(\
u_{\xi}+1)log(1+\frac{1}{\
u_{\zeta}}(\frac{e^{2 \
u_{\xi}^{2}}}{A_{\zeta}{}^{2}}))+\xi \tag{17}\]  
where \(\
u_{\xi}\) is parameter of half-t distribution, being the degree of freedom, and \(A_{\xi}\) is the value of variance.  
On the other hand, the log posterior, \(log(f(\theta))\), for logistic mixed effect would be written as:  
\[log\ f(y|X,Z,\beta,u)=(y-1)^{\intercal}(X\beta+Zu)-1\ log\ (1+exp(-(X\beta+Zu)) \tag{18}\]  
where all items were defined before.  
Now, the full likelihood after considering the transformation should include additive terms being discussed above.  
The gradients would be used in the leapfrog algorithm to direct the transitions of Markov Chain to the region of higher probability. The gradients are obtained by getting the derivatives of full log likelihood, which include log priors for transformed variables.