Data were collected through questionnaires from children attending 442 kindergartens in Longgang District, Shenzhen, between May and July 2021. The collected information encompassed basic demographic data, clinical symptoms and signs of GHRS, other symptoms and signs, and dietary behaviors. The inclusion criteria were (1) children enrolled in kindergartens in Longgang District during the survey period, and (2) guardians who were willing and able to complete the survey. The exclusion criterion was missing data on GHRS-related symptoms and signs, rendering confirmation of a GHRS diagnosis impossible.

The variable data evaluated in this study were obtained through an online questionnaire survey completed by guardians based on the children’s actual circumstances. Children’s basic information encompassed 39 items, including kindergarten name, grade level, class, relationship to the informant, gender, date of birth, feeding method during the first 6 months of life, timing of animal protein complementary food introduction, and other relevant information. Clinical symptoms and accompanying signs of GHRS were 19 symptoms and signs observed within the preceding 2 weeks, such as increased palm and sole skin temperature, heat aversion, oral malodor, warm exhaled air from the mouth or nose, and thirst accompanied by a preference for cold drinks. Other symptoms and signs were derived from manifestations occurring within the preceding year and included listlessness, anorexia, sallow complexion, flaccid muscles, a preference for quietness, weak and feeble speech, reticence, excessive sweating, cold hands and feet, aversion to cold, dislike of raw and cold foods, and 51 additional symptoms or signs.

Symptoms and signs related to GHRS diagnosis were all included in the questionnaire. The questionnaire was completed by the children’s legal guardians, who possessed basic literacy skills. Following the questionnaire instructions, guardians provided responses based on their observations and judgments regarding the children’s symptoms and signs. This study used a guardian-administered questionnaire survey and did not implement blinding of investigators or outcome assessors. GHRS diagnosis was conducted using the “Self-Rating Diagnostic Scale for Gastrointestinal Heat Retention Syndrome in Children” developed by our team in prior research, which has demonstrated good reliability and validity [21,22]. The scale comprises 10 items with different weighting coefficients. Individuals with a cumulative score of 133 points or more were classified as GHRS positive, whereas those with scores below this threshold were classified as non-GHRS, as detailed in Table 1.

Box plots were used to detect outliers in continuous variables. Outliers were defined as values above the upper quartile plus 1.5 times the IQR or below the lower quartile minus 1.5 times the IQR (Table S1 in Multimedia Appendix 1). Each outlier was replaced with the nearest boundary value to align with the main data distribution (Figure S1 in Multimedia Appendix 1). To maximize data validity, samples with missing values were excluded before model training. Categorical variables were predefined before model construction, and their data were processed using one-hot encoding. All sociodemographic subgroups were subjected to the same preprocessing procedures for outlier screening, missing value handling, and feature encoding, without differential treatment. Independent correlation analyses were performed to examine statistical associations between variables and GHRS and to assess the feasibility of model development. These analyses did not contribute to feature selection for the final model.

In this study, GHRS diagnosis required the use of assessment scales. To prevent the trained ML models from learning the scale scoring rules directly or indirectly, we initially excluded variables related to GHRS clinical symptoms and accompanying signs (including scale items), retaining 467 other variables (Table S2 in Multimedia Appendix 1). Subsequently, univariate correlation screening and recursive feature elimination (RFE) were used for feature selection to enhance model performance and stability. Univariate correlation screening utilized the Pearson correlation coefficient to calculate the absolute correlation between each feature and the GHRS diagnostic outcome. Features with an absolute correlation coefficient of 0.1 or more were retained, whereas those exhibiting weaker correlations with the target variable were excluded to minimize interference from redundant variables. RFE is a widely used feature selection method in ML that eliminates unimportant features and ultimately identifies the optimal feature combination for achieving peak model performance. During the RFE process, 10-fold stratified cross-validation was applied, and features contributing minimally to model performance were recursively eliminated to derive the optimal feature subset for each model.

According to the events per variable (EPV) rule, the minimum required sample size for the training set was calculated using the following formula:

