Bayesian Parametric Modeling of Time to Tuberculosis Co-infection of HIV/AIDS Patients Under Antiretroviral Therapy Treatment at Jimma University Medical Center, Ethiopia

Background: Tuberculosis is the most common opportunistic infection among HIV/AIDS patients, including those following Antiretroviral Therapy treatment. The risk of Tuberculosis infection is higher in people living with HIV/AIDS than in people who are free from HIV/AIDS. Many studies focused on prevalence and determinants of Tuberculosis in HIV/AIDS patients without taking into account the censoring aspects of the time to event data. Alternatively, this study used INLA methodology to build Bayesian parametric survival models for modeling time to Tuberculosis co-infection of HIV/AIDS patients following Antiretroviral Therapy treatment. Methods: A data of a retrospective cohort of HIV/AIDS patients under Antiretroviral Therapy treatment follow-up from January 2016 to December 2020 until Tuberculosis was clinically diagnosed or until the end of the study was collected from Antiretroviral Therapy treatment center of Jimma University Medical Center, Ethiopia. In order to identify the risk factors which have association with Tuberculosis free survival time, we ﬁtted Exponential, Weibull, Log-logistic, Gamma and Log-normal models to our data set.. Results: About 26.37% of the study subjects had been co-infected with tuberculosis during the study period. Among the parametric Accelerated failure time models, the Bayesian log-logistic Accelerated failure time model was found to be the best ﬁtting model for the data. Conclusions: Tuberculosis co-infection survival time was signiﬁcantly associated with place of residence, smoking, drinking alcohol, family size, WHO clinical stages, functional status, CD4 count, BMI and hemoglobin level. The ﬁnding of this study provide timely information on the risk factors associated with TB co-infection survival time for healthy policy makers and planners.


