We obtained TB surveillance data from the Jiangsu Tuberculosis Information Management System (TBIMS), spanning January 1, 2011, to December 31, 2021. The original dataset contained 354,274 infection cases reported in Jiangsu Province during this period, including patient information such as names, ages, sex, occupations, sources of patients, types of diagnosis, tracking status, types of hospital, dates of birth, onset dates, diagnostic dates, and so on. A total of 332,091 patients were analyzed in this study following exclusion criteria: (1) patients diagnosed and reported between January 1, 2011, and December 31, 2021, at health care institutions outside Jiangsu Province (n=8010); (2) patients with missing critical information (n=12,889); and (3) patients whose standardized z value for the total diagnostic delay that exceeded 3 (ie, more than 3 SD from the mean; n=1266) [22].

Total diagnostic delay, defined as the interval from the onset of TB symptoms to formal diagnosis, comprises both patient delay and health system delay [13,23]. Patient delay refers to the time between symptom onset and the first medical consultation, while health system delay spans from the time of the first health care visit to diagnosis. This study focused exclusively on total diagnostic delay, as the dataset lacked information regarding the date of the first medical visit. At the individual level, the outcome was defined as the diagnostic delay status, which indicated whether a patient’s total diagnostic delay exceeded 28 days, a commonly adopted threshold for total TB diagnostic delay based on previous studies [13]. At the county level, the outcome was the TB diagnostic delay rate, which was calculated by dividing the number of delayed patients (ie, those with more than 28 d of total diagnostic delay) by the total number of patients in each county each year.

In addition to the individual-level TB surveillance data from Jiangsu TBIMS, a total of 9 annual county-level explanatory variables were included, categorized into three distinct domains: (1) demographic factors, comprising the annual proportions of older adult patients (≥60 y), male patients, local patients, and agricultural-worker patients among reported cases; (2) socioeconomic and health care indicators, including gross domestic product (GDP) per capita (adjusted to the 2021 Consumer Price Index), resident population size, TB incidence rate (per 1000 population), and the health care technicians (professionals per 1000 population); and (3) pandemic period, a binary variable (1=2020‐2021; 0=2011‐2019) introduced to adjust for the potential impact of social isolation policies and health care resource diversion during the COVID-19 pandemic. Data regarding health care technicians, GDP, and resident population were sourced from annual county statistical yearbooks, while other variables were aggregated directly from the TB surveillance data.

To explore the spatial correlation of the TB diagnostic delay rate across 89 districts and counties in Jiangsu Province from 2011 to 2021, we calculated the global Moran I to measure the spatial correlation [24,25]. When the result for Moran I is statistically significant, a positive value for Moran I suggests spatial clustering, while a negative value suggests spatial dispersion. The closer the value of I is to 1 or −1, the stronger the spatial association, while a value near 0 indicates a random spatial distribution of TB diagnostic delay rates.

Prior to the spatiotemporal modeling, multivariable binary logistic regression was used to identify risk factors associated with individual diagnostic delay status (binary outcome: 1 if total delay >28 d, 0 otherwise). Subsequently, a logit-link Bayesian Beta regression model was then applied, incorporating fixed effects for the proportion of older adult patients, proportion of male patients, proportion of local patients, proportion of agricultural-worker patients, GDP, TB incidence rate, number of health care technicians, and resident population, along with spatial and temporal random effects.

Specifically, let y(s,t) denote the TB diagnostic delay rate in year t over district or county s, with values ranging from 0 to 1. Here, t=1, 2, ..., 11 represents the years 2011 to 2021, and s=1, 2, ..., 95 represents the 95 counties or districts in Jiangsu Province, China. We further assumed that y(s,t) follows a β distribution with mean μ(s,t)∈(0,1) varying over time across counties or districts, and a constant precision parameter ϕ>0. Namely,

where the vector x(s,t) represents the regional-level variables in Table 3. The ϵ(s,t) term is an unstructured random effect in the model. We employed the first-order Gaussian random walks (RW1) model and the Besag-York-Mollié 2 model to capture the overall temporal random effect δt and spatial random effect us [26].

The Besag-York-Mollié 2 model was used to capture spatial random effects by combining structured and unstructured spatial components through a mixing parameter [27]. The spatial effect for area i is expressed as:

