Abstract

This study compares five probabilistic distributions whose HLPKD, Weibull, Gamma, Gompertz and Lindley for modeling annual neonatal mortality rates in Burkina Faso over the period 1969-2023. The objective is to identify the most appropriate model for describing the distribution of this indicator and to evaluate estimator stability through simulation. Parameters were estimated using the maximum likelihood method based on data compiled by the United Nations Inter-agency Group for Child Mortality Estimation (UN IGME), accessible via the CEIC database. Monte Carlo simulations were conducted using R software (version 4.3.2) with $N=500$ replications for sample sizes $n=40,80,120,160,200$ . Performance criteria included bias, relative bias, mean squared error (MSE), root mean squared error (RMSE), average confidence interval lengths (AL₉₀ and AL₉₅) and coverage probabilities (CP₉₀ and CP₉₅). For the real data, goodness-of-fit was assessed using log-likelihood (lnL), Akaike information criterion (AIC), Bayesian information criterion (BIC), the Kolmogorov-Smirnov test (KS), mean absolute error (ASAE), the Cramér-von Mises test ($W$ ) and the Anderson-Darling test ($A$ ). Bootstrap-corrected p-values were also computed to account for parameter estimation uncertainty. The simulations showed that the Weibull, Gamma, Gompertz and Lindley models converge rapidly to stable estimates, whereas HLPKD requires larger sample sizes to achieve comparable precision. The application to real data identified the Gamma and HLPKD models as the best-fitting, with high p-values (0.8431 for Gamma and 0.9050 for HLPKD) and small discrepancies between theoretical distributions and observations. The Gamma model stands out for its parsimony (two parameters) and estimation stability, while HLPKD, despite its flexibility, suffers from greater uncertainty in its parameter $\theta$ . Bootstrap-corrected tests confirmed that Gamma and HLPKD are the only models not rejected at the 5% level. The Gamma model emerges as the most reliable for describing the distribution of neonatal mortality rates in Burkina Faso, confirming that the increased flexibility of HLPKD does not translate into a decisive advantage for aggregated data.

Keywords

Neonatal Mortality, Burkina Faso, HLPKD, Weibull, Gamma, Lindley, Gompertz, Maximum Likelihood, Monte Carlo Simulation

Share and Cite:

1. Introduction

Neonatal mortality, defined as the number of deaths occurring within the first 28 days of life per 1000 live births, is a fundamental health indicator. It reflects not only the quality of maternal and child health care but also the level of development of a health system. While advanced countries have managed to stabilize this rate at a low level, typically between 2 and 5 per 1000 live births, the situation remains concerning in many sub-Saharan African countries. Recent studies conducted in West Africa, particularly in Gambia and Burkina Faso, report neonatal mortality rates reaching 11.5 per 1000 live births, with an alarming concentration of deaths (82%) occurring within the first 72 hours after birth [1]. This critical period, during which the risk is highest, gives the distribution of survival times a particularly asymmetric shape.

Statistical modeling of survival durations or mortality rates relies on the choice of an appropriate probabilistic distribution. Historically, the law proposed by [2] laid the foundations for the analysis of mortality phenomena by introducing the concept of a mortality force increasing exponentially with age [3]. This distribution, initially developed for human life tables, paved the way for a rich tradition of parametric modeling in demography and epidemiology. Building on this work, more recent distributions have been proposed to address the growing complexity of survival data. The Lindley distribution, for example, has been the subject of recent methodological developments, including approaches using hybrid censoring for the analysis of medical and biological data, thereby demonstrating its relevance in contexts where information is partially observed [4].

Concurrently, sustained efforts are being made to increase the flexibility of existing models. The introduction of additional parameters aims to better capture complex behaviors such as heavy tails or pronounced skewness. The Gamma distribution, widely used in survival analysis, has recently been revisited within the framework of neutrosophic logic to address the imprecision and uncertainty frequently encountered in health statistics, with a specific application to neonatal mortality data [5]. These studies show that the development of new models capable of handling data ambiguity constitutes an active and promising research direction.

In this vein, a major recent contribution is the introduction of the Half-Logistic Power Komal Distribution (HLPKD) by [6]. By combining the half-logistic generator with the Power Komal distribution, these authors proposed a three-parameter model offering increased flexibility for the analysis of biomedical and radiation data. Their study demonstrated, on individual survival data from cancer patients, that HLPKD could outperform numerous competing distributions, including classical models such as Weibull or Gamma, as well as more recent models such as the power Lindley or half-logistic Weibull Pareto distributions. This work notably highlighted the ability of HLPKD to adapt to various shapes of hazard functions (increasing, decreasing, or J-shaped). A key methodological feature of their work lies in conducting an extensive simulation study, with varying sample sizes ($n=50,100,150,200,300$ ) and different parameter combinations, enabling the assessment of the stability and reliability of maximum likelihood estimators.

This simulation approach, which consists of generating artificial samples from known parameter values to evaluate estimator behavior, is a methodological standard in validating new statistical models. It allows, in particular, examining estimator convergence, bias reduction, and precision improvement as sample size increases. Incorporating such an approach into the present study aims to evaluate not only the goodness-of-fit of models to observed data but also their stability under varying sampling conditions.

