Revealing HIV viral load patterns using unsupervised machine learning and cluster summarization

Human subjects protection

This proposal was reviewed and approved by the University of Rochester Human Subjects Review Board (protocol number RSRB00068884). Consent was waived by the review board due to de-identification of the data set. The analysis in this paper is presented in compliance with Center for Medicare Services (CMS) current cell size suppression policy²⁰. Data were coded such that patients could not be identified directly in compliance with the Department of Health and Human Services Regulations for the Protection of Human Subjects (45 CFR 46.101(b)(4)).

Study data

We obtained medical encounter data from all patients with an HIV diagnosis in the University of Rochester Medical Center’s electronic medical record system (EMR) between 2011–2016, including age, gender, race, ethnicity, and VL measurements. There were a total of 1,892 patients with at least one VL measured, with 1,576 of these patients having at least three VL measurements.

Measurements ≤48 copies/mL, present as categorical values “NEG", “POS < 20", or “POS < 48" were transformed into numerical values of 0, 20, and 48 respectively. The deidentified study data containing only viral load measurements and relative time are available at https://doi.org/10.5281/zenodo.1313245²¹.

Hardware and software specifications

Analyses were performed on a Windows 8 server with Intel(R) Xeon(R) CPUs E5-2620 v2 @ 2.10GHz and 256GB of RAM. Python 2.7 was used for most data mining and machine learning under Spyder v.3 installed from Anaconda2 (64-bit). The default packages available in Anaconda were used for analysis, including, but not limited to: NumPy, scikit-learn, SciPy, datetime, csv, math, Matplotlib, pip, operator, copy, random, and time. Using pip we installed the webcolors and pydotplus packages for rendering a decision tree. SQLite was used to store, query, and clean ~the data. Analytic code is available for download at https://github.com/SamirRCHI/Viral_Load_Data_Categorization.

Viral load analysis methods

Since VL data is asynchronous and noisy, with variable numbers of data points for each subject, we excluded patients with ≤ 2 VL measurements as too few to accurately assess VL patterns. Based on temporal patterns of VL described in the literature, the VL pair distribution of our data (Figure S1), and a further extensive investigation into the data, we hypothesized six potential temporal VL patterns, defined in Table 1 and illustrated in Figure 1.

Figure 1. Possible HIV viral load patterns.

Examples of each type of viral load pattern. Note that actual viral load patterns are noisier and may often be more difficult to distinguish. The magnitudes of viral load values reflect those found in the dataset.

It is important to note that these definitions are pattern based, and do not explicitly select absolute VL cutoff levels or a specific temporal window, as other reports have done^8–10,13. This has the advantage of allowing the absolute VL levels and critical time windows to emerge from the analysis. It also does not preclude incorporation of absolute levels (e.g. VLM>400) at a later stage into the pattern specification.

Feature vector definition

Mathematical notations for this work are described in Table 2. We next designed a feature vector to capture characteristics that would allow us to distinguish between VL patterns. VL values at the lower limit of detection are a function of the specific assay used, and appear in our data set as 0, 20, and 48 copies/mL (Figure S1). Thus, plots of the log₁₀ transformed data have discretely spaced values at the lower level of detection, capturing the undetectable range of viral load. Additionally, we adjusted the data by log₁₀[V L + 10] to avoid log₁₀[0]. The addition of 10 to VL (instead of 1) is used to minimize the distance between the undetectable values: 0, 20, and 48 (copies/mL). Thus, in our notation, all the values related to viral load are assumed to have been adjusted to this measure. For example, min_{V L} = log₁₀[0 + 10] = 1 and max_{V L} = log₁₀[10⁷ + 10] ≈ 7.

Using the transformed VL data, we next extract several relevant features of the VL measurements over time. These features are used for machine learning classification of individual patient VL time series, and designed to distinguish patterns in VL change while minimizing the effects of noise. We do not limit feature extraction based on the total elapsed time of viral load measurements because the optimal time-point for determining viral load class is not well established. The attributes for feature extraction are: relative area of viral exposure, weighted recency reliability, adjusted maximal difference, and interquartile range. The definitions include:

Statistical analysis

Machine learning methods for cluster classification were compared by calculating F₁ scores, the harmonic mean of precision and recall²², defined by Equation 9–Equation 11.

Analytic terminology

