Bayesian estimation of log-linear Poisson models: analytic and MCMC approaches with applications to reliability and medical count data

Abstract

Introduction:

In reliability engineering and biomedical research, event information is often recorded as counts within fixed time intervals rather than exact event times. Log-linear Poisson models provide a principled framework for modeling such count data through covariate-dependent intensity functions.

Methods:

This study develops a Bayesian formulation of the Poisson regression model using weakly informative priors and compares analytic and simulation-based approaches for posterior estimation, including Laplace approximation and several Markov Chain Monte Carlo (MCMC) methods. The proposed framework is illustrated using a system reliability dataset and a medical event-count dataset. Posterior summaries, convergence diagnostics, prior sensitivity analysis, and posterior predictive checks were employed to assess estimation performance and model adequacy.

Results:

The results demonstrate that Laplace approximation produces estimates comparable to full MCMC methods for moderate sample sizes. Posterior estimates, convergence diagnostics, and model adequacy assessments showed strong agreement across Bayesian computational approaches.

Discussion:

Simulation-based methods offer greater flexibility for more complex settings and provide a useful benchmark for validating analytic approximations. The findings support the application of Bayesian log-linear Poisson models as a reliable framework for analyzing count data in reliability and medical research.

1 Introduction

In reliability engineering and biomedical research, event information is often recorded as counts observed within fixed time intervals rather than as exact event times. Examples include the number of system component failures per month, or the number of clinical events observed during follow-up. Such type of data often arise from reporting practices or monitoring systems that aggregate events over time. When events occur independently within short intervals and the probability of occurrence in a small subinterval is low, the Poisson distribution provides a natural modeling framework for failure and event counts (1, 2, 19, 28).

The primary objective of this work is to develop a Bayesian approach for log-linear Poisson regression models in which the event intensity depends on covariates through a log-link function. Within this framework, regression parameters are assigned weakly informative priors to provide regularization while avoiding overly restrictive assumptions. Weakly informative priors are particularly useful when prior knowledge is limited but complete noninformativeness may lead to instability in estimation (3, 22, 23).

Because posterior distributions for Poisson regression models do not admit closed-form solutions, Bayesian inference requires numerical approximation. This study compares analytic and simulation-based computational approaches for posterior estimation. Specifically, Laplace approximation is used as a deterministic second-order approximation to the posterior density (4), while simulation-based methods include Independent Metropolis sampling, Gibbs sampling (5, 21, 26, 27), and Hamiltonian Monte Carlo implemented through the No-U-Turn Sampler (39). These approaches represent different strategies for exploring posterior distributions and differ in computational efficiency and accuracy.

All methods are implemented within the R environment (6), using LaplacesDemon for Laplace approximation and Independent Metropolis sampling (7), R2jags for Gibbs sampling (8, 29, 33, 34), and RStan for Hamiltonian Monte Carlo (9). Prior sensitivity analysis is conducted to evaluate the influence of prior specification on posterior inference, and posterior predictive checks are used to assess model adequacy and the plausibility of the Poisson modelling assumptions (23, 25).

The methodology is illustrated using two datasets: a system component reliability dataset involving monthly failure counts and a medical event-count dataset related to metastasis occurrence. The goal is not to introduce a new probabilistic model, but to provide a structured comparison of Bayesian computational approaches for log-linear Poisson models and to assess their empirical agreement, stability, and practical suitability across distinct applied domains.

The remainder of the paper is organized as follows. Section 2 presents the Poisson regression model and Bayesian formulation. Section 3 describes the computational methods. Section 4 reports the empirical results for the two applications. Section 5 discusses diagnostic assessment and model validation, and Section 6 concludes with a discussion of findings and limitations.

2 The Poisson model

This section presents the log-linear Poisson regression model, its likelihood formulation, and the associated Bayesian framework. These elements form the basis for the computational approaches examined in this study, including Laplace approximation, Independent Metropolis sampling, Gibbs sampling, and Hamiltonian Monte Carlo via the No-U-Turn Sampler.

2.1 The model and likelihood

Let for denote observed event counts within fixed time intervals. We assume:

With probability mass function.

Here, represents the expected event count for observation . To ensure positivity and incorporate covariate effects, a log-link function is adopted:

Where is the design matrix and is the vector of regression coefficients.

The likelihood for the full dataset is:

When deriving the posterior distribution, the factorial terms do not depend on and can be absorbed into the proportionality constant (Gelman et al., 2013). The likelihood kernel is therefore.

