Prediction of First and Multiple Antiretroviral Therapy Interruptions in People Living With HIV: Comparative Survival Analysis Using Cox and Explainable Machine Learning Models

Introduction

The Cox proportional hazards (CPH) model has long been a foundational tool for time-to-event analysis, widely used in health care research to examine treatment interruptions, including among patients on antiretroviral therapy (ART). The CPH model is a favorable choice due to its straightforward interpretability; its coefficients can be easily understood as hazard ratios, which quantify the relative effect of each covariate on the likelihood of an event, offering clear insights for clinical decision-making and risk assessment [1,2]. In addition, its flexibility in handling time-dependent covariates allows it to incorporate variables that change over time, providing a more accurate reflection of patient dynamics [3-5]. Our previous work [6] focused on using time-dependent covariates within the CPH model to enhance its application in modeling the complex patterns of multiple treatment interruptions in patients on ART, aiming to capture the nuanced and changing influences on treatment adherence.

Despite its popularity, the CPH model has limitations. A key disadvantage is its handling of high-dimensional data with many covariates, as it can be prone to overfitting and may yield unreliable results due to multicollinearity [4]. This risk is particularly pronounced in time-dependent CPH models, where changing covariates can amplify the potential for erroneous inferences [4]. In addition, as a semiparametric model, it does not fully specify the baseline hazard function, potentially limiting its ability to capture complex, nonlinear relationships [7]. These limitations highlight the need for more flexible modeling approaches, especially when working with complex health care datasets with nonlinear, time-varying, or high-dimensional characteristics.

Machine learning (ML) models are increasingly used in time-to-event analysis due to their ability to capture complex, nonlinear relationships at scale. In HIV and ART research, however, most survival ML applications such as random survival forests and gradient boosting primarily rely on baseline or time-invariant covariates, limiting their ability to reflect the dynamic nature of ART care. Although some methods, including recursive partitioning (RP), have been adapted to incorporate temporal structures and time-varying effects [8-10], support for fully time-varying covariates remains limited in many survival ML frameworks. RP provides a flexible alternative by modeling nonlinear effects and interactions with time-varying covariates, yet its application to large-scale, longitudinal ART cohorts has been limited. In addition, advances in model interpretability have made ML models more explainable, helping to bridge the gap between predictive accuracy and practical applicability. Through methods such as permutation importance, partial dependence plots, and Shapley Additive Explanations (SHAP) values [11-13], ML models now provide insights into variable importance and individual predictions, making them increasingly useful for clinical decision-making and patient-centered care.

The aim of this paper is first to compare the performance of traditional CPH model and ML models in predicting treatment interruptions among patients on ART. This comparison aims to evaluate whether ML models, with their flexibility in handling nonlinear relationships and time-dependent covariates, can provide more accurate and reliable predictions than the CPH model. Second, provide clear, intuitive explanations for the predictions at both the global level to understand overall patterns and key risk factors and the local level for individualized predictions for a patient. By combining predictive accuracy with interpretability, this study aims to support health care providers in identifying patients at high risk for treatment interruption and in implementing personalized interventions to enhance treatment retention.

Methods

Study Population and Data

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.

Data Preprocessing and Definitions

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.  

Analysis and Models

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:

• Gradient boosting machine (GBM): an ensemble learning method that builds a series of decision trees, where each new tree corrects the errors of the previous ones. It is powerful for time-invariant survival analysis due to its ability to handle complex, nonlinear relationships and improve predictive accuracy by minimizing loss iteratively.
• Extreme gradient boosting (XGBoost): an optimized version of GBM, known for its speed and performance, especially with large datasets. It is well suited for our analysis due to its efficiency, scalability, and ability to manage missing data, offering robust predictions in survival tasks. The Accelerated Failure Time model was trained for a time-varying dataset.
• Regularized generalized linear models: Extends traditional regression models with regularization techniques—Ridge (L2), Lasso (L1), and Elastic Net (a combination of both)—to prevent overfitting and handle multicollinearity. These models are beneficial in time-invariant survival analysis for selecting relevant covariates while improving model generalization, especially with high-dimensional data.
• Recursive partitioning (LTRCtrees/Rpart): Uses decision trees to model survival data, splitting the dataset based on covariates to form homogeneous subgroups. It is advantageous due to its interpretability and ability to handle nonlinear relationships without making parametric assumptions.

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).