Background
Tuberculosis is one of the infectious diseases that affects the lungs and other sites [1]. Tuberculosis has been the main public health problem affecting millions worldwide and it remains the top infectious killer in the world causing close to about 4000 lives a day [2]. Around 10.0 million people estimated to have developed TB disease in 2019 worldwide, and there were around 1.2 million TB deaths among HIV-negative people and an additional 208, 000 deaths among people living with HIV [3].
Tuberculosis is the most common opportunistic infection among HIV positive people, including those following ART treatment, and it is the major cause of HIVrelated death [4]. UNAIDS report of 2018 showed that Sub-Saharan Africa is the hardest hit region of the world, as it has around 70% of all people living with HIV/AIDS and TB co-infection in the world [5].
Ethiopia is one of the highest 30 TB and HIV co-infection countries in the world, with TB incidence rate of about 164 per 100,000 TB, 112/100,00 TB among people living with HIV [3]. Research shows that the TB co-infection in Ethiopia is high in areas where HIV is highly distributed [6]. The majority of the People living with HIV were from the Amhara (30%), Oromia (26%) and Addis Ababa (18%) regions [7].
The HIV virus infects CD4 cells causing reduction of the number of immune cells which causes the body fail to control viral multiplication which increases the chance of opportunistic infection with tuberculosis being the most common opportunistic infection at HIV diagnosis [8]. The treatment outcome of HIV-positive following ART treatment has remarkably changed with a a reduction of plasma viral copies and an increase of CD4 counts [9]. It had been reported that the ART treatment has reduced the incidence of TB in HIV patients by about 70 -90% [10]. Even with the advantages of ART treatment, still HIV/AIDS patients on ART treatment develop TB with about prevalence rate of 2.5-30.1 [8].
Though, Tuberculosis can affect everyone, the risk of Tuberculosis infection is higher in people living with HIV than in people who are free from HIV [11]. Studies revealed that certain HIV-infected people develop TB, while others do not. The risk of developing active TB in people living with HIV is about 21 times higher than those free from HIV which may be due to the altered immune state of HIV positive patients [18]. There are various factors that increase the chance of TB infection among HIV/AIDS patients including CD4 cell count and the number of vi-ral loads [19,20], household family size, cigarette smoking, baseline CD4 cell counts, WHO clinical stages, having a history of diabetics [58], ans etc. However, these factors have not been studied in the context of Survival Analysis, where association between risk factors and time to TB co infection might be of interest.
Majority of the study focused on prevalence and predictors of TB in HIV patients. In order to determine the important determinants of TB co-infection in HIV patients, most of the methodologies in the literature used logistic regression with outcome being the TB's viability through follow up time of HIV/AIDS patients taking ART treatment [12,14,17]. In logistic regression, our interest is to study how risk factors were associated with the presence or absence of a disease( or an event) without taking into account the effect of time [15]. These approaches help to provide odds ratios for significant variables associated with the risk of TB infection but rejects the censoring aspects of time-to-event data.
The aim of this study was to model the predictors of time to TB co-infection in HIV/AIDS patients following ART treatment using Bayesian parametric survival analysis approach based on INLA methodology. There have been advances in computational and modeling techniques using Bayesian approach of survival data [23,44]. Due to the complex likelihood functions to accommodate censoring, survival models are generally very difficult to fit. Bayesian approach to survival analysis may overcome this by using the MCMC techniques and other numerical integration methods like INLA [34]. Additionally, the Bayesian approach can help to capture prior information and incorporate it seamlessly via a rigorous, probabilistic framework [35]. Therefore, Bayesian Survival Analysis approach is preferred over the usual frequentist technique as the power of information obtained from the approach is much better as it is the combination of likelihood of the data and prior information about the distribution of the parameter [37].
Due to the fact that Bayesian inference cannot be performed analytically, MCMC is particularly useful in Bayesian inference because of the focus on posterior distributions which are often difficult to work analytically [41].
The popular software packages used with MCMC method include BUGS, JAGS and STAN [48]. The first problem with these packages is that they are slow and some times need to wait days to get the result. The second difficulty with MCMC methods is that it requires high skill of programming language to specify the likelihood functions of different models. The integrated nested Laplace approximation method for approximate Bayesian inference was developed by Rue, Martino, and Chopin as an alternative to the MCMC method [23]. INLA is an alternative method for Bayesian inference on latent Gaussian models when the focus is on posterior marginal distributions. It substitutes MCMC methods with accurate, deterministic approximations to posterior marginal distributions. Integrated nested Laplace approximation provides a fast and exact approach to fitting latent Gaussian models which include many statistical models, including survival models [48].
Survival models can be written into a latent Gaussian model which allows us to perform Bayesian inference using integrated nested Laplace approximations [44]. Survival analysis consists of a great body of work using latent Gaussian models and it is one of the statistical models on which INLA has been successfully applied [27,28]. The main advantage of INLA over MCMC techniques is its simplicity of computation [28]. Using INLA results are generated within seconds and minutes even for models with a large dimensional latent field, where as MCMC algorithm would take hours or even days. The other advantage of INLA is that INLA treats latent Gaussian models in a unified way, thus allowing greater automation of the inference process.
Even though Bayesian approaches to the analysis of survival data can provide a number of benefits, they are less widely used than the classical approaches [29]. Therefore, the motivation to apply Bayesian Survival Analysis for this study stems from the above mentioned advantages of Bayesian survival analysis approach over the classical survival analysis approach.

Data source
The nature of the data set used for this study was survival data. In the data set, time until TB infection was clinically diagnosed in HIV/AIDS patients was investigated.

Inclusion and exclusion criteria
All adult HIV patients following ART treatment and who were 18 years old and above during the study period and TB free at the inception of the study with at least two follow up period were included in the study. Patients whose date of TB co-infection was unknown were excluded from the study. Also, Patients with insufficient information about one of the variables in the study were not included.

Study design, population and sample size
A data of a retrospective cohort of adult HIV/AIDS patients under ART follow-up from January 2016 to December 2020 until TB was clinically diagnosed or until the end of the study was collected from ART clinic of JUMC, Ethiopia. In this study, the source population was all adult HIV/AIDS patients who were 18 years old and above at Jimma University Medical Center. There were a total number of 3069 HIV/AIDS pateints. Among the total patients 421 of the patients were included in the study based on the inclusion and exclusion criteria.