This log-linear formulation accommodates both reliability failure counts and medical event counts within a unified modeling framework.

2.2 Bayesian framework

Bayesian inference proceeds by combining the likelihood with prior distributions for the regression coefficients. By Bayes’ theorem,

Weakly informative Normal priors are assigned to the regression parameters to provide regularization without imposing strong prior structure (3, 7, 23, 30, 32). For centered and scaled predictors, we consider priors of the form.

Where is chosen large enough to avoid undue influence while stabilizing computation (18, 22).

The resulting posterior distribution does not admit a closed-form expression due to the non-conjugacy introduced by the log-link. Consequently, marginal posterior distributions and posterior expectations require numerical approximation.

2.3 Computational methods

Bayesian inference for the log-linear Poisson model reduces to evaluating integrals that are analytically intractable (10). We therefore consider both analytic approximation and simulation-based methods.

2.3.1 Laplace approximation

The Laplace approximation provides a second-order Taylor expansion of the log-posterior density around its mode, yielding a Gaussian approximation centered at the maximum a posteriori estimate (4, 37). This deterministic approach is computationally efficient and often accurate when the posterior is unimodal and approximately symmetric (35, 36). Its performance depends on the curvature of the log-posterior and is most reliable in moderate to large samples (23).

2.3.2 Markov chain Monte Carlo methods

Markov Chain Monte Carlo (MCMC) methods generate dependent samples from the posterior distribution using a Markov chain whose stationary distribution is the target posterior (5, 11, 12, 38). In this study, we implement Independent-Metropolis (IM) sampling and Gibbs sampling to explore the posterior distribution under the same likelihood and prior specifications (30, 32).

Independent Metropolis relies on a fixed proposal distribution, while Gibbs sampling draws sequentially from full conditional distributions when available. These methods are flexible but may exhibit slow mixing in higher-dimensional or highly correlated parameter spaces.

2.3.3 Hamiltonian Monte Carlo and the no-U-turn sampler

Hamiltonian Monte Carlo (HMC) improves sampling efficiency by incorporating gradient information from the log-posterior to suppress random-walk behavior (13); (Gelman et al., 2013). The No-U-Turn Sampler (NUTS) adaptively selects trajectory length and step size, eliminating the need for manual tuning and improving robustness (Hoffman and Gelman, 2014). This method is particularly effective for higher-dimensional regression models.

2.4 Prior sensitivity and posterior predictive assessment

Because prior specification can influence posterior inference, especially in smaller samples, we conduct prior sensitivity analysis by varying prior scale parameters and examining the resulting posterior summaries. Stability across specifications indicates that inference is primarily data-driven and robust to reasonable prior choices.

Model adequacy is assessed using posterior predictive checks, in which replicated datasets are generated from the posterior predictive distribution and compared with observed data (Gelman et al., 2013). Discrepancies between simulated and observed summaries provide evidence of potential model misspecification.

3 Datasets: empirical applications and simulation study

The proposed Bayesian log-linear Poisson framework is illustrated using two applications: a system component reliability dataset and a medical event-count dataset. All analyses were conducted in R using the packages LaplacesDemon, R2jags, and RStan to implement the analytic and simulation-based estimation procedures described in Section 2.

3.1 System component reliability data

The first dataset, referred to as the System Component Reliability (SCR) data, is obtained from SAS Institute Inc. (14). The data record the monthly number of maintenance removals for a complex system over a four-year period (1987–1990). The system consists of numerous components that periodically fail and are replaced or repaired. During the observation period, the system operated under approximately steady conditions.

The response variable is the monthly count of component removals, modeled as a Poisson-distributed outcome. The observed counts are presented in Table 1. The discrete nature of the data distribution is illustrated in Figure 1.

Figure 1

The discrete nature and moderate variability of the counts make this dataset suitable for illustrating Bayesian Poisson regression in a reliability setting.

3.2 Medical event-count data

The second application considers medical count data motivated by lung cancer progression, where the response variable represents the number of metastasis events observed for a patient. Covariates include demographic and clinical characteristics such as age, gender, tumor size, smoking history, cancer stage, and treatment type.

Because access to real individual-level clinical data is restricted, a simulated dataset was generated using a log-linear Poisson model with clinically plausible parameter values. This controlled structure allows evaluation of Bayesian estimation procedures in a multivariable medical setting. A sample of the dataset is shown in Table 2.