These recent developments raise a fundamental question: can a model’s superiority be transferred from one type of data to another? The Weibull distribution, for instance, is frequently used in epidemiological studies in Africa to identify risk factors for neonatal mortality, whether for newborns admitted with perinatal asphyxia or for general populations [7] [8]. Its shape parameter enables characterizing the evolution of risk over time, thereby providing valuable insights for clinical intervention. This widespread use of Weibull contrasts with the scarcity of applications of more complex models such as HLPKD to aggregated data or to specific geographic contexts like West Africa.

The present article addresses this issue. It aims to evaluate, for the first time to our knowledge, the performance of the HLPKD distribution on neonatal mortality data from Burkina Faso, comparing it with a set of classical and widely validated distributions. The choice of this geographic context is justified by the persistence of high neonatal mortality rates and by the availability of time series covering a sufficiently long period to allow robust modeling. The analysis focuses on the marginal distribution of annual neonatal mortality rates, a form of aggregated data that differs from the individual survival data for which HLPKD was initially validated [6]. It is explicitly acknowledged that the annual observations exhibit a strong downward trend over time and are not independent; the study does not model the temporal dynamics or serial dependence. Rather, the objective is to identify which probability distribution best describes the overall shape and variability of the observed rates when treated as a sample from an underlying fixed distribution.

More specifically, this study pursues the following objectives:

• Fit the HLPKD distribution to neonatal mortality data from Burkina Faso over the period 1969-2023 using the maximum likelihood method;

• Estimate the parameters of competing distributions, namely Weibull, Gamma, Lindley and Gompertz, on the same dataset;

• Conduct a simulation study with different sample sizes and parameter combinations to evaluate the stability of estimators for the five considered models;

• Compare the goodness-of-fit of the five models using information criteria (AIC, BIC), goodness-of-fit tests (Kolmogorov-Smirnov) and graphical indicators (fitted density functions, P-P plots, etc.);

• Identify, based on these comparisons, the most appropriate statistical model for describing the behavior of the neonatal mortality rate in this specific context.

The remainder of this article is organized as follows. Section 2 presents the methodology employed, including data description, specification of the five candidate distributions, estimation techniques, simulation design and comparison criteria. Section 3 reports the results obtained, both for the simulation study and for the fitting to real data. Section 4 discusses these results in light of existing literature and the Burkinabe health context. The conclusion summarizes the main contributions of the study and suggests perspectives for research and action.

2. Methodology

The methodological approach adopted in this study aims to compare the performance of five probabilistic distributions for modeling annual neonatal mortality rates in Burkina Faso. This section successively presents the statistical models selected, the data used and how they were collected, the simulation procedures implemented to evaluate estimator stability, as well as the estimation methods and comparison criteria employed.

2.1. Model Specification

This subsection presents in detail the five probabilistic distributions selected to model annual neonatal mortality rates in Burkina Faso. The choice of these models is based on several criteria: their classical use in survival analysis (Weibull, Gamma, Gompertz), their parsimony and flexibility (Lindley), and their ability to capture complex behaviors such as pronounced skewness and heavy tails (HLPKD).

To facilitate reading and comparison across the different models, a standardized parameter notation has been adopted. Throughout this study, the symbol $\alpha$ consistently denotes a shape parameter, $\beta$ a scale (or rate) parameter, $\theta$ an additional flexibility parameter (for the HLPKD model), and $\eta$ the risk growth parameter for the Gompertz distribution. For the Lindley distribution, $\theta$ represents its single parameter. This convention ensures coherent presentation of estimation results and avoids confusion regarding the roles of parameters across models.

For each distribution, five fundamental functions are presented: the probability density function (PDF), the cumulative distribution function (CDF), the survival function, the hazard function, and the quantile function. These functions constitute essential tools for parameter estimation, simulated data generation, and goodness-of-fit assessment. The mathematical properties of these distributions are well documented in standard references such as [10]. The following sections successively detail the Weibull, Lindley, Gamma, Gompertz and HLPKD distributions, highlighting their mathematical properties and relevance for mortality data analysis.

2.1.1. Density Functions

Probability density functions (PDF) describe the distribution of annual neonatal mortality rates for each considered model. Table 1 presents these functions for the five selected distributions, adhering to the standardized parameter notation established above.

Table 1. Density functions (PDF).

Source: [2] [6] [9]-[11].

The parameters play a fundamental role across the different models:

• $\alpha$ : shape parameter, controls the shape of the distribution and the hazard function.

• $\beta$ : scale (or rate) parameter, adjusts dispersion or the rate of decay.

• $\theta$ : additional flexibility parameter, specific to HLPKD. It also serves as the single parameter of the Lindley distribution.

• $\eta$ : risk growth parameter, specific to the Gompertz distribution.

• $\Gamma \left(\alpha \right)$ : Gamma function, defined by $\Gamma \left(\alpha \right)={\displaystyle {\int }_0^{\infty }{y}^{\alpha -1}{\text{e}}^{-y}\text{d}y}$ , used to normalize the density of the Gamma distribution.

