Developing ML-based diagnostic tools requires carefully devised strategies for the selection, utilization, evaluation, and optimization of the ML models and tools. Accordingly, the type of data to be used, computational resources available for model building and deployment, expected model performances, undesired model behaviors such as overfitting, model optimization strategies to be used, and the nature of the disease to be diagnosed are the major parameters that determine the ML methods and strategies to be used. For small-sized datasets, support vector machines (SVMs), decision trees, and random forests (RF) are the preferred ML algorithms, while the convolutional neural network (CNN)–based DL models can also be used [14,15]. However, the tree-based model extreme gradient boosting (XGBoost) is the most suitable model for tabular data since it requires fewer hyperparameter tuning tasks compared to the DL-based models, as stated in [16]. Generally, the gradient-boosted tree models of XGBoost, light gradient boosting machine (LightGBM), and CatBoost are well-suited for tabular data, although they can be challenged by model complexity [17], higher computational resource requirements such as GPUs for large datasets [18], and possible issues during real-time deployment on edge devices [19].

There are multiple ML strategies and models for handling dataset-related issues, as it is sometimes challenging to acquire complete datasets that are sufficient to train ML and DL models. For small-sized tabular datasets, traditional ML models and ensemble approaches are the most suitable methods due to their stability and interpretability, as indicated in recent studies. Accordingly, the tree-based ensemble models of RF and extra trees are most often the top-performing models for the classification of diseases in dermatology based on tabular patient metadata, including the situation of having small-sized tabular datasets [20]. In terms of performance reliability, k-nearest neighbor (KNN) and SVM are the most stable models that resulted in superior performance compared to complex neural network–based models for small-sized structured clinical datasets [21,22]. While the boosted tree-based model XGBoost requires a larger number of samples to reach performance stability, the RF model appears to be a more reliable model for small-sized datasets [22]. Regarding missing data, studies suggested that the gradient-boosted tree models of CatBoost, XGBoost, and LightGBM are highly robust to missing data when trained on incomplete patient records. Accordingly, CatBoost (with nan_mode enabled) is the most robust model, as it has an internal mechanism to treat “NaN” values as a separate split direction [23]. The skin disease detection study by Aquil et al [20] highlighted that XGBoost and LightGBM models have the ability of learning “default direction” for missing values during training. Accordingly, if a given symptom (feature) has a null value, the model automatically assigns this feature to the branch that minimizes the loss based on other patients.

Most previous studies conducted for the diagnosis of skin NTDs heavily rely on skin images. Accordingly, the previous studies [23-25] used DL-based methods to apply diagnostic models for skin NTDs, where all these studies used only skin image-based approaches. On the other hand, previous works that implemented ML methods for the diagnosis of skin NTDs using only structured clinical data were difficult to find, except 2 studies that used tabular patient metadata along with skin image data. Accordingly, the study by Barbieri et al [11] used structured clinical data of leprosy patients to integrate with skin lesion images and used ML algorithms that include RF and XGBoost for handling the clinical patient data in their multimodal data fusion study. Similarly, Achary et al [12] used tabular clinical metadata to apply multimodal data fusion by combining the metadata with skin images of patients using generative models for text data augmentation. In both studies, the tabular clinical data of patients are not used independently for skin NTD diagnosis; instead, they are merged with image data based on multimodal data fusion methods.

Overall, the results of the literature review underscored that structured clinical metadata of patients are not fully used, let alone the independent utilization, demonstrating the development of ML-based diagnostic models for skin NTDs based on tabular clinical patient data, indicating a clear research gap. One of the biggest challenges in building AI-based diagnostic tools for skin NTDs is the lack of large-scale, complete, and properly annotated datasets, as affirmed by the study by Achary et al [12]—the study used generative models to generate textual data to fill the gap of data scarcity. Conversely, the selection and use of specific ML strategies require further investigation, given that skin NTD data are naturally scarce resources, including limited previous studies for comparative analysis of methods. All these factors clearly mark research-wide gaps that forced the authors of this study to conduct this study.

For this study, we created a new hand-crafted dataset using the data collected from patients with skin NTDs living in Gacho Baba district of the Gamo zone, southwest of Ethiopia, which is one of the remote areas representing underserved communities in terms of health care services. The data were collected by a team of researchers from the Collaborative Research and Training Center for NTDs, College of Medicine and Health Sciences of Arba Minch University, in a project-based research conducted for the assessment of skin NTD burden. For this study, the original (first-hand) data were acquired through institutional collaboration (acknowledgments section) and used to achieve the objective of this study.