The dataset includes both continuous and categorical covariates: patient age and tumor size (continuous), gender and smoking history (binary), cancer stage (I–IV), and treatment type (Chemo, Radiation, Surgery, or None). The response variable is the number of metastasis events, modelled as a Poisson-distributed count.

The presence of both continuous and categorical predictors produces a higher-dimensional regression structure than in the reliability example, allowing assessment of computational efficiency across estimation methods.

4 Bayesian implementation and results

4.1 Bayesian analysis of system’s components reliability data

The monthly number of component removals, denoted by , is modeled within the unified log-linear Poisson regression framework introduced earlier. For month , the outcome satisfies.

With conditional mean.

This specification defines a generalized linear model with canonical log link, where the intercept corresponds to the baseline category (Year 1987, Month 1), the year indicators measure inter-annual deviations, and the month indicators capture seasonal variation. The same likelihood structure (as defined in sec. 2) and prior specification are retained across all computational strategies to ensure methodological consistency.

Independent weakly informative priors are assigned as:

Which allow the data to dominate posterior inference while maintaining numerical stability (3, 16, 17, 23). Because the posterior distribution does not admit a closed-form solution, both analytic and simulation-based approximation methods are applied.

Before applying full posterior computation, preliminary fit was obtained using the bayesglm function from the arm package in R (3). This provides rapid approximate posterior summaries (Table 3) under Student- priors and serves as a diagnostic starting point rather than the primary inferential method used for comparison.

Although informative, these estimates indicate moderate seasonal effects and relatively small year-to-year differences. However, these summaries rely on asymptotic normal approximations and do not represent the full posterior inference targeted in this study.

A model matrix was constructed to implement the Bayesian estimation procedures in LaplacesDemon, JAGS, and Stan. The same likelihood and prior specification were retained across all computational approaches to ensure methodological.

4.1.1 Bayesian estimation using Laplace approximation method

The Laplace approximation provides a deterministic second-order approximation to the posterior density (4, 31).

Where the likelihood arises directly from the log-linear Poisson specification and the prior is Gaussian. A second-order Taylor expansion of around its posterior mode yields:

Where is the negative Hessian matrix of the log-posterior. Optimization was performed using the limited-memory BFGS algorithm within the LaplacesDemon package. This approach produces posterior modes and curvature-based credible intervals derived directly from the log-linear Poisson likelihood.

Table 4 presents the posterior summaries obtained via Laplace approximation. A clear graphical summary of these results illustrated in Figure 2.

Figure 2

The approximation yields stable posterior modes and symmetric credible intervals, consistent with a unimodal posterior structure.

4.1.2 Bayesian estimation using independent Metropolis (IM) algorithm

Full simulation-based inference was conducted using the Independent Metropolis (IM) algorithm implemented in LaplacesDemon. The proposal distribution was centered near the Laplace mode to improve efficiency. The IM algorithm targets the same posterior density defined by the log-linear Poisson likelihood and Normal priors.

The observed acceptance rate of 0.49 falls within the recommended efficiency range for Metropolis-type samplers (38), indicating adequate mixing.

The resulting posterior summaries (Table 5) closely align with the Laplace approximation, validating the analytic results (see Figure 3).

Figure 3

4.1.3 Bayesian estimation using Gibbs sampling (JAGS)

The model was re-expressed in JAGS using the same likelihood structure:

With and normal priors for .

Gibbs sampling iteratively draws from full conditional distributions implied by the joint posterior of the log-linear Poisson model (5, 33, 34). Convergence was assessed using the potential scale reduction factor (15, 20), with all parameters satisfying . Convergence diagnostics indicate stable mixing and posterior estimates consistent with both Laplace and Independent Metropolis results. Table 6 represents the full posterior summary via Gibbs sampling using JAGS language (see Figure 4).

Figure 4

4.1.4 Bayesian estimation using Hamiltonian Monte Carlo (NUTS) sampler (Stan)

Hamiltonian Monte Carlo leverages gradient information from the log-posterior to reduce random-walk behavior (13); (Gelman et al., 2013). The No-U-Turn Sampler (NUTS) adaptively tunes trajectory length and step size (Hoffman and Gelman, 2014), improving efficiency in higher dimensions. The gradient of the log-posterior,

Is computed directly from the log-linear Poisson likelihood and Normal priors, enabling efficient exploration of the posterior surface. Four parallel chains were run in Stan, and all parameters achieved .

Posterior summaries (Table 7) are nearly identical to those obtained from the other MCMC methods, confirming consistency across computational strategies (see Figure 5).