Here we formally define keywords appearing in the analysis: Let Feature extraction be the process of determining the values Ȧ, wRR, $\stackrel{˅}{D}$, and IQR from a set of patients (using their viral load patterns) with the formulations given above. Then a feature vector ($\overrightarrow{F}$_p) contains the values Ȧ_p, wRR_p, ${\stackrel{˅}{D}}_p$, and IQR_p extracted from patient p’s viral load pattern. The words sample or point are also used here RVL (black; n= 237) and HVLS (purple; n=316) clusters. interchangeably. The term feature (F) can be thought of as a column vector for all patients in the dataset consisting of the four attributes: F_Ȧ, F_wRR, $F_{\stackrel{˅}{D}}$, and F_IQR. Finally, the terms label assignment, VL pattern membership assignment, patient categorization, and prediction, all refer to the same principle: To assign the most appropriate label which characterizes the viral load pattern of a patient. However, while the principle is the same, the method of assigning such an appropriate label differs depending on the categorization or the learning method used.

Feature extraction and normalization

We began by transforming viral load data by min-max normalization²² to equally weight the temporal features of the VL series (Equation 12). That is, we normalize the features, F, to a range between [0, 1] using Equation 12 where F^* = f(F).

Next, we examined each of the four features for all patients with ≥ 3 viral load measurements (N = 1,576 patients), and did not find distinct bi-variate clustering (Figure S2). A feature correlation coefficient analysis (Supplementary Table S1) revealed that the Adj MD feature is linearly independent of Area and wRR. In contrast, there is modest linear dependence between IQR and Adj MD, and between Area and both wRR and IQR. As expected, the largest linear dependency is between wRR and IQR. These results suggest the separation between viral load patterns will be most noticeable between the Area and the Adj MD features - as we designed them to be. Also, although Adj MD is dependent upon wRR, we find that their correlation coefficient is very low (0.033).

Hierarchical clustering

We then performed hierarchical clustering of the individual subject VL patterns using a Euclidean distance metric and Ward’s linkage criterion²³ to minimize the total within-cluster variance. Patients showed a clear separation into 5 distinct groups, which had clinical significance (Figure 2 and Figure 3). The cluster with the lowest viral loads and the highest weighted recency reliability (n=442) corresponds to the DSVL patient group. The patients corresponding to the SHVL group (orange; n=46) exhibited the highest relative Area and very low IQR. Compared to the DSVL cluster, the blue cluster (n=535) has slightly greater area and IQR with a significant difference in the weighted recency reliability. Using this information, along with the general patterns shown by Figure 3, we identify this as the SLVL group. The algorithm also identifies the RVL (black; n=237) and HVLS (purple; n=316) clusters. The RVL cluster has a low weighted recency reliability and high IQR. In contrast, the HVLS cluster has a lower area, higher weighted recency reliability, indicating little variation in the terminal portion of the VL time series, and most importantly very low adjusted maximal differences (Figure 4).

Figure 2. Dendrogram of hierarchically clustered patients.

Clustered using the Euclidean distance along with Ward’s method. Numbers on the bottom axis show number of patients in each cluster. The corresponding viral load pattern plots can be found in Figure 3. DSVL = Durably Suppressed Viral Load, SHVL = Sustained High Viral Load, SLVL = Sustained Low Viral Load, RVL = Rebounding Viral Load, HVLS = High Viral Load Suppression

Figure 3. Extracted patient viral load patterns.

For each cluster categorization of the patient from Figure 2, the days since first viral load measurement are plotted against the viral load counts. The points on the plots indicate the last viral load measurement.

Figure 4. Feature segregation from hierarchical clustering.

Each patient is colored corresponding to the results from the hierarchical clustering in Figure 2. The artificial line of points is a result of the grounding function used in Adj MD. Area = relative area of viral load exposure, wRR= weighted recency reliability, IQR = interquartile range, AdjMD = adjusted maximal viral load difference.

VL patterns are similar within clusters, and dissimilar between clusters Figure 3. Interestingly, there are patients within each cluster whose last VL measurement occurs near 1,827 days. This is equivalent to the full span of five years of VL monitoring data set. This suggests that these clusters don’t disappear after some elapsed time, but rather each type of pattern can be found at virtually any time point.