• For the HLPKD distribution, the complexity of the PDF reflects the increased flexibility of the model, allowing it to capture different shapes of asymmetry and tail behavior [6].

2.1.2. Cumulative Distribution Functions

The cumulative distribution function (CDF) allows assessing the probability that the neonatal mortality rate is less than or equal to a given value. Table 2 presents these functions for the five considered models.

Table 2. Cumulative distribution functions (CDF).

Source: [2] [6] [9]-[11].

• $\gamma \left(\alpha ,x\right)$ : incomplete Gamma function, defined by $\gamma \left(\alpha ,x\right)={\displaystyle {\int }_0^x{y}^{\alpha -1}{\text{e}}^{-y}\text{d}y}$ , used to express the cumulative distribution function of the Gamma distribution.

2.1.3. Survival Functions

The survival function $S\left(t\right)=1-F\left(t\right)$ represents the probability that the neonatal mortality rate exceeds a given value. It serves as a complementary tool to the cumulative distribution function for mortality data analysis. Table 3 presents these functions for the five distributions.

2.1.4. Instantaneous Hazard Functions

The instantaneous hazard function, denoted $h\left(t\right)=\frac{f\left(t\right)}{S\left(t\right)}$ , measures the intensity of mortality at a given time, conditional on survival up to that time. It allows analyzing the temporal dynamics of the neonatal death risk. Table 4 presents these functions for the five models studied.

Table 3. Survival functions $S\left(t\right)$ .

Source: [2] [6] [9]-[11].

Table 4. Instantaneous hazard functions (Hazard function $h\left(t\right)$ ).

Source: [2] [6] [9]-[11].

2.1.5. Quantile or Simulation Functions

The quantile function, denoted $Q\left(p\right)=F^{-1}\left(p\right)$ , is the inverse of the cumulative distribution function. It plays a fundamental role in random number generation for simulation studies and in the calculation of empirical quantiles (median, quartiles, etc.). Table 5 presents these functions for the five distributions considered.

Table 5. Quantile function $Q\left(p\right)$ .

Source: [2] [6] [9]-[11].

• $W_{-1}(\cdot )$ : negative branch of the Lambert W function, used to express the quantile of the Lindley distribution.

Obtaining quantiles for the Gamma, Lindley and HLPKD distributions is not straightforward. For the Gamma and HLPKD distributions, the implicit equation requires numerical resolution (Newton-Raphson method or Brent’s algorithm). For the Lindley distribution, although an analytical formula exists via the Lambert function, its evaluation also requires numerical methods. These procedures are implemented in statistical software (R, Python, MATLAB) through specialized functions such as qgamma(), qlindley() or nonlinear equation solvers.

2.2. Variable and Data Collection

The main variable in this study is the neonatal mortality rate in Burkina Faso, observed annually over the period 1969 to 2023. This indicator, expressed per 1000 live births, measures the number of deaths occurring within the first 28 days of life. It constitutes an essential marker of maternal and child health, as well as a reflection of a country’s health development level.

The data used come from annual estimates compiled by the United Nations Inter-agency Group for Child Mortality Estimation (UN IGME). This organization gathers and harmonizes data from several international reference institutions, including UNICEF, the World Health Organization (WHO), the World Bank and the United Nations Population Division [12]. The CEIC database disseminates these data for analytical use. The values are model-based estimates produced by UN IGME, which combines survey, census and civil registration data within a statistical modeling framework. No additional transformation or preprocessing was applied to the extracted values.

The resulting time series covers 55 consecutive years, from 1969 to 2023, providing sufficient historical depth for robust statistical modeling. The data were extracted directly from the CEIC platform as provided by UN IGME, without any smoothing, interpolation or adjustment. Using these annual aggregated data, rather than individual survival data, directs the analysis toward modeling the distribution of the rates themselves. This approach allows studying the inter-annual variability of the indicator and identifying the probability distribution that best describes its overall statistical behavior.

2.3. Data Simulation Methods

Simulation studies using the Monte Carlo method occupy a fundamental place in the validation of statistical models, particularly when comparing several competing distributions. As noted by [13], simulation allows evaluating estimator behavior under controlled conditions, where asymptotic developments reach their limits, especially for moderate sample sizes. In the context of modeling neonatal mortality in Burkina Faso, where real data cover 55 years of observation, simulation becomes an indispensable tool for assessing the reliability and stability of the different models considered.

The main objective of this simulation study is to evaluate the performance of maximum likelihood estimators for each model under different sampling conditions. More specifically, it aims to examine estimator behavior in terms of bias, precision and convergence as sample size varies. This approach, inspired by the work of [6], helps determine whether the models produce stable and reliable estimates, a necessary condition for their application to real mortality data.

2.3.1. Simulation Design

The Monte Carlo simulation design was developed to cover a broad spectrum of realistic situations. Five sample sizes were selected: $n=40,80,120,160$ and 200. This choice allows studying estimator behavior from relatively small to larger samples, thus providing a complete view of asymptotic convergence.

To ensure result robustness, each scenario was repeated $N=500$ times, following the recommendations of [14], who advocate a sufficient number of replications to obtain stable estimates of performance measures.

