Bayesian Shared Frailty Models for Time to First Birth of Married Women in Ethiopia: Using EDHS 2016

Introduction The primary effect of the fertility process is the birth of the first child. The ages at which women establish marital union and give their first birth depend on and result in varying demographic features. This research demonstrates how to examine the effect of numerous factors on married women's delay to first birth in Ethiopia using Bayesian parametric models with gamma shared frailty distribution. Methods This study analyzed data from the 2016 EDHS on factors related to the time of married women to first birth. A sample of 8810 married women from all parts of Ethiopia participated in the study. The Akaike information criterion (AIC) and Bayesian information criterion (BIC) were used to compare several parametric models with gamma shared frailty distributions to find the best model (BIC). Finally, when the prior data was taken into account, the chosen model was proven to be accurate (Bayesian approach). Results The median survival time for the first birth after marriage is 24 years (95% CI; 23.4, 25.3). The result shows that the place of residence, the access to media, the level of education of the mother, the education level of the husband, the use of the head of the contraceptives, and the sex of the household are statistically associated with the time to first birth of married women. The Weibull-gamma shared frailty model under the Bayesian approach was found to be the best model that fit the time to first birth data in this study. The result also showed that there is heterogeneity between regions of married women. Conclusion To slow the increase in the Ethiopian population, families must be taught how to use contraception, and rural populations must be educated on the necessity of increasing the length of the first birth gap rather than encouraging early marriage. In general, attempts to reduce fertility by raising the age of the first marriage must consider the social and cultural settings in which marriage takes place. On the other hand, the campaign against early marriage should focus on the sociocultural, physiological, and psychological effects, as well as the reduction of reproduction.


Background
The primary effect of the fertility process is the birth of the first child [1]. The ages at which women establish marital union and give their first birth depend on and result in varying demographic characteristics [2]. Every year, between the ages of 15 and 19, approximately 16 million young girls give birth. Young mothers are responsible for roughly one in every ten childbirths worldwide, with developing countries accounting for 95% of this [3]. Adolescent fertility, especially among the youngest age groups, poses serious health hazards to both the mother and the child [4][5][6]. The health and social causes of early-age pregnancy and childbearing are well documented, and they become major issues for the well-being of a mother and her child [7]. In the poorest parts of the world, 20% of females give birth before 18 years old, and in Kenya, this number jumps to more than a third (35%) [8]. Furthermore, women less than 15 years are five times higher chance to die, and those between the ages of 15 and 19 are twice as likely to die during pregnancy or childbirth as those between the ages of 20 and 24 [8,9]. While sub-Saharan African countries dropped their teenage birth rate from 140 to 101 births per 1000 women aged 15-19 years over the Millennium Development Goals era , other regions such as South Asia, North Africa, and the Middle East saw far greater reductions [10]. In 2013, Nigeria's median age at the first birth was 20 years [11]. More than a third (34%) of Ethiopian women aged 20 to 49 years give birth before the age of 18 years and 54% at the age of 20 [12]. The high frequency of poor mother-and-child health in third-world countries can be attributed to several variables. Deep-seated sociocultural and spiritual practices, illiteracy, and low income are all part of this [11,13,14]. Socioeconomic characteristics were consistently identified as a predictor of age at the first birth. Many studies have shown that women with no or low levels of education had a higher chance of having their first child at a younger age than women with higher levels of education [15][16][17]. A study carried out in Nigeria on Bayesian semiparametric multilevel survival modelling of age at the first birth showed that variation in age at the first birth in Nigeria is determined more by individual households than by community and that substantial geographical variations exist in timing of the first birth [18]. Even though the age at which a child is born has a significant impact on maternal and infant survival, both individual and cumulative levels of fertility, and comprehensive consequences for women's roles and social changes in general, there are few studies on this topic in Ethiopia on adolescent pregnancy [19] and age at first birth after marriage [20]. However, those studies did not consider unobserved heterogeneity (shared frailty). If the data come from different groups or the nature of the data has repeated measures, heterogeneity between individuals should be taken into account [21]. Numerous issues could arise if heterogeneity is ignored, including an overestimation of the relative hazard rate, inaccurate estimates of the regression coefficients, and a tendency for the regression parameter estimate to approach zero [21]. Frailty provides a more precise estimate of the parameters compared to standard AFT models [22].
In addition, Bayesian approaches have recently been employed as the best method over classical approaches in numerous research investigations, particularly in the field of medicine. One of the problems is that traditional methods rely on asymptotic considerations that are usually only true for large datasets [23]. In the Bayesian approach, no assumption is made as to the shape of the percentile distribution; rather, the data themselves specify the distribution and the Bayesian approach has the possibility of improving the precision of the results by introducing external information in terms of the a priori distribution [23]. Understanding the time and factors that influence the first birth in the country would aid in the development of effective measures to improve mother and child health. Thus, taking into account the aforementioned limitations, this study is meant to estimate the time to first birth and find predictors among all married women in Ethiopia using a Bayesian approach while accounting for random impact (shared frailty) between regions of Ethiopia.

