Materials and methods

Data source

The study data were mainly extracted from the GBD 2021 database, which evaluates the disease burden caused by 369 diseases and 87 risk factors in 204 countries, with detailed information on age, gender, region, and other dimensions. Data on OA attributed to smoking-induced digestive system diseases in the GBD 2021 database are mainly derived from death surveillance systems, vital registration systems, and the Chinese Center for Disease Control and Prevention [[5], [6], [7], [8]]. For this study, we extracted data on crude incidence rate, age-standardized incidence rate (ASIR), and disability-adjusted life years (DALYs) of OA in China from 1990 to 2021 to analyze the disease burden over this period. The data for this study were accessed on September 10, 2025. During the data extraction process, the authors did not have access to any personally identifiable images or data that could identify individual participants. The GBD database has released its 2023 updated version recently; however, official full public access to China-specific 2023 OA stratified data (age, gender, joint subtypes) is currently unavailable.

Statistical methods

Time trend analysis

The Joinpoint regression model (Joinpoint 4.2.0.1) was used to assess the time trend of OA disease burden in China. This model divides the entire study period into segments and identifies time points where significant trend changes occur. The annual percentage change (APC) and its 95% confidence interval were estimated for each segment, and the average annual percentage change (AAPC) and its 95% confidence interval were calculated for the entire period using a log-linear model. A permutation test was used to determine whether the APC of each segment was statistically significant, and the optimal model recommended by the Joinpoint regression program was adopted [[9], [10], [11]]. Set maximum allowed joinpoint numbers per regression based on sample size following standard epidemiological practice. An APC or AAPC >0 indicates an increasing trend, while a value < 0 indicates a decreasing trend. The significance level was set at α = 0.05 [12].

Autoregressive Integrated Moving Average (ARIMA) model

The ARIMA model is a popular and powerful method for time series prediction, combining the “memory” feature of autoregressive (AR) models, the stationarization capability of differencing (I), and the error term processing capability of moving average models. Before constructing ARIMA models, Augmented Dickey–Fuller (ADF) stationarity tests were performed on original time series to judge whether differencing transformation was required. The NINIC procedure in SAS software was used for optimal model order selection [13,14]. The upper limit of autoregressive order p and moving average order q was set to 5 to avoid overfitting. The values of p and q were determined based on autocorrelation and partial autocorrelation plots, and the optimal model was selected according to the minimum Akaike Information Criterion [15]. To further evaluate model adequacy, Ljung-Box white noise tests were conducted to detect residual autocorrelation; ten-fold cross-validation was also adopted to calculate out-of-sample mean absolute percentage error (MAPE) and assess generalization ability. Only models with Ljung-Box P > 0.05 and low cross-validation MAPE were retained for long-term extrapolation. The significance level was set at α = 0.05.

Brown's linear exponential smoothing

The Brown model was selected for projecting the Age-Standardized Incidence Rate (ASIR). Unlike the DALYs rate, the incidence trend exhibited a stable linear pattern without significant fluctuations. Brown's linear exponential smoothing is statistically optimal for forecasting data with a linear trend. Given the long-term projection horizon (20 years), this method provides more stable and robust estimates for such a consistent trend compared to ARIMA models, which are more sensitive to short-term stochastic variations.