We found large VL spikes within the time series of the HVLS group. We hypothesize that this may be due to the asynchronous timing of measurements between subjects, the natural variation in biological responses, or patient variability in adherence to therapy. This observation also reflects one limitation of asynchronous outcomes data sampling, which lacks a “completion" endpoint characteristic of most prospective, randomized clinical trials. If measurements ended at a spike, the adjusted maximal difference feature may be weighted in the favor of the patient being classified as RVL. This may indicate that some patients classified as suppressing their viral loads should have been classified as having rebounding viral loads. Alternatively, may indicate that these features do not restrict a patient to forever to one category, but allow for dynamic classification as a function of biological or therapeutic responses.

Comparison of categorization methods

Using the same data set, we next compared our VL pattern categorization method to those previously published in the literature. Visually, we find that the SLVL group detected by our method is very similar to the LLVR group defined by Greub et al. (Figure 5). Furthermore, it appears that the methods trying to capture SHVL, viral rebound, and viral failure patients did not succeed as well as the identification of SHVL and RVL patients in our method. RMVL repeat continuous visually appears to have performed very well in identifying patients whom have suppressed their viral loads. However, the results suggest that our analysis performs slightly better in identifying the suppression group (HVLS), as we find that the last VL measurements (black dots in Figure 5) are consistently low using our method.

Figure 5. Comparison of patient categories with existing methods.

A 2D binning of VLM counts for every patient category. Each row uses a different categorization method, and the method name is located to the right of the row, and the title of each subplot is the category assigned by the indicated method. The columns of each 2D bin are normalized based on the maximum number of logged viral load measurement (VLM) counts in the column: log₁₀[1 + V L M ]. Bin color for a count of 0 is copper, and other bin colors range from white to teal (the maximum of the log₁₀[V L M counts] in the column of the bin). The black dots represent the last viral load measurement for the patient (opacity ≥ 0.3; 2D bins have variable opacity for the dots). The bottom row is our analysis is the same as Figure 3, but represented as a 2D bin. DSVL = durably suppressed viral load; LLVR = low level viral rebound; SLVL = sustained low viral load; SHVL = sustained high viral load; HVLS = high viral load suppression, RVL = rebounding viral load.

The other methods may not have performed as well as they rely on a window or a consecutive pair measure, which may be too subjective for assigning VL pattern membership. Furthermore, notice that patients with baseline VL<200 (Figure 5) contain VL patterns which can reach as high as 10⁶ copies/ml, which is in contrast to Rose et al.’s assumption that these patients have consistently low viral variation. Lastly we wish to emphasize that while some of these categorization methods are successful in identifying a specific group of patients, our method is unique as it attempts to associate each VL pattern to a specific category, without using categories such as “Not Suppressed”, “Unspecified”, or “Omitted”.

Supervised learning of VL patterns

We next used the classes identified by hierarchical clustering to compare several machine learning models, with the goal of identifying methods that could be trained to prospectively assign HIV patients to VL categories (i.e. SHVL, SVL, SLVL, DSVL, and RVL). Unsupervised learning (e.g. hierarchical clustering) is useful for establishing the data structure of VL categories and their locations in the feature vector space. Once the model is established (e.g. cluster boundaries), supervised learning methods are better suited for prospective cluster assignment, given a robust "ground truth" for model training, as they do not depend on re-analysis of the entire population.

To this end, we compared the predictive power of several supervised learning methods for HIV cluster assignment, including: k-nearest neighbors (kNN), decision tree, support vector machine (SVM), Adaboost, and random forests. Models were trained on the original data set, and we then ranked their prediction power by their average F₁ score derived by leave-one-out cross-validation (LOOCV) on the clustered results (Table 3). We compared the ability of these methods to reconstruct the originally identified clusters, even when allowing for variability in cluster numbers (e.g. kNN with k={7, 9}, or DT without a maximum depth specification). All methods performed comparably well, with the notable exception of Adaboost. This generally high performance was expected because the VL pattern categories are well-separated as a result of the clustering. k-Nearest Neighbors and k=5, was computationally efficient and yielded the best results in Table 3.