Model Evaluation

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.

Model Interpretability

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.

Global-Level Explanations

The global-level explanations are listed as follows:

1. Variable importance: Using a permutation-based approach, we assessed how removing the effect of selected variables through multiple permutations impacts model performance.
2. Partial dependence profiles: We explored the effects of selected variables by examining how the model’s predicted values change as a function of individual variables, based on randomly selected observations.

Individual-Level Explanations

The individual-level explanations are listed as follows:

1. Variable attributions: This explains how a model’s prediction for a single observation differs from the average prediction and distributes the difference among explanatory variables. We applied two methods: (1) Breakdown additive attribution, a procedure that decomposes the model’s prediction into contributions from each explanatory variable. (2) SHAP, the averaging of variable attribution values over all (or several) possible orderings to show their individual impacts.
2. CP profiles: This explores how the model’s prediction changes when the value of a single explanatory variable is altered, while keeping all other variables constant.

Ethical Considerations

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.

Results

Our dataset comprised 621,115 unique patients on ART, with 432,041 first treatment interruption events and 959,170 multiple treatment interruption events. Data were split into an 80:20 ratio for the training and test sets, respectively. The hyperparameters for each model, along with the final parameter combinations used in the final model, are shown in Table 1. The predictive performance of ML models and their comparison with the CPH are shown in Table 2 and Figure 1, summarized by the estimated C-index scores and 95% CIs for both the train and test sets.

‎

In the time-invariant analysis, the glmnet models—Lasso, Ridge, and Elastic-Net—produced comparable results to the CPH model, with train C-index values ranging from 0.7799 to 0.7801 and test C-index values ranging from 0.7794 to 0.7796. XGBoost performed slightly worse than the CPH model, while the GBM model had the lowest performance. The LTRCtrees model outperformed all other ML models and the CPH model with a C-index value of 0.8115 for the train set and 0.8105 for the test set. The LTRCtrees model remained the top performer in the time-varying analysis with a C-index of 0.8918 for the train set and 0.8901 for the test set. XGBoost was the next best performing ML model and slightly comparable with the CPH model. We further evaluated the best-performing model (ie, LTRCtrees), using a 200-bootstrap resamples; the C-index interquartile range remained stable across resamples.

Figures 2 and 3 show the permutation-based variable importance for the top 10 variables in both the time-invariant and time-varying models. MMD emerged as the most important variable in the time-invariant model, while the duration of treatment interruption was most important in the time-varying model. In addition, variables such as viral load, ART start year, stability assessment, baseline regimen, and occupation were consistently ranked among the top 10 most influential variables in both models. The permutation-based variable importance results were similar to the ranking of the absolute values of the coefficients of the CPH models (Table 3). MMD, missed and defaulted appointments, year of ART initiation, occupation, and stability assessment consistently ranked among the top 10 variables in both the time-invariant LTRCtree and CPH models. In the time-varying models, duration of interruption, MMD, year of ART initiation, baseline regimen, and occupation were consistently top influencers in both LTRCtree and CPH models.

‎

‎

The full results of the time-invariant and time-varying CPH models are included in Table S2 in Multimedia Appendix 1. The partial dependence plots for the top 4 variables from both models are displayed in Figures 4 and 5. These plots show the average predicted risk for each of these variables. Short dispensing interval (MMD of 1 month), initiating ART in 2021, and having not missed or defaulted on an appointment were associated with higher predicted risks of first treatment interruption. Conversely, longer dispensing intervals (6+ months MMD), initiating ART in 2023, and having missed or defaulted on an appointment had a lower predicted risk of first treatment interruption. Short interruption durations (<30 days), unknown viral load, separated marital status, and unspecified alcohol use were all associated with higher predicted risk of multiple treatment interruption, while having no prior interruptions (ie, interruption duration of zero), lower viral load (<200), and no alcohol consumption showed lower predicted risk of multiple treatment interruption.

‎

‎

To examine individual-level predictions for the time-invariant model, we randomly selected 2 observations from the dataset: one representing a patient (labeled “patient A”) who experienced a first treatment interruption, and another representing a patient (labeled “patient B”) who did not experience an interruption by the end of the follow-up period.

