This study, conducted in September 2023 in Anhui, China, involved students aged 14–17 from nine pilot high schools. These pilot schools included both general senior high schools and vocational schools, and were located in central and non-central cities, ensuring the inclusion of students with diverse socioeconomic backgrounds. The inclusion criteria required participants to meet the following conditions: (1) no history of major illnesses; (2) no plans to transfer schools or relocate during the upcoming year; (3) the ability to participate in follow-up surveys for 1 year.

At baseline, participants underwent anthropometric assessments, including body mass and height measurements with digital scales and wall-mounted stadiometers, conducted by trained researchers. BMI was calculated as weight in kilograms divided by height in meters squared ( BMI = weight ( kg ) heigh t 2 ( m 2 ) ) . Basal metabolic rate (BMR) was estimated using the FAO/WHO/UNU adolescent (10–18 y) predictive equations: for male, BMR (kcal/day) = 16.6 × Weight (kg) + 77 × Height (m) + 572; for female, BMR (kcal/day) = 7.4 × Weight (kg) + 482 × Height (m) + 217 (23, 24). Overweight and obesity status was classified based on age- and sex-specific BMI reference standards established by the Working Group on Obesity in China for school-aged children and adolescents (25). In addition, students and their parents completed a structured electronic questionnaire through a publicly accessible online system available on both mobile phones and computers. Developed based on prior studies (7–9, 26, 27), the questionnaire covered a range of factors, including genetic predispositions, socio-demographic characteristics, daily habits, physical activity patterns, self-perception of body status, and health literacy. The system required participants to complete all items before submission, ensuring no missing data at the individual question level. After 1 year, a follow-up anthropometric assessment was conducted to evaluate BMI changes in the cohort, with assessors blinded to baseline measurements.

This study conducted a power analysis based on a medium effect size (Cohen’s f² = 0.15) and a significance level of 0.05 to assess whether the sample size was adequate for detecting meaningful effects in the statistical analyses (28). To assess potential clustering by school, we fit a two-level random-intercept linear mixed-effects model and computed the intraclass correlation coefficient (ICC) as the ratio of between-school variance to total variance (29). The impact of loss to follow-up was examined by comparing baseline characteristics of included vs. excluded participants using the standardized mean difference (SMD), reported as absolute values (|SMD|) (30). Differences were considered negligible when |SMD| < 0.10, and larger values were regarded as imbalanced (31). In addition, we quantified differential attrition in key subgroups (age, gender group, and baseline BMI category) by reporting attrition rates and risk differences (RDs) with 95% CIs relative to a prespecified reference level. Descriptive statistics were reported as means ± standard deviation for continuous variables and as counts for categorical variables. To assess regression to the mean (RTM), the change score was computed (Δ = Y₂ − Y₁, with Y₁ = baseline BMI and Y₂ = follow-up BMI), and a simple linear model was fitted: (Δ = α + ꞵY₁ + ε). A negative β indicates that higher baseline values are associated with greater negative change (shrinkage toward the mean). The point estimate of β, its 95% confidence interval (CI), and the p-value were reported. R² quantified the proportion of variance in Δ attributable to RTM, with a 95% CI obtained via nonparametric bootstrap.

This study applied the same data preprocessing to all algorithms to keep an identical feature space. Z-score standardization was applied to numerical features, ordered categorical variables were encoded with prespecified ordinal levels, and unordered categorical variables were represented using one-hot encoding. All preprocessing components were fit on training data (or training folds) only and applied to validation/test sets to prevent leakage.

In the development of our prediction model, we adhered to the TRIPOD+AI checklist, and the completed checklist is provided as Supplementary Table 1.