Study
Area. This research was conducted in Ethiopia, which is Africa's second most populated country, after Nigeria, and is located in the Horn of Africa. Ethiopia has nine regional states (Tigray, Afar, Amhara, Oromiya, Somali, Benishangul-Gumuz, Southern Nations Nationalities and People (SNNP), Gambela, and Harari), as well as two city governments (Addis Ababa and Dire Dawa) [24].

Study
Design. Enumeration areas (EAs) were the sampling units for the first stage of the 2016 EDHS sample, which was chosen using a stratified two-stage cluster design. During the 2007 census, each kebele (ward) was partitioned into census enumeration areas (EAs) to make the census easier to administer. There were 645 EAs in the sample, with 187 in urban regions and 437 in rural areas. The second stage of the sampling involved households. The interviews were open to all women aged 15 to 49 years. Individual interviews were done with 16,583 eligible women from the interviewed households; complete interviews were conducted with 15,683 [25]. In the current study, 8810 married women from nine regions and two city administrations were included ( Figure 1).

Study Population.
The current study used married women data from the Ethiopian Demographic and Health Surveys (EDHS) conducted in 2016.

Variables Included in the Study
The response variable is the time to first birth of married women in Ethiopia. It is measured as the length of time from birth to the age at the first birth, which is measured in years.

Explanatory
Variables. The expected explanatory variables included socioeconomic, demographic, health, and environmental factors (Table 1).

Statistical Analysis.
After obtaining the consent letter for use of the measured DHS, the dataset was acquired through the website https://dhsprogram.com. A data extraction technique was used to extract variables from the EDHS 2016 dataset for children and individual women. The study population was characterized using descriptive measures like graphs and frequency tables after editing and coding. The Kaplan-Meier (K-M) and the log-rank test were calculated to show the time of the first birth and to compare the survival time between covariates, respectively. We first analyzed our data using Cox proportional hazard [26], accelerated failure time [26], and parametric shared frailty models. Then, using AIC and BIC, the best model was selected. Finally, the Bayesian parametric gamma shared frailty model was fitted, and the best model was selected using DIC criteria. Data were entered and cleaned using SPSS-22 and analyzed using WinBugs1.4.
2.6. Survival Analysis. Survival analysis is a set of statistical processes for data analysis for which the outcome variable of interest is the time until an event occurs. By the time, this means year, month, week, or days from the beginning of follow-up of an individual until an event occurs. By event, it means death, disease incidence, relapse from remission, recovery (e.g., return to work), or any designated experience 2 Computational and Mathematical Methods in Medicine of interest that may happen to an individual. Despite these, survival models have a long history in the biostatistical and medical literature [27].