Minimum required sample size for the training set = number of factors included × EPV/incidence of positive events in the training set

where the number of included factors refers to the number of features in the optimal feature subset obtained after univariate correlation screening and RFE. The positive event incidence rate in the training set represents the probability of a positive GHRS diagnosis. The EPV value was set at 10. According to existing research findings, the prevalence of GHRS among children is approximately 51.88% [17]. Therefore, the minimum required sample size for the training set should be at least 19.27 times the number of features in the optimal feature subset.

Using stratified random sampling, the preprocessed dataset was divided into training and test sets at an 8:2 ratio. Final feature selection was conducted exclusively within the training set and cross-validation process to avoid information leakage from the test set. The test set was used solely for final model performance evaluation and did not participate in feature screening, resampling, hyperparameter tuning, or model selection. This study included samples only from multiple kindergartens in Longgang District, Shenzhen, and did not contain multiregional subsets; therefore, cross-subset heterogeneity was not assessed. After retaining complete cases, and considering that complete-case screening may alter the distribution of the training set, SMOTETomek (Synthetic Minority Oversampling Technique Tomek Links) was prespecified as a conservative resampling step in the training pipeline. Given the nearly balanced ratio of GHRS-positive to GHRS-negative cases, this step was not considered necessary to correct severe class imbalance but was retained as part of the primary analysis pipeline. SMOTETomek, feature standardization, and model training were integrated into an ImbPipeline (Imbalanced Learning Pipeline) and fitted exclusively on the training set and within cross-validation training folds. Validation folds and the test set were used solely for performance evaluation and were not involved in resampling, feature screening, hyperparameter tuning, or model selection, thereby minimizing the risk of information leakage. Seven ML models—Decision Tree, K-Nearest Neighbors, Linear Support Vector Machine (LinearSVC), LR, naive Bayes (NB), Random Forest (RF), and XGBoost—were used to estimate the probability of GHRS positivity. To optimize model performance, hyperparameter tuning was performed using grid search based on 3 rounds of 3-fold cross-validation to determine the final hyperparameters for each model using the optimal feature subset. Finally, because the outputs of the Decision Tree and LinearSVC models do not represent reliable probability estimates, Platt scaling (sigmoid) was applied for probability calibration to all models. Specifically, a calibration function was fitted using 3-fold cross-validation to convert model outputs into calibrated probability values ranging from 0 to 1.

After model construction, internal validation was performed. We conducted 10-fold stratified cross-validation on the training set. In each cross-validation iteration, feature standardization, SMOTETomek resampling, model training, and probability calibration were performed exclusively within the training folds, whereas the validation folds were used solely for performance evaluation. Performance metrics, including the area under the receiver operating characteristic curve (AUC), sensitivity, specificity, and Brier score, were calculated for each fold. The optimal threshold for each fold was determined using the Youden Index. Results were then aggregated across all folds to estimate the average performance on the training set and to derive the global optimal threshold. Using the optimal feature subset, optimal hyperparameters, and global optimal threshold obtained from the training set, the final model was trained. The model was subsequently validated on the test set, and performance metrics were calculated.

Model performance was compared across 3 dimensions: discriminatory power, calibration, and clinical utility. To evaluate discriminatory power, the AUC was primarily used to assess the model’s ability to distinguish between GHRS-positive and GHRS-negative samples. The AUC ranges from 0 to 1, with values closer to 1 indicating better discriminatory performance. Sensitivity and specificity were also calculated at both the fixed threshold (0.5) and the optimal threshold. Sensitivity reflects the model’s ability to identify positive samples, whereas specificity reflects its ability to identify negative samples. For calibration assessment, the Brier score was used to evaluate the accuracy of the model’s probability estimates. The Brier score ranges from 0 to 1, with values closer to 0 indicating better agreement between the predicted probabilities and the observed outcomes. Additionally, a calibration curve was plotted with the model-predicted probability on the x-axis and the observed positive rate on the y-axis to visually assess the model’s calibration performance. To further quantify calibration, the calibration slope and calibration intercept were calculated. The calibration slope typically approaches 1, with values closer to 1 indicating better agreement between the trend of the predicted probabilities and the observed positive rate. The calibration intercept typically approaches 0, with values closer to 0 indicating lower overall prediction bias. For clinical utility evaluation, decision curve analysis (DCA) was used to calculate the net benefit across different threshold probabilities. A higher net benefit indicates greater practical value of the model at a given clinical decision threshold. Decision curves were then plotted to visually compare the net benefits of each model with those of the “intervention for all” and “intervention for none” strategies. Finally, based on the performance of the above evaluation metrics in both the training and test sets, the optimal model was selected according to its overall practical utility.