To establish a BMI prediction model, the following regression algorithms were utilized, including CatBoost (CB), LightGBM (LGBM), Neural Network (MLP), Decision Tree (DT), Support Vector Regressor (SVR), and K-Nearest Neighbors (KNN). These models represent a range of learning paradigms, allowing for a comprehensive comparison across different modeling strategies. Among them, tree-based models such as CB and LGBM are especially well-suited for handling high-dimensional, noisy, and heterogeneous health data. First, the dataset was split into training (80%) and independent test (20%) sets, with the test set held out throughout model development and tuning. On the 80% training set, all algorithms underwent 5 × 5 nested cross-validation (CV), with the inner loop performing feature selection and hyperparameter tuning and the outer loop providing unbiased performance estimates (32). Within each inner loop, predictors were prescreened on the inner-training folds using the univariable analyses (Pearson correlation for continuous variables and analysis of covariance (ANCOVA) for categorical variables; p < 0.05), followed by five-fold CV with grid search to select the hyperparameter set with the best mean validation score. Features selected across the five inner folds were aggregated by selection frequency to form that outer fold’s consensus feature set (33). In the outer loop, models were retrained on the outer-training folds using the consensus features and inner-optimal hyperparameters, and then evaluated on the outer test folds for unbiased assessment. Model performance was assessed based on four widely used indicators: root mean squared error (RMSE), mean squared error (MSE), mean absolute error (MAE) and coefficient of determination (R²). MSE reflects the average squared difference between predicted and actual values, capturing overall model fit. RMSE, the square root of MSE, reports the typical prediction error in the outcome’s original units. MAE provides a direct measure of the average prediction error. R² quantifies the proportion of variance in the outcome that is explained by the model, indicating its explanatory power. After completing all outer folds, performance on the outer test folds was summarized as mean ± SD. The model with the lowest RMSE/MSE/MAE and highest R² was deemed optimal. All outer-fold feature sets were then combined using the same frequency rule to obtain the final feature set. Final hyperparameters were selected via five-fold CV with grid search, then the model was retrained on the entire 80% training set. Complete model-tuning details (hyperparameter search grids, CV folds, seeds, early-stopping settings) are provided in Supplementary Table 2.

Heteroscedasticity was assessed on the 80% training set using out-of-fold residuals via a Breusch–Pagan test (α = 0.05) to examine whether prediction errors varied with baseline BMI (34). Results were reported without calibration when the test was not significant, and weighted least-squares (WLS) calibration was applied otherwise (35, 36). Residuals were computed from training out-of-fold predictions and used to fit an empirical variance model to derive sample weights. Using these weights, a linear recalibration of observed versus predicted values was fit on the training data to obtain fixed intercept and slope. The learned weighting function and coefficients were then applied once to the test set without any refitting, avoiding information leakage.

Overall model performance was compared on the 20% independent test set. First, generalization was assessed with 1,000 bootstrap resamples, reporting RMSE, MSE, MAE, and R² with 95% confidence intervals (CIs) for each model. Paired tests were conducted by bootstrapping the paired differences in RMSE using identical resamples across models (defined as comparator minus best model). This quantified the incremental benefit of the best-performing model, reporting ΔRMSE with its 95% CIs. Second, for the best-performing model, incremental benefit over a trivial baseline (predicting follow-up BMI = baseline BMI) was quantified on the same independent test set using paired bootstrap with 1,000 resamples. Paired differences were defined as ΔRMSE/ΔMSE/ΔMAE = baseline − best and ΔR² = best − baseline, and Δ values were reported with 95% CIs. Third, to assess robustness to a dominant predictor, a sensitivity analysis excluding baseline BMI was performed. Under an identical modeling pipeline to the primary analysis, we retrained and evaluated the model without baseline BMI using the best-performing algorithmic framework, and reported RMSE, MSE, MAE, and R² on the same independent test set.

Stratified performance was evaluated across clinically relevant subgroups, including gender (male vs. female), age groups (14–15 vs. 16–17 years), and baseline BMI category (normal vs. overweight/obesity). For each subgroup, we reported RMSE, MSE, MAE, and R² with 95% CIs. We also computed and reported between-group differences (Δ) in each metric with 95% CIs. Subgroup heterogeneity was assessed using permutation tests, and the significance level was set at 0.05.

