## 1 Random Intercept and Slope Model:  
\[logit(p_{i})=\alpha_{i}+\beta X_{i}\]  
\[\begin{bmatrix}\alpha_{i}\\ \beta_{i}\end{bmatrix}\sim MVN\left(\begin{bmatrix}\mu_{\alpha}\\ \mu_{\beta}\end{bmatrix},\Sigma\right)\]  
The Leave-One-Out (LOO) cross-validation method was used to determine the best models. LOO uses log-likelihoods from posterior simulations of the parameter values to estimate pointwise out-of-sample prediction accuracy. This is used to determine the relative predictive performance of each model based on existing posterior simulation draws, which are analysed using Pareto-smoothed importance sampling [29]. This method provides approximate standard errors as a by-product for estimating predictive errors and comparing predictive errors between two models. The lower the LOO information criterion (LOOIC), the better the out-of-sample prediction accuracy.  
Posterior predictive checks were conducted on the best models (as determined by LOO) to assess the fit of the model to the data [30]. These checks included visual comparisons of the mean, empirical cumulative distribution function (ECDF) and density of the observed data and the posterior predictive distributions. There was no systematic discrepancy between any of the models checked and the observed data.  
The posterior distributions of each of the regression coefficients for the raw coefficients and the exponentiated coefficients or ORs were calculated from the best fitting models. These coefficients were used to calculate Bayes Factors, which were calculated to quantify the evidence for the null hypothesis (that the effect of the predictor has no effect on time loss injury) versus the alternative hypothesis. The null interval was specified as \(-0.01\) to \(0.01\).  
Finally, the direction of the effect of predictor on time loss injury was assessed using the probability of direction. This calculates the proportion of the posterior distribution that is greater or less than zero. This provides a measure of the certainty that the effect is positive or negative.  
All models were fitted using Bayesian Regression Models using the Stan (brms) package [31; 32] with MCMC sampling via Stan [33]. Probability of direction and Bayes Factors were calculated using the bayestestR package [34].  
A simulation study was conducted to evaluate the impact of varying sample sizes on the precision of coefficient estimates in a Bayesian mixed-effect logistic regression model to determine if the proposed sampling methodology would result in estimates that would be fit for purpose. As such, the simulation varied the number of participants and the number of observations per participant to determine precision. Precision was assessed based on the width of the credible intervals for coefficient estimates. This analysis was conducted using integrated nested Laplace approximation (INLA) [35]. While the highest precision was observed in scenarios with larger sample sizes, a configuration involving 34 participants and an average of 70 measurements per individual provides credible precision. This provides a balance between feasible data collection efforts and the achievement of statistically reliable estimates. Table 5 outlines selected results from the sampling simulation.