The analysis used longitudinal, deidentified electronic medical record (EMR) data from patients initiated on ART between January 1, 2017, and November 30, 2023, across 2156 facilities in all 47 counties of Kenya. Eligible patients had at least 2 follow-up visits or drug pickups after ART initiation, were followed for more than 28 days postinitiation, and had returned to care if they had experienced an initial treatment interruption.

The preprocessing, cleaning, feature engineering, and covariate definitions are detailed in our previous work [6]. Treatment interruption was defined as no ART drug pickup and no clinical contact for greater than 28 days after the last expected contact [14]. Time-invariant covariates included age, gender, education, marital and employment status, noncommunicable disease, baseline regimen, and ART initiation year, while time-varying covariates were adherence, alcohol intake, World Health Organization stage, viral load, clinical stability assessment, regimen line, multimonth dispensing (MMD), differentiated service delivery (DSD) model, and prevention with positives package. Full list of covariates and description is included in Table S1 in Multimedia Appendix 1. A key preprocessing step included creating an “unknown” category for categorical variables with missing data. Engineered features included the COVID-19 pandemic period—dummy variable to indicate the COVID-19 pandemic; missed appointment—failure to attend a scheduled visit within 1-3 days of the scheduled date; defaulted appointment—failure to attend a scheduled visit within 4-28 days; and the duration of treatment interruption—average duration of any treatment interruption before a patient returned to care.  

We performed 2 sets of modeling, that is, time-invariant model, for first-ever treatment interruption, and a time-varying model, for multiple interruptions. The data structure differs between the models: in a time-invariant model, covariates were fixed, with 1 row per patient containing the time to interruption and interruption status. In a time-varying model, covariates change over time, with multiple rows per patient representing different time intervals and multiple interruption statuses.

Our base reference model was the traditional CPH model. The time-invariant CPH assessed the effect of fixed covariates on time to first interruption, with an assumption that the hazard ratio remains constant over time, while the time-varying allowed covariates to change over time, capturing the dynamic effects of these covariates on the hazard at different intervals.

Our choice of ML models was based on algorithms with the ability to fit time-varying covariates and their computational feasibility at scale. We evaluated 3 classes of algorithms: boosting methods for their efficiency and predictive performance, regularized shrinkage models for handling multicollinearity and improving generalization, and RP for its ability to capture nonlinear relationships and interactions while remaining interpretable. Although random survival forests are widely used in survival analysis, their computational demands and limited support for fully time-varying covariates constrained their applicability in this study. The trained models are as described:

We initially trained each model using default hyperparameters as a baseline, followed by hyperparameter tuning through a 10-fold cross-validation over a grid of parameter combinations. The optimal set of hyperparameters was selected based on the highest Harrell’s concordance index (C-index) score, reflecting the study’s primary objective of ranking patients by treatment interruption risk. The optimal parameter sets were then used for final model training. All analysis and modeling were performed using R Statistical Software (version 4.4.0; R Foundation for Statistical Computing).

We evaluated model performance using an 80:20 train-test split (including the Cox models to ensure comparability). Likewise, a 10-fold cross-validation was performed on the train set, and models were trained on k-1 folds and validated on the remaining fold, repeating this process 10 times, with each fold serving as the validation set once.

Final model evaluation and comparison were performed using the C-index, a widely used metric for assessing the performance of survival models. The C-index measures the model’s ability to correctly rank pairs of observations based on survival times, with values ranging from 0.5 (random prediction) to 1.0 (perfect ranking). We calculated the C-index for each model on both the training and test sets to evaluate predictive accuracy. This approach ensured that the models performed well during training and on validation data, providing a reliable measure of generalizability. Furthermore, to assess the sensitivity and uncertainty of the predictive performance measures, we conducted 200 bootstrap resampling iterations to estimate the IQR distribution of the scores for our best performing model.