Using the same independent test set, overall and stratified error analyses were performed for the calibrated best-performing model. For overall error visualization, we generated predicted-versus-observed scatterplots with a smoothing line and Bland–Altman plot. To characterize the error distribution, we calculated mean error ± SD, mean absolute error (MAD), the interquartile range of |error| (IQR|e|), and the 90th/95th percentiles of |error| (P90|e|/P95|e|), and reported these metrics both overall and stratified by gender, age group, and baseline BMI category.

This study applied SHAP to provide both global and local interpretation of the best-performing model. We quantified the contribution of each feature to model predictions by examining feature interactions, and decomposed individual predictions into additive feature contributions, using visualizations to convey overall patterns and individual differences. Non-modifiable features were excluded from SHAP visualizations. For global interpretation, we generated a SHAP summary plot based on mean absolute SHAP values across all samples, ranking modifiable predictors for BMI. In addition, we computed SHAP interaction values and displayed a heatmap of their mean strength across samples. SHAP dependence plots were used to examine the overall effect shapes of specific features across the cohort and to highlight potential interactions. For local interpretation, SHAP waterfall plots were generated to break each individual’s predicted BMI into the model’s base value (the average prediction) plus the additive contributions of features, thereby visualizing how specific factors influence the prediction at the individual level.

Analyses were run in Python 3.11.3 (scikit-learn 1.5.2, CatBoost 1.2.7, LightGBM 4.5.0, SHAP 0.46.0). Randomness was controlled by fixing seeds for dataset splitting, outer/inner nested CV splitters, final 5-fold CV, algorithm random seeds. All seed values are listed in Supplementary Table 2. Development was performed in PyCharm 2023.1.3.

At baseline, a total of 2,006 students aged 14–17 were enrolled in the study. During the one-year follow-up period, 98 students were excluded due to incomplete baseline anthropometrics data, 21 students withdrew due to lack of interest or parental refusal, 45 students were excluded due to incomplete follow-up measurements, and 15 students experienced health complications that hindered their continued participation. Consequently, 1,827 students (1,009 males and 818 females) with complete data at both time points were included in the final analysis. Notably, the electronic questionnaire system ensured all items were completed, guaranteeing no missing values in baseline data. Anthropometric measurements were 100% complete in the final cohort. The participant eligibility, follow-up, and analysis process is illustrated in Figure 1. At baseline, the mean BMI was 21.18 ± 3.63 kg/m², and after 1 year, it increased to 21.54 ± 3.59 kg/m². The mean ΔBMI over the one-year period, calculated as the individual-level difference between follow-up and baseline BMI, was 0.36 ± 1.40 kg/m², reflecting a slight upward trend within the cohort. The power analysis confirmed that the final sample size was sufficient to detect medium-sized effects with adequate statistical power, supporting the validity of the univariate analyses. In addition, the school-level ICC for baseline BMI was approximately zero (95% CI: 0.000–0.0037; p > 0.05), indicating negligible clustering. Baseline characteristics are presented in Table 1, comparing included and excluded participants. All |SMD| values were <0.10, indicating good balance of baseline characteristics with negligible differences. The overall loss to follow-up was 4.2%. In prespecified subgroup analyses of attrition rates and RDs versus the reference level, differences were small and imprecise. Age: 14–15 years 4.2% (reference) vs. 16–17 years 4.4%, RD = +0.3 percentage points (95% CI: −1.7–2.2). Gender: female 4.0% (reference) vs. male 4.5%, RD = +0.5 percentage points (95% CI: −1.3–2.3). Baseline BMI category: non-overweight/obese 4.1% (reference) vs. overweight/obese 4.7%, RD = +0.6 percentage points (95% CI: −1.5–2.7). Collectively, subgroup RDs were close to zero and all 95% CIs included zero, providing no evidence of differential attrition. Furthermore, the RTM analysis yielded β = −0.084 (95% CI: −0.122 to −0.047; p < 0.01), indicating evidence of regression to the mean. The coefficient of determination was R² = 0.051 (95% CI: 0.008–0.119), implying that baseline BMI accounts for about 5% (95% CI: 0.8–11.9%) of the variance in the observed change.