Study variables
The response variable for this study was time to TB infection in HIV/AIDS patients under ART treatment. Time is measured in months.The predictor variables which were assumed to have effect on time to TB infection in HIV/AIDS patients were Age, sex, place of residence, family size, alcohol usage status, smoking status, marital status, education level, WHO clinical stages, functional status, CD4 count, body mass index and hemoglobin level.

Survival data analysis
Survival analysis is a collection of statistical techniques for data analysis for which the outcome variable of interest is the time until an event occurs. By time, we mean years, months, weeks, or days from the beginning of follow-up of an individual until an event occurs. By event, we mean death, disease incidence, relapse from remission, recovery from disease or any designated experience of interest that may happen to an individual. Censoring is one of the common features that makes survival analysis unique from another statistical analysis. Censoring is present when we have some information about a subject's event time, but we don't know the exact event time [30]. The general reasons why censoring might occur are: A subject does not experience the event before the study ends, the patient is lost to follow-up during the study period, or the patient withdraws from the study.

Survival function
Assume that the survival time, T, is a continuous random variable. The distribution of T can be described by the usual cumulative distribution function which is the probability that a subject from the population will die (or a specific event of interest for a subject has occurred) before time t [31]. The corresponding density function of T is In survival analysis, it is common to use the survival function The relationship between f(t) and S(t) is given as follows;

Hazard function
It is also of interest, in analyzing survival data, to assess which periods having high or low chances of the event among those still active at the certain time. A suitable method to characterize such risks is the hazard function., h(t), defined by the following equation [31].
It is the instantaneous rate of failure (experiencing the event) at the time t given that a subject is alive at the time t.
The definition of the hazard function implies that A related quantity is the cumulative hazard function, H(t), defined by And thus, The Kaplan-Meier estimator of survival function The Kaplan-Meier (KM) method is a non-parametric(distributionfree) survival analysis method which is used to make comparisons of the survival rates between two or more groups [32]. To estimate the survivor function, S(t), without covariates, we can use the Kaplan Meier estimator. Let there be n individuals with observed survival times t 1 , ...,t n and r be death times amongst the individuals, where r ≥ n, j = 1,...,r. The r ordered death times are t (1) < t (2) < ... < t (r) . Let n j denotes the number of individual who are alive just before time t ( j) , including those who are about to die at this time, and let d j denotes the number who die at this time. The Kaplan-Meier estimator of the survival function at any time in the k th time interval from t (k) to t (k+1) , k = 1,..., r is given by [36].

Parametric Survival Models
In a parametric survival models, survival time is assumed to follow a known distribution [33]. Parametric models play an important role in Bayesian survival analysis, since many Bayesian analyses in practice are carried out using parametric models and parametric modeling offers straightforward modeling and analysis techniques [34].

Parametric proportional hazard models
Let h(t/x)) be the hazard function at time t for a subject given the covariate vector x = (x 1 , ..., x p ) T . The basic model proposed by Cox is as follows: where h 0 (t) is the baseline hazard function and β i 's are the unknown regression parameters to be estimated. In Parametric proportional hazard model, the baseline hazard function h 0(t) is assumed to follow a specific distribution when a fully parametric PH model is fitted to the data. The hazard ratio is hence given by HR = exp (β 1 X 1 + β 2 X 2 + ... + β p X p ).

Accelarated failre time models
Although parametric PH models are very useful to analyze survival data, there are relatively few probability distributions for the survival time that can be used with these models [39]. In these situations, the accelerated failure time model (AFT) is an alternative to the PH model for the analysis of survival time data. Under AFT models we measure the direct effect of the explanatory variables on the survival time instead of hazard, as we do in the PH model. This characteristic allows for an easier interpretation of the results because the parameters measure the effect of the covariates on the survival time.
In accelerated failure time (AFT) models, the natural logarithm of the survival time, log(t), is expressed as a linear function of the covariates, which yields therefore a linear model: We interpret the effect of the AFT model as the change in the time scale by a factor of exp(xjβ ). Based on whether this factor is greater or less than 1, survival time is interpreted to either accelerate or decelerate. Accelerated failure time does not imply a positive acceleration of time with the increase of a covariates but rather a deceleration, or, in other words, an increase in the expected waiting time until failure. AFT models have the opposite sign from similar estimates in proportional hazard models, due to the fact that the PH models predict the hazard and the AFT model predicts time.
An advantage of the AFT approach is that the effect of the covariates is described in absolute terms (i.e. number of months or years) instead of in relative terms (i.e. a hazard ratio). The acceleration factor is the central measure of association obtained in AFT models and allows you to evaluate the effect of covariates on the survival time.