To interpret our model predictions, we used the DALEX package in R [15] and applied 2 levels of model-agnostic explanation, that is, model level global explanations—using permutation variable importance and partial dependence plots, and individual-level explanations—using variable attribution—breakdown and SHAP values and ceteris paribus (CP) analysis.

The global-level explanations are listed as follows:

The individual-level explanations are listed as follows:

This study was reviewed and approved as nonhuman subject research by the institutional review boards of the University of Maryland, Baltimore (HP-00108126), and the Aga Khan University Nairobi (2023/ISERC-94 (v2)). The analysis used secondary EMR data that were fully deidentified prior to receipt by the study team in accordance with Kenya’s Data Protection Act and institutional guidelines; therefore, written informed consent was not required.

We compared the traditional CPH model with 6 different ML models to predict the risk of treatment interruption in 621,115 patients receiving ART in 256 facilities across all 47 counties in Kenya. The analysis was conducted in 2 parts: time-invariant models for predicting time-to-first treatment interruptions and time-varying models for predicting time-to-multiple treatment interruptions. Model evaluation was performed using an 80:20 train-test split, and Harrell’s C-index was applied to compare the performance of the models on both the train and test sets. Prediction from the best-performing model was then interpreted at both global and individual levels.

The RP model outperformed all other models in both the time-invariant and time-varying analyses. It achieved C-index values of 0.8115 for the train set and 0.8105 for the test set in the time-invariant model, and 0.8918 for the train set and 0.8901 for the test set in the time-varying model. These results align with findings from previous studies that have compared CPH models with ML models, demonstrating the effectiveness of RP models in survival analysis. Previous studies have reported C-index values ranging from 0.63 to 0.96, primarily from ML models such as GBM, random survival forest, and survival support vector machines in predicting survival for patients with cancer [16-19]. Other studies that focused on predicting treatment interruptions in patients on ART had typically reported metrics such as sensitivity, specificity, and area under the curve, with values ranging between 61.9% and 76.0% [20-24]. Few of these studies have incorporated time-varying covariates [16,25]. Even when compared with prior studies that evaluated ML models against CPH [17], our RP approach achieved higher discrimination, with C-index values of 0.8105 and 0.8901 compared with reported XGBoost performance of approximately 0.73. Notably, our study stands out as one of the few to apply RP within a time-varying survival framework, effectively constructing a survival tree that accommodates time-varying covariates.

The RP method offers flexibility in detecting complex interactions without assuming proportional hazards. They automatically select covariates and split points, allowing more adaptive modeling of survival times [26]. This method is advantageous for identifying nonlinear relationships and interactions, making it particularly useful in our case where covariate effects changed over time. The RP model also had the advantage of simple hyperparameter tuning and faster computational times compared with the other ML models we evaluated. We attempted to train a random survival forest model but were unsuccessful due to the large computational resources and extensive time required to handle our large dataset. Likewise, the widely used randomForestSRC R package does not currently support time-varying covariates directly. While random forests are known for their high predictive performance, their significant computational requirements make them less practical for real-world applications [16].

To ensure that our results were fully comparable with the traditional CPH, we implemented a cohesive framework for explaining the ML predictions with meaningful insights. Using a model-agnostic [15] explainer framework, we gained detailed insights into how specific variables influenced individual-level predictions and overall global outcomes, similar to the interpretations provided by hazard ratios in the CPH model. The explainer framework offers additional advantages by capturing complex, nonlinear relationships that traditional CPH may overlook.

On a global level, the permutation-based variable importance from the RP model was similar to the ranking of the absolute values of the coefficients of the CPH model for both time-invariant and time-varying models. In addition, the effects of the selected covariates we considered in the partial dependence plots were consistent with the direction of hazard ratios in the CPH model. This combination of techniques offers diverse insights into our model predictive performance. They reveal the most significant effects and risks associated with treatment interruption across various covariates without the assumptions of linearity, hence capturing potential nonlinear effects that traditional CPH models might overlook. Overall, these methods provide more comprehensive global explanations, particularly for complex variable interactions.