The Cox Proportional Hazard
Model. The Cox proportional hazards regression model is a model that assumes that the log of hazard rate is additively related to a function of covariates or the hazard rate is related multiplicatively to a function of covariates [26]. This model presents an equation for the hazard at time t for an individual with a particular collection of explanatory factors indicated by X, and it is usually expressed as where h 0 ðtÞ is the baseline hazard function at time t, X is the vector of values of the explanatory variables, and β = ðβ 1 , β 2 , ⋯, β k Þ is the vector of unknown regression parameters that are assumed to be the same for all individuals in the study, which measures the influence of the covariate on the survival experience.
2.8. Parametric Survival Model. If our survival time follows a specified probability distribution, we apply parametric survival models, and the parameters of that distribution depend on covariates. Among the popular parametric survival models, some of them are exponential, Gompertz, Weibull, log-normal, and log-logistic.
2.9. Exponential Distribution. The simplest model for the hazard function is to assume that it is constant over time.
The hazard of death at any time after the time origin of the study is then the same irrespective of the time elapsed. This famous property of the exponential distribution is known as "loss of memory" which requires that the age of the person does not affect future survival. Let t be the survival time that follows exponential distribution with parameter λ. Then, the pdf of t is and hðtÞ = λ.
2.11. Weibull Distribution. Weibull distributions are parameterized as both proportional hazard (pH) and accelerated failure time (AFT) models. The Weibull distribution is suitable for modeling data with monotone hazard rates that increase or decrease exponentially over time. For Weibull regression, λ is the scale parameter and γ is a shape parameter, and its probability distribution function is given by SðtÞ = e −λtγ and hðtÞ = λγt γ−1 .
2.12. The Log-Logistic Distribution. One limitation of the Weibull hazard is that it is a monotonic function of time. However, situations can arise in which the hazard function changes direction. In this situation, log-logistic is a preferable model.
2.13. The Log-Normal Distribution. The log-normal distribution is also defined for random variables that take positive values and so may be used as a model for survival data. A random variable T is said to have a log-normal distribution with parameters μ and δ log T, and it has a normal distribution with μ and variance σ. The probability density function of T is given by The survivor function of the log-normal distribution is SðtÞ = 1 − Φððlog t-μÞ/δÞ where Φð:Þ is the standard normal distribution function given by ΦðzÞ = ð1Þð2πÞ 1/2 Ð z −∞ exp ð−u 2 /2Þdu and hðtÞ = f ðtÞ/SðtÞ.