This study used the SHAP method to interpret the optimal ML model. SHAP is a game theory–based technique that ranks the importance of input features and explains model outputs, thereby addressing the limited interpretability of ML models. By calculating each feature’s contribution to the model output, SHAP provides both global and local interpretability, enhancing model transparency and explainability. In addition, a predicted probability distribution plot for the optimal model was generated to illustrate the agreement between model-predicted probabilities and actual outcomes, thereby supporting the reliability of the model evaluation results.

To facilitate exploratory self-assessment at home, we integrated the optimized final model into a web application developed using the Streamlit framework. After users complete the online questionnaire and enter values for all predictive features, the application calculates and returns the probability that the child meets the study-defined GHRS scale criteria.

The aforementioned ML models were developed using the scikit-learn and Keras libraries in Python (version 3.11; Python Foundation). Continuous variables with normal distributions are presented as mean (SD) and were compared using the unpaired (2-tailed) t test. Categorical variables were represented using one-hot encoding. Independent correlates in the overall cohort were identified through univariate LR analysis, collinearity testing, stepwise Wald regression analysis, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis, and multivariable LR analysis.

This study protocol adhered to the principles of the Declaration of Helsinki and received approval from the Ethics Review Committee of Beijing University of Chinese Medicine, Shenzhen Hospital (Longgang; approval number SZLDH2021LSYM-03). As a cross-sectional study, this research was not registered in a clinical trial registry before data collection. This study was conducted with informed consent obtained from kindergarten administrators and the guardians of participating children. Data collection complied with the Ethical Review Measures for Life Science and Medical Research Involving Humans and the Personal Information Protection Law of the People’s Republic of China. The dataset was rigorously anonymized and deidentified to ensure that no participants could be identified, and no participant identities could be discerned from any images or multimedia appendices associated with this article. Participants received standardized compensation determined by the research team. Detailed compensation terms are not disclosed beyond the participants and members of the research team.

This study was reported in accordance with the TRIPOD+AI (Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis + Artificial Intelligence) guidelines. The complete TRIPOD+AI abstract and full reporting checklists are provided in Multimedia Appendices 2 and 3, respectively.

Following questionnaire distribution, 120,198 questionnaires were returned. Of these, 11,751 were excluded because of missing data on GHRS symptoms and signs, leaving 108,447 questionnaires for analysis. The demographic characteristics of the children are presented in Table 2. As shown in Table 2, significant differences (P<.001) were observed between the GHRS-positive and GHRS-negative groups with respect to gender, grade level, age, height, and body weight, indicating that these demographic variables may contribute to estimating the probability that a child meets the GHRS scale criteria. Details of the study design are presented in Figure 1.

To determine the association between GHRS prevalence and questionnaire-collected variables while preliminarily assessing the feasibility of constructing a GHRS prevalence prediction model, we investigated independent correlates of GHRS in children using a dataset comprising 108,447 samples. After creating dummy variables and applying one-hot encoding to categorical variables, 469 variables were obtained. Variables such as “unclear,” which were designed to improve questionnaire completeness but lacked analytical significance, were excluded, leaving 467 variables. Univariate LR analysis identified 433 potential correlates associated with GHRS prevalence (P<.01; see Table S2 in Multimedia Appendix 1 for the complete list of P values). After collinearity testing, 335 potential correlates were retained (variance inflation factor [VIF]<10). Following stepwise regression with Wald tests, 170 potential correlates were identified (P<.001), and LASSO regression analysis identified 113 potential correlates (coefficient≠0). Finally, multivariable LR analysis identified 59 independent correlates associated with GHRS prevalence, including 10 protective factors and 49 risk factors (Table S2 in Multimedia Appendix 1 and Tables 3 and 4).