For each model, the optimization was performed using the L-BFGS-B method implemented in the optim() function of R [15]. The starting values for the optimization were set to the true parameter values used for data generation. The parameter bounds were set as follows: for HLPKD, $\alpha \in \left[{10}^{-4},20\right]$ , $\beta \in \left[{10}^{-4},20\right]$ , $\theta \in \left[{10}^{-4},20\right]$ ; for Weibull, $\alpha \in \left[{10}^{-4},100\right]$ , $\beta \in \left[{10}^{-4},100\right]$ ; for Gamma, $\alpha \in \left[{10}^{-4},50\right]$ , $\beta \in \left[{10}^{-4},50\right]$ ; for Gompertz, $\alpha \in \left[{10}^{-4},30\right]$ , $\eta \in \left[{10}^{-4},30\right]$ ; for Lindley, $\theta \in \left[{10}^{-4},50\right]$ . Multiple starting points (the true values, and 0.7 and 1.3 times the true values) were used to avoid local optima. Convergence was declared when the optimization algorithm returned a convergence code of 0 (indicating successful convergence).

2.3.2. Generation of Simulated Data

For each model and each parameter combination, simulated samples were generated using the inverse transformation method [16], based on the quantile function $Q\left(p\right)=F^{-1}\left(p\right)$ presented in Table 5. For distributions lacking an explicit quantile function, numerical solving algorithms were employed (Brent’s method implemented in the uniroot function of R).

2.3.3. Estimation and Evaluation Criteria

For each simulated sample, the parameters of each model were estimated using the maximum likelihood method. This method, widely preferred for its optimal asymptotic properties (consistency, efficiency, asymptotic normality), was implemented using the optimization algorithm provided in the optim() function of R software [15].

Estimator performance was assessed using the following indicators, commonly used in the statistical literature [17]:

• Mean: evaluates the central tendency of estimates across all replications.

$\widehat{\theta }=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{\widehat{\theta }}_i}$

• Bias: measures the systematic deviation between the estimator and the true parameter value.

$\text{Bias}\left(\widehat{\theta }\right)=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N\left({\widehat{\theta }}_i-{\theta }_0\right)}$

where ${\widehat{\theta }}_i$ is the parameter estimate for the $i$ -th replication and ${\theta }_0$ is the true value. A positive bias indicates systematic overestimation of the parameter, while a negative bias indicates underestimation. A value close to zero indicates unbiasedness of the estimator.

• Relative bias (RB): expresses the bias as a proportion of the true value, facilitating comparisons across parameters of different scales.

$\text{RB}\left(\widehat{\theta }\right)=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N\left|\frac{{\widehat{\theta }}_i-{\theta }_0}{{\theta }_0}\right|}$

Relative bias, always positive due to the absolute value, measures the average magnitude of deviations as a proportion of the true value, regardless of their direction.

• Mean squared error (MSE): measures the overall precision of the estimator by combining variance and bias.

$\text{MSE}\left(\widehat{\theta }\right)=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{\left({\widehat{\theta }}_i-{\theta }_0\right)}^2}$

• Root mean squared error (RMSE): expresses the MSE in the same unit as the parameter, facilitating interpretation.

$\text{RMSE}\left(\widehat{\theta }\right)=\sqrt{\text{MSE}\left(\widehat{\theta }\right)}$

• Average confidence interval length (AL): for confidence levels of 90% and 95%, computed from the observed information matrix.

${\text{AL}}_{90\text{\%}}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^{N}2\times 1.645\times \widehat{se}\left({\widehat{\theta }}_i\right)}$ and ${\text{AL}}_{95\text{\%}}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^{N}2\times 1.960\times \widehat{se}\left({\widehat{\theta }}_i\right)}$

• Coverage probabilities (CP): percentages of times the constructed confidence interval contains the true parameter value. These probabilities should be close to the nominal levels (90% and 95%) to validate the use of asymptotic confidence intervals.

All simulations were performed using R software (version 4.3.2), a statistical computing environment widely recognized for its flexibility and rich set of libraries [15].

2.4. Estimation Methods for Real Data

After evaluating estimator performance through simulation, we now apply the different distribution models to real neonatal mortality data from Burkina Faso. The objective is to identify the distribution that best fits the observed data.

2.4.1. Parameter Estimation

For each considered distribution, parameters are estimated using the maximum likelihood method. Denoting $x_1,x_2,\cdots ,x_n$ the observed sample (neonatal mortality rates), the likelihood function $L\left(\theta \right)$ is maximized numerically using the optim () function in R software [15].

2.4.2. Selection and Performance Criteria

The goodness-of-fit of the different models is evaluated using the following indicators, consistent with the literature [6]:

• Standard error (SE): The standard error of the estimator $\widehat{\theta }$ is obtained from the square root of the estimated parameter variance. Within the maximum likelihood framework, it is computed using the information matrix and is directly provided by R software during estimation.

• Log-likelihood (lnL): The log-likelihood value at the convergence point of the optimization algorithm. Higher values indicate better fit.

• Akaike information criterion (AIC): A measure of relative model quality that penalizes complexity (number of parameters $k$ ):

$\text{AIC}=-2\ln L+2k$

