Machine learning-based classification of HIV viral load suppression in low-resource settings

Study design and setting

This study adheres to the TRIPOD (Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis) statement for prediction model development and validation. A quantitative research approach has been adopted in this study, where machine learning techniques will be employed to predict the viral load status among PLHIV (People Living with Human Immunodeficiency Virus). This research was conducted at the University of Gondar Comprehensive and Specialized Hospital, one of the leading healthcare institutions in Gondar City, Amhara National Regional State, Ethiopia. Gondar City, located about 748 km northwest of Addis Ababa, has a population of 457,938 and is one of the largest hubs for health service delivery and research. This ART clinic was established in March 2005 and provides free ART services to over 7 million people in Gondar province and surrounding regions. At the time of this study, the clinic had enrolled 15,933 patients, of whom 5,481 were currently on treatment. This retrospective cohort study was conducted at the University of Gondar’s comprehensive and specialized hospital, which was deliberately chosen not only due to data availability but also because of its distinct demographic and epidemiological profile. This region has a high burden of HIV, with a substantial proportion of PLHIV receiving ART. Additionally, the area serves as a referral center for surrounding rural districts, providing a heterogeneous mix of urban and rural populations, which enhances the generalizability of the findings. The healthcare infrastructure in the area also allows for long-term follow-up and relatively comprehensive medical records, which are essential for a mortality prediction study.

The study included all eligible patients from the electronic medical records from March 2005 to December 2024. For this study, total population sampling was employed, rather than random sampling. The cohort size for this study is 4,152, representing all the eligible patients.

Data source and study population

The secondary data in this study were obtained from the EMR-ART at the University of Gondar Comprehensive and Specialized Hospital. The data were extracted from the electronic medical records of PLHIV who received ART at the University of Gondar’s comprehensive and specialized hospital between March 2005 and December 2024. Inclusion criteria were: (1) age ≥15 years at ART initiation, (2) confirmed HIV diagnosis, and (3) complete baseline clinical and laboratory data available at the time of ART initiation. Patients were excluded if they had incomplete baseline records, transferred in from another facility without a full clinical history, or were lost to follow-up within one month of ART initiation. The final cohort included 4,152 patients who met all eligibility criteria and were followed until death, loss to follow-up, or the end of the study period.

The extraction of the relevant information from the EMR database was done following a standardized procedure in which a query was run on the database in a structured manner to obtain information regarding all the patients with a confirmed diagnosis of HIV and initiation of ART between March 2005 and December 2024. The baseline information was obtained from the first clinical encounter during which the initiation of ART occurred, while the outcome and follow-up information were obtained from the last clinical encounter. After extraction, the information was exported in CSV format, replacing the names of the patients with unique study identifiers for confidentiality purposes. To ensure the accuracy of the extracted information, 5% of the extracted information was randomly checked against the source document by two independent researchers.

All persons living with HIV on ART who visited the hospital for care in its ART clinic formed the source population. This dataset covered various variables in demography, clinical setups, and immunological and treatment-related factors that would be necessary to perform a model predictive. The outcome variable was the viral load status, categorized into two classes: Suppressed and Unsuppressed. The independent variables were sociodemographic characteristics: age, sex, marital status, occupational status, religion, and place of residence. Other characteristics included clinical factors: duration on ART, WHO clinical stage, duration with HIV, and TB co-infection. Further hematologic and immunological factors, including baseline and current CD4 counts and viral load, were considered. Treatment-related factors such as adherence to ART, regimen line, initiation and discontinuation of TPT (Tuberculosis Preventive Therapy), and CPT (Cotrimoxazole Preventive Therapy) usage were added to the model for improved interpretability.

Regarding the handling of longitudinal data, this study employed a cross-sectional snapshot design using the most recent clinical encounter for each patient as the outcome assessment point, while incorporating key historical variables (e.g., duration on ART in months, baseline CD4 count) that capture cumulative treatment exposure and immunological history. This approach was deliberately chosen for three reasons: (1) clinical decision-making in resource-limited settings typically relies on current patient status and treatment history rather than full temporal trajectories; (2) the electronic medical record system at the study site does not consistently capture all intermediate visit data with complete temporal granularity, making comprehensive longitudinal modeling potentially unreliable; and (3) the primary objective was to develop a practical screening tool that could be deployed using routinely available cross-sectional data, maximizing real-world applicability.

Data preprocessing