The entire data collection was initially conducted in a large-scale community-based mass drug administration (MDA) campaign in the specified community based on on-site clinical diagnostic procedures. At the MDA site, however, extended medical resources like laboratory facilities, specialized medical imaging, or any advanced medical infrastructure were not available. In such settings, the initial clinical diagnostic procedures were performed by trained frontline health care workers (mostly middle-level clinicians and nurses) who were actively involved in the day-to-day clinical diagnostics of patients around the area, and a professional dermatologist from Arba Minch General Hospital who participated during the MDA campaign. Therefore, the on-site diagnostics were performed entirely based on visual clinical examination procedures by the frontline clinicians, with strict diagnostic confirmation by the professional dermatologist.

As the data collection involved actual field-based diagnoses, the establishment of initial ground truth was challenged by limitations related to field logistics. To overcome this challenge, a hierarchical diagnostic verification protocol was applied by the data collection team to clinically confirm the diagnostics. Hence, each of the initial diagnostic procedures underwent a mandatory clinical review by a professional dermatologist from Arba Minch General Hospital to confirm the initial diagnostics. While this procedure confirms the correct classification of the skin NTD diagnosed based on the signs and symptoms presented by the patients, it also helped in creating higher-level diagnostic fidelity. Overall, the ground truth we used for the classification of the skin NTDs for this study complies with the WHO-standardized clinical case definitions [26]. After all these clinical and ethical procedures were completed, the authors of this study used the clinically validated data for the ML-based research in this study after acquiring the proper authorization through ethical clearance.

Overall, the data collection was conducted based on a community-based screening approach for assessing the burden of skin NTDs, which included on-site examination, treatment, and registration of patient information. To achieve this, the data collection team used a well-prepared data collection and patient registration form, which was prepared by consulting professional dermatologists from Arba Minch General Hospital. The data collection involved 4029 participants during the initial screening and diagnosis, targeting the identification and diagnosis of 8 skin NTDs that are endemic to Ethiopia, which include mycetoma, leprosy, CL, lymphatic filariasis, podoconiosis (podo), onchocerciasis, tungiasis, and head lice. Finally, according to the data collection report, the specified community was highly affected by 3 of the skin NTDs, namely scabies, tungiasis, and podoconiosis. Thus, we selected these 3 diseases for our study based on the prevalence and availability of sufficient data. The other 5 mentioned diseases have very little (or no data at all) in the specified community at the time of data collection; hence, they are not used in the final dataset.

Generally, collecting large-scale skin NTD data from wider geographic areas (like Ethiopia) with complete representation of all endemic diseases is a complex and highly resource-intensive process, which is almost impossible to achieve in a single research project. Hence, our study is limited to a single geographic location in the Gacho Baba district, southwest of Ethiopia. This confinement to the specific location limited the number of endemic diseases to only 3 skin NTDs (scabies, tungiasis, and podoconiosis) out of the 8 overall endemic diseases, which forced us to exclude other skin NTDs like buruli ulcer, leishmaniasis, leprosy, and mycetoma. However, as a preliminary pilot study, we used the acquired skin NTD data to build a benchmark skin NTD diagnostic model, with the potential of extending the disease representation through further efforts.

During initial collection, the patient data collected for each of the 3 diseases were recorded separately using individual Microsoft Excel files for each disease. Hence, we initially created 3 separate datasets, which represent (1) scabies, which consists of 955 instances having 10 features (6 demographic and 4 disease-specific); (2) tungiasis, consisting of 474 instances having 8 features (6 demographic and 2 disease-specific); and (3) podoconiosis, comprising 66 instances having 11 features (6 demographic and 5 disease-specific), all before data preprocessing. In each of these 3 separate datasets, all disease-specific features (symptoms) are represented in binary categorical form using either “1” or “0” to signify the presence or absence of a given disease-specific symptom in a patient, which were applied after clinical confirmation of all diagnoses. For instance, “slowenlarging_swelling(feet&legs)” is a symptom specifically assessed while diagnosing podoconiosis, and all 66 podo patients were assessed if they manifest this particular symptom. Hence, after final clinical confirmation of the diagnosis, this feature is encoded as either “1” if a given patient showed this symptom or “0” otherwise, where a similar pattern was used for all symptoms of the 3 diseases.