2.
14. Parametric Shared Frailty Model. The multivariate or shared frailty model is a conditional independence model in which frailty is common to all subjects in a cluster. In this study, the clusters are regions. The concept of frailty pro-vides a suitable way to introduce random effects into the model to account for association and unobserved heterogeneity. In its simplest form, frailty is an unobserved random factor that modifies multiplicatively the hazard function of an individual or cluster of individuals. [28] introduced the term frailty, and [29] promoted the model by its application to the multivariate situation on the incidence of chronic diseases in families. There are different frailty distributions such as gamma [22], inverse Gaussian [30], log-normal [31], and positive stable [32]. However, the gamma distribution is the most common and widely used in the literature for determining the frailty effect, which acts multiplicatively on the baseline hazard [30,33].
Let us have k observations and i subgroups (regions). Each region consists of n i observations and ∑ r i=1 ni = n, where n is the total sample size. The hazard rate for the k th individual in the i th region is given by Here, frailty Z is a random variable varying over the population decrease (Z < 1) or increases (Z > 1) the individual risk. The most important point here is that the frailty is unobservable. The respective survival function S, describing the fraction of surviving individuals in the study population, is given by Sðt, X, ZÞ = exp ð−Z Ð t 0 hoðuÞdu exp ðXβÞÞ. The cumulative hazard function is given by HðtÞ = Ð t 0 h 0 ðuÞdu. The main assumption of a shared frailty model is that all individuals in region i share the same value of frailty Z i ði = 1, 2, 3, ⋯, mÞ, and this is why the model is called the shared frailty model. In our study, the women assumed to share some common frailty in the region. The shared frailty (random) effect Z i follows a gamma distribution with mean one and variance θ, as defined in the density function in equation (6).
where θ > 0 indicates the presence of heterogeneity. Therefore, the high values of θ reflect a greater degree of heterogeneity among regions of pregnant women and a stronger association within regions. The associations within group members (regions) are measured by Kendall's, and gamma frailty distribution is given by 2.15. Bayesian Parametric Gamma Shared Frailty Models. In the Bayesian approach, the critical issue is identifying the prior distribution for each parameter in the model to get a best-fit posterior distribution. The prior distribution is a probability distribution that represents the prior information associated with the parameter of interest. In this study, 4 Computational and Mathematical Methods in Medicine we used noninformative priors for all parameters of interest. Let δ ij denote the censoring indicator variable, taking value 1 if the j th subject (j = 1, 2, ⋯, m i ) of the i th cluster (i = 1, 2, ⋯, n) fails and 0 otherwise; T = ðt 11 , t 12,: ::, t nm Þ′ and X = ðx 1 , x 2 , ⋯, x n Þ, where x i is the m × n matrix of covariates. Let D = ðX, T, δ ij , ZÞ denote the complete data, and let D obs = ðX, T, δ ij Þ denote the observed data. Here, we only allow for the right-censored survival data and assume that the censoring is noninformative. The complete data likelihood is given: where β, ∧ 0 , and θ are the regression coefficient, the baseline distribution parameters, and the gamma frailty distribution parameters, respectively.
Because the observed data likelihood is simply too intricate to work with, evaluating the joint posterior distribution analytically is challenging. We employ the Gibbs sampler to obtain samples from the joint posterior distribution to avoid this problem.
2.16. Prior Distributions. In Bayesian inference, prior elicitation may be the most important factor. The prior distributions for each parameter in the model must be specified first before we can do data analysis from a Bayesian perspective. We want our data information to dominate the prior distribution by assuming suitably noninformative priors for all parameters in this model because we have little prior information for all parameters to be estimated. We used independent imprecise normal priors with a mean of 0 and a variance of 1 * 10 5 for all regression coefficients. A gamma distribution with shape parameter 1 and scale parameter 0.001 (with mean 1000 and variance 1 * 10 6 ) is used to supply noninformative priors to the scale and shape parameters in the model. We use the hazy proper gamma prior distribution with shape parameter 0.001 and scale parameter 0.001 for their reciprocals for the shared frailty parameter (precision parameters for the random effects).
2.17. Posterior Distributions. By multiplying the prior distribution across all parameters φ, by the entire likelihood function Lðϕ/yÞ, the posterior distribution is derived. The posterior distribution of the model created informs all Bayesian inferential judgments. The inference is carried out by sampling from a posterior distribution until the posterior distribution converges. The main drawback of the Bayesian technique is that, in most circumstances, the whole form of the posterior distribution cannot be derived in a closed form, implying that the posterior density may not belong to the standard distribution. This is a difficult problem to solve. MCMC simulations will be used to solve these challenges.
We start with the joint density function of observable information Y and latent information Z to obtain the posterior densities.
Then, we will pretend that φ is a random variable with a prior distribution defined by πðϕÞ. The posterior distribution, which is derived using Bayes' theorem, is then used to make inferences. The posterior distribution ϕ is then calculated as follows: f ðϕ/yÞ = Lðϕ/yÞπðϕÞ/f ðyÞ = f ðy/ϕÞπðϕÞ/ f ðyÞ, where f ðyÞ = Ð Lðϕ/yÞπðϕÞdφ is a normalizing factor (constant). Thus, we have f ðφ/yÞ ∝ Lðϕ/yÞπðϕÞ.
The joint posterior density function of a parameter at a certain failure time is derived by using this expression and assuming independence between the prior density functions of the parameters.
5 Computational and Mathematical Methods in Medicine where g i ð:Þ indicates the prior density function with known hyperparameters of the corresponding argument for baseline parameters and frailty variance. π i ðβ i Þ is the prior density function for the regression coefficients β i for i = 1, 2, ⋯, r. Also, the equation gives the likelihood function: The graphical representation of multivariate distributions well known in connection with for models high dimensional contingency tables and Bayesian inference in expert system [34]. Graph theory is also useful in the study of Markov random fields, which are important in statistical mechanics and, more recently, spatial statistics and image analysis. The general layout of this research is presented in Figure 2.
2.18. MCMC Estimation Methods. The Bayesian approach, often known as "full probability modeling," adds probability theory to a model formed from substantive information and can, in theory, deal with truly complex circumstances. However, it must be emphasized that the computations may be challenging, with the specific challenge being the integration required to determine the posterior distributions of the quantities of interest when nonstandard prior distributions are utilized in the model. For many years, these integration issues limited Bayesian applications to rather simple cases. However, there has been a lot of improvement recently in Bayesian computation methods, which generally take advantage of modern computing power to do simulations known as Markov Chain Monte Carlo (MCMC) approaches. Although the shape of the posterior distribution of interest has no known algebraic form, the MCMC simulation performs the integration numerically rather than analytically by sampling from the posterior distribution of interest [34]. All posterior summary statistics will be generated as a result of this procedure (approximately). The most commonly used algorithms in MCMC applications are of two types, and they are the Metropolis algorithm and the Gibbs sampler.
The standard deviations offer measures of precision, whereas the means of the posterior samples provide point estimates for the model parameters. Along with estimating precision, the 95 percent intervals provide an alternative indicator of the covariates' influence. The difference between the mean of the sampled values (which we use as our posterior mean estimate for each parameter) and the true posterior mean is calculated as the MC error. As a general rule, the simulation should be run until the Monte Carlo error for each parameter of interest is less than 5% of the sample standard deviation.