To assess the mechanism of missingness, we conducted systematic analyses comparing patients with complete versus missing data across key variables. Using Little’s Missing Completely at Random (MCAR) test, we found evidence that data were not MCAR (p < 0.001), suggesting a combination of Missing at Random (MAR) and Missing Not at Random (MNAR) mechanisms. For example, missing CD4 counts were more common among patients with poor adherence documentation, indicating MAR. To address potential bias from mode imputation, we implemented multiple sensitivity analyses: (1) complete-case analysis restricted to patients with no missing data (n=2,847) yielded similar model performance (Gradient Boosting accuracy 75% vs. 76% with imputation), suggesting minimal imputation-induced bias; (2) we compared mode imputation against multiple imputation by chained equations (MICE) on a random subset, finding no significant differences in coefficient estimates or model performance; and (3) we created missingness indicator variables for features with >5% missingness to allow the model to learn potential informative missingness patterns. These analyses support that our imputation approach did not substantially bias the final model.

The data was preprocessed, following a structured data preprocessing pipeline to ensure the quality and reliability of the data. Machine learning requires a high-quality dataset for prediction. Due to this, handling the missing data during the pre-processing of the dataset is a crucial phase. To assess the extent of data quality issues in the dataset, we examined the proportion of missing values across all features. As illustrated in Figure 1, several variables demonstrated varying degrees of missingness. Notably, ‘Baseline CD4 Count’ and ‘Recent CD4 Count’ exhibited the highest levels of missing data, with 15.2% and 15.0% of values missing, respectively. Other variables, such as ‘Duration on ART in Months’ and ‘Functional Status’, had missing data ranging from 1% to 5%, while the majority of features had less than 2% missing values. Overall, the dataset contained approximately 1.7% missing data. To address these missing values, we applied imputation techniques tailored to the nature of each variable. For continuous variables, such as CD4 counts and ART duration, mean imputation was employed based on the assumption of approximately symmetric distributions with minimal outliers. For categorical variables, including marital status, education level, and functional status, mode imputation was used, assuming that the most frequent category reasonably represents the underlying distribution in the absence of systematic bias. This imputation strategy was selected to preserve the integrity of the dataset while minimizing information loss and ensuring the inclusion of a maximal number of observations in subsequent analyses. The Simple Imputer class of the scikit-learn module was used to fill in the missing values in the dataset.OPEN IN VIEWER

Data pre-processing also includes encoding data, which is a crucial step. This study used one-hot and label encoding to encode categorical variables. Values with two or more category values are considered categorical if they are discrete and not continuous. In one-hot encoding, the categorical values are replaced by a number between 0 and 1.

Feature selection

From the EMR, we extract all relevant features and do correlation analysis to determine if there is a strongly correlated feature that will affect the model performance and introduce biases. The correlation heatmap illustrates the strength and direction of linear relationships among the variables in the dataset, with a specific focus on identifying factors associated with the viral load status, the primary outcome variable in this study. As shown in the heatmap, viral load status demonstrates modest positive correlations with Recent CD4 Count (r = 0.35), Recent_CD4_Category (r = 0.35), and Baseline CD4 Count (r = 0.24), suggesting that higher CD4 levels may be associated with better viral suppression outcomes. Additionally, weak positive correlations are observed with variables such as CPT Use (r = 0.05) and Duration on ART in Months (r = 0.11), indicating potential clinical relevance. Conversely, a slight negative correlation is noted between Age and viral load status (r = -0.07), implying that younger individuals may have more favorable viral outcomes (Figure 2). While most other features exhibit weak or negligible correlations with the outcome, this analysis provides initial insights into potentially influential predictors of viral load suppression, which will be further examined through multivariate modeling in subsequent stages of the research. Theoretically, the suppression status derives from the recent CD4 count and is directly related to the recent CD4 category, which is stated as low and high, derived directly from the count; we remove the recent CD4 count and the recent CD4 category from the final analysis. To prevent data leakage, we rigorously ensured that all predictor variables were based on information available before or at the time of viral load measurement, not after.OPEN IN VIEWER

Data balancing

In this dataset, the dependent variable (viral load status) was almost evenly distributed, with suppressed cases representing 50.4% and unsuppressed cases 49.6% (Figure 3). This natural balance between the classes eliminates the need for synthetic oversampling techniques, allowing models to be trained directly on the authentic data distribution. The equitable representation of both outcomes supports the development of a robust classifier without inherent bias toward a majority class. This is particularly advantageous in clinical settings, as it enables the creation of predictive models that are reliable and generalizable, ensuring accurate assessment of viral load suppression for critical decision-making.OPEN IN VIEWER

Machine learning models