Bayesian Modeling Approach for survival data
The Bayesian paradigm is based on specifying a probability model for the observed data D, given a vector of unknown parameters θ , leading to the likelihood function L(θ /D). Then we assume that θ is random and has a prior distribution denoted by π(θ ). Inference concerning θ is then based on the posterior distribution [34] , which is obtained by Bayes' theorem. The posterior distribution of θ is given by: where Θ denotes the parameter space of θ The quantity m(D) = Θ L(θ /D)π(θ )dθ is the normalizing constant of π(θ /D) , and is often called the marginal distribution of the data or the prior predictive distribution. In most models and applications, m(D) does not have an analytic closed form, and therefore π(θ /D) does not have a closed form. The Bayesian survival analysis approach considers the parameters of the model as random variables and requires that prior distributions specified for them and data are considered as fixed [40].

Likelihood Function in Bayesian
Survival Analysis Suppose we observe n independent vectors of (T i , δ i ), where T i is time to the event and δ i is indicator variable telling us whether T i is censored or not, i.e, T i = 0 for censored observation(δ i = 0) and T i = 1 for uncensored observation(δ i = 1). The likelihood function of the set of unknown parameters θ in the presence of right censoring is given as The Integrated Nested Laplace Approximation Methodology for

Bayesian Inference
For long time, Bayesian statistical inference has relied on MCMC methods to compute the joint posterior distribution of the model parameters which is usually computationally very expensive [41]. An alternative approach and fast estimation methods to MCMC which allows user to easily perform approximate Bayesian inference using Integrated Nested Laplace Approximation was proposed by Havard Rue, Martino, and Chopin [44]. INLA computes posterior marginals for each component in model, from which posterior expectation and standard deviations can easily be found.
Latent Gaussian models are subset of all Bayesian additive models with a structured additive predictor say η i . In these models, the observed variable y i is assumed to belong to an exponential family, where the mean µ i is linked to this structured additive predictor η i through a link function g(.), so that g(µ i ) = η i . The structured additive predictor η i accounts for effects of various covariates in an additive way: Here, the f j (.) are unknown functions of covariates u, the {β k } represent the linear effect of covariates z and the ε i 's are unstructured terms. A Gaussion prior is assigned to α, f j (.) , {β k } and ε i . We denote π(./.) as the conditional density of its arguments, and let x denote the vector of all n Gaussian variables η i , α, f j (.) and {β k }, and θ denotes the vector of hyper-parameters, which are not necessarily Gaussian. The density π(x/θ 1 ) is Gaussian with(assumed) zero mean and precision matrix Q(θ 1 ) with hyperparameter θ 1 .
The distribution for the n d observational variables y = {y i : iεI} is denoted by π(y/x, θ 2 ) and we assume that {y i : iεI} are conditionally independent given x and θ 2 . For simplicity, we denote by The posterior then reads( for non singular Q(θ )), The imposed linear constraints(if any) are denoted by Ax = e for a kxk matrix A of rank k. The main aim is to approximate the posterior marginals of the latent field, π(xi/y) and the posterior marginals of the hyperparameters , π(θ /y) and π(θ j /y). We can write the posterior marginal of interest as The importance of INLA is to use the above form to construct nested approximations, as this approach makes Laplace approximations very accurate when applied to latent Gaussian models.
Here,π(./.) is an approximated( conditional) density of its arguments. Approximations to π(x i /y) are computed by π(θ /y) and π(x i /θ , y) and using numerical integration to integrate out θ . The approximation of π(θ j /y) is computed by integrating out θ − j fromπ(θ /y). The posterior marginal π(θ ) of the hyperparameters θ is approximated using a

