This retrospective cohort study used data from the Jiangsu Pneumoconiosis Follow-up Online Reporting System, which initially comprised 27,239 registered patients. Stringent exclusion criteria were implemented to confirm data quality and reliability. Individuals were excluded if they were lost to follow-up (n = 403), had missing cause of death (n = 8,537), or presented with logical irregularities in their records (e.g., years of dust exposure exceeding chronological age or death dates preceding diagnosis dates) (n = 235). Consequently, a total of 18,064 patients first diagnosed in Jiangsu Province between 1960 and 2024 were included in the final analytical cohort. The detailed selection process is shown in Supplementary Figure S1. All patients in the final cohort (n = 18,064) had complete data for variables used in regression models. Patients with missing covariates were excluded during data cleaning. Therefore, complete case analysis was used without imputation. Characteristics of excluded patients are presented in Supplementary Table S1.

Diagnoses and staging were performed at certified medical institutions in Jiangsu Province authorized to diagnose occupational diseases. Diagnostic criteria complied with the national occupational health standards in effect at the time, specifically the Diagnostic Criteria of Occupational Pneumoconiosis (GBZ 70–2002, GBZ 70–2009, GBZ 70–2015) and the Diagnostic Criteria for Pneumoconiosis on X-ray (GB 5906–1997). The diagnostic procedure included a consensus-based evaluation by an expert panel of at least three qualified radiologists. This evaluation incorporated occupational dust exposure history, clinical symptoms, workplace epidemiological data, and high-kilovolt chest X-rays. Based on the severity of pulmonary fibrosis, cases were classified into Stages I, II, or III. The study was approved by the Ethics Committee of the Jiangsu Provincial Center for Disease Control and Prevention (Approval No. JSJK2022-B002-01).

Follow-up duration was calculated from the date of initial pneumoconiosis diagnosis to either the date of death or the last follow-up. Follow-up was closed on December 31, 2024. Vital status was ascertained through a combination of active telephone follow-up and linkage with the provincial death registration system. Outcomes were ascertained based on underlying cause-of-death codes recorded in death certificates, classified according to ICD-10. Death information was collected through registry linkage, CDC follow-up, and institutional reporting. Underlying cause of death was determined by trained coders following ICD-10 rules. ICD-10 was used because it remains the standard for China’s medical documentation and death registration during the study period. The follow-up system does not systematically track migration status. Patients who could not be contacted through telephone follow-up and had no death records in the provincial registry were classified as lost to follow-up and excluded from the analysis. The final analytical cohort consisted only of patients with verifiable vital status through either confirmed death records or successful telephone contact.

The primary endpoint was pneumoconiosis-related death, defined as deaths where the underlying cause was coded as J60–J65. Non-pneumoconiosis-related death, comprising mortality from any cause other than these specific codes, was classified as competing death (13). This definition was intentionally conservative to maximize specificity for pneumoconiosis as the underlying cause. Patients who survived to the end of follow-up were censored.

Disease types were classified as silicosis, coal workers’ pneumoconiosis (CWP), welder’s pneumoconiosis (WP), or other pneumoconiosis. The “other pneumoconiosis” category included asbestosis, talc pneumoconiosis, and other dust-related lung diseases not classified as silicosis, CWP, or WP. Regions were classified as southern Jiangsu (Nanjing, Suzhou, Wuxi, Changzhou, Zhenjiang), central Jiangsu (Nantong, Yangzhou, Taizhou), or northern Jiangsu (Xuzhou, Huai’an, Yancheng, Lianyungang, Suqian) based on patients’ residential locations.