Figure 5

4.1.5 Comparative assessment

Across all four methods—Laplace approximation, Independent Metropolis, Gibbs sampling, and Hamiltonian Monte Carlo—the posterior summaries are highly consistent. The Laplace method provides a fast deterministic approximation (4), while IM and Gibbs sampling confirm its accuracy through stochastic simulation (Robert and Casella, 2010). Hamiltonian Monte Carlo offers improved efficiency and scalability, particularly relevant for more complex multivariable extensions.

The strong agreement across analytic and simulation-based methods directly supports the study objective by demonstrating empirical equivalence across Bayesian computational strategies for well-behaved Poisson posteriors.

4.2 Bayesian analysis of simulated lung cancer data

The lung cancer dataset was generated under a log-linear Poisson regression structure to evaluate Bayesian estimation performance in a multivariable medical setting.

Let denote the number of metastasis events observed for patient . The outcome is modeled as:

With log-link specification.

Baseline categories are Female, Stage I, and No Treatment. Under this parameterization, each coefficient represents a multiplicative change in the expected metastasis rate through the relationship , holding other covariates constant. Continuous predictors such as age and tumor size act linearly on the log scale, while categorical predictors induce relative rate shifts.

The simulated dataset was generated with a prespecified sample size of and with covariates drawn from distributions chosen to reflect a clinically plausible lung cancer setting. Specifically, age was generated from a normal distribution with mean 65 and standard deviation 10, tumor size was generated from a uniform distribution on [2, 8], whereas gender and smoking status were generated from Bernoulli distributions, stage was generated from a four-category multinomial distribution, and treatment was generated from a categorical distribution with four levels. The true regression coefficients used in the simulation were fixed at β = (−1.5, 0.02, 0.10, 0.15, 0.20, 0.30, 0.50, 0.70, 0.20, 0.40, 0.10), which were selected to produce interpretable metastasis rates and moderate separation among the predictor effects. The data were generated under the standard Poisson assumption without additional overdispersion, so that the simulation reflects the intended Poisson regression structure rather than a mis specified count model.

To assess the stability of the Bayesian estimators, the single simulated dataset was analyzed using the four estimation strategies considered in this study. Because the purpose of this section is to provide an objective benchmark with known parameter values, the simulated data allow direct assessment of estimator accuracy, uncertainty estimation, and computational efficiency.

The same likelihood,

And prior specifications.

Are maintained to ensure comparability with the reliability analysis. The simulation design was intended to provide an objective benchmark because the true parameter values are known in advance. This allows comparison of posterior estimates, uncertainty measures, and computational efficiency across methods under controlled conditions.

As the posterior distribution does not admit a closed-form solution, the same four estimation strategies were also implemented for the lung cancer data set: Laplace approximation, Independent Metropolis sampling, Gibbs sampling, and Hamiltonian Monte Carlo (NUTS). All four methods target the same posterior distribution but differ in approximation mechanism, mixing behavior, and computational efficiency.

For each simulated dataset, we recorded posterior means, posterior standard deviations, credible intervals, convergence diagnostics, and wall-clock computation time.

The simulation results were summarized using bias, root mean squared error, empirical coverage probability, average interval width, maximum absolute difference in posterior means, relative error in standard deviations, and effective sample size for the MCMC methods.

These summaries provide a reproducible basis for comparing the accuracy and efficiency of the four Bayesian estimation approaches under a controlled Poisson regression setting with known parameter values.

4.2.1 Bayesian estimation using Laplace approximation method

Posterior modes indicate that tumor size, advanced cancer stage, and chemotherapy are strong positive determinants of metastasis counts. Age shows a modest positive association, while gender and smoking effects are comparatively smaller. The intercept is negative, reflecting a low baseline metastasis intensity under reference conditions. Table 8 reports full posterior summaries.

4.2.2 Bayesian estimation using independent-Metropolis algorithm

Simulation-based posterior mean closely match the Laplace modes. Uncertainty estimates are moderately larger, particularly for stage and treatment effects, reflecting full posterior exploration rather than local curvature approximation (see Table 9, Figure 6).

Figure 6

4.2.3 Bayesian estimation using Gibbs sampler via JAGS

The Gibbs-based posterior summaries remain highly consistent with the Independent Metropolis results. Convergence diagnostics indicate stable chains, and credible intervals confirm the dominant effects of tumor size, stage progression, and chemotherapy (see Table 10, Figure 7).