For the estimation of viral load status, we evaluated eight ML algorithms, comprising both linear and ensemble-based models. Logistic regression was used as a baseline due to its interpretability and wide acceptance in clinical research. Ensemble methods such as random forest and gradient boosting (including its optimized variants like XGBoost) were selected for their ability to handle non-linear relationships and reduce overfitting through averaging or boosting techniques. Support vector machines, Naïve Bayes, and k-nearest neighbors were included to assess performance under different assumptions of class boundaries and data similarity. Decision trees were also used due to their intuitive rule-based structure, which aligns with the need for model interpretability in medical settings. These algorithms were chosen for their balance between predictive performance and interpretability, which is essential for clinical applicability and stakeholder trust. Each algorithm has its own assumptions: Logistic Regression assumes a logit linear relationship between predictors and outcome variable and independence of errors; Naive Bayes assumes conditional independence between features and outcome variable; SVM assumes separability and correct kernel choice (RBF kernel used after testing linear and polynomial kernel on data); tree-based algorithms (Random Forest, Gradient Boosting, Decision Tree, XGBoost) make no assumptions and can automatically deal with non-linear relationships. Deep learning models, such as neural networks, were not included in this study due to the relatively small dataset size, which could lead to overfitting and reduce generalizability. Furthermore, deep learning models often require extensive computational resources and hyperparameter tuning, which were beyond the current study’s scope. Before modeling, we performed feature selection using correlation analysis and mutual information to retain the most informative independent variables and reduce noise, thereby enhancing model performance and interpretability.

Model training and optimization

Model training and hyperparameter tuning were conducted in a standard computational environment using Python’s scikit-learn library. Training times ranged from approximately 3-8 minutes per algorithm, with simpler models (logistic regression, decision tree) converging faster than ensemble methods (Gradient Boosting, XGBoost) which required longer due to iterative boosting. Five-fold stratified cross-validation was employed during the training phase to prevent overfitting and ensure generalizability. A fixed random seed (42) was used for all models to enable exact reproducibility. The testing set (20% holdout) was never accessed during hyperparameter tuning or feature selection; it was used only for final model evaluation.