The Fine–Gray proportional subdistribution hazards model was applied to account for the risk of competing death. Univariate analyses were initially performed to evaluate crude associations between baseline characteristics and pneumoconiosis-related mortality. Subsequently, to adjust for potential confounders and identify independent predictors, age at diagnosis, gender, industry category, region, disease type, era of diagnosis, baseline stage, and years of dust exposure were entered into the multivariable model. These variables were selected based on clinical relevance and prior evidence of their associations with pneumoconiosis outcomes. Information on smoking status, comorbidities, and workplace-level exposure controls was not systematically available in the registry data. Effect estimates were expressed as subdistribution hazard ratios (SHR) with 95% confidence intervals (CI). Multicollinearity was assessed using the variance inflation factor (VIF); all VIF values were less than 5, indicating the absence of severe multicollinearity (14). Additionally, the proportional hazards assumption was evaluated using log-minus-log survival curves and scaled Schoenfeld residual plots. Detailed results are provided in Supplementary Figures S2, S3.

To develop and validate a tool for individualized prognostic assessment, the entire cohort was randomly partitioned into training and validation sets at a 7:3 ratio. Guided by the multivariable Fine–Gray model results and the principle of parsimony, variables demonstrating independent statistical significance (p < 0.05) and substantial prognostic contribution were selected. These variables included age at diagnosis, disease type, era of diagnosis, and baseline stage. A nomogram was subsequently established to predict the 3-, 5-, 10-, and 15-year cumulative incidence of pneumoconiosis-related death. Model performance was evaluated based on discrimination and calibration. Discrimination was quantified using the Concordance Index (C-index), with higher values indicating better predictive accuracy. Calibration was assessed by calibration curves to visually demonstrate the agreement between predicted and observed probabilities. Internal validation was performed using bootstrap resampling with 1,000 repetitions to calculate standard bootstrap confidence intervals for the C-index and ensure model robustness.

To verify the robustness of the primary findings among different subpopulations, subgroup analyses were conducted, stratified by age at diagnosis (<50 vs. ≥50 years), era of diagnosis (before 2000 vs. 2000 or later), and disease type (silicosis vs. coal workers’ pneumoconiosis). Within each subgroup, the prognostic impact of baseline stage (Stage II and Stage III vs. Stage I) on pneumoconiosis-related death was evaluated. Interaction tests were performed to assess potential effect modification by these stratifying factors. Furthermore, to demonstrate the necessity of employing a competing risk framework, a traditional Cox proportional hazards model was fitted for methodological comparison. In the Cox model, competing deaths were treated as censored observations, with covariates identical to those in the Fine–Gray model. Discrepancies between hazard ratios (HR) from the Cox model and SHRs from the Fine–Gray model were quantified using relative bias, calculated as (|HR−SHR|/SHR) × 100%.

The normality of continuous variables was evaluated using the Shapiro–Wilk or Kolmogorov–Smirnov tests. Normally distributed variables were expressed as mean ± standard deviation (SD) and compared using one-way analysis of variance (ANOVA). Non-normally distributed variables were reported as medians (interquartile range [IQR]) and compared using the Kruskal-Wallis test. Categorical variables were summarized as frequencies (percentages) and analyzed using Pearson’s chi-square test. All statistical analyses were conducted using R software (version 4.3.1). The cmprsk package (version 2.2–11) was used for cumulative incidence estimation and Fine–Gray modeling; the survival package (version 3.5–5) for Cox regression; and the rms (version 6.7–0) and regplot (version 1.1) packages for nomogram construction. All tests were two-sided, with statistical significance defined as p < 0.05.

As presented in Table 1, the final analysis included 18,064 patients with pneumoconiosis. By the end of follow-up, 15,441 (85.5%) patients remained alive, while 626 (3.5%) had died from pneumoconiosis-related causes and 1,997 (11.1%) from competing causes. The cohort was predominantly male (93.4%), with a median age at diagnosis of 52.0 years (IQR: 45.0–62.0 years). The median age at diagnosis was significantly higher in the pneumoconiosis-related death (59.0 years) and competing death groups (60.0 years) compared with the survival group (51.0 years; p < 0.001). Regarding occupational exposure, the mining industry predominated (49.8%), followed by manufacturing (25.6%). Geographically, 60.1% of patients were from southern Jiangsu, 36.2% from northern Jiangsu, and 3.7% from central Jiangsu.