Prior Distributions in INLA
Bayesian statistical inference depends on the posterior distribution which is obtained by updating the prior beliefs by new evidence. Prior distribution can be broadly classified into non-informative, weakly informative and informative prior distributions. Non-informative prior distributions, also known as objective prior distributions, are designed to have minimal impact on the posterior distribution so that the data alone can be the source of inference [45]. The non-informative prior distribution often produce the same results as maximum likelihood estimates. On the other hand, the informative prior distributions that aim to construct a prior distribution that reflect the current knowledge on the values of the parameters and the uncertainties that surround the knowledge about the parameters in question [46]. In INLA, it is assumed that fixed effects follow Gaussian distribution with mean zero and small number of precision matrix Q(θ 1 ) and only the parameters in the precision matrix of the random effect need a prior which was considered as a hyper-parameter [47]. For this study, Gaussian prior distribution (non-informative) with mean zero and variance equal to 1000(precision equal to 0.001) was used for the fixed effects and the intercept [23]. And for hyper-parameters a non-informative prior of Gamma distribution prior is a common non -informative prior to be assigned [48]. In INLA, the Latent component of the model, η i = β 0 + β 1 z 1 + ... + β p z p must follow a Gaussion distribution [48]. In this study, it was assumed that fixed effect(coefficients) associated with covariates have a Normal distribution with mean 0 and variance 10 2 , i.e, β p , p = 0,.....,, i.e, β p ∼ N(0, 10 2 ) [34]. Then for this study, to complete the model we have assigned a non-informative Gamma prior for for the hyperparameter of the model τ i ∼ Γ(a, b) and α ∼ Γ(a, b) with a =1 and b = 0.001 which is similar with prior distribution used by many researchers worked on Bayesian survival analysis [23,34,44,49].

Bayesian parametric survival models Exponential Model
The exponential model is the most fundamental parametric model in survival analysis [34]. Suppose we have independent identically distributed (i.i.d.) survival times t = (t 1 ,t 2 , ...,t n )', each having an exponential distribution with parameter λ . Denote the censoring indicators by δ We build a regression model by introducing covariates through λ , and write x i 'is a p x 1 vector of covariates, β is a p x 1 vector of regression coefficients, ϕ(.) and is a known functionand D = (n,t,X, δ ) denotes the observed data for regression model Using these, we get the likelihood function [44].
Suppose we specify a normal prior for β with mean µ 0 and variance σ 2 0 . Then the posterior distribution of β is given by where π(β /µ 0 , σ 0 ) is the normal density with mean µ 0 and variance σ 2 0 .

Weibull model
The Weibull model is perhaps the most widely used parametric survival model [34]. Suppose we have independent identically distributed survival times t = (t 1 ,t 2 , ...,t n )', each having a Weibull distribution, denoted by ω(α, γ). It is often more convenient to write the model in terms of the parameterization λ = log(γ), leading to We can write the likelihood function of (α, λ ) as To build the Weibull regression model, we introduce covariates through λ and write λ i = x ′ i β . Where x i is a px1 vector of covariates, β is a px1 vector of regression coefficients. Assuming a normal prior with parameters µ 0 , σ 2 0 for λ and gamma prior with parameters (α 0 , κ 0 ), the joint posterior distribution of (α, λ ) is given by Where D = (n,t, x, δ ) denote the observed data for regression model.

Log-logistic model
The log-logistic model possesses a rather supple functional form [50]. The Log-logistic distribution is among the parametric survival models where the hazard rate initially increases and then decreases. If we have independent identically distributed survival times t =(t 1 ,t 2 , ...,t n )', each having an log-logistic distribution, denoted by T ∼ LL(α, λ ), with density for α > 0, λ > 0 and t ≥ 0. And, the survival function is given by for t > 0. We can write the likelihood function of (α, λ ) as To build the regression model, we introduce covariates through λ , and write λ i = x ′ i β . Where x ′ i is px1 vector of covariates, and β is px1 regression coefficients. If we assume gamma prior with parameters (α 0 , κ 0 ) for α, we will have the following joint posterior