Finally, we handcrafted our initial dataset (named it IDS) by merging the 3 individual disease-specific datasets. Originally, IDS comprises a total of 1495 instances with 17 trainable features (6 demographic and 11 diagnostic) and the target variable “disease_diagnosed” defining the class labels. After initial data preprocessing, IDS is transformed into a new structure having 1495 instances with 28 trainable features (17 demographic and 11 diagnostic) plus the target variable. However, the merger of the 3 separate disease-specific datasets introduced a different challenge due to the disease-specific symptoms, as not all symptoms are applicable to all diseases—a specific symptom that can be used to diagnose scabies is not applicable to the other 2 diseases. This created 3 blocks of missing values in IDS, showing a structural missingness issue [27]. Hence, using ML-based operations, we created 2 additional datasets through transformation: first, we created a final dataset (FDS), having 28 trainable features (17 demographic and 11 disease-specific features); next, we created a different dataset by applying feature engineering on the final dataset (FEFDS), having 38 trainable features (17 demographic and 11 disease-specific features, including 10 symptom applicability indicator flags introduced after feature engineering). We created all these 3 different datasets in 3 separate CSV files to address structural missingness, ensure reproducibility, and maintain dataset integrity across different experimental settings, thereby allowing us to establish a consistent training and evaluation pipeline.

Multimedia Appendix 1 presents the overall distribution of instances (n=1495), disease classes, and missing values in IDS. As shown in Figure S1 of Multimedia Appendix 1, the scabies disease class represents the largest proportion with a total of 955 (63.9%) instances, while tungiasis represents the second largest proportion with 474 (31.7%) instances, and podoconiosis has only 66 (4.4%) instances. Gender-wise, 66.6% (995/1495) of the patients are male, while 33.4% (500/1495) are female. These statistics clearly highlight that the dataset exhibits a severe class imbalance among the 3 disease classes, having an extremely distorted distribution with an overall imbalance ratio of 14.5 between scabies (majority class) and podoconiosis (minority class).

Correspondingly, Figure S2 in Multimedia Appendix 1 depicts the statistical summary of missing values in the dataset. As stated, our new dataset is created by merging 3 separate disease-specific datasets, each of which has a unique set of clinically relevant symptoms that are not applicable to the other diseases. For instance, the symptom “PPSEs” (which represents “papule, pustule, scratches, and excoriations”) is assessed only for scabies patients, where this field is structurally empty for podoconiosis and tungiasis patients. Ultimately, this created 3 blocks of missing values that define the case of structural missingness or missing by design [27], where our new dataset contains a total of 11,347 such structural null values. The missing values in our dataset are not random omissions created as a result of the random absence or exclusion of symptoms that are called “missing completely at random (MCAR)” or “missing at random (MAR)” [28]. Overall, as depicted in Multimedia Appendix 1, the scabies class has the highest number of null values, with 58.9% (6685/11,347) of the total number of null values (having 7 null values per each of the 955 rows), while tungiasis has 37.6% (4266/11,347), and podoconiosis has 3.5% (396/11,347) since it has the least number of instances.

Table 1 presents the overall summary of our first dataset, IDS, including its overall statistical summary, including its overall statistical summary of class distribution, feature categories, and the initial feature sets of IDS.

As shown, our new dataset (IDS) primarily contains two major categories of patient data: (1) demographics, which include 2 numerical features (age and weight) and 3 categorical features (gender, marital status, educational level, and occupation) of the patient information, and (2) clinical data, which comprises 11 disease-specific symptoms from the 3 classes. Class-wise, podoconiosis has 5 specific features (disease-specific symptoms) that are not applicable to the other 2 diseases, and the scabies class has 4 unique symptoms applicable only to scabies, while there are only 2 tungiasis-specific symptoms (features). Overall, our dataset has 2 features representing numerical data (age and weight), while the remaining features are categorical, which in turn include 3 ordinal categorical features, representing sex, educational level, and one disease-specific feature (the jigger infestation classification feature of tungiasis); 2 nominal categorical features (marital status and occupation), representing the unordered categorical data that include marital status occupation; and 10 binary categorical features, represented as “1” or “0” indicating the presence or absence of a given disease-specific symptom. Table S1 in Multimedia Appendix 2 presents the complete dataset information.