Patients from northern Jiangsu comprised 53.8% of competing deaths versus 36.6% of survivors (p < 0.001). Notably, the proportion of miners was highest in the pneumoconiosis-related death group (61.0%), significantly exceeding that in the survival (50.0%) and competing death groups (44.7%). Additionally, the median duration of dust exposure was longest in the pneumoconiosis-related death group (18.0 years) compared with the survival (15.0 years) and competing death groups (14.0 years) (p < 0.001). Clinically, silicosis was the most common disease type (68.0%). While 84.9% of the overall population presented with Stage I disease at baseline, the prevalence of Stage II and Stage III disease was markedly higher in the pneumoconiosis-related death group (27.8 and 11.7%, respectively) than in the survival group (10.9 and 2.6%; p < 0.001). Disease progression was most frequent in the pneumoconiosis-related death group (15.0%) compared with the survival (6.9%) and competing death groups (7.5%), with statistically significant differences across the three groups (p < 0.001). Significant disparities were also observed in the era of diagnosis and baseline stage among the outcome groups. Median follow-up time varied by era [30.0 years (IQR: 26.0–33.0) for before 2000, 17.0 years (15.0–21.0) for 2000–2010, and 8.0 years (5.0–11.0) for after 2010] and by baseline stage [17.0 years (11.0–27.0) for Stage I, 17.0 years (10.0–23.0) for Stage II, and 11.0 years (6.0–17.0) for Stage III] (Supplementary Table S2).

The overall cumulative incidence curves illustrate long-term trends in cause-specific mortality within the study cohort (Figure 1). The analysis indicated that the cumulative incidence of competing death surpassed that of pneumoconiosis-related death throughout the follow-up period. Furthermore, the divergence between the two curves progressively widened over time. Specifically, at 10 years of follow-up, the cumulative incidence of competing death and the event of interest was 4.6 and 1.3%, respectively. By 30 years, these figures rose to 16.3 and 5.0%, respectively. At the 50-year endpoint, the cumulative incidence of competing death reached 39.0%, substantially surpassing the 13.4% observed for pneumoconiosis-related death. These data underscore the substantial burden of competing risks on the prognostic outcomes of this study population.

Stratified analyses of pneumoconiosis-related death revealed significant associations between baseline clinical characteristics and cumulative mortality (Figure 2). In the analysis stratified by baseline stage (Figure 2A), patients with Stage III disease exhibited a steep increase in cumulative mortality, significantly exceeding that of patients with Stage I and Stage II (Gray’s test p < 0.001). Distinct temporal disparities were observed for the era of diagnosis (Figure 2B). Compared with patients diagnosed before 2000, those diagnosed after 2010 demonstrated a significantly reduced cumulative incidence of pneumoconiosis-related death (p < 0.001). Furthermore, marked prognostic differences existed among disease types (Figure 2C); cumulative mortality for patients with silicosis rose rapidly during early follow-up and was significantly higher than for patients with coal workers’ pneumoconiosis. Age at diagnosis also showed a distinct gradient (Figure 2D), with the older group (≥50 years) showing markedly higher cumulative mortality than the younger group (<50 years) (p < 0.001). Supplementary Figure S4 presents the corresponding cumulative incidence curves for competing death stratified by baseline characteristics. Stratified analyses by disease progression and industry categories are provided in Supplementary Figure S5.

Stratified analysis by region revealed contrasting patterns (Supplementary Figure S6). For pneumoconiosis-related death, cumulative incidence curves overlapped across regions (Gray’s test p = 0.168). For competing death, northern Jiangsu showed significantly higher cumulative incidence than southern regions (Gray’s test p < 0.001), indicating that regional disparities primarily affect non-pneumoconiosis mortality.