Log-Normal Model
Another commonly used parametric survival model is the log-normal model [34]. For this model, we assume that the logarithms of the survival times are normally distributed.
If t i has a log-normal distribution with parameters (µ, σ ) , denoted by ιN(µ, σ ), then The survival function is given by We can thus write the likelihood function of (µ, σ ) as Then, To construct the regression model, we introduce covariates through µ, and write µ i = x ′ i β . Assuming β /τ ∼ N p (µ 0 , τ −1 ς 0 ), the joint posterior for β , τ is given by

Gamma model
The gamma model is a generalization of the exponential model [34]. For this model, t i ∼ ζ (α, λ ) and its density function is given by: The survival function is given by We can thus write the likelihood function of (α, λ ) as To construct the regression model, we introduce covariates through λ , and write λ i = x ′ i β . Assuming β ∼ N(µ 0 , σ 2 0 ), we are lead to the joint posterior

Model Comparison Methods
Integrated Nested Laplace Approximation computes a number of Bayesian criteria for model assessment and model selection [51]. Model selection criteria will be of help when selecting among different models. The following methods of model selection techniques was used in this study.

Marginal likelihood
The marginal likelihood of a model is the probability of the observed data under a given model [42]. The marginal likelihood approximation provided by INLA is computed as.π (y) = π(θ , x, y) π G (x/θ , y) | x=x * (θ ) dθ

Information-based criteria (DIC and WAIC)
The deviance information criterion (DIC) is a popular criterion for model choice [52]. It takes into account goodnessof-fit and a penalty term that is based on the complexity of the model via the estimated effective number of parameters. The DIC is defined as where, D(.) is the deviance,x andθ the posterieor expectations of the latent effects and hyperparameters, respectively, and pD is the effective number of parameters. The effective number of parameters pD can be computed as pD = E[D(.)] -D(x,θ ) The Watanabe-Akaike information criterion, also known as widely applicable Bayesian information criterion, is similar to the DIC but the effective number of parameters is computed in a different way. The final formula to calculate WAIC is.
Where, ∑ n i=1 logp post(y i ) is the sum of predictive density for each data point and pD is the effective number of parameters.

Model Diagnosis
Diagnosis for the accuracy of INLA approximation for the models The Kullback-Leibler divergence (kld): This value describes the difference between the normal approximation and the simplified Laplace approximation. Small values indicate that the posterior distribution is well-approximated by a normal. Effective number of parameters(pD): The posterior summary results from INLA also contain, effective number of parameters which is another measure of the accuracy of approximation. In particular, if the effective number of parameters is low compared to the sample size, then one expects the approximation to be good. A two types of "Goodness of fit" reported by INLA are: Conditional predictive ordinates (CPO): Conditional predictive ordinates are a cross-validatory criterion for model assessment [53]. It is computed for each observation as