Figure 7

4.2.4 Bayesian estimation using (NUTS) sampler using RStan

The NUTS results show high effective sample sizes and stable posterior summaries. Parameter estimates are virtually identical to the other MCMC approaches, reinforcing the robustness of the simulated regression structure (see Table 11, Figure 8).

Figure 8

4.2.5 Comparative assessment

All four methods converge to nearly identical posterior means, confirming posterior stability under the log-linear Poisson structure. Laplace approximation slightly underestimates uncertainty relative to simulation-based methods, as expected from its quadratic assumption. Independent Metropolis and Gibbs sampling yield comparable inference but require larger iteration counts to achieve similar precision. NUTS demonstrates the highest computational efficiency and most stable effective sample sizes.

The small maximum absolute differences between the Laplace estimates and the MCMC posterior means indicate that the fitted posterior surface is regular and approximately Gaussian in this simulated setting. The relative standard deviation errors are also small, suggesting that the different methods provide closely aligned uncertainty estimates, although Laplace approximation remains slightly more local in nature.

Among the simulation-based methods, NUTS provides the best balance of sampling efficiency and numerical stability, as reflected by its higher effective sample size per second and more efficient exploration of the posterior distribution. Independent Metropolis and Gibbs sampling produce broadly similar posterior summaries, but their lower effective sample size per second values indicate slower posterior exploration and greater computational burden.

To complement the graphical comparison, Table 12 reports quantitative agreement across methods using maximum absolute differences from the Laplace estimates, relative standard deviation error, wall-clock run time, and effective sample size per second for the MCMC methods. These summaries provide a concise measure of both estimation accuracy and computational efficiency.

For this simulated dataset, methodological agreement is strong because the posterior surface is regular and approximately Gaussian. In more complex or sparse clinical count data, divergence between Laplace and MCMC methods may signal non-normal posterior geometry and motivate full simulation-based inference.

The table shows that the posterior means are highly consistent across all four methods, with only small method-specific differences. The MCMC samplers provide comparable point estimates, while their relative efficiency is better summarized by ESS/s and run time than by posterior means alone (see Figure 9).

Figure 9

5 Model adequacy: prior sensitivity analysis and posterior predictive assessment

Bayesian inference depends jointly on the likelihood specification and prior assumptions. Even when computational convergence is achieved, inference may be unreliable if the assumed log-linear Poisson structure fails to represent key data features. Model adequacy is therefore evaluated through prior sensitivity analysis and posterior predictive assessment, following standard Bayesian practice (20, 23).

5.1 Prior sensitivity analysis

Throughout Sections 4.1 and 4.2, regression coefficients were assigned independent weakly informative Normal priors to stabilize computation while allowing the likelihood to dominate inference (3); (Gelman et al., 2013). Because the two applications differ in structure and complexity, different prior scales were considered for sensitivity analysis.

5.1.1 System component reliability data (SCRD)

For the system component reliability data, prior sensitivity was examined using three increasingly concentrated Normal priors,

Representing weakly to moderately informative specifications.

Posterior means and credible intervals under these alternative priors were nearly identical to those reported in Section 4.1. Year effects remained centered near zero, while seasonal month effects—particularly February, October, and December—consistently retained positive posterior mass away from zero. Posterior uncertainty exhibited only minor changes across prior scales, indicating that inference is driven primarily by the likelihood rather than prior specification (see Figures 10, 11).

Figure 10

Figure 11

5.1.2 Simulated lung cancer data

For the simulated lung cancer dataset, broader prior scales were adopted to reflect higher parameter uncertainty in a multivariable clinical setting. Prior sensitivity was assessed using.

Posterior summaries for tumor size, cancer stage indicators, and chemotherapy effects remained stable across prior specifications. Minor widening of credible intervals was observed for weaker predictors such as smoking status and radiation therapy; however, effect direction and substantive interpretation were unchanged. The intercept consistently reflected a low baseline metastasis intensity under reference conditions (see Figure 12).

Figure 12

5.2 Posterior predictive assessment

Posterior predictive checks were used to assess whether data generated from the fitted models resemble the observed counts (Gelman et al., 2013). Under the log-linear Poisson specification,

Replicated observations were drawn from the posterior predictive distribution.

5.2.1 System component reliability data (SCRD)