Evaluation of the proportional hazards assumption using log-minus-log survival curves (Supplementary Figure S2) and scaled Schoenfeld residual tests (Supplementary Figure S3) showed that although the global test was statistically significant (p < 0.001), this primarily reflects the large sample size. Individual covariates including era 2000–2010 (p = 0.655), welder’s pneumoconiosis (p = 0.736), and central Jiangsu (p = 0.047) showed no significant time-varying effects. Log-minus-log curves for stage, era, disease type, region, and age group demonstrated generally parallel patterns, supporting the adequacy of the proportional hazards assumption for the core variables.

Univariate Fine–Gray competing risk regression analysis indicated that age at diagnosis, employment in the manufacturing industry, specific disease types (CWP and WP), era of diagnosis, baseline stage, and years of dust exposure were all significantly associated with the risk of pneumoconiosis-related death (Table 2). In contrast, gender and other industry categories did not show notable correlations at the univariate level.

In the multivariable model adjusting for all candidate variables, age at diagnosis, regional economic level, disease type, era of diagnosis, and baseline stage emerged as independent predictors of pneumoconiosis-related death. Baseline stage exhibited the strongest prognostic impact, showing a distinct gradient. Compared with Stage I, the adjusted SHRs for pneumoconiosis-related death were 2.94 (95% CI 2.41–3.60) for Stage II and 5.18 (95% CI 3.91–6.87) for Stage III. Age at diagnosis was positively associated with mortality risk; each 1-year increase corresponded to a 9% increase in risk (SHR 1.09, 95% CI 1.08–1.10). Region demonstrated a significant independent effect. Compared with southern Jiangsu, patients from central Jiangsu (SHR 1.92, 95% CI 1.24–2.98) and northern Jiangsu (SHR 1.44, 95% CI 1.16–1.78) exhibited significantly elevated pneumoconiosis-related mortality risk. Marked variations in mortality risk were also identified across disease types. Relative to silicosis, patients with CWP had a 38% lower risk (SHR 0.62, 95% CI 0.49–0.79).

Notably, after multivariable adjustment, the association with the era of diagnosis reversed from a positive to a negative direction compared with the univariate analysis, revealing a significant independent protective effect. Compared with patients diagnosed before 2000, the risk was 34% lower for those diagnosed between 2000 and 2010 (SHR 0.66, 95% CI 0.54–0.80) and 62% lower for those diagnosed after 2010 (SHR 0.38, 95% CI 0.28–0.52). Additionally, employment in public administration and social organizations was associated with lower mortality risk compared with mining (SHR 0.66, 95% CI 0.47–0.92). In contrast, gender, employment in manufacturing or other industries, WP, other pneumoconiosis types, and years of dust exposure did not retain statistical significance in the multivariable model after adjusting for other covariates.

Univariate Fine–Gray regression analysis indicated that age at diagnosis, gender, industry category, region, disease type, era of diagnosis, baseline stage, and years of dust exposure were all significantly associated with the risk of competing death (Supplementary Table S3). In the multivariable model, age at diagnosis, gender, region, disease type, baseline stage, and years of dust exposure emerged as independent predictors. Age at diagnosis was positively associated with mortality risk; each 1-year increase corresponded to a 9% rise in the risk of competing death (SHR 1.09, 95% CI 1.09–1.10). Females exhibited a lower risk than males (SHR 0.77, 95% CI 0.63–0.94). Region demonstrated a strong independent effect on competing mortality. Compared with southern Jiangsu, patients from central Jiangsu (SHR 1.44, 95% CI 1.09–1.91) and northern Jiangsu (SHR 1.90, 95% CI 1.69–2.14) exhibited significantly elevated competing death risk, with northern regions showing a 90% higher risk.

Regarding disease type, relative to silicosis, patients with CWP demonstrated a 46% lower risk (SHR 0.54, 95% CI 0.46–0.64), and those with WP showed a 51% lower risk (SHR 0.49, 95% CI 0.27–0.90). Baseline stage exhibited a significant but modest effect on competing death. Compared with Stage I, both Stage II (SHR 1.49, 95% CI 1.29–1.72) and Stage III (SHR 1.41, 95% CI 1.10–1.80) were associated with increased competing risk. Years of dust exposure showed a weak but statistically significant negative correlation with the risk (SHR 0.99, 95% CI 0.98–1.00). Notably, although industry category and era of diagnosis showed strong associations in the univariate analysis, these associations attenuated and failed to maintain statistical significance in the multivariable model after adjusting for covariates such as age, region, and disease type.