Unusually small or large values of CPO i indicate a surprising observation. Predictive integral transform (PIT):
The predictive integral transform (PIT) measures the probability of a new value to be lower than the actual observed value for each observation [54]. It is computed as An unusual large or small value indicates possible outliers.
Due to howπ(x i /y −1 , θ ) are computed there may be cases where this computation fails due to inaccurate tail behavior ofπ(x i /y i , θ j ). To monitor the reliability of the CPO and PIT values computed, failure variable computed for each i (or y i ) is defined as follows.
• Ifπ(x i /y i , θ j ) is monotone increasing or decreasing, then failure is set to 1 and thenπ(x i /y i , θ j ) is set to the 0-function. In this case,π(x i /y i , θ j ) is known to be just wrong. • Ifπ(x i /y i , θ j ) is has a (local) maximum either at minx i or at maxxi, thenπ(x i /y i , θ j ) is set to zero in that part whereπ(x i /y i , θ j ) is decreasing (starting from min(x i ) or increasing (starting from maxx i . When the expected failure is 0 then the computed value of CPO and PIT seems to be reliable, and when the expected failure is 1 then the computed value of CPO and PIT is known to be completely unreliable.

Summary of descriptive statistics results
A total of 421 adult HIV patients following ART treatment at Jimma University Medical Center, Ethiopia were included in the analysis. During the follow-up period, 111(26.37%) of the study subjects had experienced the event(had been co-infected with TB).
Descriptive results of demographic and clinical characteristics of patients were presented in Table 1 and

Kaplan-Meier estimate of survival functions
From the plot of the overall Kaplan-Meier survival curve given in the figure 1 below, it can be seen that, a large number of TB co-infection recorded at the earlier time of the initiation of ART treatment and there is a decreasing pattern of TB co-infection through the follow up period. In order to explore differences between TB co-infection free survival time between or among groups, separate KM survival function curves are constructed for categorical covariates and results are given in figure 2,3,4 and 5. In general, if the pattern of one survivor-ship function is above the other,

Results of Bayesian Log-logistic AFT model
The table 4 and Table 5 below shows the posterior summary results of Bayesian Log-logistic AFT model. The decision about the significance of the variables is based on the 95% Credible Interval for the posterior mean of the coefficients.
In our study, it appeared that residence, smoking status, alcohol consumption status, WHO clinical stages, functional status, family size, CD4 count, Body Mass Index and hemoglobin level of the patients were significant risk factors associated with TB co-infection free survival time of HIV patients following ART treatment. The interpretation of the estimated posterior mean of parameters of the model was done using estimated acceleration factor(γ = exp(β j ). In order to decide the significance of the covariates in the model, the 95% credible interval was used. The factors whose credible intervals for posterior mean of parameters contained 0, or whose credible intervals for acceleration factor contained 1, implied that these factors were not sig-nificant. The final model can be written as: Where, T represents time to TB co-infection for each subject. I is an indicator variable for categories of variables where I (.=1) is considered as a reference category.
In Log logistic AFT model, the positive estimated posterior β coefficients indicate a longer TB co-infection free survival time, where as the negative estimated posterior β coefficients indicate shorter TB co-infection free survival time for the patients.
The estimated acceleration factor for patients who reside in urban wasγ = exp(-1.258) = 0.2842 with 95% CI of [0.1571, 0.5035] which does not include one. This means that, keeping all other factors constant the expected TB coinfection free survival time of patients who reside in urban area decreases by a factor 0.2842 as compared to patients residing in rural area.
The second significant predictor of TB co-infection free survival time of HIV/AIDS patient was smoking status, in which acceleration factor of smoker patients was estimated to 0.5137 with 95% CI of [0.2814, 0.9361]. Thus, keeping all other factors constant, the expected TB co-infection free survival time of smoker patients decreases by a factor of 0.5137 as compared to non-smoker patients.
If we look at the alcohol usage status of the patients, the estimated acceleration factor for posterior mean of its coefficient was estimated to be 0.4151 with 95% CI of [0.2324, 0.7408] which can be interpreted as, keeping all other factors constant, the expected TB co-infection free survival time of alcoholics patients decreases by a factor of 0.4151 as compared to non-alcoholics patients.      It can also be seen that patients with low body mass index and patients with severe anemic status was found to be the significant risk factors for TB co-infection in HIV patients. The estimated acceleration factor for underweight patients was 0.3667 with 95% CI of [0.2141, 0.6955], and the estimated acceleration factor of severe anemic patients was estimated to be 0.3036 with CI of [0.1408, 0.6564]. Thus, keeping all other factors constant, the TB co-infection free survival time of underweight HIV patients decreases by a factor of 0.3667 as compared to patients with normal weight,and, the TB co-infection free survival time of severe anemic HIV patients decreases by a factor of 0.3036 as compared to patients with normal anemic status.

Model Diagnosis results
The Kullback-Leibler divergence (kld) From the tables 4 and Table 5, we can see that the values of kld is zero which means the marginal posterior densities of regression coefficients were well approximated by the Normal distribution. Effective number of parameters(pD) In this study, the ratio of sample size (421) and effective number of parameters (28.93) was found to be 14.55, suggesting a reasonably good approximation. The ratio can be interpreted as the number of equivalent replicates corresponding to the number of observations for each expected number of effective parameters.

Discussions
The finding of this study revealed that being urban resident, smoker, alcoholics, in WHO clinical stages III and IV, in bedridden functional status, in larger family size, in lower CD4 count, underweight and anemic patients shortens TB co-infection free survival time of HIV/AIDS patients under ART treatment at Jimma University Medical Center.
Among the study participants who fulfilled the inclusion criteria, about 26.37% of subjects had been co-infected with Tuberculosis. The proportion of of TB co-infection in this study cohort is smaller compared to other study settings in different parts of Ethiopia in which it was found to be 62.3% and 40.1% in a retrospective study conducted in seven ART clinics located at Addis Ababa and in Northeast Ethiopia with respectively. The proprtion of TB coinfection observed in this study is consistent with the findings from Noth-west Ethiopia(26.4%) [55], Amhara re-   [63]. However, The proportion of TB/HIV infection reported in this study is larger than the reported prevalence of TB among HIV patients in Albert Luthuli municipality of South Africa(18.2%) [56].
The findings of this study showed that patients' residence place were significantly associated with TB co-infection free survival time of HIV/AIDS patients. Patients who reside in urban areas are more likely to be infected with TB as compared to patients residing in rural areas. This result is consistent with the report of the retrospective study conducted by Beshir et al. at Adama Referral Hospital and Medical College, Oromia, Ethiopia [57]. They reported residence place as one of the significant risk factors of TB coinfection in HIV/AIDS patients. However, the finding of this study is inconsistent with those of other studies undertaken in Ethiopia [55,58,59].
The finding of this study have shown that smoking was an important predictor of TB co-infection free survival time of HIV/AIDS patients following ART treatment. Smoker patients are more likely to be infected with TB compared to non-smoker patients. This result agrees with the result reported by Anye et al. based on four year retrospective data of 1077 HIV patients in the Bameda regional hospital of Cameroon [60]. Our result also agrees with report of studies undertaken in Ethiopia [55,59,61]. Their results suggested that being smoker is significantly associated with TB co-infection free survival time in HIV/AIDS patients.
Our result also revealed that being alcoholics is associated with TB-co-infection free survival time of HIV/AIDS patients. Being alcoholics leads to shorter TB co-infection free survival time of HIV/AIDS patients attending ART treatment. Our findings showed alcoholic patients were more likely to co-infected with TB compared non-alcoholic patients. This finding is consistent with the reports of Ahmed et al. and Abdu et al. [55,61].
The result of our study suggested that baseline WHO clinical stages were one of the clinical factors associated with TB co-infection free survival time of HIV/AIDS patients. Patients with WHO clinical stages III and IV are more likely to be co-infected with TB than those with clinical stage I. This finding supports the findings of the study undertaken by Chang et al. [62]. Their study suggested that being in advanced WHO clinical stages is associated with higher risk of developing TB compared WHO clinical stages I and II [62]. Our finding also coincides with the study conducted in Amhara region of Ethiopia by Aweke et al. [63].
Similarly, in this study, patients' functional status at baseline was found to be the predictor of TB co-infection free survival time of HIV/AIDS patients. This result is consistent with the report of the study conducted by Aemro et al. at Debra Markos referral hospital, Northwest Ethiopia [64]. Accorging to our study, patient with bedridden functional status at baseline had shorter TB co-infection free survival time compared to patients with working functional status at baseline.

Conclusion
Bayesian survival analysis approach with INLA method was applied to fit the parametric survival models to our data set. Among the parametric AFT survival models, Bayesian Log-logistic AFT model was found to be the best fitting model for our data set.
The TB-free survival time of HIV/AIDS patients was significantly affected by residence place, smoking, drinking alcohol, family size, WHO clinical stages, functional status, CD4 count, body mass index and hemoglobin level of the patients. Thus, the finding of this study suggests that all stake holders are expected to work on prevention of TB by giving awareness on the risk factors of TB co-infection in HIV/AIDS patients.

Competing interests
The authors declare that they have no competing interests.

Funding
The study was funded by Jimma University.

Availability of data and materials
The data set used in this study can be obtained from Jimma University Medical Center.

Ethics approval and consent to participate
Letter of ethical clearance was obtained from Department of Statistics of Jimma University and submitted to Jimma University Medical Center to get permission to undertake the research. This study was developed in accordance with established legislation and complies with the norms of good clinical practice, and informed consent was being not necessary as personal identifying information was kept separate from the research data.

Consent for publication
Not applicable.