Model Selection Criterion by the Bayesian Approach.
Computing posterior model probabilities is a typical method for comparing Bayesian models. We employ the Akaike information criterion (AIC), the Bayesian information criterion (BIC), the deviance information criterion (DIC), and the Bayes factor to compare the proposed models. These are the most prevalent Bayesian model evaluation techniques.
Given a class of competing models for a dataset, Akaike (1973) proposed that the model that minimizes be chosen.
where P is the number of model parameters. D ð b θÞ is a deviation estimate based on the posterior mean, which is equal to b θ = E ðθ/dataÞ. The deviation is equal to b θ = −2 Hyper parameter Parameters Data z Figure 2: A graphical illustration of the Bayesian frailty model. Note: ∧ 0 is the baseline distribution parameter, β is the regression parameter of each covariate, Z is the frailty distribution, and θ is the hyperparameter which follows a gamma distribution. 6 Computational and Mathematical Methods in Medicine log LðθÞ, where LðθÞ is the likelihood function of the model and θ is a vector of unknown parameters of the model. Schwarz proposed the Bayesian information criterion (BIC) (1978). [35,36] have [33] shown that, even asymptotically, AIC tends to overstate the number of parameters required. According to the Schwarz criterion, the model that minimizes has the best posterior probability.
The number of observations, or sample size, is denoted by the letter n.  Table 2 show that a total of 8810 women who got their first marriage were included in this study from nine regional states and two administrative cities. Among the  More than half of the women (68%) have not had current work. 61.5% of the total women were uneducated while 38.5% of them were attaining primary and above education level. Furthermore, 4459 (50.6%) of the women had the experience using contraceptive methods while 4351 (49.4%) of them had no experience of using a contraceptive. Regarding exposure to the media, 6280 (71.3%) of the women had no access to mass media and 2530 (28.7%) of them had access to mass media.

The Kaplan-Meier Estimate of Time to First Birth of Married Women.
To depict the survival of Ethiopian women's time to first birth at varying degrees of the covariate, a nonparametric survival analysis is necessary. It also explains the structure of the survival and hazard functions in the first birth interval dataset. A visualization of the KM curves for survival and hazard experience from time to first birth is shown in Figure 3. At first, the survival plot drops at a rising rate and then decreases at a decreasing rate. This indicates that the majority of women gave birth to their first child soon after marriage. Ethiopian women have a median survival period of 24 years after marriage for their first child (95% CI; 23.4, 25.3).

Comparison of the Different Covariates in terms of Survival Time.
To compare the differences between each categorical variable, a formal test was performed using the log-rank test. According to the results of the log-rank test (Table 3), there is a statistically significant difference in the survival experience between the levels of each covariate.

Parametric Shared Frailty Model
Results. As a baseline distribution, we used several parametric models with and without gamma shared frailty distributions to fit the data. The model with the lowest AIC and BIC is the best. The results showed that the Weibull-gamma shared frailty model was the best fit for the data (Table 4).