Under nested cross-validation, the CB model performed the best, achieving the lowest RMSE, MSE, and MAE, and the highest R², as shown in Table 2. The final hyperparameters selected for all models are reported in Supplementary Table 2. The final CB model included the following predictors: baseline BMI, baseline BMR, level of health literacy, recognize self-weight status correctly, sedentariness duration on weekends, participation in professional sports training, frequency of staying up late, daily sleep duration, frequency of high-calorie food intake, physical activities duration on weekends, post-exercise sensations, satisfaction with body size, family residence location, and on-campus residence. Table 3 summarizes the performance of all final models on the independent test set, reporting RMSE, MSE, MAE, and R² with 95% CIs. The CB model demonstrated the best generalization. In Addition, paired bootstrap RMSE difference (comparator—CB) for MLP 0.080 (95% CI: 0.012–0.151), LGBM 0.124 (95% CI: 0.065–0.185), SVR 0.076 (95% CI: 0.011–0.142), KNN 0.453 (95% CI: 0.340–0.578), and DT 0.171 (95% CI: 0.074–0.268), with all intervals strictly positive, confirming lower RMSE for CB and its significant superiority over the comparator models.

The Breusch–Pagan test on training out-of-fold residuals indicated heteroscedasticity of errors with baseline BMI (p < 0.05), so we applied WLS calibration learned on the training data and then applied to the test set without refitting. After calibration, test-set performance was: RMSE 1.200 (95% CI: 1.101–1.303), MSE 1.440 (95% CI: 1.211–1.697), MAE 0.895 (95% CI: 0.818–0.981) and R² 0.902 (95% CI: 0.882–0.918).

As a trivial baseline model (predicting follow-up BMI equals baseline BMI), performance on the independent test set was MSE 2.212 (95% CI: 1.766–2.767), RMSE 1.487 (95% CI: 1.329–1.663), MAE 1.065 (95% CI: 0.966–1.181) and R² 0.850 (95% CI: 0.808–0.883). Using bootstrap of paired differences, the incremental performance benefit of the CB-based model over the trivial baseline model was: ΔRMSE (baseline − CB) 0.276 (95% CI: 0.159–0.389), ΔMSE 0.752 (95% CI: 0.415–1.115), ΔMAE 0.169 (95% CI: 0.093–0.242) and ΔR² (CB − baseline) 0.051 (95% CI: 0.027–0.078). All intervals exclude zero, indicating that the CB model provides a statistically significant improvement over the baseline.

After excluding baseline BMI, the CB model achieved RMSE 2.497 (95%: CI 2.270–2.728), MSE 6.233 (95% CI: 5.154–7.440), MAE 1.877 (95% CI: 1.709–2.045), and R² 0.576 (95% CI: 0.502–0.644) on the independent test set. Compared with the primary model, performance declined, suggesting that baseline BMI is likely a dominant predictor.

In subgroup performance evaluation, Table 4 presents the stratified metrics (RMSE, MSE, MAE, R²) of the calibrated CB model, with paired-bootstrap between-group differences (Δ with 95% CIs). In permutation tests, p-values for all performance metrics by gender and by age group were >0.05. For baseline BMI categories, p < 0.05. These results indicate that performance differences across gender and age groups were small, with no statistically detectable heterogeneity. Statistically significant differences were observed across baseline BMI categories.

In error distribution analysis, the calibrated CB model showed near-zero overall bias, mean error = 0.03 ± 1.21, with central and tail dispersion MAD = 0.66, IQR|e| = 1.05, P90|e| = 1.98, P95|e| = 2.71. By gender, males (−0.03 ± 1.20) and females (0.10 ± 1.22) were similar, indicating no material heterogeneity by gender. By age, 14–15 years (−0.01 ± 1.23) and 16–17 years (0.08 ± 1.17) were only mildly different, indicating no material heterogeneity by age. By baseline BMI category, errors in the overweight/obesity group (0.24 ± 1.65) were larger and more right-shifted than in the normal group (−0.06 ± 0.98).