In this study, we used IDS in 3 different settings. First, IDS is used with minimal preprocessing operations for developing baseline models. This dataset is used to mimic the original real-world scenario, which is sometimes characterized by incomplete datasets having multiple missing values. As a base dataset, the majority of the essential data preprocessing operations are applied to IDS, while we applied only normalization for FDS and FEFDS, as these datasets were derived from the feature-encoded IDS. Overall, for data preprocessing, we applied data encoding, which mainly included data encoding using label encoding and one-hot encoding methods, including normalization.

For encoding the ordinal categorical features (“sex,” “educational_level,” “JI_classification(#Jig-Inf),” including the target variable “disease_diagnosed”), the Label Encoder method is used. As the remaining 12 features represent nominal categorical data, which include “marital_status,” “occupation,” and 10 disease-specific symptoms, we applied the one-hot encoding method to transform the data under these features. Data normalization was the other preprocessing operation we performed. Accordingly, the numerical features of age (with minimum 2 and maximum 77) and weight (with minimum 8 and maximum 80) are normalized using the “standard scaler” technique to bring these features on the same scale with the other features. This allowed us to maintain a standard normal distribution with zero mean and unit variance within the dataset [29].

As stated, our new dataset (IDS) contains a total of 11,347 structural null values, which are not random omissions created as a result of the random absence or exclusion of symptoms that are called MCAR or MAR [28]. Over reliance on conventional imputation techniques such as simple imputation or using other methods like mean, median, or multivariate imputation by chained equation (MICE) would introduce biases that could potentially cause both statistical and clinical issues [27,28], having contradictory implications: positively, imputing missing values with “0” values is logically correct, easy to apply, and does not cause any loss of information; while negatively, imputing the missing values with “0” values will introduce hidden danger if a patient develops a secondary (additional) symptom due to a different cause. In our case, for instance, if our models are trained using the fact that “if the values of the ‘Slowenlarging_swelling(feet&legs)’ column are non-zero then, the disease cannot be scabies or tungiasis.” However, such restrictions will affect the generalizability of our models when dealing with patients having mixed symptoms if a patient also has swelling on the feet due to different causes. Therefore, this requires careful analysis and selection of strategies to be used for the imputation, which adds robustness to our models.

Therefore, we created the final dataset by applying a hybrid data encoding strategy: (1) we encoded all symptoms as binary, putting 1 if a given symptom is present (applicable) to a particular disease and putting 0 otherwise; and (2) we introduced a “symptom missingness indicator” parameters, which flags 1 if a given symptom applies to a particular disease and 0 if that symptom does not apply. Therefore, we created a new column called “DISEASE-SYMPTOM_applicable” as a missingness indicator, indicating if that “SYMPTOM” is applicable to the “DISEASE” or not. For instance, since “PPSEs” is a symptom applicable only to scabies and not applicable to podoconiosis or tungiasis, the “PPSEs_applicable” column for all scabies patients is encoded as “1” whether a patient has “PPSEs” or not, while this column is encoded as “0” for podoconiosis and tungiasis patients. The repeated application of this method across 10 disease-specific features adds 10 additional confirmatory columns, which expanded the number of trainable columns of the dataset to 38.

This study is conducted based on the stratified k-fold cross-validation (CV) method to evaluate possible performance variations resulting from changes in the dataset and ML models and methods. The hold-out method is excluded from the experimental setup due to high-level overfitting recorded during preliminary experiments. Therefore, we conducted extensive experiments using our new dataset in 3 structures (IDS, FDS, and FEFDS) with an 80/20 train-test split based on the stratified 5-fold CV method. We designed 4 separate experimental settings that include (1) baseline training, where IDS (the original dataset) is used with minimal data preprocessing (only data encoding); (2) handling structural missingness, which includes 2 separate training experiments conducted for handling the structural missing values—the first experiment includes training using FDS that applies the simple imputer method, while the next experiment was conducted on FEFDS that applies the symptom applicability indicator flags created after feature engineering—each of these two experiments exactly reflects the other except for the datasets used, aiming to comparatively evaluate performance impacts between the two methods; (3) conditional class weighting; and (4) the hybrid approach, the experiment that harmoniously combines class weighting and the robust dual validation methods, all applied on the dataset created using feature engineering (experiment 3).

To address the severe class imbalance, we applied a class weighting method. Due to the difference in the internal nature of one of our selected models in accepting the “sample_weight” parameter of the scikit-learn library [30], we conditionally applied the class weighting method. Therefore, for all 7 models that inherently accept the “sample_weight” parameter, these sample weights are automatically computed using the built-in method with the “class_weight” option set to “balanced,” then directly passed to their fit methods. Conversely, as the multilayer perceptron (MLP)–based model does not accept the “sample_weight” parameter, we implemented a fallback strategy based on resampling to replicate the instances of the minority class. In this setting, instances from the training set are resampled with replacement based on probabilities equal to the sample weights computed at the beginning of the class weighting method. We applied this method to facilitate the consistent application of the class weighting method across all the selected models for a complete and fair performance comparison.

Given the size and distribution of our dataset, using a single k-fold (5-fold) split might not help to critically assess the stability and reliability of the models, where the repeated stratified k-fold CV method plays vital roles in overcoming performance instability issues [29,31]. Therefore, to establish robust and unbiased performance estimates, we harmoniously applied the repeated and nested stratified CV methods.

This study uses datasets having 1495 instances with only 38 trainable features, where all the features are mostly linear, as shown by the principal component analysis (PCA) projection in Figure 1.

As shown, the features of the 3 disease classes are highly linear, which shows the tendency of the features to be linearly separable. Hence, given the PCA projection analysis and the datasets to be used, 6 parameters are used to select the models. The parameters used include the model’s ability in handling small-sized datasets, the nature of the data used (mostly linear), the robustness of the model in handling higher-level class imbalances, the robustness of the model in handling data with missing values, higher model performance expected, and model simplicity and interpretability for future deployment of the models on portable devices. Therefore, based on the extended analysis of size, complexity, and overall nature of our datasets, ML models, and the intent of the study, we selected 8 ML models: 3 nontree models, namely naïve Bayes (NB), SVM (with linear kernel), and a simple CNN-based model using MLP; and 5 tree-based models, which include RF (ensemble tree-based model) and 4 boosted tree models of CatBoost, AdaBoost, LightGBM, and XGBoost.

For this study, we applied a multifaceted evaluation strategy, given the nature of our new dataset. First, the study entirely relies on the stratified k-fold CV method, as it allowed us to achieve more robust and lower-variance estimates. Hence, we applied a 5-fold (k=5) stratified k-fold CV method, as the hold-out method resulted in total overfitting during preliminary experimental trainings. For this study, using the standard accuracy might be misleading given the small-sized dataset used, the simplicity or complexity of models with respect to the dataset used, and the tendency of models to favor the majority class due to the extreme class imbalance [34]. To address this, we used “balanced accuracy,” which is a metric used to reliably evaluate models by giving equal weights for all classes by averaging sensitivity (recall) and specificity (ie, rate of correctly predicted negative classes) of the models [28,29].

The study highly focuses on the F₁-score metric to balance between the precision and recall. The G-mean score, which represents the geometric mean of TPR (true positive rate or recall) and TNR (true negative rate or specificity) [35], is used to evaluate the class-wise accuracy of the models on the 3 disease classes. We used the G-mean along with AUPRC (area under the precision-recall curve) as these metrics are not affected by class imbalances [36], where the AUPRC metric is used to robustly assess the models’ predictive performance on the minority class. To evaluate the models’ class-wise separation performance, the ROC AUC metrics are computed based on the one-vs-one and one-vs-rest classification approaches. Apart from the macro average scores, the class-specific performance of each model is also evaluated to assess how well each model performs on each class, given the severe class imbalance. Hence, the class-specific precision (positive predictive value), recall (sensitivity), specificity (true negative rate), F₁-score, and NPV (negative predictive value) are used to individually assess the models on each class. Confusion matrices are heavily used to compute and visualize the class-specific scores for each model across all experimental settings. Further, we used the feature importance (FI) scores of each model to evaluate model performance and stability, monitor overfitting due to imbalance, include clinical interpretability, and identify noise and redundant features in the dataset. To apply the FI evaluation, we implemented 3 different but complementary techniques that include using built-in importance for the tree-based models, absolute coefficient magnitude for the linear model (SVM), and permutation importance for MLP and NB. Finally, based on the detailed evaluation of the models’ performance, we conducted deeper analysis using descriptive statistical metrics, which include computing average or mean (μ) values, computing and evaluating numerical variations or differences (Δ) of models’ performance, and analyzing performance stability of the models using SD or σ, across the different experimental settings.

The second experiment that applied the simple imputation method for the structural missing values introduced the first performance variation, allowing the 4 models (NB, SVM, MLP, and AdaBoost) to have full performance scores. The third experiment brought the next variations, creating 3 divergent performance score groups after the application of disease-specific symptom applicability indicator flags. This experiment resulted in substantial performance changes, particularly to MLP, as evidenced by the model’s absolute improvements of Δ+0.097 in balanced accuracy and macro recall, Δ+0.560 (σ=0.396) in podo precision, Δ+0.460 (σ=0.325) in podo F₁-score, Δ+0.295 (σ=0.209) in tungiasis precision, and Δ+0.173 (σ=0.122) in tungiasis F₁-score. According to the results, all improvements, particularly the class-specific scores, were achieved at a slight expense of the scabies class recall Δ−0.005 (σ=−0.004) and F₁-score Δ−0.002 (σ=−0.001). All these improved scores of MLP underscore that the feature engineering technique was significantly effective in helping the model identify the minority classes. The confusion matrices in Figures 2A and 2B depict all these changes shown by MLP across the 2 experiments. The other 6 models showed all perfect per-class scores with zero misclassification.

Unlike MLP, the generative model NB exhibited overall performance declines of Δ−0.026 in balanced accuracy and macro recall, Δ−0.015 in macro F₁-score, including precision (Δ−0.004) and G-mean (Δ−0.014), including the per-class scores. The remaining 6 models maintained the perfect scores across all metrics. The sample weighting method further improved the performance of MLP with an absolute increase of Δ+0.003 in balanced accuracy, recall, and G-mean, Δ+0.026 in macro F₁-score, and Δ+0.044 in precision, while the remaining 7 models remained consistent.

Overall, across all the experimental phases, 4 models (SVM, RF, CatBoost, and AdaBoost) showed absolute stability of performance, with no performance variation (σ=0.0), while the other 4 models (NB, MLP, XGBoost, and LightGBM) exhibited slight performance variations. Table 12 presents an analysis of models’ performance variance across the 6 major evaluation metrics throughout the 5 experiments. As shown, MLP has the highest overall performance variance in its balanced accuracy and macro recall (σ=0.043), macro F₁-score and precision (σ=0.099), including G-mean (σ=0.03) and macro AUPRC (σ=0.001). Conversely, NB exhibited dropping patterns in balanced accuracy and macro recall (σ=−0.014), macro F₁-score (σ=−0.008), and precision, G-mean, and AUPRC (σ=−0.002, −0.007, and −0.015, respectively). These overall declines were introduced due to the introduction of symptom applicability indicator flags, which created a set of highly correlated features, leading to incompatibility with the model’s probabilistic assumption [39].

According to the results of the nested CV method (Table 12), LightGBM and XGBoost exceptionally exhibited marginal performance drops. These 2 models dropped their performance by Δ−0.007 (σ=−0.003) in balanced accuracy and macro recall, Δ−0.004 (σ=−0.002) in F₁-score and G-mean, while showing no change in precision and AUPRC.

Similar patterns were seen in the class-wise performance of these 4 models, as shown in Table 13.

MLP has the highest per-class performance variance with σ=0.33, σ=0.26, σ=0.49, and σ=0.39 in podo-precision, F₁-score, specificity, and recall, respectively. These higher variations are caused by the lowest per-class scores achieved during the third experimental training that demonstrates the structural missing values using missingness indicator features. NB exhibited smaller variance in tungiasis precision (σ=−0.006), podo recall (σ=−0.04), podo-F₁-score (σ=−0.023), tungiasis-F₁-score, and specificity (σ=−0.003). However, the 6 models maintained a totally stable per-class performance (σ=0).

While the complete FI scores of the 38 features across 8 models are presented in Table S3 of Multimedia Appendix 2, in this section, we present an overall summary of the FI scores. Accordingly, the FI evaluation results revealed that 4 models (NB, SVM, MLP, and RF) used 37 features out of the 38 trainable features, as shown in Figure 3A.

Regarding individual features used, Figure 3B shows the top 10 features according to their use by the 8 models across all experiments. Accordingly, the jigger infestation classification feature (JI_classification(#jig-inf)) appears to be the top feature having the highest mean FI value. This feature appears in 6 models in their top 5 models list. Model-wise, we conducted deeper analyses of the feature relevance evaluation scores, presented in Table 14, that reveal a diverse set of feature distribution, utilization, and selectivity patterns among models.

As shown, all the 5 tree-based models exhibited a right-skewed FI distribution, which signifies that smaller subsets of the entire feature set are highly relevant, with multiple remaining features being least or nonrelevant [40]. This demonstrated that these 5 tree-based models are highly selective in using features, as evidenced by the feature preference and sparsity ratio scores. Conversely, NB and MLP exhibited a left-skewed FI distribution pattern, highlighting that these models equally designated the majority of the features as highly relevant—MLP assigned the same FI score of 0.665 for 27 features, whereas NB assigned a uniform 1.0 FI score for 36 features. These models have higher nondiscriminative tendencies to identify and ignore less relevant features, as shown by the higher feature preference and lower sparsity ratio scores (Table 14). Finally, SVM is the only model that demonstrated an approximately balanced (symmetric) FI score distribution, as shown by its lower sparsity and information decay scores. Overall, the model-specific FI distributions displayed 3 different distribution trends, which is confirmed by the aggregate summary of FI distribution results (mean 0.311, SD 0.390)—verifying that the overall distribution of the skin NTD diagnostic features is right-skewed. Given the higher mean (μ=0.311) having an SD that exceeds the mean FI score (σ=0.390), which creates higher aggregate variance factors, the right-skewed distribution is naturally expected to be reflected [40].

In this study, we crafted a new structured skin NTD clinical dataset, confronted by multifaceted data-related challenges (sample size limitation, severe class imbalance, and lack of full disease representation). We used our novel dataset in 3 different structures with the aim of experimentally analyzing, comparing, and identifying the best-performing ML model for the diagnosis of skin NTDs based on tabular patient data. In achieving the objectives, the study recorded multiple findings with different implications. The use of the initial dataset (having the structural missing values) clearly demonstrated the resiliency of the 4 tree-based models (RF, CatBoost, LightGBM, and XGBoost) to null values due to their built-in mechanisms of handling missing values, which include considering the missing values as a separate group [41,42]. These results underlined the effectiveness of these tree-based architectures to establish robust real-time diagnostic tools, including the scenarios of having incomplete diagnostic data.

Given the size and distribution of the new dataset, our study underlined that the selection and utilization of sound ML methods is the key to beat expected generalizability challenges, as demonstrated by the progressive performance gains across experiments. As an initial experimental finding, the Δ+0.097 (10%) and Δ+0.204 (24%) macro recall and F₁-score boosts of MLP validate the adoption of feature engineering as one big initial strategy in addressing and demonstrating our methodological success. Another important finding highlighted that some models, specifically LightGBM and XGBoost, exhibited unreliable performance due to predictive bias, a situation where the standard validation methods (the test CV performance in our case) were pretentiously presenting biased prediction estimates [38,39]. Such performance traps were exposed by the nested CV scores, proving the robustness of the dual CV-based approach in establishing unbiased disease prediction estimates, given the highly constrained dataset used in this study.

Overall, this study highlighted that stable model performances with higher predictive accuracies are achievable through a synergistic integration of multiple ML strategies working in harmony, instead of individual and isolated methods—first, establishing a reliable data preprocessing pipeline followed by a robust validation pipeline incorporating a synchronized implementation of class, all applied on top of the robust data preprocessing pipeline. However, as all these results are achieved against all the odds of data-related constraints, it is mandatory to have further investigations with extended resources to maintain the high achievements.

Based on the overall test performance scores, the majority of the models (except NB and MLP) exhibited consistent performance scores across all metrics. While desirable, perfect scores are not always good signs, given the dataset constraints of this study. Specifically, the perfect scores achieved by the highly feature-selective tree-based models, such as XGBoost, that used only 3 features for disease prediction, raise critical performance concerns. Positively, the perfect scores indicate the maximum disease discriminative power of the model. However, given the fact that the model makes disease predictions based on only 3 features out of 38, and these scores showed declines during the nested CV training, these perfect scores have negative implications, indicating that these scores are fragile and misleading. Such scores indicate that the models are facing the problem of over-optimism in their overall test performance. The use of such models for the diagnosis of diseases using only 3 diagnostic features (3 different symptoms of 3 different diseases) is highly risky and clinically infeasible, as the model would fail if one symptom is missing or misrepresented.

According to the FI analysis reports, the diagnostic features have the highest importance for the classification of the skin NTDs and are most widely used by all models, while the demographic features are the least frequently used set of features for the classification of the skin NTDs. These results highlight that models using almost the whole feature set, such as random forest, achieved 100% performance scores across all metrics, while LightGBM (8 features) and XGBoost using only 3 features scored the same 100% performance. These tendencies of the models in using the least proportions of the trainable features could suggest that the models are applying complex architectures and are suffering from overfitting, given the limited dimensionality of our dataset used. Conversely, the tendencies of models to use the highest proportion of the trainable features could also point to the same problem, though such challenges can be alleviated through the use of ML-based strategies, as in the case of this study. Therefore, all these factors suggest that tradeoffs need to be made for the final model selection. Overall performance analysis and FI scores revealed that the 8 models exhibited different behaviors, forming 3 groups of different model characteristics.

Four models (NB, SVM, MLP, and RF) generally fall under this group due to the behavior they commonly share by equally ranking 37 features, while each model has specific determining characteristics. Accordingly, NB showed a different (outlier) behavior due to its invariable selection character that considers all features are equally important, having no interdependency or highly correlated relationship. This high-level feature-inseparability raises an issue, as it is not feasible for 2 or more features to be equally important in diagnosing skin NTDs. For instance, according to the results of our experimental results, a given value of the “educational level” of all patients has never been equally important as the diagnostic features, as confirmed by the dominant diagnostic feature “JI_classification(#Jig-Inf)” and its applicability flag. Further, the introduction of the symptom applicability flags created performance distortions. Both these factors diluted the reliability of the NB models, making them unable to differentiate between symptoms. This, in turn, results in a model inability to differentiate between diseases, which led to perfect scores, where such models cannot be reliable candidates for clinical diagnostics. MLP is the other all-inclusive model with a slightly different feature ranking behavior compared to NB, where 27 features (14 diagnostic and 13 demographic) showed an equal feature rank of 0.665. This depicts invariable ranks (σ=0) among the 27 features, highlighting that the model failed to effectively differentiate disease signal from noise, which also proved that this neural network-based model showed overfitting due to the model’s complex representational ability, given our small-sized tabular dataset. The other 2 models (SVM and RF) showed similar patterns of feature ranking scores with increasing feature sparsity, where the variance in FI scores increased among features.

Results showed that 3 of the tree-based models (XGBoost, LightGBM, and AdaBoost) exhibited a “lazy learning” pattern by using a limited number of features (3, 8, and 11 features), while all these 3 models achieved perfect scores. Especially with XGBoost showing absolute feature parsimony (showing high feature selectivity), the perfect scores using only 3 diagnostic features highlight a potential tendency of following shortcut paths by memorizing the feature patterns with the highest performance. Given the dataset constraints and the complexity of these models, this can underscore the plausibility of this shortcut-based (lazy) learning tendency by these 3 models, as exceptionally exposed by the overall performance declines recorded during the nested CV experiment.

After the extensive analysis of the models’ performance, the study concludes by identifying the model that achieved the optimal performance as the final skin NTD diagnostic model based on structured clinical data. According to the analytic results of the feature preference scores, CatBoost demonstrated the average performance standard compared to all models in our study. Statistically, based on the number of features identified as important (27 features, 71%), the distribution of the features (17 diagnostic and 10 demographic), and performance stability across experimental phases including nested CV, CatBoost showed an optimal separation pattern, unlike the other 2 groups of models. Additionally, its overall internal working logic of ordered boosting (to reduce prediction variances) and inherent categorical feature handling [40], including its fast and stable inference capabilities [43], could also play vital roles in helping the model achieve higher and more stable performance based on an optimal set of features. Therefore, we identified and proposed the CatBoost model as the potential benchmark skin NTD diagnostic model.