The dataset was split into training and testing subsets, with 80% allocated for training and 20% for testing. To optimize model performance, hyperparameter tuning was conducted using grid search for each model. For logistic regression, tuning focused on regularization strength (C), explored over a logarithmic scale from 10^-4 to 10³ and penalty type (both L1 and L2, which are for Lasso and Ridge regularization, respectively), while for Random Forest, these were several trees (between 50 to 500), maximum depth between 5 to 50, minimum samples per leaf from 1 to 10. In Gradient Boosting-based models, such as XGBoost and Gradient Boosting, important tunings were the learning rate, ranging from 0.001 to 0.3; estimators ranging from 50 to 500; and maximum depth between 3 to 15. For support vector machines (SVM), tuning was performed on kernel type (linear, polynomial, and RBF (Radial Basis Function), regularization parameter (C) ranging from 10^-3 to 10², and gamma values between 10^-4 to 10³ and decision trees were fine-tuned by varying max depth (3 to 20), minimum samples split (2 to 10), and evaluating both Gini impurity and entropy criteria.

KNN (k-Nearest Neighbors) tuned the number of neighbors from 1 to 50 and the distance metric by Euclidean, Manhattan, and Minkowski. Cross-validation was performed during model training to reduce overfitting. The use of k-fold cross-validation with k=5 was used for the training to ensure the model generalizes to unseen data. The study developed each model to have optimal predictive performance with generalizability and computational efficiency by systematically tuning hyperparameters and incorporating strategies to prevent overfitting.

Evaluation metrics

Model performance was evaluated using a suite of metrics selected for clinical relevance and robustness. The Area Under the Receiver Operating Characteristic Curve (AUC-ROC) quantified overall discriminative power. Primary emphasis was placed on class-specific precision and recall to ensure high sensitivity and reliability for all patient groups. The F1-score was used as a primary summary metric to effectively balance these competing clinical priorities. To compare the performance of the eight machine learning models statistically, paired evaluation metrics were used across cross-validation folds. This approach ensures not only robust and fair comparison but also aligns the model selection process with both statistical validity and clinical relevance in HIV management under resource-constrained conditions.

Overall, the process involves collecting EMR data, preprocessing it, selecting and training models, optimizing performance, and evaluating results to identify the best model (Figure 4).OPEN IN VIEWER

Implementation

All models were implemented in Python using machine learning libraries. The preprocessing, training, and testing of models are done on a Jupiter Notebook environment to ensure the reproducibility and transparency of the analysis. Results and methods have been duly recorded for future research and verification purposes.

Descriptive statistics

The dataset consisted of 4,152 patients, with 50.4% having suppressed viral load and 49.6 % having unsuppressed viral load. The mean age of participants was 46.72 years (SD = 11.38), with the majority being female (59.51%) (Table 1).

OPEN IN VIEWER

Table 1.

Socio-demographic variables.

Variable	Category	Frequency	Percentage
Sex	F	2471	59.51
	M	1681	40.49
Marital Status	Married	2700	65.03
	Divorced	698	16.81
	Never Married	437	10.53
	Widowed	317	7.63
Education level	Secondary Education	2145	51.66
	No Education	805	19.39
	Primary Education	789	19.00
	Higher Education	413	9.95
Residence	Urban	3177	76.52
	Rural	975	23.48
Religion	Orthodox	3910	94.17
	Muslim	206	4.96
	Protestant	26	0.63
	Catholic	9	0.22
	Other	1	0.02

Regarding clinical characteristics, the median duration of ART was approximately 12.7 years (152.5 months), and most patients were in WHO clinical stage 4. The dataset exhibited a balanced class, with suppressed viral load being the predominant class (50.4%), and unsuppressed (49.6 %) viral load categories. The near-equal distribution of suppressed (50.4%) and unsuppressed (49.6%) viral load in this cohort differs from many routine HIV care settings where suppression rates typically exceed 80-90%. This unique distribution reflects the study’s sampling strategy and setting characteristics: (1) The cohort includes patients enrolled over 19 years (2005-2024), encompassing eras with less effective early ART regimens and higher failure rates; (2) as a tertiary referral center, the University of Gondar Hospital manages a disproportionate number of complex, treatment-experienced patients transferred from peripheral clinics after initial treatment failure; (3) the retrospective design captured all available patients rather than only those with optimal outcomes; and (4) the study period predates and overlaps with the global scale-up of viral load monitoring, meaning some patients may have had undetected viremia for extended periods. While this balanced distribution is advantageous for model training (avoiding class imbalance challenges), it may not reflect the current suppression rates in optimally managed cohorts, which should be considered when generalizing findings (Table 2).

OPEN IN VIEWER

Table 2.

Clinical-related variables.

Variable	Category	Frequency	Percentage
BMI	Normal	2607	62.79
	Overweight	733	17.65
	Underweight	630	15.17
	Obese	182	4.38
TB screening Result	No	4070	98.03
	Yes	82	1.97
Functional Status	Working	3861	92.99
	Ambulatory	163	3.93
	Bedridden	128	3.08
Viral Load Status	Suppressed	2092	50.4
	Unsuppressed	2060	49.6
Adherence	Good	3611	86.97
	Poor	355	8.55
	Fair	186	4.48
Regimen Line	First Line	3811	91.79
	Second Line	319	7.68
	Third Line	22	0.53
TPT Started	Yes	2958	71.24
	No	1194	28.76
CPT Use	Yes	3313	79.79
CPT Use	No	839	20.21

Model training results

In this study, the comparative evaluation of machine learning models revealed clear distinctions in performance across the selected algorithms. Gradient Boosting and Random Forest emerged as the top-performing models, achieving accuracies of 0.73 and 0.72, respectively, with consistently balanced precision, recall, and F1 scores (all around 0.70). This highlights their robustness and reliability in classifying the data effectively. Naive Bayes stood out for its exceptionally high recall (0.84), suggesting a strong ability to identify positive cases, albeit at the expense of lower precision (0.59), which may increase false positives. Decision Tree and XGBoost offered stable mid-range performance, with accuracies of 0.67 and 0.69 and balanced metrics, making them reasonable alternatives. Logistic Regression yielded moderate results (accuracy = 0.62), while SVM and KNN underperformed, both with accuracies of 0.54 and comparatively lower overall scores (Table 3). Collectively, these findings indicate that ensemble-based methods, particularly Gradient Boosting and Random Forest, are the most effective approaches for this dataset, providing a promising direction for building accurate and reliable predictive models. Beyond discriminative performance, we assessed the calibration of the optimized Gradient Boosting model and the agreement between predicted probabilities and observed frequencies of viral suppression. The Brier score, which measures the mean squared difference between predicted probabilities and actual outcomes, was 0.18 (range 0-1, with 0 indicating perfect calibration).

OPEN IN VIEWER

Table 3.

Model training results.

Model	Accuracy	Precision	Recall	F1 score
Logistic Regression	0.62	0.60	0.55	0.58
Random Forest	0.72	0.70	0.70	0.70
Gradient Boosting	0.73	0.70	0.70	0.70
Naive Bayes	0.66	0.59	0.84	0.69
SVM	0.54	0.51	0.55	0.53
KNN	0.54	0.50	0.58	0.54
Decision Tree	0.67	0.65	0.62	0.64
XGBoost	0.69	0.66	0.67	0.66

The ROC curves shown in Figure 5 illustrate the trade-off between the true positive rate and false positive rate for all models evaluated. The Gradient Boosting model achieved the highest area under the curve (AUC = 0.79), closely followed by Random Forest (AUC = 0.78) and XGBoost (AUC = 0.76), confirming their strong discriminative ability. Naive Bayes also performed well, with an AUC of 0.73, demonstrating competitive sensitivity in distinguishing between classes. In contrast, Logistic Regression (AUC = 0.66) and Decision Tree (AUC = 0.67) provided moderate results, while SVM (AUC = 0.57) and KNN (AUC = 0.56) performed poorly, with curves lying closer to the diagonal reference line, indicating limited predictive power. Overall, the ROC analysis reinforces the superiority of ensemble methods, particularly Gradient Boosting and Random Forest, as they consistently achieved higher AUC values and demonstrated better classification performance across thresholds compared to traditional algorithms.OPEN IN VIEWER

Examining class-specific performance for the optimized Gradient Boosting model, the confusion matrix (Figure 6) reveals a sensitivity (recall) of 0.65 for the unsuppressed class (211 correctly identified of 324 total unsuppressed cases) and 0.76 for the suppressed class (339 correctly identified of 447 suppressed cases). This indicates that the model correctly identifies approximately two-thirds of patients with unsuppressed viral load, those at highest risk of treatment failure and clinical progression, while misclassifying 35% as suppressed (false negatives). The higher sensitivity for suppressed cases reflects the model’s better performance on the majority class, even in a balanced dataset.OPEN IN VIEWER

Feature importance

The feature importance analysis from the tuned Gradient Boosting model (Figure 7) highlights the variables that most strongly contributed to predicting viral load suppression status. Feature importance was calculated using the built-in feature importance attribute of the optimized Gradient Boosting model, which averages the reduction in impurity (Gini importance) across all trees for each feature. Baseline CD4 Category emerged as the most influential predictor, contributing more than half of the total importance score, underscoring the critical role of immunological status at treatment initiation. The second most significant predictor was Duration on ART in Months, indicating that longer treatment duration is strongly associated with viral load suppression outcomes. Other features, such as sex and age, contributed moderately, while factors like regimen line, CPT use, and marital status had relatively smaller impacts. Variables including WHO stage, adherence, BMI, educational level, TPT initiation, and TB screening result showed minimal importance in the model. These findings suggest that immunological markers and treatment duration are the most decisive factors in determining viral load suppression, while sociodemographic and clinical background factors play a comparatively limited role in prediction.OPEN IN VIEWER

Model explainability

SHAP (SHapley Additive exPlanations) analysis was performed using the KernelExplainer with 1,000 background samples from the training set, as implemented in the SHAP Python library (version 0.44.0). The interpretation of the model using SHAP revealed that the most important features driving predictions were the duration on ART, baseline CD4 category, and patient age. The direction of impact demonstrated that longer ART duration and a higher baseline CD4 count were strongly associated with a decreased risk of the outcome, whereas advanced WHO clinical stage and male sex were associated with an increased risk. The analysis also highlighted the protective effect of initiating TPT. This feature importance ranking and directional analysis align with established clinical understanding, thereby enhancing the credibility and interpretability of the model’s decision-making process (Figure 8).OPEN IN VIEWER

The model’s prediction for a specific individual was further elucidated using a SHAP waterfall plot. The analysis begins with the baseline model output, E [f(X)] = -0.201, representing the average prediction across the dataset. For this patient, the most influential factors increasing their predicted risk were a WHO Stage of 4 (+1.3) and being male (Sex=1, +1.27). Their age of 53 also contributed significantly to a higher risk (+0.86). These factors were partially offset by a protective effect from a longer duration on ART (162 months), which decreased the risk (-0.41) (Figure 9). Contributions from other features, such as CPT Use and Baseline CD4 Category, had smaller positive effects. Ultimately, the cumulative effect of all these feature contributions shifted the model’s prediction from the population average to a final, higher value for this individual, clearly illustrating how their specific clinical and demographic profile led to the elevated risk score.OPEN IN VIEWER