A model with a lower AIC is preferred.

• Corrected Akaike information criterion (CAIC): A corrected version of AIC that accounts for sample size $n$ :

$\text{CAIC}=-2\ln L+k\left(\ln n+1\right)$

• Bayesian information criterion (BIC): Also known as the Schwarz information criterion, BIC penalizes complexity more heavily than AIC:

$\text{BIC}=-2\ln L+k\ln \left(n\right)$ where $k$ is the number of model parameters and $n$ is the sample size. As with AIC, a model with a lower BIC is preferred.

• Kolmogorov-Smirnov test (KS): The KS test statistic measures the maximum distance between the empirical cumulative distribution function $F_n\left(x\right)$ and the theoretical cumulative distribution function $F\left(x|\widehat{\theta }\right)$ of the estimated model:

$D={\text{sup}}_x\left|F_n\left(x\right)-F\left(x|\widehat{\theta }\right)\right|$

The associated p-value indicates whether the null hypothesis that the data follow the theoretical distribution is rejected ($p<0.05$ ) or not.

• Mean absolute error (ASAE): The mean absolute error between observed and theoretical probabilities, computed as:

$\text{ASAE}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\left|\widehat{F}\left(x_i\right)-F\left(x_i|\widehat{\theta }\right)\right|}$ where $\widehat{F}\left(x_i\right)=\frac{i}{n}$ is the empirical cumulative distribution function.

• Cramér-von Mises test (W): The Cramér-von Mises test statistic measures the quadratic distance between empirical and theoretical distributions:

$W={\displaystyle {\sum }_{i=1}^n{\left[F\left(x_i|\widehat{\theta }\right)-\frac{2i-1}{2n}\right]}^2}+\frac{1}{12n}$

The associated p-value ($W$ p-value) is provided to assess model adequacy.

• Anderson-Darling test (A): The Anderson-Darling test statistic gives more weight to distribution tails:

$A=-n-\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\left(2i-1\right)\left[\ln F\left(x_i|\widehat{\theta }\right)+\ln \left(1-F\left(x_{n-i+1}|\widehat{\theta }\right)\right)\right]}$

The associated p-value ($A$ p-value) complements the goodness-of-fit assessment.

To identify the best model, we follow a hierarchical approach:

• Models with a KS test p-value greater than 0.05 are considered adequate (null hypothesis not rejected).

Among these models, those with the lowest AIC and CAIC values are preferred;

• The $W$ and $A$ statistics (with their p-values) along with ASAE confirm the ranking.

3. Results

This section presents the study results, organized into two parts. The first part presents the Monte Carlo simulation results conducted to evaluate the stability and performance of maximum likelihood estimators for each of the five considered models. The second part presents the results of fitting these models to real neonatal mortality rate data from Burkina Faso, along with the various goodness-of-fit measures and corresponding graphical visualizations.

3.1. Simulation Results

The simulation study enabled evaluating the behavior of maximum likelihood estimators for the five distributions under different sampling conditions. For each model, reference parameter values were set to generate random samples of increasing sizes ($n=40,80,120,160,200$ ). For each combination, $N=500$ replications were performed. All indicators presented in the methodology were computed to assess model stability. Results are presented in Tables 6-10. Figures 1-5 illustrate the density and instantaneous hazard functions for each model, plotted using the parameters employed in the simulation.

Table 6. Simulation results for the HLPKD model with $\alpha =2.0$ , $\beta =1.0$ and $\theta =1.5$ .

Table 7. Simulation results for the Weibull model with $\alpha =2.0$ and $\beta =1.50$ .

Table 8. Simulation results for the Gamma model with $\alpha =2.2$ and $\beta =1.4$ .

Table 9. Simulation results for the Gompertz model with $\alpha =0.5$ and $\eta =1.3$ .

Table 10. Simulation results for the Lindley model with $\theta =2.0$ .

The examination of the simulation tables reveals several common trends across all models studied. For each distribution, the mean of the estimates approaches the theoretical parameter value as sample size increases. This observation is particularly clear for the Weibull model, where the mean of $\widehat{\alpha }$ decreases from 2.0782 ($n=40$ ) to 2.0061 ($n=200$ ), thus approaching the reference value $\alpha =2.0$ . Similarly, the Lindley model shows marked convergence, with the mean of $\widehat{\theta }$ evolving from 2.0589 to 1.9990 while the theoretical value is $\theta =2.0$ .

Bias and relative bias systematically decrease as $n$ increases. For the HLPKD model, the bias of $\widehat{\alpha }$ drops from 0.2775 ($n=40$ ) to 0.0463 ($n=200$ ), while the bias of $\widehat{\theta }$ decreases from 0.9885 to 0.0762. This bias reduction confirms the asymptotic consistency property of maximum likelihood estimators. The two-parameter models (Weibull, Gamma, Gompertz) as well as the one-parameter model (Lindley) generally exhibit lower bias than the three-parameter HLPKD model, which is attributable to the greater estimation complexity of the latter.