We next considered the trade-off of predictive precision versus model interpretability. Critical clinical evaluation of machine learning results is important to protect against mis-categorization and clinical error. For this reason, many have advocated using models that are more clinically interpretable. kNN is dependent upon the entire training set for prediction, as it does not inherently “learn” patterns²⁴, hence it does not meet our interpretability criteria. In comparison, SVM offers a simpler model, but it’s results could be non-intuitive for clinicians. And although Decision Trees offer the best interpretability, overly complex trees may be generated, as occurred in our study (Figure S3).

We found that pruned decision tree rules, with a maximum depth of 5 levels, met this interpretable criteria, however at a slight cost to the predictive power (Table 3). The extracted decision rules are shown in Figure 6. Each category has a rule with a high proportion of true positive samples following the rule relative to all samples for the category (support). Similarly, a high proportion of the predicted class was found in the rule (precision), indicating that the rules can be summarized into a majority rule. Note that the sum of the support does not necessarily add up to one for each class because some samples belonging to that class may have been otherwise placed into a different rule, making the precision of that rule weaker.

Figure 6. Extracted rules from pruned decision tree and polyhedral CM rule region.

Support is the fraction of true positives satisfying the rule relative to all samples of the class. Precision is the proportion of true positives versus all positives in the rule. Rules are sorted in order of application, first by the level of the decision tree depth (Depth), and then by descending precision. The colored regions represent the the values for which the rule holds (rule feature space). For the centroid method (CM; shaded gray) bounds were calculated by the polyhedron method, where the rectangular bar is the center and the radius is the area inside the parentheses. Area = relative area of viral load exposure, wRR= weighted recency reliability, IQR = interquartile range, AdjMD = adjusted maximal viral load difference.

As an alternative interpretable model we explored the use of centroid cluster summarization, which is often used in clustering algorithms, and is flexible enough to accommodate different centroid determination methods^22,25. To determine the effects of different centroid determination algorithms, we compared seven different methods: multidimensional mean, multidimensional median, best representative center, bounding box method, smallest disk method, polyhedral center, and a novel “push and pull" (PnP) method inspired by force-directed graph drawing such as the Fruchterman-Reingold’s algorithm^26,27 (see Supplementary File 1). Force directed clustering methods maximize inter-cluster center distances, while minimizing intra-cluster distance, and are the basis for modularity clustering in graph theory²⁸.

We then combined the centroid cluster summarization approach with a radius-based classification prediction algorithm. Let c_i be the ith cluster center with corresponding radius r_i, where r_i is calculated as the distance to the farthest intra-cluster sample from c_i, then for a new sample s choose its predicted cluster membership j such that $\frac{{\Vert s-c_j\Vert }_2}{r_j}$ is a minimum. We refer to this method as radial normalization classification.

Comparing the representative F₁ power of the centroid radial normalization methods (italicized in Table 3) to common machine learning algorithms, we find that the centroid interpretation loses some predictive power. However, the centroid summary is highly interpretable because the entire model can be expressed concisely (Supplementary Table S2), and understood clearly. For example, a clinician classifying a patient by VL time series values would compare observed feature values with the ranges given in Figure 6, and find which classification the patient’s data fits best within. In the case of the centroid method, if an observed value appears to fall in multiple categories, then they should be assigned to the one closest to the center (this allows a clinician to cross-check model predictions).

Temporal state variation

HIV patient viral load states are often fluid, with class changes (e.g. SHVL → HVLS) occurring due to therapy, viral genetics, social and other factors. To examine this aspect of classes, we use the k-Nearest Neighbors (k=5) model, fit to the original clusters, to predict the class state of each patient with ≥ 3 VL measurements using only partially retained VL data. For example if a patient has 6 viral load measurements, then we predict the class state at 3, 4, 5, and 6 VLM, which may yield SHVL → SHVL → RVL → HVLS as its prediction. We then constructed a state-transfer network using the trace-route method²⁹, revealing several interesting relationships:

Figure 7. Class state variation.

A) Classification using kNN, with k=5, trained on the original five clusters to predict on partially retained viral load for patients ≥ 900 days of data. The number of patients in one class between 0–900 days are shown relative to the first state classification (i.e. third viral load measurement). B) A trace-route map of class state transfers (class₁ → class₂) as a function of partially retained viral load derived from model. Nodes represent viral load classification and arrows reflect the volume of state transitions between successive VL measurements (e.g. SHVL→DSVL). Self-loops (e.g. RVL→RVL) indicate no change in state reflecting stable classification.