The row labeled “intercept” represents the overall average prediction value (1.516) for the model—baseline prediction (Figure 6A). The subsequent rows show how the mean prediction changes when a specific covariate is fixed. Green bars indicate an increase in the risk of treatment interruption, while red bars show a decrease in the risk. The final row (prediction), shown with a blue bar, displays the overall mean value and change, giving the final predicted risk of treatment interruption for the patient. Patient A had a high predicted risk of treatment interruption, with a risk score of 2.051, which is 0.535 points higher than the model’s average prediction (Figure 6). Covariates influencing this prediction are short MMD intervals (1 month), viral load of 200‐999, did not receive a prevention with positive package, not defaulted on appointment, non–COVID-19 pandemic period, aged 25‐34 years at the start of ART, and baseline regimen of 2NRTIs + NNRTI. Patient B, on the other hand, had a low predicted risk (0) of interruption, largely driven by the long MMD interval (6+ months), enrolled in the DSD fast-track model, and had prior missed and defaulted appointments contributing to a reduced risk (Figure 7).

‎

‎

The right plots of both Figures 6 and 7 show the SHAP plots for the patients, where the SHAP values represent the computed mean attribution across 25 random orderings of the explanatory variables. The red bars indicate an average increase in risk, while the green bars show a decrease in risk. The violet box plots summarize the distribution of the attributions over the 25 orderings. The plots reveal similar influential variables as in the breakdown attributions, with a slight reordering of importance and the inclusion of additional variables—such as education level and stability assessment for patient A, and COVID-19 pandemic period for patient B.

To examine the time-varying model’s prediction for a single time-interval instance for a patient, we randomly selected 2 patients: patient C, who had experienced 1 or more interruptions and subsequently returned to care; patient D, who had no history of treatment interruption. Figure 6 shows the breakdown and SHAP attribution plots for a single time interval for patient C. After considering all variable contributions, the final predicted risk of subsequent treatment interruption for this patient was 2.43, which is 1.614 points above the average prediction. The key covariate contributing to this risk prediction is the patient’s previous interruption duration (duration of 30‐179 days). In contrast, being on a 6-month MMD interval and having not missed or defaulted on any appointments were associated with a lower risk of subsequent interruptions. The SHAP plot reflects similar variable contributions (Figure 8B). The predicted risk of treatment interruption for patient D was 0.181, which is 0.635 lower than the average model prediction (Figure 9A). Key covariates contributing to this reduced risk include no prior treatment interruptions (interruption duration of zero), being on a 3‐ to 5-month MMD interval, having a viral load between 200 and 999 cp/mL, the COVID-19 pandemic period, no defaulted appointments, good adherence, and starting ART between the ages of 35 and 44 years. SHAP plots (Figure 9B) further highlighted covariates associated with a lower risk for this patient.

‎

‎

Using CP profile plots, we further examined how the model’s predictions for these patients would shift if key covariates—such as the duration of the previous interruption, viral load, DSD model, MMD, missed appointments, and defaulted appointments—were altered. This allowed us to explore how changes in time-varying covariates impact the risk of multiple treatment interruptions. The blue dot and dotted line in the CP plot (Figures 10 and 11) represent the values of the covariates and the corresponding prediction for the patient. Each green bar in the plot shows the predicted change in value and direction for a specific covariate, with bars pointing to the right of the dotted line indicating an increased risk of interruption and bars pointing to the left indicating a decreased risk.

‎

‎

The CP plot (Figure 10) for patient C indicates that their risk of future treatment interruptions would increase significantly if they were to miss or default on an appointment or be placed on a shorter MMD interval. However, the predicted risk remains unchanged if they were to transition across different DSD models. The future risk of treatment interruption would increase for patient D if they miss or default on an appointment, had an increase in their viral load, were placed on a 1-month MMD interval, or experienced a treatment interruption, irrespective of how long it took to return to care (Figure 11). Transitioning from their current standard care DSD model to fast track or community ART distribution health care worker–led models would slightly reduce the risk of interruption.

Discussion

Principal Findings

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.

Limitations

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.

Conclusions

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.