Using the dataset comprising 108,447 samples, records with missing values were excluded, yielding a complete-case sample of 49,798. Significant differences were observed in 304 out of 467 (65.1%) variables between the included and excluded pediatric groups (Table S3 in Multimedia Appendix 1). After stratified random sampling at an 8:2 ratio to create the training and test sets, the training set contained 39,838 samples, and the test set contained 9960 samples. Subsequently, the feature selection pipeline was reestablished. In the training set, after univariate correlation screening using the Pearson correlation coefficient, a total of 162 features were retained (Table S4 in Multimedia Appendix 1). Following RFE and a 10×10 cross-validation process, the optimal feature subset for each model was obtained (Table 5).

Using the maximum number of optimal features among all models, the minimum required sample size for the training set was calculated to be 3122. With the training and test sets divided at an 8:2 ratio, the corresponding minimum required test set sample size was 625 cases. Therefore, the total minimum sample size required for model development and validation was 3747 cases.

On the training set, we performed 10×10-fold internal cross-validation and hyperparameter grid search (Table S5 in Multimedia Appendix 1) to construct 7 ML models (Table 6, Figure 2). The results showed that XGBoost achieved the best AUC performance (mean 0.7417, SD 0.0086; 95% CI 0.7352-0.7482). Using 0.5 as the classification threshold, RF demonstrated the highest sensitivity (mean 0.6752, SD 0.0115; 95% CI 0.6666-0.6839), whereas NB exhibited the highest specificity (mean 0.7387, SD 0.0086; 95% CI 0.7323-0.7452). At the optimal threshold determined by the maximum Youden index, NB achieved the highest sensitivity (mean 0.6977, SD 0.0384; 95% CI 0.6688-0.7267), whereas LR demonstrated the highest specificity (mean 0.6783, SD 0.0304; 95% CI 0.6554-0.7013). For the Brier score, XGBoost exhibited the best performance (mean 0.206, SD 0.0031; 95% CI 0.2037-0.2084), with LR ranking second (mean 0.206, SD 0.0030; 95% CI 0.2037-0.2083).

On the test set, we evaluated the discriminatory performance, calibration, and clinical utility of each model (Table 7, Figure 3). Among the discrimination metrics, the highest AUC was achieved by XGBoost; the highest specificity at the 0.5 threshold by RF; the highest sensitivity at both the 0.5 threshold and the optimal threshold by LinearSVC; and the highest specificity at the optimal threshold by NB. DeLong test results indicated statistically significant differences in AUC among the models (Table 8). Among the calibration metrics, the mean Brier scores for XGBoost, LR, and RF were closest to 0.20, with RF exhibiting the narrowest CI for the Brier score. This finding suggests smaller fluctuations in calibration performance and slightly greater stability across different sampling scenarios. Regarding the calibration slope and calibration intercept, LR and XGBoost performed slightly better than RF, indicating that RF exhibited mildly conservative calibration and a tendency to underestimate risk relative to these models. For clinical utility assessment, bootstrapped CIs for the net benefit derived from DCA (Table S6 in Multimedia Appendix 1), together with the decision curves, showed that within the threshold probability range of 0.01-0.90, the net benefit of the RF model exceeded those of both the “treat-all” and “treat-none” strategies. Within the threshold probability range of 0.13-0.85, the net benefit of the RF model was substantially higher than that of the other models, demonstrating superior clinical decision-making utility across the study’s threshold spectrum.

