Our study complies with the Guideline for Accurate and Transparent Health Estimates Reporting (GATHER) recommendations (Supplementary Materials S3).

Data sources and extraction

This study is part of the Global Burden of Disease Study 2021, which provides comprehensive estimates of disease burden for 371 diseases and injuries across 204 countries and territories. GBD 2021 quantifies incidence, prevalence, mortality, years of life lost (YLLs), years lived with disability (YLDs), and disability-adjusted life years (DALYs) for 25 distinct age groups, stratified by sex and Socio-Demographic Index (SDI) [8, 9].

Data for this study were manually retrieved from the Global Health Data Exchange (GHDx) query tool (available at https://ghdx.healthdata.org/gbd-2021). We extracted pre-processed estimates of incidence, prevalence, and disability-adjusted life years (DALYs) for pediatric urolithiasis, stratified by age (< 5 years, 5–9 years, 10–14 years), sex, year (1990–2021), and location (5 socio-demographic index super regions, 21 subregions, and 204 countries/territories). Data were independently extracted by two authors and subsequently verified by two additional authors. Manual cross-checking was conducted to compare key variables including age groups, diseases, locations, and years between datasets. In addition, consistency checks were performed using custom scripts to ensure the reliability of the extracted data. The dataset included absolute counts and rates (per 100,000 population), along with 95% uncertainty intervals (UIs) and SDI values. No additional data cleaning or handling of missing data was performed, as the GBD 2021 database provides quality adjustments and imputations applied by the GBD collaborative team[11].

Definition of the urolithiasis

In GBD 2021, urolithiasis was identified according to the International Classification of Diseases (ICD), Tenth Revision (ICD- 10), and ICD- 9. Urolithiasis was defined as the formation of stone anywhere along the genitourinary tract (coded as N20–N23.0, 592–592.9, 594–594.9)[12].

Data calculation

For the raw GBD data, uncertainty intervals (UIs) were provided to account for variability and measurement uncertainty in the original estimates. The 95% UIs were derived from the 25 th and 975 th percentiles of 1,000 draws, representing the range of plausible values. GBD studies draw on a wide array of data sources and models, so uncertainty stems not only from sampling error (which is addressed by confidence intervals, or CIs) but also from factors such as data quality, model assumptions, and parameter estimates.

The 0–14 age group age-standardized rates (ASRs) were obtained using data from < 5 years, 5–9 years, and 10–14 years, combined with the ageadjust.direct function from the epitools package in R[13, 14]. The 2021 GBD global standard population was used as the reference population[12]. The 95% confidence intervals (CIs) for the ASRs were calculated using the gamma distribution method to account for the uncertainty in the estimates[15]. Although the original GBD estimates include 95% UIs that incorporate a broader range of uncertainty sources such as model-based and data input variability, the ASRs in this study were derived from GBD point estimates using standard epidemiological procedures. As a result, the reported CIs reflect only sampling variability and may not capture the full uncertainty inherent in the original GBD data.

Joinpoint regression

Joinpoint regression analysis was conducted to estimate the average annual percent change (AAPC) and assess long-term trends at both the global and Socio-demographic Index (SDI) levels. The maximum number of joinpoints was set at five (six-line segments) to ensure sufficient model flexibility while preventing overfitting. This upper limit was chosen to balance trend detection sensitivity and statistical robustness, avoiding excessive responsiveness to short-term fluctuations[16].

The Monte Carlo permutation test was applied to determine the actual number of joinpoints, minimizing the sum of squared errors while adjusting for multiple testing. To ensure stable trend estimation, a parametric method was used to compute the annual percent change (APC) and AAPC[17]. The sensitivity analyses and model fit statistics are provided in supplementary Sect. 4.

Cross-country inequalities analysis

SDI is a measure of socio-economic development based on income per capita, average years of schooling for females under 25, and fertility rates. Countries are ranked by SDI from lowest to highest on a scale from 0 to 1. The SDI values are categorized into five levels: low SDI (0–0.455), low-middle SDI (0.455–0.608), middle SDI (0.608–0.690), high-middle SDI (0.690–0.805), and high SDI (0.805–1.0)[12].

To quantify inequalities in the distribution of disability-adjusted life years across countries ranked by their SDI, we employed the Slope Index of Inequality (SII) and the Relative Concentration Index (RCI).

SII measures absolute inequality in disability-adjusted life years by comparing the rates between the most and least advantaged countries along the SDI spectrum. It assumes a linear relationship between disability-adjusted life years and SDI, meaning health burdens increase or decrease in a predictable way as SDI changes. Additionally, SII incorporates population size, ensuring that the results reflect true health disparities across countries of different sizes. RCI focuses on relative inequality, measuring how DALYs are distributed across SDI groups in proportion to their population size. It normalizes the data to show whether DALYs are disproportionately concentrated in either high-SDI or low-SDI countries. RCI assumes that SDI rankings represent socio-economic hierarchies and adjusts for population size to make comparisons more accurate, accounting for the differences in country sizes[18].

The SII was calculated using a repeated iterative weighted least squares model in R. The model was fitted with the rlm function from the MASS package, which accounts for heteroscedasticity in the data. The RCI was calculated using the conindex command in Stata[19, 20]. Detailed explanations of SII, RCI, and SDI ranking can be found in the supplement materials Sect. 5.

Spearman correlation and regression analysis

To examine the relationship between disability-adjusted life years and SDI at the country level, we calculated Spearman correlation coefficients to determine the strength and direction of the non-linear relationship between these variables. Additionally, we applied locally estimated scatterplot smoothing (LOESS) regression to visualize the trend in disability-adjusted life years across the SDI spectrum.

Predictive analysis

To forecast the age-standardized incidence rate from 2022 to 2045, we employed a Bayesian generalized age-period-cohort (BAPC) power model. Model parameters were estimated within a Bayesian inference framework, incorporating the Iterative Leave-out-Node Approach (ILNA) to enhance prediction accuracy[21]. We specified a second-order random walk (RW2) for both period and cohort effects to allow flexible temporal trends, and used Penalized Complexity (PC) priors on the precision parameters to control smoothness. The age effect was modeled as a fixed categorical factor. Sensitivity analyses were conducted to assess the robustness of model results to different prior distributions and smoothness assumptions. The RW2 model with a PC prior (0.5, 0.05) was adopted for final projections.

The estimated parameters were applied to Global Burden of Disease (GBD) Study population projections to predict age-standardized rates through 2045. Credible intervals (CIs) were derived using Integrated Nested Laplace Approximations (INLA), with the 2.5 th and 97.5 th percentiles defining the 95% CI.

Statistical analysis and data visualization

Most statistical analyses and data visualization were conducted using R (version 4.3.3; R Foundation for Statistical Computing, Vienna, Austria). Joinpoint Trend Analysis software (version 5.2.0.0; National Cancer Institute, Rockville, MD, USA) was employed for trend analysis. The epitools (version 0.5–10.1) was used for age-standardized rate calculations, while INLA (version 23.09.09) and BAPC (version 0.0.36) were applied for BAPC projections. MASS package (version 7.3–60.0.1) was used to calculate SII, and Stata Statistical Software: Release 16 (StataCorp, 2019) was used for calculating RCI.A statistical significance threshold of p < 0.05 was applied across all analyses.