Mean squared errors (MSE) and their square roots (RMSE) also decrease with sample size. For the Gamma model, the MSE of $\widehat{\alpha }$ falls from 0.3479 ($n=40$ ) to 0.0457 ($n=200$ ). This reduction indicates improved estimation precision when more observations are available. The Weibull and Lindley models show the smallest MSE values, reflecting better estimator stability for these distributions.

Average confidence interval lengths (AL₉₀ and AL₉₅) shrink as $n$ increases, reflecting improved estimation precision. For the Weibull model, the AL₉₀ of $\widehat{\alpha }$ decreases from 0.8476 ($n=40$ ) to 0.3646 ($n=200$ ). Coverage probabilities (CP₉₀ and CP₉₅) stabilize around the nominal levels of 90% and 95%, particularly for sample sizes $n$ greater than or equal to 120. This finding confirms the validity of asymptotic confidence intervals for moderate sample sizes.

Comparison among the five models shows that two-parameter distributions (Weibull, Gamma, Gompertz) and the one-parameter distribution (Lindley) exhibit better estimator stability than the three-parameter HLPKD model, for which MSE values remain higher even at $n=200$ . The Weibull model stands out with the lowest MSE and the best coverage probabilities, followed by the Lindley and Gamma models. The Gompertz model shows intermediate performance, with notable estimator stabilization starting at $n=120$ .

Figure 1. Density functions (PDF) and mortality intensity function $h\left(t\right)$ of the HLPKD model.

Figure 2. Density functions (PDF) and mortality intensity function $h\left(t\right)$ of the Weibull model.

Figure 3. Density functions (PDF) and mortality intensity function $h\left(t\right)$ of the Gamma model.

The HLPKD model plots illustrate a highly flexible density, capable of capturing pronounced skewness and marked peaks. The mortality intensity functions show diverse behaviors: rapid growth, linear progression or sustained increase depending on the parameters. This variety confirms that HLPKD can adapt to contexts where the distribution of mortality rates exhibits complex and non-monotonic patterns.

Figure 4. Density functions (PDF) and mortality intensity function $h\left(t\right)$ of the Gompertz model.

Figure 5. Density functions (PDF) and mortality intensity function $h\left(t\right)$ of the Lindley model.

Weibull densities display shapes ranging from simple decay to bell-shaped distributions, depending on the shape and scale parameters. The mortality intensity functions illustrate varied dynamics: decreasing, constant or increasing. This versatility explains Weibull’s popularity in survival analysis, as it can represent situations where the intensity of mortality rates decreases, remains stable or increases.

Density curves show that the Gamma distribution can take various forms depending on the parameters, ranging from rapid decay to a unimodal distribution. The associated mortality intensity functions reflect increasing intensity, which corresponds to a model suitable for describing phenomena where the intensity of mortality rates gradually increases. This confirms Gamma’s flexibility for modeling distributional shapes of mortality rates.

Gompertz model densities highlight strongly asymmetric distributions, with high probability concentrated near zero. The mortality intensity functions reveal exponential growth, sometimes very steep, depending on the values of $\alpha$ and $\eta$ . This behavior is consistent with the nature of the model, designed to represent mortality intensity that increases strongly with the values of the rates.

Lindley model densities show a high concentration of probability near zero, with rapid decay. The mortality intensity functions exhibit an increase that tends toward a plateau, reflecting stable asymptotic behavior. This simple, single-parameter model stands out for its ability to provide robust fits in contexts where mortality rates are concentrated at low values.

3.2. Application to Real Data (Neonatal Mortality Rate in Burkina Faso, 1969-2023)

Analysis of the real data began with a descriptive examination of the neonatal mortality rate time series. Figure 6 presents the evolution of this indicator over the period 1969-2023. Parameter estimates and their standard errors for the five models are reported in Table 11. Selection and performance criteria, including log-likelihood, AIC, CAIC, BIC, the Kolmogorov-Smirnov (KS) statistic and its p-value, mean absolute error (ASAE), the Cramér-von Mises ($W$ ) and Anderson-Darling ($A$ ) statistics along with their respective p-values, are presented in Table 12. Figures 7-11 show the profile likelihood functions for the parameters of each model. Finally, Figures 12-14 present the fitted density functions, cumulative distribution functions and P-P plots for all models compared against the observed data.

Figure 6. Evolution of the neonatal mortality rate in Burkina Faso from 1969 to 2023.

The graph (Figure 6) illustrates the evolution of the neonatal mortality rate in Burkina Faso over a long period, from the late 1960s to 2023. A general downward trend is observed, with initially high values (around 70 deaths per 1000 live births) progressively decreasing to below 30 deaths by the end of the period. This trajectory reflects a notable improvement in maternal and child health conditions, as well as a progressive strengthening of health systems.

Table 11. Parameter estimates and standard errors for the five models.

Note: For the Weibull model, the scale parameter estimate remained at the boundary of the optimization domain ($\widehat{\beta }=30$ ) even after increasing the upper bound. This reflects the inability of the two-parameter Weibull distribution to capture the empirical mean (44.92) and is consistent with its rejection by all goodness-of-fit tests.