The above results indicate that, among the 7 ML models, the RF model did not outperform all other models across every performance metric. XGBoost and LR demonstrated slight advantages in certain discrimination or calibration metrics. According to the prespecified comprehensive utility-oriented selection criterion, RF was chosen as the primary analysis model because of its narrower CI for the Brier score, acceptable calibration performance, and relatively greater net benefit in DCA, indicating superior overall utility.

To improve the interpretability and credibility of the established prediction model, we used the SHAP method to analyze the contribution of each variable to the model output. We also generated a probability distribution plot for the RF model on the test set to further support the reliability of its prevalence probability predictions.

The SHAP method provides model interpretability through both global and local explanations. Global interpretability characterizes overall model behavior and is illustrated by SHAP summary bar plots, SHAP summary dot plots, and SHAP decision plots, whereas local interpretability is represented by SHAP dependence plots. Given the large number of features and their lengthy names in this study, we created a feature ID-feature name mapping table (Table S7 in Multimedia Appendix 1) to facilitate simplified feature representation. For global interpretation, SHAP summary bar plots were used to evaluate feature contributions to the model based on mean absolute SHAP values and to display the top 20 contributing features in descending order (Figure 4A). SHAP summary dot plots were used to visually demonstrate the direction and magnitude of each feature’s influence on model outputs (Figure 4B). SHAP decision plots illustrated the decision-making process for the first 100 samples (Figure 4C). The 5 most influential features were Feature19 (always quiet: none), Feature20 (speak in a low voice: none), Feature17 (soft muscles, not as strong as other children: rarely), Feature18 (soft muscles, not as strong as other children: sometimes), and Feature13 (yellow complexion, not as ruddy as other children: none). Higher values for Feature19 (always quiet: none) and Feature20 (speak in a low voice: none), together with lower values for Feature17 (soft muscles, not as strong as other children: rarely), Feature18 (soft muscles, not as strong as other children: sometimes), and Feature13 (yellow complexion, not as ruddy as other children: none), were associated with a higher predicted probability of GHRS in the model output. For local interpretability, SHAP dependence plots were generated for these 5 features to further illustrate their individual effects on model predictions (Figure 4D).

The probability distribution plot was used to assess the agreement between the model-predicted sample distribution and the observed positive rate (Figure 5). Test-set samples were stratified into 3 subgroups according to their predicted probabilities of meeting the study-specific GHRS scale criteria: low probability (0%-30%), moderate probability (30%-60%), and high probability (60%-100%). The low-probability group (n=9960) consisted of 1898 (19.06%) test-set samples, including 311 positive cases, corresponding to an observed positive rate of 16.39%. This observed rate was substantially lower than the upper threshold of 30%, suggesting that the model was conservative within this probability range and could preliminarily identify children with a low likelihood of meeting the study-specific GHRS scale criteria. The moderate-probability group (n=9960) comprised 5302 (53.23%) test-set samples, including 2301 positive cases, corresponding to an observed positive rate of 43.41%. This finding suggests relatively good calibration performance within this probability range. The high-probability group (n=9960) included 2760 (27.71%) samples, of which 2008 were positive cases, corresponding to an observed positive rate of 72.75%. This result indicates that the RF model was able to appropriately stratify children with a high likelihood of meeting the study-specific GHRS scale criteria. Overall, the predicted probabilities of GHRS positivity generated by the RF model were consistent with the observed trend in positive rates, demonstrating satisfactory risk-stratification performance.

To facilitate public use, we integrated the RF model into a web application developed using the Streamlit framework for use in everyday settings. We established a definitive mapping between the 158 features included in the RF model and the questionnaire items used in the Streamlit application. The questionnaire content is consistent with the survey instrument used for data collection in this study, and all items can be understood and completed by a child’s guardian for the purpose of information collection (Table S8 in Multimedia Appendix 1). After the child’s guardian completes the electronic questionnaire consisting of 75 single-choice questions, the application automatically estimates the probability that the child currently meets the GHRS diagnostic criteria, outputs the corresponding disease probability category, and provides recommendations based on the results. This application can be accessed online [23].