For overall error visualization, Figure 2 shows the predicted-versus-observed scatter with a smoothing line, indicating overall fit with only slight departures at the extremes of the prediction range. Figure 3 displays the Bland–Altman plot with a near-zero mean bias and approximately symmetric limits of agreement (−2.33 to 2.38), and the data points show no systematic drift with the mean, suggesting negligible bias.

Population-level feature importance was assessed using SHAP global interpretation applied to the CB-based model. To focus on modifiable features for intervention insights, non-modifiable features such as baseline BMI, gender, family residence location, among others, were excluded from SHAP visual analysis. Figure 4 displays the SHAP summary plot, ranking modifiable features by their mean absolute SHAP values, which represent their average contribution to BMI prediction across the entire population. The plot also visualizes the distribution of SHAP values, with feature color indicating feature values (red: high, blue: low). Features with positive SHAP values contribute positively to BMI prediction, while negative values indicate a decreasing effect. The modifiable features ranked in descending order were: level of health literacy, recognize self-weight status correctly, sedentariness duration on weekends, participation in professional sports training, frequency of staying up late, daily sleep duration, frequency of high-calorie foods intake, physical activities duration on weekends.

The SHAP interaction analyses are presented in Figure 5. The heatmap shows the mean absolute SHAP interaction value for each feature pair. Color intensity encodes interaction strength, with darker/warmer colors indicating stronger interactions and lighter/cooler colors indicating weaker ones. Diagonal cells approximate main effects, while off-diagonal cells reflect pairwise interactions. As shown in Figure 5, most feature pairs exhibit near-zero interaction values, indicating predominance of main effects and no strong pairwise interactions.

Figure 6 shows SHAP dependence plots for the top four features ranked by global importance (mean |SHAP|). The x-axis shows feature values and the y-axis shows SHAP values. Point color indicates the strongest interacting feature, and vertical dispersion reflects potential interactions. Level of health literacy was positive for “lowest/lower-middle” and negative for “upper-middle/highest” (Figure 6A). Recognize self-weight status correctly was negative for “Yes” and positive for “No” (Figure 6B). Sedentariness duration on weekends >9 h/day showed the largest positive SHAP values, whereas 3–9 h/day was negative to mildly negative (Figure 6C). Participation in professional sports training was negative for “Yes” and near-zero to mildly positive for “No” (Figure 6D). Overall, the plots show smooth directional trends with limited vertical spread, indicating predominance of main effects and modest interactions.

SHAP values quantify each feature’s association with the model’s prediction, providing a detailed view of how individual modifiable features relate to the model’s predicted BMI values for specific samples. SHAP waterfall plots for two specific samples are shown in Figure 7. The red bars represent positive contributions, indicating an increase in the predicted outcome, while the blue bars represent negative contributions, indicating a decrease in the predicted outcome. The red bars indicate positive contributions (increasing the predicted BMI) and the blue bars indicate negative contributions (decreasing the predicted BMI). The bar length and its numeric label reflect the effect magnitude in kg/m², and longer bars denote larger increases or decreases in the predicted BMI. As shown in Figure 7A, the predicted BMI for this individual is 21.00 kg/m². Negative contributions dominate and reduce the overall prediction. The main negative features and their numerical contributions to the predicted BMI are: PPST (Participation in professional sports training) = Yes (−0.45 kg/m²), SDOW (Sedentariness duration on weekends) = 3–5 h/day (−0.15 kg/m²), and RSWSC (Recognize self-weight status correctly) = Yes (−0.13 kg/m²). Positive contributions are smaller, including DSD (Daily sleep duration) = > 8 h (+0.26 kg/m²) and FHCFI (Frequency of high-calorie foods intake) = Always (+0.14 kg/m²). In Figure 7B, positive contributions predominate, leading to a higher predicted BMI for this individual.