Further exploration using breakdown attribution and SHAP plots reveals how individual variables influence each patient’s predicted risk within specific time intervals. On a global level, the time-invariant model was primarily influenced by MMD. Considering this at the patient level, the breakdown and SHAP plots provide detailed insights into how MMD increases or decreases risk, explaining why patient A was at a higher risk than patient B based on MMD. In addition, these plots highlight the most influential covariates for predicting the risk of a first treatment interruption for each patient, emphasizing individual variability in risk factors. For the time-varying model, we furthered the instance-level exploration using CP profiles. These plots provided insights into how a patient’s predicted risk would change if certain variables were adjusted. For example, with patient C—who has no history of treatment interruption—we explored the predicted risk of an interruption if the patient were to miss or default on an appointment, transition to a less intensive DSD model, or receive a different drug-dispensing interval. In contrast, for patient D—a patient with a history of interruptions—the focus was on understanding how the predicted risk of another interruption might change if they transitioned to a less intensive care model, adjusted their dispensing interval, or missed an appointment. The changes highlighted for each patient provide actionable insights to guide transitions between DSD models and targeted interventions to prevent treatment interruptions. Following ART initiation, CP profiles can identify patients whose interruption risk under standard care would decrease with earlier DSD enrollment or longer dispensing intervals. For patients returning to care after an interruption, CP profiles quantify how changes in follow-up intensity or MMD duration influence recurrent interruption risk, supporting tailored rather than routine care. These scenario-based insights align with recent World Health Organization and global guidance [27,28] emphasizing person-centered, targeted care packages for patients reengaging in treatment.

Bringing the various components together, our model could be integrated into an EMR system to predict a patient’s risk of treatment interruption in real time. This integration would allow the model to continuously process each patient’s historical data, current visit date, and next scheduled visit date, creating a new data point for each interaction. With each visit, the EMR system would automatically update risk predictions, enabling health care providers to proactively identify high-risk patients before interruptions occur. The Kenya EMR currently includes an ML module that generates interruption risk scores based on a limited set of patient-centric covariates. By leveraging our model’s broader range of covariates, the current system can be optimized for improved prediction accuracy and comprehensive risk assessment.

Our study had some limitations. First is the issue of missing data, as highlighted in our previous work [6]. This was addressed by creating an “unknown” category for affected variables, given that the data were not missing at random. Variables with more than 60% missingness, such as baseline and follow-up CD4 counts, were excluded, which may have introduced bias. Sensitivity analyses comparing complete-case models with alternative “unknown” category specifications were conducted in prior work [6] and showed consistent results. Multiple imputation was not considered, given data were not missing at random. Second, while the explainer we used is model-agnostic, it lacks built-in support for survival-specific metrics such as Brier score and area under the curve. Survex explainer is specifically designed for survival models, but it currently does not support the LTRCtrees package or time-varying survival models. A potential work-around could be treating the model as a classification task, focusing on whether there was an interruption or not within a specified fixed time, although this approach will sacrifice the temporal dimension of survival analysis or future studies can consider implementing other Survex-supported algorithms with the ability to fit time-varying covariates. Finally, generating breakdown and SHAP attribution profiles for individual predictions was computationally intensive, particularly for SHAP, which may limit real-time EMR use. However, subsampling across observations or feature permutations can substantially reduce computational cost while preserving accurate SHAP approximations.

Our study compared the predictive performance of the traditional Cox model with 6 different ML models for predicting the risk of treatment interruptions in patients on ART. Overall, our results demonstrate that ML models outperform the traditional Cox model, with RP, in particular, effectively capturing nonlinear effects and complex interactions without relying on proportional hazards assumptions, making it well suited for time-varying survival analyses. Our model-agnostic explanation aligns well with the Cox model results while also offering additional, patient-specific insights for targeted interventions. The combination of enhanced predictive performance and interpretability demonstrates how ML models can effectively support patient-centered care strategies to reduce the likelihood of treatment interruptions in ART.