where τb is the overall precision, ui∗ is a standardized intrinsic conditional auto-regressive component capturing structured spatial dependence, and vi∗ is a standardized Gaussian noise term representing unstructured spatial variability. The mixing parameter ϕ∈[0,1], is a spatial smoothing parameter, measuring the proportion of the marginal variance explained by the structured random effect.

As a latent effect implemented in the R-INLA package, the first-order Gaussian random walk was used to model temporal dependence [26]. For a latent Gaussian field u=(u1,…,un)T, it is a random walk of order 1 if the increments Δui=ui−ui−1 are independent and identically distributed Gaussian random variables with zero mean and precision τ>0 (inverse variance).

Based on the delay rates estimated from the Bayesian spatiotemporal Beta model, we dichotomized at the median and then fitted a Bayesian spatiotemporal binomial model to identify factors associated with a higher likelihood of diagnostic delay. Furthermore, to explore the temporal dynamics and potential drivers of delay rate, we conducted a panel Granger causality analysis [28].

All statistical analyses were conducted using R software (version 4.4.1; R Foundation for Statistical Computing). The INLA, stats, and lmtest packages were used to conduct the Bayesian spatiotemporal modeling, logistic regression, and the Granger causality analysis, respectively, with default prior distributions specified for the INLA hyperparameters.

Anonymized data were obtained from the Jiangsu TBIMS, with all personal identifiers (eg, name and ID number) removed prior to analysis. The study protocol was reviewed by the ethical review board of the Jiangsu Provincial Center for Disease Control and Prevention (Jiangsu CDC) and granted an official exemption (acceptance: SL2025-B030-01), as the study was deemed retrospective and the data were deidentified. Data access and usage were strictly governed by a formal Data Usage Agreement between the Jiangsu CDC and the study authors. Informed consent was not required for this retrospective study using anonymized data.

Table 1 summarizes the demographic and clinical features of patients with TB with delayed (≥28 d) and nondelayed (<28 d) diagnoses, respectively, in Jiangsu Province (2011-2021). Among all patients, the majority were male participants (239,692/332,091, 72.18%), working in agriculture (204,983/332,091, 61.72%), and reported by general hospitals (148,960/332,091, 44.86%), with a high proportion of local residents (239,440/332,091, 72.10%).

Figure 1A shows the spatial distribution of average TB diagnostic delay rates in Jiangsu Province over the period from 2011 to 2021. The cities of Yancheng and Huaian generally exhibited the highest delay rates, and the Binhai district in Yancheng recorded a peak delay rate of 93.35% in 2016. Figure 1B illustrates the temporal trend across all 95 counties. An increasing trend was evident from 2011 to 2014, followed by a plateau. Notable distinct declines in the average delay rate were observed in specific months, including December 2018, November 2019, and November 2020.

Figure 2 presents the results of a logistic regression analysis used to investigate potential risk factors associated with individual TB diagnosis delay status. Specifically, males had significantly lower odds of experiencing delayed TB diagnosis compared to females (odds ratio [OR] 0.891, 95% CI 0.874-0.906), while patients older than 60 years had slightly higher odds of delayed TB diagnosis compared to those who were younger than 60 years (OR 1.070, 95% CI 1.053-1.088). We also found that patients who worked in agriculture had higher odds of delayed TB diagnosis than patients in other occupations. Particularly, the differences between patients working in agriculture and patients of occupations in education (OR 0.769, 95% CI 0.740-0.799), housekeeping (OR 0.807, 95% CI 0.791-0.824), service (OR 0.794, 95% CI 0.749-0.843), and worker (OR 0.778, 95% CI 0.757-0.800) were statistically significant. Local patients (from the same county/district) had significantly lower odds of delayed TB diagnosis than patients from different counties (OR 1.598, 95% CI 1.525-1.675) or cities (OR 1.508, 95% CI 1.480-1.536) in Jiangsu. Regarding types of diagnosis, the suspected cases had significantly lower odds of delayed TB diagnosis compared to the confirmed cases (OR 0.242, 95% CI 0.221-0.266), while the clinically diagnosed cases had significantly higher odds of delayed TB diagnosis (OR 1.021, 95% CI 1.003-1.037) compared to the confirmed cases. Patients diagnosed at community health centers (CHCs) had higher odds of delayed diagnosis compared to those diagnosed at designated TB hospitals (OR 1.103, 95% CI 1.075-1.133). However, patients diagnosed at local CDCs had considerably lower odds of delayed diagnosis than those diagnosed at designated TB hospitals (OR 0.687, 95% CI 0.668-0.707). Finally, the COVID-19 period was associated with significantly lower odds of diagnostic delay (OR 0.954, 95% CI 0.933-0.975).

