Study setting
The study was conducted in the South Gondar Zone, Amhara Region, in northwest Ethiopia. The South Gondar zone comprises town administrations and rural districts. There are over 500 public health facilities at various levels, including health posts, that are directly or indirectly involved in the TB control program. The routine practice of the existing TB control program in the South Gondar Zone is passive case detection. Direct Observation Therapy (DOT) is implemented in all public health facilities, including health posts. In addition to modern healthcare facilities, there are several informal healthcare providers, such as traditional healers and holy water centers. Therefore, we specifically selected this zone due to the widespread traditional healthcare practices.
Study design
A cluster randomized controlled trial was conducted from April 1 to January 30, 2024, across four districts and two town administrations (totaling six districts). These districts were randomly selected from a total of 13 districts and eight town administrations. The six districts were assigned to either the intervention or control groups in a 1:1 ratio. The randomization process was conducted by experts who were not involved in the study as researchers. The number and size of the clusters involved in the intervention and control groups were not considered. All health facilities (clusters) that offer TB diagnostic and DOT programs were included in the study. Out of all the health facilities in the study area, 23 facilities were randomly selected. Of these, thirteen health facilities were located at the intervention site, and ten health facilities were found at the control site.
Intervention group
The intervention group received integrated TB care, which integrated traditional and modern TB services. This approach included: 1) training for health professionals and traditional care providers; 2) TB screening at traditional healthcare sites; and 3) referral linkage from traditional to modern healthcare. The training was conducted in three rounds, aiming to increase knowledge, foster a favorable attitude, and enhance skills in TB screening and patient referral activities. The intervention was designed in four phases. In the initial phase, investigators prepared comprehensive training manuals that underwent review and standardization by invited experts, including physicians, public health experts, nurses, and language specialists. A workshop was conducted to enhance the training manuals, and experts provided valuable feedback on the content, depth, readability, and comprehensibility. The manuals covered various aspects of TB, including causes, symptoms, signs, transmission, screening and referral procedures, diagnostic approaches, case detection techniques, treatment outcomes, advantages of early detection, challenges of late diagnosis, and TB control and prevention mechanisms. The training also included models for integrating traditional and modern healthcare systems, which were approved by senior experts.
During the training provision phase, both traditional and modern care practitioners underwent training in three rounds. The first round focused on traditional practitioners (e.g., traditional healers, religious leaders) for five days, while healthcare providers working at DOT clinics and TB focal persons received training for two days. Subsequent one-day training sessions were conducted three and six months after the initial training. The training was facilitated by researchers and TB experts who held trainers' certificates. Traditional healers and religious leaders who demonstrated proficiency in knowledge, attitude, and skills through a post-test and practical assessment were included in the intervention to screen and refer presumptive cases. Details of the operational procedures for the intervention packages are provided in a published protocol(18) and are available in the supplementary materials (S1-Table 1).
The intervention was fully implemented when patients were screened and referred to the health facilities near them. Traditional healers and religious leaders used standardized screening tools to identify patients showing symptoms of TB, referring all suspected cases to nearby health facilities in the intervention districts. Trained TB focal persons re-screened and diagnosed patients according to national TB treatment guidelines. To ensure effective intervention implementation, regular monthly supervision was conducted by experts. In the final phase, the end-line outcome was assessed by comparing the difference between the end-line and baseline results in both the intervention and control groups. The details of screening and referral formats are provided in the supplementary materials (S2: screening form and S3: referral format).
Control group
In the control group, routine TB care continued without any additional intervention from the research team. Routine TB care in Ethiopia involves identifying individuals with TB when they visit healthcare services on their own due to symptoms or other health concerns. The control group served as a reference comparator for measuring the effectiveness of integrated interventions. Baseline information was collected simultaneously in both the control and intervention groups.
The outcome of interest and its definition
This study focused on the delay in diagnosing TB, which includes patient, health system, and treatment delays. Diagnosis delay is defined as the time interval between the onset of symptoms and the confirmation of TB in the patient. Patient delay refers to the period from the onset of the first symptom to the first medical consultation. Health system delay is also defined as the time from the first consultation to the date of TB diagnosis. Treatment delay refers to the time from diagnosis to the start of anti-TB medication. Detailed definitions of the outcome of interest and other variables are available in the supplementary materials (S4-Table 2). To measure the association between the dependent and independent variables, the median diagnosis time was used as the cutoff time. The time of measurement for this study was during the day. The time exceeding the median was considered a “delayed diagnosis,” while the time from the onset of symptoms to the median was considered "not delayed."
Participant recruitment and randomization
Two districts and one town administration were assigned to the intervention group, while the other two districts and one town administration were assigned to the control group. Twenty-three public health facilities located within the intervention and control districts were included in the study. Additionally, 29 traditional care centers located at the intervention site were selected for the implementation of the intervention. Using a random sampling approach, Dera woreda, Libokemekem woreda, and Worta town were chosen as intervention districts, while Farta, Gunabeyemeder, and Debre-Tabor town were selected as control districts. The districts and town administrations allocated to the intervention and control groups had similar baseline characteristics. In the study area, the level of training for providers, Directly Observed Treatment (DOT) programs, laboratory supplies, diagnostic techniques, and guidelines were consistent across all health facilities. Furthermore, study participants were selected randomly. A buffering zone between the intervention and control groups was implemented to minimize information contamination.
Data safety and adverse effects monitoring
The trial examined the integration of traditional care with modern care without involving invasive procedures or the administration of any drugs. Participant adherence in both intervention and control groups was assessed through self-reports and direct observation by trained field supervisors, with regular communication and feedback maintained between supervisors and researchers. Ethical approval was obtained from the institutional review board (IRB) of Bahir Dar University, College of Medicine and Health Sciences, with ethical review Ref No. 353/2021. Written consent for adult participants and written assent for pediatric participants were obtained from each participant.
Sample size determination
The sample size was determined using a two-sample comparison of the mean, incorporating data from a previous study on diagnosis and treatment delays in southwestern Ethiopia with a total median of 55 (interquartile range (IQR): 32–100) days (median: m1 = 55)(19). Assuming a 14-day reduction in diagnosis delay in the intervention group compared to the control group (m2 = 41), considering a type 1 error probability of 0.05, a 95% confidence interval, 80% power, and accounting for a 10% non-response rate, the sample size was 537. In consideration of the study being a cluster randomized control trial, the design effect was considered with a determined value of 0.95 based on the recommended intraclass correlation coefficient (ICC) of less than 0.052 for studies with more than 10 clusters. Then, by multiplying the original sample size by the design effect, the final sample size was 537*0.95 = 510(255 for each group).
Data collection
Structured interviewing questionnaires were used to collect data. The data were collected using interviewing questionnaires that contained sociodemographic, clinical, and behavioral variables, the onset of symptoms, health-seeking behavior, date of diagnosis, and treatment commencement. Experienced BSc nurses and public health officers collected the data.
Statistical analysis
Data were entered into EpiData version 4.6 and analyzed using Stata software version 17.0. Descriptive statistics were utilized to analyze the characteristics of the baseline data. Principal component analysis was used to calculate the household wealth index, considering factors such as land ownership and livestock. Initially, descriptive statistics were employed to determine the frequency, percentage, mean, and chi-square for comparing baseline and end-line data. The t-test was utilized to analyze the mean, mean difference, and standardized group mean difference (effect size) of the diagnosis delay. Comparisons of diagnosis and treatment delays within and between the intervention and control groups were conducted using t-tests.
The time from the onset of illness to the diagnosis of TB was used to calculate the incidence rate per person-day. The Kaplan-Meier survival function was estimated to determine the probability of time to diagnosis delay. The log-rank test was used to compare survival curve probabilities between the intervention and control groups. Subsequently, the Cox proportional hazard model was applied for semi-parametric multivariable analysis. Additionally, a parametric approach incorporating completely parametric survival models was used to address cluster variations more effectively. Univariable analysis was conducted to assess each explanatory variable, with variables showing significance at a level of 0.25 considered for multivariable analysis(20). In multivariable analysis, variables with a p-value of less than 0.05 were considered statistically significant. The Cox proportional hazard assumption was assessed using both graphical and statistical methods. The graphical method indicated that the assumption was met, but it relies on subjective judgment. We plotted the Weibull distributions of selected variables against their respective Kaplan-Meier curves to display the Cox-Snell residual and Nelson-Aalen cumulative hazard Cox-Snell residual, demonstrating how well the estimated Weibull survival plots fit the data (see Fig. 1). The Schoenfeld residual proportional hazard test confirmed that the assumption was met (statistically insignificant, with a p-value of 0.611). A multivariable accelerated failure time shared frailty model was utilized to determine the predictors of time to failed diagnosis. While the Cox proportional hazard model assumes a constant hazard ratio between individuals over time, our study involved 23 clusters (facilities) that may exhibit variations among them. The intra-cluster correlation was taken into account to address the unexplained covariates of the clusters. To address this correlation among the clusters, a shared frailty model was employed to identify the variance within each cluster by introducing a random effect model where individuals in a cluster are assumed to share the same frailty value. The hazard function is expected to follow a certain distribution and is influenced by an unobservable random frailty effect shared by participants within a cluster. In a shared frailty model assuming a Weibull distribution, the hazard function at time "t" for the "jth" individual, where "j = 1, 2,..., ni," in the "ith" group, where "i = 1, 2,..., g," is expressed as: hij(t) = Ziexp(¬β′xij) ρt ρ-1. Here, xij represents a vector of explanatory variables for the "jth" individual in the "ith" group, β is the vector of regression coefficients, ρt ρ−1 is the baseline hazard function, ρ is a shape parameter, and the zi are frailty effects shared by all "ni" individuals within the "ith" group.
When considering a parametric survival model characterized by its hazard function, h(t), all functions are affected by any covariates. Whether we parameterize the model as having proportional hazards (PH) regarding changes in covariate values or accelerated failure time (AFT) due to the covariates, the hazard function at time t for individual i with covariate Xi is given by: hi(t) = exp(xiβ)pt p−1.
A frailty model in the univariate case introduces an unobservable multiplicative effect, denoted by α, on the hazard. This means that conditional on the frailty, h(t|α) = αh(t), where α is a random positive quantity assumed to have a mean of one for the purpose of model identifiability and a variance of θ. Individuals with α > 1 are considered more frail, leading to an increased risk of failure, though the reasons for this frailty are not explained by the covariates. On the other hand, individuals with α < 1 are less frail and tend to survive longer, all else being equal (i.e., given a certain covariate pattern). Since α is a multiplicative effect that accounts for the cumulative impact of one or more omitted covariates, the relationship between the hazard and survival functions can be shown as the individual survival function conditional on frailty, S(t|α) = {S(t)}α, where S(t) is the survival function from a standard survival model that may include ancillary parameters and covariate effects.
Shared frailty distribution and parametrization: A parametric survival model follows a known distribution. We fitted the Weibull, exponential, log-logistic, lognormal, and generalized gamma distributions by considering both the gamma and inverse-Gaussian frailty distributions. The Weibull AFT inverse-Gaussian shared frailty model, which had the smallest AIC and BIC values, was selected to analyze the data. Finally, the variance of the random effect (θ), Kendall's Tau (τ), the regression coefficients, and the acceleration factor (δ) with a 95% confidence interval were estimated. The estimated variability (heterogeneity) in the population of clusters (facilities) was determined using the Weibull inverse-Gaussian shared model. In addition, the goodness of fit was checked using the Cox-Snell residuals plot. The model closely followed the 45-degree straight line, with a slight deviation on the left tail. This indicated that the model was well-fitted to the time-to-diagnosis delay (Fig. 1).