Bayesian Weibull-Gamma Shared Frailty Model Results.
The Gibbs sampler algorithm was used to construct parameter inferences, and it was implemented using 20,000 iterations in two separate chains; then, 1,000 terms were removed due to the burn-in state to avoid autocorrelation, and 42,000 samples were collected for the whole posterior distribution. We used a baseline distribution of the Weibull parametric model with and without gamma shared frailty distribution to fit the data. The model with the lowest DIC value is the most suitable. The results revealed that the best model to represent the data was Weibull-gamma shared frailty (Table 5). In this model, frailty is assumed to follow a gamma distribution with a mean of one and a variance of σ 2 . Table 5 shows the results of the Weibull-gamma shared frailty model, which reveals that σ 2 = 0:95 suggests regional heterogeneity. Table 6 shows that all parameter estimates have a Monte Carlo error (MC error) value of less than 5% of the standard deviation, indicating that the parameters are converging. The researcher uses this posterior summary as the final results because of this rationale and convergence graphs. The hazard ratio and 95 percent credible interval of Bayesian technique calculated values were used to analyze the final model findings. At the 5% level of significance, the confidence intervals of the mean for covariates that do not include 0 are significant.
Women who have never used contraceptives have e 0:1134 = 1:12 times increased hazard of first birth (95% CI: 0.06222, 0.1643) than women who have ever used any contraceptive method. Similarly, a researcher can say that the credible interval for the Bayesian hazard ratio did not include the one by exponentiation of the Bayesian credible interval (HR = e 0:1134 = 1:12, 95% CI: e 0:06222 , e 0:1643 ). The . And those who have married greater than 34 years have e −2:864 = 0:0561 times decreased hazard of first birth at an earlier age than individuals who married at age less than18 years keeping other variables constant. The husband's education level also affects the timing of the first birth. The hazard of the first birth at an earlier age among married women whose husbands have a primary level of education was increased by 11% (HR = e 0:09979 = 1:11, 95% CI: 1.047, 1.17), and that among married women whose husbands have secondary or above education level was increased by 73% (HR = e 0:1148 = 1:123, 95% CI: 1.04, 1.21) when compared to uneducated husbands keeping other covariates constant. The hazard of the first birth among educated married women was increased by 20% (HR = e 0:1826 = 1:20, 95% CI: 1.15, 1.27) when compared to uneducated women keeping other variables constant.
The estimated coefficient of the parameter for married women who get access to the media was -0.1168. The sign of the coefficient was negative, implying a decreased hazard rate for time to the first birth than those who did not (reference group). The other important variable in this study is the sex of the head of the household. That is, the risk of the first birth among married women whose household head is 3.6. Assessment of Convergence 3.6.1. History Plot. The posterior history plots (Figure 4) showed that for 20,000 iterations for each covariate, the model is converging, because the history plot seems tight and able to respond to all parameters.
3.6.2. Density Plot. Figure 5 shows the density plots associated with the coefficient of some selected covariates. Estimates for all parameters revealed good results because the density plot tends to smooth shape.
3.6.3. The Brooks-Gelman-Rubin (BGR) Statistics. The BGR convergence diagnostic graphs in Figure 6 show the line converted into one for stability indicating the convergence of the algorithm.