For the system component reliability data, replicated monthly counts closely matched the empirical distribution. Posterior predictive means aligned with observed averages, and dispersion characteristics were comparable. No systematic underprediction of high-count months or overprediction of low-count months was detected. Discrepancy measures based on total counts and maximum monthly removals yielded Bayesian p-values near 0.5, providing no evidence of lack of fit.

Figure 13 illustrates the posterior predictive density comparisons between observed and replicated monthly removal counts for the SCRD across the Bayesian estimation methods.

Figure 13

Simulated lung cancer data

For the simulated lung cancer data, posterior predictive distributions successfully reproduced the observed variability and right-skewed count structure induced by heterogeneous covariate effects. Because the data were generated under the assumed model, close predictive agreement is expected; nevertheless, these results confirm correct implementation of the Bayesian computational procedures.

Figure 14 presents the posterior predictive checks for the lung cancer model. The top panels compare observed and replicated count distributions, while the bottom panels compare posterior predictive mean and variance across Bayesian estimation methods.

Figure 14

Convergence diagnostics further support adequacy. Across all MCMC implementations, potential scale reduction factors satisfied (15), effective sample sizes were sufficient, and no pathological sampling behavior was detected under NUTS sampler. The close agreement between Laplace approximation and simulation-based methods additionally indicates a regular and approximately unimodal posterior surface.

Overall, both prior sensitivity analysis and posterior predictive assessment support the adequacy of the log-linear Poisson framework for the applications considered.

6 Discussion and conclusion

This study examined Bayesian estimation of log-linear Poisson regression models with applications to reliability engineering and medical event-count data. The primary objective was to provide a structured comparison of analytic and simulation-based Bayesian approximation techniques under a unified likelihood and prior specification.

For both applications, the log-linear Poisson model.

Offered a coherent representation of count data, allowing covariate effects to act multiplicatively on the intensity scale. The use of weakly informative Normal priors ensured computational stability while permitting the likelihood to dominate posterior inference.

Across both data sets, strong agreement was observed among the Laplace approximation, Independent Metropolis sampling, Gibbs sampling, and Hamiltonian Monte Carlo. The small maximum absolute differences in posterior means and the limited relative error in posterior standard deviations indicate that all four methods target essentially the same posterior structure, with only minor numerical variation. When the posterior distribution is approximately Gaussian and unimodal, the Laplace approximation provides accurate posterior summaries with minimal computational cost (4, 31, 37). Simulation-based methods validate these results while offering full posterior uncertainty (5, 38). The efficiency metrics further clarify the practical differences among the methods. In particular, effective sample size per second and total wall-clock time show that Hamiltonian Monte Carlo via the No-U-Turn Sampler delivers the most efficient posterior exploration among the MCMC methods, combining high sampling efficiency with stable convergence behavior (13); (Hoffman and Gelman, 2014). Gibbs sampling and Independent Metropolis produced comparable posterior estimates, but their lower effective sample sizes per unit time indicate slower exploration of the posterior surface and greater computational burden. By contrast, Laplace approximation was the fastest deterministic method and produced results close to the MCMC summaries, but its local quadratic approximation may slightly understate posterior uncertainty relative to full simulation-based inference.

Prior sensitivity analysis confirmed that posterior inference remained stable under reasonable variations in prior scale for both applications. Posterior predictive checks further indicated that the fitted models adequately reproduced key features of the observed count distributions in reliability and medical contexts, supporting the adequacy of the assumed log-linear Poisson structure.

A key limitation of the present framework is the Poisson assumption of equi-dispersion. In settings with overdispersion, zero inflation, or unobserved heterogeneity, extensions such as negative binomial or hierarchical Poisson models may be more appropriate (2). Moreover, in sparse or highly correlated designs, analytic approximations may understate posterior uncertainty, reinforcing the importance of full simulation-based inference.

In conclusion, Bayesian log-linear Poisson regression offers a flexible and computationally efficient framework for modeling reliability and medical count data. The close agreement between analytic and MCMC-based approaches observed in this study supports the use of Laplace approximation for rapid inference in moderate settings, while Hamiltonian Monte Carlo provides a robust and scalable solution for more complex models. These findings establish a practical methodological foundation for future extensions to hierarchical reliability systems and multilevel clinical count data analyses.

References

• 1.

  GrantGMRoesenerWSHallDGAtwoodCLGentillonCDWolfTR. Reliability Study: High-Pressure Coolant Injection (HPCI) System, 1987–1993. Idaho Falls: Idaho National Engineering and Environmental Laboratory (1999).

  • Google Scholar
  • View reference in article