Subgroup analyses demonstrated that the prognostic impact of baseline stage on pneumoconiosis-related death remained robust across most stratifications (Figure 3). Regarding Stage II (overall SHR 2.82) and Stage III (overall SHR 4.83), no notable interactions were observed with age at diagnosis or era of diagnosis (P for interaction > 0.05). However, disease type significantly modified the effect of stage (P for interaction = 0.011 for Stage II; p = 0.035 for Stage III). Among patients with silicosis, Stage II (SHR 3.45, 95% CI 2.74–4.34) and Stage III (SHR 6.25, 95% CI 4.59–8.52) were associated with substantially elevated mortality risks. In contrast, among patients with CWP, the associations of Stage II (SHR 1.41, 95% CI 0.81–2.43) and Stage III (SHR 1.44, 95% CI 0.47–4.45) with mortality were not statistically significant.

Subgroup analyses regarding competing death showed distinct patterns (Supplementary Figure S7), indicating that the prognostic impact of baseline stage on competing death was generally modest (SHR 1.25–1.37). Notably, significant effect modification by disease type was observed for the comparison between Stage II and Stage I (P for interaction = 0.015). Specifically, the risk of competing death was significantly elevated among patients with Stage II silicosis (SHR 1.40, 95% CI 1.20–1.64), whereas no significant association was found for those with CWP (SHR 0.82, 95% CI 0.53–1.25). Furthermore, no significant interaction by disease type was detected for the comparison between Stage III and Stage I.

The study population was randomly partitioned into a training cohort (N = 12,645) and a validation cohort (N = 5,419) at a 7:3 ratio. Guided by the multivariable Fine–Gray model, variables exhibiting considerable prognostic contribution, including age at diagnosis, disease type, era of diagnosis, and baseline stage, were selected to construct a nomogram predicting the 3-, 5-, 10-, and 15-year cumulative incidence of pneumoconiosis-related death (Figure 4). In the nomogram, age at diagnosis and baseline stage spanned the widest point ranges, indicating that these two variables contributed most markedly to individualized prognostic assessment. Validation analysis showed favorable discrimination in both cohorts. The C-indices were 0.769 (95% CI 0.742–0.797) in the training cohort and 0.753 (95% CI 0.714–0.791) in the validation cohort. The consistency between these values suggested the absence of significant overfitting. Calibration analysis showed favorable agreement between predicted and observed probabilities. At the 3-, 5-, 10-, and 15-year time points, the intercepts of the calibration curves were close to 0 (range −0.001 to −0.010) while slopes ranged from 1.10 to 1.33 (Supplementary Figure S8). For calibration assessment, patients in the training cohort were divided into five groups based on quintiles of predicted risk, with each quintile containing approximately 1,069 to 1,099 patients. The wider confidence intervals observed at higher predicted probabilities reflect the lower absolute number of events in these groups. The calibration curves generally aligned with the 45-degree ideal reference line, indicating consistency between model-predicted probabilities and actual results.

Supplementary Table S3 compares parameter estimates between the Fine–Gray and Cox models. The magnitude and direction of discrepancy varied across covariates. The Cox model showed substantial differences for era of diagnosis (relative bias +36.0% for 2000–2010 and +65.6% for after 2010) and moderate differences for baseline stage (+12.5% for Stage II and +5.2% for Stage III). However, several covariates showed minimal discrepancy: age (+3.5%), dust exposure (−0.4%), and regional economic level (+1.9% for central Jiangsu and −2.3% for northern Jiangsu). Some covariates showed negative bias, with Cox producing smaller estimates than Fine–Gray for CWP (−15.8%) and female gender (−9.6%).