Table 11 presents the parameter estimates and their standard errors for the five models. The HLPKD model provides estimates but with substantial uncertainty, particularly for parameter $\theta$ . The Weibull model appears more stable, with precise estimates for $\alpha$ and $\beta$ . The Gamma distribution yields relatively high parameters with moderate standard errors, reflecting good reliability. The Gompertz model shows a very precise estimate for $\alpha$ , but greater uncertainty for $\eta$ . Finally, the Lindley model, with its single parameter, stands out with a very precise and robust estimate.

Table 12. Selection and performance criteria for the five models.

Table 12 presents the selection and performance criteria for the five models studied. The log-likelihood (lnL), AIC, CAIC and BIC values show that the HLPKD and Gamma models achieve the lowest scores, reflecting a better fit to the data compared to the other models. In contrast, the Weibull and Lindley models display substantially higher values, indicating a less satisfactory fit.

Goodness-of-fit tests (Kolmogorov-Smirnov, Anderson-Darling and Cramér-von Mises) confirm this trend. The HLPKD and Gamma models exhibit low test statistics and high p-values (above 0.80 for the KS test), demonstrating that the data do not deviate significantly from these distributions. The Gompertz model achieves intermediate results, with acceptable but less convincing p-values. Conversely, the Weibull and Lindley models are rejected by the tests, with p-values close to zero.

Overall, the convergence of information criteria and goodness-of-fit tests indicates that the HLPKD and Gamma models offer the best performance for describing the neonatal mortality data. Between the two, Gamma stands out for greater parameter stability, while HLPKD provides additional flexibility. These results show that Gamma constitutes the most reliable model, with HLPKD as a competitive alternative.

Table 13. Bootstrap-corrected goodness-of-fit p-values ($B=1000$ parametric bootstrap replications).

Figure 7. Profile likelihood functions for the parameters of the HLPKD model.

Table 13 reports the bootstrap-corrected p-values. The bias correction lowers the p-values for all models relative to the classical asymptotic results, as expected when the tested parameters are estimated from the same sample. Gamma and HLPKD remain the only two models not rejected by the KS and $W$ tests at the 5% level, confirming their superiority over the three other distributions. The Anderson-Darling test, which is more sensitive to tail departures, yields p-values slightly below 0.05 for both models ($p=0.021$ for HLPKD and $p=0.039$ for Gamma). Given the bootstrap variability and the fact that the KS and $W$ tests do not reject either model at the 5% level, these results do not constitute strong evidence against the adequacy of HLPKD and Gamma. The overall conclusion remains that both models provide a satisfactory description of the data. Gompertz is rejected by the bootstrap-corrected tests at the 5% level ($p\le 0.017$ ), confirming that it does not provide as good a fit as HLPKD or Gamma. Weibull and Lindley are clearly rejected, with bootstrap p-values of zero.

Figure 8. Profile likelihood functions for the parameters of the Weibull model.

Figure 9. Profile likelihood functions for the parameters of the Lindley model.

Figure 10. Profile likelihood functions for the parameters of the Gamma model.

Figure 11. Profile likelihood functions for the parameters of the Gompertz model.

These graphs show the profile log-likelihoods for the different parameters of the models studied. Each curve illustrates how the likelihood changes as one parameter varies while the others are held at their optimal values. The peak of the curve, marked by a red point and a vertical line, corresponds to the maximum likelihood estimate of the parameter under consideration. The presence of reference horizontal lines allows identifying the associated confidence intervals.

Figures 12-14 respectively present the density functions (PDF), cumulative distribution functions (CDF) and P-P plots for the five models fitted to the real data. These representations allow visually comparing the goodness-of-fit of each distribution: the PDFs show the correspondence between theoretical curves and the empirical histogram, the CDFs illustrate the proximity between observed and theoretical cumulative distributions, and the P-P plots indicate the agreement between expected and observed values. Overall, the HLPKD and Gamma models appear most consistent with the data, Gompertz occupies an intermediate position, while Weibull and Lindley show substantially less satisfactory fits.

Figure 12. Estimated PDFs for the study models.

Figure 13. Estimated CDFs for the study models.

Figure 14. P-P plots for the study models.

4. Discussion

The results obtained in this study allow assessing the comparative performance of five probabilistic distributions applied to annual neonatal mortality rates in Burkina Faso. Monte Carlo simulations showed that maximum likelihood estimators converge to stable values as sample size increases, but with notable differences across models. Weibull and Lindley stood out for their robustness even with modest sample sizes, while Gamma and HLPKD required larger samples to reduce bias. It is important to clarify that the intensity functions presented in this study are mathematical properties of the fitted distributions applied to mortality rates; they are not estimates of individual neonatal survival risks over time. They describe the intensity of the distribution of annual rates. Gompertz, for its part, showed precise estimation of the shape parameter but greater variability for the growth parameter.

The application to real data confirmed these trends. Information criteria and goodness-of-fit tests indicated that HLPKD and Gamma models offer the most satisfactory fits, with high p-values and small discrepancies between theoretical distributions and observed data. Gompertz showed intermediate performance, while Weibull and Lindley were rejected by the tests, their fits being substantially less consistent with the data. These findings align with those of [6], who demonstrated HLPKD’s flexibility in biomedical contexts, and confirm that this model can also be relevant for aggregated neonatal mortality data. These conclusions were confirmed by a parametric bootstrap correction of the goodness-of-fit p-values, which additionally revealed that the Gompertz model, initially classified as intermediate, is in fact rejected at the 5% level once the bias induced by parameter estimation is accounted for, and that both Gamma and HLPKD show a mild but statistically detectable departure in the tails of the distribution under the Anderson-Darling test.