Discussion
The major objective of this study was to use Bayesian parametric shared frailty models to find determinant factors for married women's time to first birth in Ethiopia. Ethiopian women had a median survival time of 24 years after marriage for their first child (95 percent CI; 23.4, 25.3). The AIC and BIC criteria were used to compare the model distributions, with the model with the lowest AIC and BIC being accepted [37]. MCMC iteration is used to start the Bayesian approach parametric survival analysis until all parameters have converged. Because there is no known way for choosing an acceptable number of iterations and burn-in size, the MCMC iterations were generated by setting the initial values and burn-in state with no criterion. Rather, the researcher uses a trial-and-error method with the ultimate goal of obtaining stable parameter estimates with the lowest possible simulation error [38]. MCMC simulation improved the accuracy of the results by reducing the credibility interval and lowering the standard error but had no effect on the direction of the results [39,40]. The Gibbs sampler algo-rithm of the MCMC iteration method was used to create 42,000 samples in this investigation. 1,000 samples were utilized for the burn-in state, and 20,000 samples with two chains were used for posterior inference utilizing the Win-BUGS program for iteration and parameter convergence checks. Weibull distributions with and without the frailty term were compared using the DIC value after the parameters had converged, and the distribution with the smaller DIC value was preferred, as stated by Spiegelhalter [34].
In this study, the Bayesian Weibull-gamma shared frailty model revealed that residence, media exposure, women's education level, husband's education level, contraceptive use, and sex of the household head are statistically significant for married women's survival time to first birth in Ethiopia. This is in line with the findings of a study published in [41,42].
One of the factors of time to first birth, according to our findings, was one's place of residency. Married women in rural areas have a worse chance of surviving than married women in cities. This is in line with studies [1,20]. One explanation could be that metropolitan women are more likely to be educated and well informed about contraception use and the consequences of early childbearing [43,44]. Women who used the contraceptive had a long time to first birth than the nonusers [42,45,46]. This is due to the contraceptive service, which helped them in preventing early and unplanned pregnancy throughout their marriage. Age at the first childbirth was also linked to age at marriage. Females who married young had their first child sooner than women who married later. This conclusion is supported by research conducted both domestically and internationally [41,47,48]. The reason for this could be that older women need to have a baby shortly after marriage to have the necessary number of children before their reproductive lives come to an end. A woman who marries young can use a contraceptive to delay her first child until she is physically and psychologically mature.
A higher level of education shortens the timing of the first birth after marriage, which is consistent with the findings of [20,[49][50][51]. Women with a greater level of education had a shorter first birth interval. The link between schooling and the length of the first hiatus from childbearing appears to be indirect. In Ethiopia, increasing women's educational levels is linked to increased labor force participation, media access, and social standing [52]. High levels of education limit parents' traditional roles in deciding on their daughters' marriage partners and encourage self-selection of spouses, which sometimes takes longer. Highly empowered women are, indeed, more likely to have control over not only who and when they marry but also when they have children. This social empowerment of women and the lengthy process of spousal selection are believed to strengthen the intimacy of marriage partners and thus allow them to build strong confidence in each other, which, in turn, increases their desire to have a child to maximize marital satisfaction, resulting in shorter first birth intervals [53]. Furthermore, education improves marital stability by providing stable financial resources. This is also thought to minimize the first birth interval, as they are emotionally ready, biologically 10 Computational and Mathematical Methods in Medicine mature, and financially comfortable to have a child when they enter married life. In this study, media exposure was found to be inversely related. This is in line with previous research [16,54,55]. That is, women who have access to the media had a longer survival period for their first child after marriage. Different advertisements and education about the dangers of early marriage and early childbirth can reduce early marriage and sexual experience and improve understanding of reproductive health issues, which could explain the inverse relationship between media exposure and motherhood at a young age [56]. Women who do not have access to the media, on the other hand, lack adequate awareness of the high-risk period of becoming pregnant and are unaware of family planning options and the costs of early childbearing on the health of mothers and children [57].

Strengths and Limitations of This Study.
Although the EDHS is a sizable, nationally representative dataset, this research may have some drawbacks. First, recall bias might be present since survey participants were asked to recall things that happened five years ago, and they might have forgotten some specifics. A second drawback of the data is its cross-sectional design, which makes it challenging to establish causal links between exposure and outcome variables. Third, this study relies on self-reported data from life histories obtained from a nationally representative survey, which is prone to a number of sources of error. As a result,

12
Computational and Mathematical Methods in Medicine estimates for particular countries may be impacted by these limitations over time, and these should be treated with caution.
While acknowledging these limitations, the EDHS is a nationally representative dataset and has been rigorously designed and deployed by the Centers for Disease Control 13 Computational and Mathematical Methods in Medicine and Prevention using a global framework, and therefore, the findings can easily be generalized throughout the country. International comparisons of the findings will also be possible because DHSs adopt similar instruments across countries. We anticipate that the findings of this study will have strong policy implications for Ethiopia at the national level, as this is a study that has identified determinant factors of time from birth to first birth.

Conclusion
The Bayesian method parametric survival model was used to show the determinants of time to first birth among married women in Ethiopia. Based on the DIC value, the Weibull baseline distribution with the gamma shared frailty model was chosen as the best model and the Weibull gamma shared frailty distribution was chosen as the final model to   14 Computational and Mathematical Methods in Medicine suit the time to first birth dataset. As a result, the results of the Bayesian approach and Weibull AFT model analysis revealed that residence, media exposure, women's education level, husband's education level, contraceptive use, and household head's sex are statistically significant for married women's survival time to first birth in Ethiopia. To curb Ethiopia's rapid population growth, it is necessary to teach families how to use contraception and to raise awareness among rural populations about the importance of increasing the length of the first birth gap and not encouraging early marriage. In general, attempts to reduce fertility by extending the age at first marriage must take into account the socioeconomic and cultural contexts in which marriage takes place. On the other hand, the fight against early marriage should not just focus on reducing reproduction but also on the sociocultural, physiological, and psychological consequences. The impacts of early marriage on the contribution of late marriage to completed family size, on the other hand, need to be investigated further.
Abbreviations DIC: Deviance information criterion AIC: Akaike information criterion BIC: Bayesian information criterion CSA: Central Statistical Agency EDHS: Ethiopia Demographic and Health Survey HR: Hazard ratio WHO: World Health Organization.

Data Availability
The dataset and supporting materials will be obtained from the corresponding author on a reasonable request.

Consent
Consent is not applicable to this study.

Conflicts of Interest
There are no competing interests declared by any of the authors.

Authors' Contributions
KDF was involved in the study design, performed data extraction, and analyzed and drafted the manuscript. SMF, HBB, SBA, and MWM were involved in the study design and reviewed the manuscript. All authors have read critically and approved the final manuscript.