• 2.

  HamadaMSWilsonAGReeseCSMartzHF. Bayesian Reliability. New York, NY: Springer (2008).

  • Google Scholar
  • View reference in article

• 3.

  GelmanAJakulinAPittauMGSuY. A weekly informative default prior distribution for logistic and other regression models. Ann Appl Stat. (2008) 2:1360–83. doi: 10.1214/08-AOAS191

  • CrossRef
  • Google Scholar
  • View reference in article

• 4.

  TierneyLKadaneJB. Accurate approximations for posterior moments and marginal densities. J Am Stat Assoc. (1986) 81:82–6. doi: 10.1080/01621459.1986.10478240

  • CrossRef
  • Google Scholar
  • View reference in article

• 5.

  GemanS.GemanD. (1984). Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. In: IEEE Transactions on Pattern Analysis and Machine Intelligence

  • Google Scholar
  • View reference in article

• 6.

  R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna: R Core Team (2020).

  • Google Scholar
  • View reference in article

• 7.

  Statisticat, LLC. LaplacesDemon: Complete Environment for Bayesian Inference. Bayesian-Inference.Com. R Package Version 16.1.4. London: Statisticat, LLC (2020).

  • Google Scholar
  • View reference in article

• 8.

  SuY.YajimaM. (2014) R2jags: a package for running JAGS from R. R package version 0.04-01, Available online at: http://CRAN.R-project.org/package=R2jags (accessed on February 19, 2026).

  • Google Scholar
  • View reference in article

• 9.

  Stan Development Team (2025) RStan: the R interface to Stan (R package version 2.32.7). Available online at: https://mc-stan.org/ (accessed on February 19, 2026).

  • Google Scholar
  • View reference in article

• 10.

  BernadoJSmithA. Bayesian Theory. West Sussex: John Wiley & Sons (2000).

  • Google Scholar
  • View reference in article

• 11.

  HastingsW. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. (1970) 57:97–109. doi: 10.1093/biomet/57.1.97

  • CrossRef
  • Google Scholar
  • View reference in article

• 12.

  MetropolisNRosenbluthARosenbluthMTellerATellerA. Equation of state calculations by fast computing machine. J Chem Phys. (1953) 21:1087–92. doi: 10.1063/1.1699114

  • CrossRef
  • Google Scholar
  • View reference in article

• 13.

  GordanNJSalmondDJSmithAFM. Novel approach to non-linear/non-Gaussian Bayesian state estimation. IEEE Proc F Radar Sig Proc. (1993) 140:107–17. doi: 10.1049/ip-f-2.1993.0015

  • CrossRef
  • Google Scholar
  • View reference in article

• 14.

  SAS Institute Inc.SAS/STATr 9.2 User’s Guide. Cary, NC: SAS Institute Inc (2008).

  • Google Scholar
  • View reference in article

• 15.

  GelmanARubinDB. Inference from iterative simulation using multiple sequences. Stat Sci. (1992) 7:457–511. doi: 10.1214/ss/1177011136

  • CrossRef
  • Google Scholar
  • View reference in article

• 16.

  AkhtarMTKhanAA. Bayesian Analysis of Generalized Log-Burr Family with R. Berlin: SpingerPlus (2014a).

  • Google Scholar
  • View reference in article

• 17.

  AkhtarMTKhanAA. Log-logistic distribution as a reliability model: a Bayesian analysis. Am J Math Stat. (2014b) 4:162–70. doi: 10.5923/j.ajms.20140403.05

  • CrossRef
  • Google Scholar
  • View reference in article

• 18.

  AkhtarMTKhanAA. JMASM 51: Bayesian reliability analysis of binomial model—application to success/failure data. J Mod Appl Stat Methods. (2018a) 17:eP2623. doi: 10.22237/jmasm/1553803862

  • CrossRef
  • Google Scholar
  • View reference in article

• 19.

  AkhtarMTKhanAA. Bayesian analysis of Poisson reliability model with R and JAGS. Int J Recent Sci Res. (2017) 8:21837–41. doi: 10.24327/ijrsr.2017.0811.1151

  • CrossRef
  • Google Scholar
  • View reference in article

• 20.

  BrooksSPGelmanA. General method for monitoring convergence of iterative simulations. J Comput Graph Stat. (1998) 7:434–55. doi: 10.1080/10618600.1998.10474787

  • CrossRef
  • Google Scholar
  • View reference in article