The comparison between our results and those of [6] also reveals an interesting contrast. In their study, HLPKD clearly outperformed other distributions on individual survival data from cancer patients, due to its ability to capture complex hazard functions. In our work, HLPKD remains competitive but does not significantly improve fit compared to the simpler Gamma model. This difference can be explained by the aggregated nature of our data, which likely smooths individual variations and reduces the advantage of high parametric flexibility. Previous work on other health indicators [5] has already emphasized that model choice depends as much on data type as on the research question.

Comparison with previous studies highlights the value of diversifying models used in survival analysis. Studies by [7] and [8] extensively employed the Weibull distribution to characterize neonatal risk in Africa, but our results show that this model is not always best suited to aggregated data. Similarly, recent developments around the Gamma distribution applied to health statistics [5] resonate with our findings, since this model proved both performant and stable. The introduction of the Lindley distribution in medical contexts [4] illustrates its simplicity, but our analyses indicate that this parsimony does not guarantee satisfactory fit in the case studied.

This study highlights that the Gamma distribution emerges as the most reliable model for describing annual neonatal mortality rates in Burkina Faso, with HLPKD as a competitive alternative offering additional flexibility. The relevance of these two models lies in their ability to balance goodness-of-fit with estimation stability, which is essential for robust statistical interpretation. These results open the way for broader use of flexible models in mortality data analysis, particularly in contexts where time series are long and inter-annual variability must be properly represented.

The limitations of this study must be clearly stated. The real sample size ($n=55$ ) remains modest for fitting three-parameter models, which may explain the uncertainty observed in HLPKD estimates. Moreover, the use of annual aggregated data does not capture intra-annual variability or individual heterogeneity in risks. A strong downward trend is observed over the study period (from 71.1 to 24.5 deaths per 1000 live births), and the annual observations exhibit serial dependence. The present analysis does not model this temporal structure; instead, it focuses on the marginal distribution of the rates, treating the 55 annual values as a sample from a fixed probability distribution. The conclusions regarding model fit are therefore specific to the distributional properties of these aggregated rates and do not imply that the annual observations are independent or identically distributed over time. Additionally, the reference parameter values used for simulation, while chosen to be representative of observed magnitudes, do not cover all possible configurations.

Despite these limitations, this study brings several new contributions. It constitutes the first comparative application of HLPKD to aggregated neonatal mortality data in sub-Saharan Africa. It shows that model performance depends closely on data structure, and that a parsimonious model like Gamma can be as effective as a more complex one when observations are smoothed. From a public health perspective, the results indicate that the decline in neonatal mortality in Burkina Faso follows a dynamics well described by a Gamma distribution, which can help anticipate future trends.

Further research could explore alternative two-parameter models, such as the log-normal distribution, or approaches accounting for the temporal structure of the data. Incorporating covariates (year, socioeconomic factors) within a parametric regression framework could also enrich the analysis. An extension to individual neonatal survival data, when available, would allow verifying whether HLPKD’s flexibility regains its full advantage.

5. Conclusion

This study aimed to evaluate the ability of five probabilistic distributions namely HLPKD, Weibull, Gamma, Gompertz and Lindley to model annual neonatal mortality rates in Burkina Faso over the period 1969-2023, combining a simulation-based assessment of maximum likelihood estimator stability with fitting to real data. Simulations showed that the Weibull, Gamma, Gompertz and Lindley models converge rapidly to stable estimates, while the more complex HLPKD requires larger samples to achieve comparable precision. Application to real data identified the Gamma and HLPKD models as the best-fitting, with high p-values and small discrepancies between theoretical distributions and observations; the Gamma model standing out for its parsimony and estimation stability, while HLPKD offers additional flexibility but suffers from greater uncertainty in its parameter $\theta$ . The parametric bootstrap goodness-of-fit results confirm that Gamma and HLPKD are the only models not rejected by the KS and Cramér-von Mises tests at the 5% level, further supporting their adequacy. Whereas HLPKD had proven superior on individual survival data in the work of [6], its advantage diminishes on aggregated data where a simpler model like Gamma offers a more favorable parsimony-to-fit ratio. From a public health perspective, the results indicate that the decline in neonatal mortality in Burkina Faso follows a dynamics well described by a Gamma distribution, confirming a continuous improvement in health conditions. Study limitations relate to the modest size of the real sample and the use of aggregated data, which smooth individual variations. This study therefore identifies the Gamma model as the most reliable in this context and confirms that the increased flexibility of HLPKD does not translate into a decisive gain on aggregated data, underscoring the importance of aligning model complexity with the structure of observed data.

Author Contributions

Daouda Traoré conceived the idea of the article and wrote the entire manuscript. Issouf Traoré contributed to verifying the objectives and assessing whether they were achieved. Ibrahim Traoré verified the mathematical foundations.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References