Table 2 presents the global Moran I indices for TB diagnostic delay rates from 2011 to 2021. From 2011 to 2014, the Moran I values were close to zero (all P>.05), indicating no significant spatial autocorrelation. However, a marked increase in Moran I was observed from 2015 onwards, suggesting the emergence of substantial spatial clustering. Consequently, the Bayesian spatiotemporal Beta model was used to fully account for the spatial autocorrelation found in TB diagnostic delay rates across the years.

Table 3 presents the fixed effect estimates for the proportion of older adult patients, proportion of male patients, proportion of local patients, proportion of agricultural-worker patients, GDP, TB incidence rate, number of health care technicians, and resident population from the Bayesian spatiotemporal Beta model. For each 1-unit increase in the proportion of local patients, the TB diagnostic delay rate decreases by 33.9% (1−exp[−0.415], 95% CI 0.128-0.498). Additionally, each 100,000-person increase in resident population is associated with a 2% (95% CI 0.005-0.033) decrease in TB diagnostic delay. For random effects, the precision parameter τb was estimated as 3.411, indicating a moderate degree of spatial variation in the diagnostic delay rates across counties. In addition, the mixing parameter ϕ was estimated as 0.318, meaning that the structured spatial component accounts for 31.8% of the total spatial variation in TB diagnostic delay rate.

Figure S1 in Multimedia Appendix 1 illustrates the spatiotemporal distribution of the estimated TB diagnostic delay rate obtained from the Bayesian spatiotemporal Beta model. Higher estimated delay rates were concentrated in the northern coastal cities such as Huaian, Yancheng, and Lianyungang, while lower rates were observed in Xuzhou and Suzhou in the southern region. Among the 95 districts and counties over the 11-year period, Qingjiangpu of Huaian ranked among the top 3 with the highest estimated delay rates, with values of 0.648 in 2021, 0.645 in 2014, and 0.644 in 2018. Meanwhile, from 2011 to 2013, the delay rate exhibited an increasing temporal trend, which was particularly evident in Suqian. In contrast, Figure S2 in Multimedia Appendix 1 presents the observed delay rates, which show much larger fluctuations, reflecting random noise and potential instability due to small-area sample variation. The differences between Figures S1 and S2 in Multimedia Appendix 1 arise because the Bayesian spatiotemporal Beta model smooths random noise and incorporates both spatial and temporal dependencies. By integrating relevant covariates and modeling residual correlation through random effects, it effectively adjusts for unobserved dependence and confounding, producing spatially coherent and statistically reliable estimates.

The median of the estimated TB diagnostic delay rate was 0.352. Based on the obtained median, the results from the Bayesian spatiotemporal binomial model (Table S1 in Multimedia Appendix 1) show that counties with a higher proportion of agricultural workers and higher GDP were more likely to experience high diagnostic delays. Conversely, counties with larger shares of local patients, larger resident populations, and the COVID-19 period were associated with a lower risk of high delays. Furthermore, Table S2 in Multimedia Appendix 1 indicated that GDP, TB incidence, and health care technicians had significant Granger causal effects on the temporal changes in the delay rate (P<.05).

The sensitivity analysis confirmed the consistency of our main results. We applied penalized complexity priors to spatial and temporal model parameters to constrain model complexity and prevent overfitting (Tables S3-S6 in Multimedia Appendix 1) [29]. Across these analyses, the direction and statistical significance of the main variables remained consistent, with a slight change in the magnitude of some coefficients. To examine the impact of risk factors associated with TB diagnostic delays over shorter periods, we divided the study period into 2 subperiods (2011‐2015 vs 2016‐2021), representing China’s 12th and 13th Five-Year Plans for Tuberculosis Prevention and Control (Table S7 in Multimedia Appendix 1) [30]. We found that the significant associations between diagnostic delay and residency or occupation observed in the first subperiod disappeared in the second subperiod, likely due to expanded health care access and continuous public health developments in Jiangsu.