• 21.

  CasellaGGeorgeE. Explaining the Gibbs sampler. Am Stat. (1992) 46:167–74. doi: 10.1080/00031305.1992.10475878

  • CrossRef
  • Google Scholar
  • View reference in article

• 22.

  GelmanA. Scaling regression inputs by dividing by two standard deviations. Stat Med. (2008) 27:2865–73. doi: 10.1002/sim.3107

  • CrossRef
  • Google Scholar
  • View reference in article

• 23.

  GelmanACarlinJBSternHSDunsonDBVehtariARubinDB. Bayesian Data Analysis (3rd ed.). Boca Raton, FL: Chapman & Hall/CRC Press (2013).

  • Google Scholar
  • View reference in article

• 24.

  GelmanAHillJ. Data Analysis Using Regression and Multilevel/Hierarchical Models. New York, NY: Cambridge University Press (2007).

  • Google Scholar
  • View reference in article

• 25.

  GelmanASuY-S. Arm: Data Analysis Using Regression and Multilevel/Hierarchical Models. Vienna: R Core Team (2020).

  • Google Scholar
  • View reference in article

• 26.

  GelfandAHillsSRacine-PoonASmithA. Illustration of Bayesian inference in Normal data models using Gibbs sampling. J Am Stat Assoc. (1990) 85:972–85. doi: 10.1080/01621459.1990.10474968

  • CrossRef
  • Google Scholar
  • View reference in article

• 27.

  GelfandASmithA. Sampling-based approaches to calculating marginal densities. J Am Stat Assoc. (1990) 85:398–409. doi: 10.1080/01621459.1990.10476213

  • CrossRef
  • Google Scholar
  • View reference in article

• 28.

  JohnsonVEMoosmanACotterP. A hierarchical model for estimation the early reliability of complex systems. IEEE Trans Reliab. (2005) 54:224–31. doi: 10.1109/TR.2005.847262

  • CrossRef
  • Google Scholar
  • View reference in article

• 29.

  KeryM. Introduction to WinBUGS for Ecologist: A Bayesian Approach to Regression, ANOVA, Mixed Models, and Related Analyses. Cambridge, MA: Academic Press (2010).

  • Google Scholar
  • View reference in article

• 30.

  AkhtarMT. Robustness of Bayesian Methods in Healthcare System Assessment: A Comprehensive Review. Int. J. Stat. Med. Res. (2025) 14:662–675. doi: 10.6000/1929-6029.2025.14.62

  • CrossRef
  • Google Scholar
  • View reference in article

• 31.

  KhanNAkhtarMTKhanAA. A Bayesian approach to survival analysis of inverse Gaussian model with Laplace approximation. Int J Stat Appl. (2017) 6:391–8. doi: 10.5923/j.statistics.20160606.08

  • CrossRef
  • Google Scholar
  • View reference in article

• 32.

  AkhtarMTKhanNKhanAA. Bayesian Estimation and Model Assessment of the Exponentiated Rayleigh Survival Model Using Laplace and MCMC Techniques: Applications with Right-Censored Medical Data. Statistics, Optimization & Information Computing. (2026) 15:3485–3511. doi: 10.19139/soic-2310-5070-2986

  • CrossRef
  • Google Scholar
  • View reference in article

• 33.

  NtzoufrasI. Bayesian Modeling Using WinBUGS. New York, NY: John Wiley & Sons (2009).

  • Google Scholar
  • View reference in article

• 34.

  PlummerM. (2003) JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. In: Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003), March 2022

  • Google Scholar
  • View reference in article

• 35.

  TannerMA. Tools for Statistical Inference. New York: Springer-Verlag (1996).

  • Google Scholar
  • View reference in article

• 36.

  TannerMAWongW. The calculation of the posterior distributions by data augmentation. J Am Stat Assoc. (1987) 82:82–6. doi: 10.2307/2289457

  • CrossRef
  • Google Scholar
  • View reference in article

• 37.

  TierneyLKassRKadaneJ. Fully exponential Laplace approximations to expectations and variances of nonpositive functions. J Am Stat Assoc. (1989) 84:710–6. doi: 10.1080/01621459.1989.10478824

  • CrossRef
  • Google Scholar
  • View reference in article

• 38.

  RobertCPCasellaG. Introducing Monte Carlo Methods with R. New York, NY: Springer (2010).

  • Google Scholar
  • View reference in article

• 39.

  HoffmanMDGelmanA. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. J Mach Learn Res. (2014) 15:1593–1623.

  • Google Scholar
  • View reference in article