Title: PREDICTING MEDICATIONS FROM DIAGNOSTIC CODES WITH RECURRENT NEURAL NETWORKS

Abstract: It is a surprising fact that electronic medical records are failing at one of their primary purposes, that of tracking the set of medications that the patient is actively taking. Studies estimate that up to 50% of such lists omit active drugs, and that up to 25% of all active medications do not appear on the appropriate patient list. Manual efforts to maintain these lists involve a great deal of tedious human labor, which could be reduced by computational tools to suggest likely missing or incorrect medications on a patient's list. We report here an application of recurrent neural networks to predict the likely therapeutic classes of medications that a patient is taking, given a sequence of the last 100 billing codes in their record. Our best model was a GRU that achieved high prediction accuracy (micro-averaged AUC 0.93, Label Ranking Loss 0.076), limited by hardware constraints on model size. Additionally, examining individual cases revealed that many of the predictions marked incorrect were likely to be examples of either omitted medications or omitted billing codes, supporting our assertion of a substantial number of errors and omissions in the data, and the likelihood of models such as these to help correct them.

1. INTRODUCTION:

The idea of exploiting the large amounts of data captured in electronic medical records for both clinical care and secondary research holds great promise, but its potential is weakened by errors and omissions in those records (Safran et al., 2007; de Lusignan & van Weel, 2006) . Among many other problems, accurately capturing the list of medications currently taken by a given patient is extremely challenging (Velo & Minuz, 2009) . In one study, over 50% of electronic medication lists contained omissions (Caglar et al., 2011) , and in another, 25% of all medications taken by patients were not recorded (Kaboli et al., 2004) . Even medication lists provided by the patients themselves contain multiple errors and omissions (Green et al., 2010) .

Many efforts have been made to ensure the correctness of medication lists, most of them involving improved communication between patients and providers (Keogh et al., 2016) , but these efforts have not yet been successful, and incorrect or incomplete medication documentation continues to be a source of error in computational medical research. In this work we attempt to identify likely errors and omissions in the record, predicting the set of active medications from the sequence of most recent disease-based billing codes in the record. Predictions from such a model could be used either in manual medication reconciliation (a common process undertaken to correct the medication record) or to provide a prior to other models, such as an NLP model attempting to extract medication use from the narrative clinical text.

Given the sequential nature of clinical data, we suspected that recurrent neural networks would be a good architecture for making these predictions. In this work we investigate this potential, comparing the performance of recurrent networks to that of similarly-configured feed forward networks.

The input for each case is a sequence of ICD-9 billing codes (Section 2.1), for which the model produces a single, multi-label prediction of the therapeutic classes (Section 3.1) of medications taken by the patient during the period of time covered by the billing code sequence. This work is designed to test how well the complete set of medications a patient is actively taking at a given moment can be predicted by the sequence of diagnostic billing codes leading up to that moment, in the context of non-trivial label noise. It also explores whether sequence-oriented recursive neural nets can do a better job of that prediction than standard feed-forward networks.

2. BACKGROUND:

2.1 MEDICAL BILLING CODES Each time a patient has billable contact with the healthcare system, one or more date-stamped billing codes are attached to the patient record, indicating the medical conditions that are associated (or suspected to be associated) with the reason for the visit. While these codes are notoriously unreliable because they are only used for billing and not actual clinical practice (O'Malley et al., 2005) , they are nevertheless useful in a research context (Bastarache & Denny, 2011; Denny et al., 2010) , especially if they are used probabilistically (Lasko, 2014) . In our institution, codes from the International Classification of Diseases, Ninth Revision (ICD-9) have historically been used, although we have recently transitioned to the tenth revision (ICD-10). For this project, we used ICD-9 codes.

The ICD-9 hierarchy consists of 21 chapters roughly corresponding to a single organ system or pathologic class (Appendix B). Leaf-level codes in that tree represent single diseases or disease subtypes. For this project, we used a subset of the two thousand most common leaf-level codes as our input data.

2.2. RECURRENT NEURAL NETWORKS AND VARIATIONS:

Most of the ICLR community are very familiar with recurrent neural networks and their variations, but we include a conceptual description of them here for readers coming from other fields. More thorough descriptions are available elsewhere (Graves, 2012; Olah, 2015) .

A recurrent neural network is a variation in which the output of one node on input x t loops around to become an input to another node on input x t+1 , allowing information to be preserved as it iterates over an input data sequence (Figure 1 ). They were introduced in the 1980s (Rumelhart et al., 1986) , but achieved explosive popularity only recently, after the development of methods to more reliably capture long-term dependencies, which significantly improved their performance on sequence-tosequence mapping (Hochreiter & Schmidhuber, 1997; Sutskever et al., 2014) .

The basic RNN unit has a simple internal structure (Figure 2a ). Output from the previous iteration h t-1 and the next input in a sequence x t are both fed to the network on the next iteration. The Long Short-Term Memory configuration (LSTM) introduces new, more complex internal structure (Figure 2b ) consisting of four neural network layers and a cell state (c t ), which is carried from one iteration to another. The additional layers form forget, input and output gates, which allow for the information to be forgotten (reset) or passed on to varying degrees.

The LSTM model and its variations are commonly used in applications where sequence and temporal data are involved, such as in image captioning (Vinyals et al., 2014) , language translation (Sutskever et al., 2014) , and speech recognition (Graves et al., 2013) . In many cases LSTM models define the state of the art, such as with a recent conversational speech recognizer that (slightly) outperforms professional transcriptionists (Xiong et al., 2016) .

A recent variation on the LSTM architecture is the Gated Recurrent Unit (GRU) (Cho et al., 2014) , which introduces a single update gate in place of input and forget gates (Figure 2c ). GRUs perform as well as or better than LSTMs in many cases (Chung et al., 2014; Jozefowicz et al., 2015) , and have the additional advantage of a simpler structure.

In this work we try both an LSTM and a GRU on our learning problem.

2.3. RELATED WORK:

Little research in the computational medical domain has used recurrent neural networks. The earliest example we are aware of is the use of an LSTM model that produced reasonable accuracy

x recurrent unit h x 1 recurrent unit h 1 x 2 recurrent unit h 2 x 3 recurrent unit h 3 x t recurrent unit h t h 1 h 2 h 3

Figure 1 : Simplified representation of a recurrent neural network (left) and an unrolled recurrent neural network (right). x i is a single element in an input sequence x, h i is an output after a single pass through the recurrent unit. Adapted from Olah (2015) .

⨯ ⨯ ⨯ forget tanh input output tanh x t + h t-1 c t-1 h t c t h t tanh x t h t-1 h t h t ⨯ ⨯ output update tanh 1- x t + h t-1 h t h t ⨯ a) b) c)

Figure 2 : Architectures of (a) Simple RNN, (b) LSTM, and (c) GRU units. x t : a single element in an input sequence being considered in the current iteration, h t-1 , h t : the output from the previous and current iterations, c t-1 , c t : the cell states of the previous and current iterations. Adapted from Olah (2015) .

(micro-AUC 0.86) in a 128-dimensional multi-label prediction of diagnoses from regularly sampled, continuously-monitored, real-valued physiologic variables in an Intensive Care Unit setting. This was an interesting initial application, but it turned out to be only 0.001 better than the baseline classifier, which was a multi-layer perceptron with expert-designed features (Lipton et al., 2016) .

Given the dataset size (10,401 patient records) the lack of improvement may have been due to insufficient data to power accurate feature learning in the recurrent network.

Very recent work, contemporary with ours, used a GRU model with a semantic embedding in 32,787 patient records to predict the development of heart failure 3 -6 months in the future, from medication orders and billing codes in an 18-month window. The model achieved respectable accuracy (0.88 AUC), and demonstrated a meaningful 0.05 AUC improvement over a deep feedforward network (Choi et al., 2016b) .

Other recent work from the same group used a GRU model in a multi-label context to predict the medications, billing codes, and time of the next patient visit from a sequence of that same information for previous visits, using 263,706 patient records. It achieved a recall@30 of 72.4 for the task, an improvement of 20 over a single-hidden-layer MLP with 2000 units (Choi et al., 2016a ). This is an example of using one of the strengths of a recurrent network -predicting the next element in a sequence. It contrasts with our work that exploits a different strength of recurrent networkspredicting a sequence or class that is semantically distinct from but parallel to the elements of the input sequence.

The closest work to ours from a medical domain perspective is a series of collaborative filter models (including co-occurrence counting, k-nearest neighbors, and logistic regression) that predict missing medications using a leave-one-drug-out evaluation design, with predictions based on the rest of the medications, ICD-9 billing codes, and demographic data. The models were trained and tested on data from 419 patients in three different clinics, with accuracy varying by clinic, as expected, but not appreciably by model. Most models ranked the missing drug in the top 10 results between 40 and 50% of the time, and ranked the therapeutic class of the drug in the top 10 results between 50 and 65% of the time.

Many aspects of our work can be found in these prior efforts, but none addresses our particular problem in the same way. Our work is unique in its learning problem of identifying all drugs a patient is likely to be taking, based only on the billing codes in the record. Like most others cited, we use recurrent neural networks in a multi-label predictive context, but in contrast to them we compare to the most similar non-recurrent model we can construct, in order to evaluate the contribution of the temporal sequence information to the solution. Finally, we use one to four orders of magnitude more data (3.3 million instances, see Section 3.1) than these prior efforts, which we hope will give us a more realistic assessment of the various deep architectures we use on our problem.

3. EXPERIMENTS:

3.1 DATA Our source database was the deidentified mirror of Vanderbilt's Electronic Medical Record, which contains billing codes, medication histories, laboratory test results, narrative text and medical imaging data for over 2 million patients, reaching back nearly 30 years (Roden et al., 2008) . We obtained IRB approval to use this data in this research.

For this experiment we filtered all records in our database to include only the top 1,000 most common medications and the top m = 2000 most common billing codes, which cover 99.5% of all medication occurrences and 85.1% of all billing code occurrences. We then included all records from the filtered data that had at least one medication occurrence and at least ten billing code occurrences. This resulted in 610,076 complete patient records, which we divided 80/5/15 into training, validation, and final test sets.

A data instance d = {E, T, y} consisted of a sequence E = {e 1 , . . . , e n }, of one-hot billing code vectors e i ∈ {0, 1} m and their associated times T = {t 1 , . . . , t n }, t i ∈ R as input, and a multi-label vector y ∈ {0, 1} k of medication classes as the output target. The most recent n = 100 billing codes to a selected reference time point in a given patient record were collected into the input sequence E, and their occurrence times into T , zero padding if necessary. All medications that occurred during the time span of T were then collected into the output vector y. Practice patterns change over time, so simply taking the most recent 100 codes for each patient could produce a biased result. To avoid this, we chose random reference points, stratified by medication. In other words, the reference points were randomly chosen from the occurrences of each medication in the entire dataset, up to 10,000 points per medication. This resulted in 3.3 million data instances, an average of 5.4 instances per patient record. Each patient's data was included in at most one of the training, validation, or test sets.

Because there are often many approximately equivalent medication choices for a given therapeutic purpose, we converted medication names to their therapeutic class (beta blocker, immunosuppressant, corticosteroid, etc.) as a synonym reduction step. This step also aggregated generic with brand names, as well as different formulations of the same active ingredient. For this task we used the Anatomical Chemical Classification System (ATC)1 , which is a multi-level ontology of medications, organized by both anatomic and therapeutic class. The top level is a broad categorization of medications (Appendix B), the bottom (fifth) level is individual medications, and we used the third level, which contains 287 therapeutic classes of the approximately appropriate abstraction level for our purpose. We used a publicly available mapping2 to translate between our medication names and ATC codes, with manual mapping for the minority of medications that had no mapping entry. Our set of medications used k = 182 third-level ATC codes, rendering our output label a 182-elementlong multi-label vector, in which an element is set y i = 1 if a medication in that class appeared in the set of medications identified for that instance, y i = 0 otherwise. Some medications mapped to more than one class, and we set y i = 1 for all of them.

Our medication data was collected from structured order entry records and extracted using NLP (Xu et al., 2010) from mentions in the narrative text of a patient record that included the medication name, dose, route and frequency. As discussed above, we assumed (and our results demonstrate) that the medication data is incomplete, and our hope was that a model learned from a sufficiently large dataset will be robust to the missing data.

This configuration represents the input billing codes in a sequence, but the output medications as a multi-label vector. This is because ICD-9 codes are represented sequentially in our source data, but medications are not. They are represented as a list that changes over time in the record. The usual goal of clinicians is to verify the list of medications at each visit, and if omissions or additions are indicated by the patient, to change the list to reflect that. But in the time-constrained reality of clinical practice, this reconciliation happens sporadically, and many clinicians are hesitant to change an entry on the medication list for which they were not the original prescriber, so the timing of the changes in the documentation do not reflect the timing of changes in reality. Therefore we are reduced to predicting a single multi-label vector, representing the medications that the patient probably took during the span of time represented by the input codes. (We actually did attempt some full sequence-to-sequence mappings, with various orderings of the medication sequences, but we did not achieve any promising results in that direction.)

3.2. CLASSIFIERS:

Our main technical goal was to test the performance of recurrent neural networks on this sequencecentric prediction problem. To evaluate the specific gains provided by the recurrent architectures, we compare performance against a fully connected feed-forward network configured as similarly as possible to the recurrent networks, and (as baselines) a random forest and a constant-prevalence model. We discuss the specific configurations of these classifiers in this section.

3.2.1. RECURRENT NEURAL NETWORKS:

We tested both LSTMs and GRUs in this experiment. We configured both architectures to first compute a semantic embedding x i ∈ R b of each input e i vector, before appending the times t i (Figure 3 ) and feeding the result to three layers of recurrent units. The final output from the last pass of recurrent unit is as a multi-label prediction for each candidate medication.

The optimal hyperparameters for the model were selected in the randomized parameter optimization (Bergstra & Bengio, 2012) , with the embedding dimension b = 32, number of layers, and number of nodes optimized by a few trials of human-guided search. Other optimized parameters included the fraction of dropout (between layers, input gates and recurrent connections), and L1 and L2 regularization coefficients (final values are presented in Appendix A).

Both models were implemented using Keras (Chollet, 2015) and trained for 300 iterations using cross-entropy under the Adadelta optimizer (Zeiler, 2012) .

e A code observation e i passes through an embedding layer, producing an embedding vector x i , which is then appended with time t. The processed matrix then passes through either recurrent layers or feed-forward layers. The output in both cases is a single vector y of label probabilities.

3.2.2. FULLY CONNECTED NEURAL NETWORK:

The fully connected network used as similar an architecture as possible to the recurrent networks, in an attempt to isolate the gain achieved from the recurrence property. Specifically, we used the same architecture for embedding and timestamp appending (Figure 3 ).

Hyperparameters were optimized using random search over the number of layers, number of nodes, dropout, activation function between layers, L1 and L2 regularization coefficients (Appendix A).

(Surprisingly, the optimizer chose tanh over ReLU as the optimal activation function.)

The models were also implemented using Keras, and were trained using cross-entropy for 500 iterations under the Adadelta optimizer.

3.2.3. RANDOM FOREST:

Because the random forest model is not easily structured to operate on sequences, we represented the input data as either binary occurrence vectors v ∈ {0, 1} m , or bag-of-codes vectors w ∈ N m (counts of each code value in the sequence) rather than as sequences of codes with associated times.

No embedding was used, because random forest code was not able to cope with the large size of the data in the (dense) embedded space.

Even in the (sparse) original space, the full dataset was too large for the random forest code, so we implemented it as an ensemble of ten independent forests, each trained on one tenth of the training data, and their average score used for test predictions.

Models were implemented using scikit-learn (Pedregosa et al., 2011) with parameters optimized under random search (Appendix A).

While other models could reasonably serve as a baseline for this work, we chose a random forest because they tend to perform well on widely varying datasets (Fernández-Delgado et al., 2014) , they are efficient to train and test, and they don't require a huge effort to optimize (in order to produce a fair comparison).

3.3. CONSTANT-PREVALENCE MODEL:

This minimum baseline model simply predicts the prevalence of each label for all instances. For example, if there were three possible medications, with prevalences of 0.3, 0.9, and 0.2, then the prediction of this model would be a constant [0.3, 0.9, 0.2] for each instance. We include this model in order to mitigate the fact that while all of our evaluation measures are suitable for comparing models on the same data, some are not well suited for external comparison because they depend, for example, on the prevalence of positive labels (Section 3.4). By including this model we can at least establish a true minimum baseline for reference.

3.4. EVALUATION:

Our main evaluation focused on the models, although we also performed a separate evaluation of the embedding.

3.4.1. MODELS:

There are several possibilities for evaluation in a multi-label classification context (Sechidis et al., 2011; Zhang & Zhou, 2014) . We chose micro-averaged area under the ROC curve (AUC) and label ranking loss as the primary methods of evaluation, because they treat each instance with equal weight, regardless of the nature of the positive labels for that instance. In other words, we wanted primary measures that did not give a scoring advantage to instances with either very many or very few positive labels, or that included very rare or very prevalent labels. Additionally, both of these measures appeal to us as intuitive extensions of the usual binary AUC, when seen from the perspective of a single instance. However, because these two measures don't reflect all aspects of multi-label prediction performance, we also include macro-averaged AUC, label ranking average precision and coverage error measures.

Micro-averaged AUC considers each of the multiple label predictions in each instance as either true or false, and then computes the binary AUC as if they all belonged to the same 2-class problem (Zhang & Zhou, 2014) . In other words, micro-averaged AUC A µ is:

EQUATION

where S = {(x, l) : l ∈ Y } is the set of (instance, label) pairs with a positive label, and Y = {y i : y i = 1, i = 1 . . . k} is the set of positive labels for input x.

Label ranking loss L R gives the average fraction of all possible (positive, negative) label pairs for each instance in which the negative label has a higher score than the positive label (Tsoumakas et al., 2010) :

EQUATION

where the superscript (j) refers to the jth test instance (of N total instances) and r(l) is the predicted rank of a label l.

Macro-averaged AUC can be thought of as averaging the AUC performance of several one-vs-all classifiers, one model for each label. It treats each model equally, regardless of the prevalence of positive labels for that model. This gives a score of 0.5 to the constant-prevalence model, at the cost of weighting instances differently in order to achieve that. This is in contrast to micro-averaged AUC, which can be thought of as averaging across instances rather than labels. It weighs each instance equally, at the cost of a 0.5 score no longer being the random-guessing baseline.

Label ranking average precision gives the mean fraction of correct positive labels among all positive labels with lower scores for each label. The coverage error function calculates the mean number of labels on the ranked list that are needed to cover all the positive labels of the sample. Both of these depend on the prevalence of positive labels in a test instance.

3.4.2. EMBEDDING:

We evaluated the embedding based on how strongly related in a clinical semantic sense the nearest neighbor to each code is (in the embedding space).

A licensed physician manually annotated the list of all 2000 codes with its match category m ∈ {strongly related, loosely related, unrelated}, and we computed the empirical marginal probability P (m) of each category, the empirical conditional probability P (m|d) of the match category given the nearest neighbor (Manhattan) distance d and the empirical marginal probability P (d). For comparison, we computed P (m) under 100 random code pairings.

4. RESULTS AND DISCUSSION:

The GRU model had the top performance by all measures, although the LSTM was a close second (Table 1 ), a performance pattern consistent with previous reports (Chung et al., 2014) . The deep neural net performance was about 0.01 worse in both measures, suggesting that the recurrent models were able to use the sequence information, but only to a small advantage over the most similar nontemporal architecture. However, we note that both RNNs' performance peaked at the top end of our tractable range for model size, while the feed-forward network peaked using a model about one third that size (Appendix A). Experimenting with the architecture, we found that increasing the number of nodes or layers for the feed-forward network increased training time but not performance. This suggests that the RNN performance was limited by the hardware available, and increasing the size of the model may further increase performance, and that the feed-forward network was limited by something else.

Both random forest models were weaker than the deep neural net, as might be expected from the need to resort to binary and bag-of-codes representations of the input data. This 0.9/0.01 rule of thumb may loosely translate to our AUC variants, but it can directly translate to Label Ranking Loss L R (2). If we think of a single output prediction ŷ ∈ [0, 1] k as a set of predictions for k binary labels, then 1 -AUC for that set of predictions is equivalent to L R for the original instance ŷ. Therefore, values of L R 0.1 may be clinically useful, and L R 0.05 may be very useful.

Subjectively examining performance on 20 randomly selected cases, we find very good detailed predictions, but also evidence of both missing medications and missing billing codes. An example of a good set of detailed predictions is from a complex patient suffering from multiple myeloma (a type of cancer) with various complications. This patient was taking 26 medications, 24 of which had moderate to high probability predictions (Figure 4 ). (We have found by eyeball that a prediction cutoff of 0.2 gives a reasonable balance between sensitivity and specificity for our model.) In the other direction, only two of the high-prediction classes were not actually being taken, but those classes, along with several of the other moderately-predicted classes, are commonly used for cancer and are clinically reasonable for the case. (Details of this and the two cases below are in Appendix C).

A good example of missing medications is a case in which the record has multiple billing codes for both osteoporosis (which is very commonly treated with medication) and postablative hypothyroidism (a deliberately induced condition that is always treated with medication), but no medications of the appropriate classes were in the record. The GRU model predicted both of these classes, which the patient was almost surely taking.

A good example of either missing billing codes or discontinued medications that remain documented as active is a case in which the record has at least five years of data consisting only of codes for Parkinson's disease, but which lists medications for high cholesterol, hypertension, and other heart disease. The GRU model predicted a reasonable set of medications for Parkinson's disease and its complications, but did not predict the other medications that are not suggested by the record.

Given how easy it was to find cases with apparently missing codes and medications, we conclude that there is indeed a substantial amount of label noise in our data, and we therefore interpret our models' performance as lower bounds on the actual performance. We are encouraged that this kind of a model may actually be useful for identifying missing medications in the record, but of course a more thorough validation, and possibly a more accurate model, would be necessary before using in a clinical scenario. A definitive experiment would use off-line research, including reconciling information from various electronic and human sources to establish the ground truth of which medications were being taken on a particular day, but such efforts are labor intensive and expensive, and can only be conducted on a very small scale.

An interesting byproduct of these models is the semantic embedding of ICD-9 codes used in the recurrent networks (Figure 5 ). Transforming input to a semantic embedding is a common pre- processing step to improve performance, but clearly the semantic understanding it provides to an algorithm can be useful beyond the immediate learning problem (Mikolov et al., 2013) . Investigating the embedding learned in this experiment shows some generalizable potential, but it also reveals the need for further refinement before it can be truly useful. Specifically, while it's easy to find tight groups of ICD-9 codes that are strongly clinically related in our embedding, we also find groups for which we cannot see a meaningful clinical relationship.

For example, we see two groups of codes relating to kidney failure and diabetes mellitus, two classes of very prevalent disease (Figure 5 , insets). In other iterations with different parameter settings, the kidney failure codes were even embedded in a sequence reflecting the natural progression of the disease, with the code for dialysis (an intensive treatment for end-stage kidney failure) embedded at the appropriate place. Interestingly, these were not the parameter settings that optimized overall prediction performance. In other settings, such as our performance-optimal setting, the sequence is close to the natural progression of the disease, but not quite identical. Nevertheless, this is an exciting result that suggests great potential.

Further evaluation of the embedding found that 49% of codes were strongly related semantically to their nearest neighbor, 10% were loosely related, and 41% unrelated. This fraction of strongly related nearest neighbors was lower than we had hoped, but much higher than expected by chance (Figure 6 ), and it definitely improved classification performance. Furthermore, it was obvious by inspection that in general, codes closer in the embedding were more semantically related than distant codes, but interestingly, the distance to the nearest such neighbor showed the opposite relationship -nearest neighbors that were very close were less likely to be semantically related than nearest neighbors that were far, and this trend is roughly linear across the full range of d (Figure 6 ). So the sparser the points are in the embedded space, the more semantically related they are to their nearest neighbor, but the causal direction of that effect and the technical reason for it are beyond the scope of this initial work.

For this prediction problem, we settled on predicting the medications that occurred in the record during the same time span as the billing codes used. Originally, we intended to predict only the medications listed on the day of the reference point, but that turned out to greatly exacerbate the missing medication problem. After trying medications that fell on the reference day only, the week prior to the reference day, and the six months prior, our best performance both subjectively and objectively was achieved using the full time range of the input data.

While the performance of the recurrent networks was quite good, we believe it could be improved by including additional input data, such as laboratory test results, demographics, and perhaps vital signs. We also suspect that if we can devise a way to convert our medication data into reliablyordered sequences, we can more fully exploit the strengths of recurrent networks for medication prediction. We look forward to trying these and other variations in future work.

APPENDIX C.:

This appendix presents results from three illustrative cases from the dozen cases randomly selected for individual evaluation.

CASE 1. Medication predictions for a simpler patient. Note that the high-prediction medications are clinically reasonable given the billing codes in the sequence. Medication predictions for a patient with only one ICD-9 code, repeated many times over five years. The medications listed under true labels are not indicated for paralysis agitans (Parkinson's disease), but the patient was surely taking them for reasons not documented in the ICD-9 sequence. The model predicted mostly reasonable medications for a patient with Parkinson's disease, especially Dopaminergic agents, which is the primary treatment for the disease. Figure representation as in case 1, above.

:

http://www.whocc.no/atc/structure and principles/

https://www.nlm.nih.gov/research/umls/rxnorm/

Figures and Tables:

Figure 3: Recurrent (left) and feed-forward (right) neural network architectures.Arrows indicate the flow of information. Input for both models is sequence of billing code observations e and sequence of corresponding timestamps t. A code observation e i passes through an embedding layer, producing an embedding vector x i , which is then appended with time t. The processed matrix then passes through either recurrent layers or feed-forward layers. The output in both cases is a single vector y of label probabilities.

Figure 4: Medication predictions for a complicated patient. Each vertical bar represents the prediction for a single medication class, with the height of the bar representing the confidence of the prediction. Black labels with arrows indicate ATC therapeutic classes for medications the patient was actually taking. Colors and letters below the axis indicate organ system groups. More detail in Appendix C.

Figure6: Semantic relatedness of nearest neighbors vs. the distance between them. Solid lines are the conditional probabilities P (m|d) for the three values of m, dashed line is the marginal probability P (d) of nearest neighbor distances d. Surprisingly, nearest neighbors that are farther away (but still the nearest neighbor) are more strongly related than nearest neighbors that are closer in the embedding space. Shaded regions, colored to correspond to the three values of m, are the 95% CI for empirically estimated P (m) under random pairings, and represent the expected null result.

actual medication classes for Case 3.

Results of multi-label classification for each model. Baseline is the constant-prevalence model. Perfect is the best possible performance for our data under the given measure. is what performance is good enough for clinical use. While there is little clinical experience with multi-label classifiers, we would generally expect clinicians using a binary classifier in an advisory role to find an AUC 0.9 to be useful, and AUC 0.95 to be very useful. An AUC difference of 0.01, and perhaps 0.005 are potentially noticeable in clinical use.

Table content: <table><tr><td/><td/><td>Label Ranking</td><td/><td colspan="2">Label Ranking Coverage</td></tr><tr><td>Model</td><td>Micro-AUC</td><td>Loss</td><td colspan="2">Macro-AUC Avg. Precision</td><td>Error</td></tr><tr><td>GRU</td><td>0.927</td><td>0.076</td><td>0.861</td><td>0.603</td><td>62.6</td></tr><tr><td>LSTM</td><td>0.926</td><td>0.077</td><td>0.859</td><td>0.600</td><td>63.0</td></tr><tr><td>NN</td><td>0.916</td><td>0.086</td><td>0.835</td><td>0.570</td><td>67.3</td></tr><tr><td>RF (binary)</td><td>0.903</td><td>0.102</td><td>0.804</td><td>0.523</td><td>73.7</td></tr><tr><td>RF (counts)</td><td>0.894</td><td>0.111</td><td>0.787</td><td>0.497</td><td>77.3</td></tr><tr><td>Baseline</td><td>0.828</td><td>0.172</td><td>0.500</td><td>0.355</td><td>97.2</td></tr><tr><td>Perfect</td><td>1.0</td><td>0.0</td><td>1.0</td><td>1.0</td><td>15.0</td></tr><tr><td colspan="2">A natural question</td><td/><td/><td/><td/></tr></table>

Figure5: A t-SNE representation of our final embedding. The insets highlight two groups of codes (diabetes mellitus and kidney failure) that are strongly related clinically, and a third group that is not. Codes are colored by whether their nearest neighbor in the embedding space (which may be different from the nearest neighbor in this t-SNE space) is strongly related (blue), loosely related (orange), or unrelated (gray) from a clinical perspective.

Table content: <table><tr><td colspan="2">585.9 Chronic kidney disease, unspecified</td></tr><tr><td colspan="2">585.3 Chronic kidney disease, Stage III (moderate)</td></tr><tr><td colspan="2">585.6 End stage renal disease</td></tr><tr><td colspan="2">585.4 Chronic kidney disease, Stage IV (severe)</td></tr><tr><td>V45.11 Renal dialysis status</td><td/></tr><tr><td colspan="2">585.5 Chronic kidney disease, Stage V</td></tr><tr><td colspan="2">V45.1 Postsurgical renal dialysis status</td></tr><tr><td colspan="2">285.21 Anemia in chronic kidney disease</td></tr><tr><td colspan="2">787.1 Heartburn</td></tr><tr><td colspan="2">727.00 Synovitis and tenosynovitis</td></tr><tr><td colspan="2">309.24 Adjustment disorder with anxiety</td></tr><tr><td colspan="2">831.00 Closed dislocation of shoulder</td></tr><tr><td>724.3 Sciatica</td><td>701.4 Keloid scar</td></tr><tr><td colspan="2">250.81 Diabetes with other specified manifestations, type I …</td></tr><tr><td colspan="2">362.01 Background diabetic retinopathy</td></tr><tr><td colspan="2">250.80 Diabetes with other specified manifestations, type II or unspecified, … 250.40 Diabetes with renal manifestations, type II or unspecified, …</td></tr><tr><td colspan="2">250.50 Diabetes with ophthalmic manifestations, type II or unspecified, …</td></tr><tr><td colspan="2">250.42 Diabetes with renal manifestations, type II or unspecified, uncontrolled</td></tr><tr><td colspan="2">250.01 Diabetes without complication, type I</td></tr><tr><td colspan="2">250.62 Diabetes with neurological manifestations, type II or unspecified,</td></tr><tr><td colspan="2">uncontrolled</td></tr><tr><td colspan="2">250.60 Diabetes with neurological manifestations, type II or unspecified …</td></tr><tr><td colspan="2">357.2 Polyneuropathy in diabetes</td></tr></table>

Medication predictions for a complicated patient. Each vertical bar represents the prediction for a single medication class, with the height of the bar representing the confidence of the prediction. Black labels above arrows indicate ATC therapeutic classes for medications the patient was actually taking. Colors and letters below the axis indicate high-level therapeutic class groups.Predicted vs. actual medication classes for the patient in Case 1. The four-character sequence in the first and fourth columns is the ATC code for the medication therapeutic class, and an asterisk in the first column indicates that the predicted medication is in the actual medication list. Probabilities listed are the model predictions for the listed therapeutic class. In the predicted medications column, all predictions with probability at least 0.2 are listed.Predicted vs. actual medication classes for Case 2. Table structure as in case 1.

Table content: <table><tr><td>CASE 2.</td><td/><td/><td/><td/><td/></tr><tr><td/><td/><td>ICD-9 code</td><td colspan="3">Code description</td><td>Time estimate (ago)</td></tr><tr><td/><td/><td>735.4</td><td colspan="2">Other hammer toe (acquired)</td><td/><td>2.4 years ago</td></tr><tr><td/><td/><td>729.5</td><td>Pain in limb</td><td/><td/><td>2.4 years ago</td></tr><tr><td/><td/><td>244.1</td><td colspan="2">Other postablative hypothyroidism</td><td/><td>1.5 years ago</td></tr><tr><td/><td/><td>285.9</td><td>Anemia, unspecified</td><td/><td/><td>1.5 years ago</td></tr><tr><td/><td/><td>244.1</td><td colspan="2">Other postablative hypothyroidism</td><td/><td>1.2 years ago</td></tr><tr><td/><td/><td>244.1</td><td colspan="2">Other postablative hypothyroidism</td><td/><td>11.5 months ago</td></tr><tr><td colspan="2">ICD-9 code</td><td>733.00 733.01</td><td colspan="3">Code description Osteoporosis, unspecified Senile osteoporosis</td><td>11.5 months ago Time estimate (ago) 7.7 months ago</td></tr><tr><td/><td>203.00 273.1 285.9 276.50 733.00 203.00 203.00 203.01 273.1 273.1 279.3 203.00 781.2</td><td colspan="5">Multiple myeloma, without mention of having achieved remission 268.9 Unspecified vitamin D deficiency Monoclonal paraproteinemia 729.5 Pain in limb Anemia, unspecified 174.9 Malignant neoplasm of breast (female), unspecified Volume depletion, unspecified 722.52 Degeneration of lumbar or lumbosacral intervertebral disc Osteoporosis, unspecified 279.3 Unspecified immunity deficiency myeloma, without mention of having achieved remission 733.01 Senile osteoporosis Multiple myeloma, without mention of having achieved remission 733.01 Senile osteoporosis Multiple myeloma, in remission 244.1 Other postablative hypothyroidism Monoclonal paraproteinemia 401.1 Benign essential hypertension Monoclonal paraproteinemia V58.69 Long-term (current) use of other medications Unspecified immunity deficiency 733.01 Senile osteoporosis Multiple myeloma, without mention of having achieved remission 244.1 Other postablative hypothyroidism Abnormality of gait V58.69 Long-term (current) use of other medications</td><td>7.7 months ago 7.7 months ago 7.7 months ago 7.7 months ago 7.7 months ago 6.4 months ago 6.2 months ago 6.0 months ago 6.0 months ago 1.9 weeks ago now now now</td><td>4.8 months ago 4.8 months ago 4.8 months ago 4.8 months ago 4.8 months ago 4.8 months ago 2.9 months ago 2.9 months ago 2.9 months ago 1.6 months ago 1.6 months ago 1.6 months ago 3.7 weeks ago</td></tr><tr><td/><td>203.00</td><td colspan="4">Multiple myeloma, without mention of having achieved remission</td><td>3.7 weeks ago</td></tr><tr><td/><td>401.9 V12.54</td><td colspan="5">Unspecified essential hypertension Personal history of transient ischemic attack (TIA), and cerebral infarction without residual deficits Top predictions Prob. True labels</td><td>3.7 weeks ago 3.7 weeks ago</td><td>Prob.</td></tr><tr><td>S03B* S01C* S02B* L01A D07X* H02A* D07A* S01B* M05B D10A* H03A H05A C05A* A11C* N02B A04A A01A R01A* A12A N06A* N02A* A01A* J05A* C07A</td><td colspan="5">794.31 786.09 Corticosteroids Nonspecific abnormal electrocardiogram [ECG] [EKG] Other respiratory abnormalities 97.01% 273.1 Monoclonal paraproteinemia Antiinflammatory agents and antiinfectives in combi-95.54% 203.00 Multiple myeloma, without mention of having achieved remission S01C nation V58.69 Long-term (current) use of other medications Corticosteroids 95.54% S02B 794.31 Nonspecific abnormal electrocardiogram [ECG] [EKG] Alkylating agents 94.00% D07X 203.00 Multiple myeloma, without mention of having achieved remission Corticosteroids, other combinations 93.37% H02A V42.82 Peripheral stem cells replaced by transplant Corticosteroids for systemic use, plain 91.06% D07A Top predictions Prob. 203.01 Multiple myeloma, in remission Corticosteroids, plain 90.83% S01B 38.97 Central venous catheter placement with guidance Antiinflammatory agents 90.79% D10A Drugs affecting bone structure and mineralization 88.18% A11C V42.82 Peripheral stem cells replaced by transplant Anti-acne preparations for topical use 88.56% C05A Thyroid preparations 84.82% N06A V58.81 Fitting and adjustment of vascular catheter Parathyroid hormones and analogues 66.33% C10A 203.00 Multiple myeloma, without mention of having achieved remission Agents for treatment of hemorrhoids and anal fissures 88.52% R01A Vitamin a and d, incl. combinations of the two 39.42% N03A V42.82 Peripheral stem cells replaced by transplant for topical use Other analgesics and antipyretics 37.58% C09C 203.01 Multiple myeloma, in remission Antiemetics and antinauseants 87.95% J05A Stomatological preparations 23.05% L02B 203.00 Multiple myeloma, without mention of having achieved remission Decongestants and other nasal preparations for topi-87.02% A01A Calcium 21.59% V42.82 Peripheral stem cells replaced by transplant Antidepressants 20.88% cal use 203.00 Multiple myeloma, without mention of having achieved remission V42.82 Peripheral stem cells replaced by transplant Opioids 84.86% A12C Stomatological preparations 86.11% B05C Direct acting antivirals 86.8% N02A Beta blocking agents 20.81%</td><td>3.7 weeks ago 3.7 weeks ago Antiinflammatory agents and antiinfectives in combi-Corticosteroids 3.7 weeks ago nation 3.6 weeks ago Corticosteroids 3.6 weeks ago Corticosteroids, other combinations 3.4 weeks ago Corticosteroids for systemic use, plain 4 days ago Corticosteroids, plain True labels 4 days ago Antiinflammatory agents 3 days ago Anti-acne preparations for topical use Vitamin a and d, incl. combinations of the two 3 days ago Agents for treatment of hemorrhoids and anal fissures Antidepressants 3 days ago Lipid modifying agents, plain for topical use 3 days ago Decongestants and other nasal preparations for topi-Antiepileptics 3 days ago Angiotensin ii antagonists, plain cal use 2 days ago Direct acting antivirals Hormone antagonists and related agents 2 days ago Stomatological preparations 1 day ago 1 day ago Other mineral supplements now Irrigating solutions now Opioids</td><td>97.01% 95.54% 95.54% 93.37% 91.06% Prob. 90.83% 90.79% 39.42% 88.56% 20.88% 88.52% 17.05% 15.61% 87.02% 10.38% 4.22% 86.83% 86.11% 79.50% 82.56% 84.86%</td></tr><tr><td>B05C*</td><td colspan="2">Irrigating solutions</td><td/><td>82.56%</td><td>B05X</td><td>I.V. solution additives</td><td>74.84%</td></tr><tr><td>A12C*</td><td colspan="2">Other mineral supplements</td><td/><td>79.50%</td><td>L04A</td><td>Immunosuppressants</td><td>68.76%</td></tr><tr><td>B05X*</td><td colspan="2">I.v. solution additives</td><td/><td>74.84%</td><td>N02B</td><td>Other analgesics and antipyretics</td><td>57.24%</td></tr><tr><td>L04A*</td><td colspan="2">Immunosuppressants</td><td/><td>68.76%</td><td>S01A</td><td>Antiinfectives</td><td>54.59%</td></tr><tr><td>N05A</td><td colspan="2">Antipsychotics</td><td/><td>58.64%</td><td>J01D</td><td>Other beta-lactam antibacterials</td><td>43.40%</td></tr><tr><td>N02B*</td><td colspan="2">Other analgesics and antipyretics</td><td/><td>57.24%</td><td>C03C</td><td>High-ceiling diuretics</td><td>39.88%</td></tr><tr><td>S01A*</td><td colspan="2">Antiinfectives</td><td/><td>54.59%</td><td>J01M</td><td>Quinolone antibacterials</td><td>29.78%</td></tr><tr><td>L03A</td><td colspan="2">Immunostimulants</td><td/><td>45.96%</td><td>C07A</td><td>Beta blocking agents</td><td>27.08%</td></tr><tr><td>A02B</td><td colspan="3">Drugs for peptic ulcer and gastro-oesophageal reflux</td><td>44.56%</td><td/></tr><tr><td/><td>disease</td><td/><td/><td/><td/></tr><tr><td>J01D*</td><td colspan="2">Other beta-lactam antibacterials</td><td/><td>43.40%</td><td>N03A</td><td>Antiepileptics</td><td>20.00%</td></tr><tr><td>C03C*</td><td colspan="2">High-ceiling diuretics</td><td/><td>39.88%</td><td>J01X</td><td>Other antibacterials</td><td>5.88%</td></tr><tr><td>B01A</td><td colspan="2">Antithrombotic agents</td><td/><td>37.80%</td><td>M03B</td><td>Muscle relaxants, centrally acting agents</td><td>5.09%</td></tr><tr><td>V03A</td><td colspan="2">All other therapeutic products</td><td/><td>34.18%</td><td/></tr><tr><td>R06A</td><td colspan="2">Antihistamines for systemic use</td><td/><td>31.78%</td><td/></tr><tr><td>A06A</td><td colspan="2">Drugs for constipation</td><td/><td>31.57%</td><td/></tr><tr><td>J01M*</td><td colspan="2">Quinolone antibacterials</td><td/><td>29.78%</td><td/></tr><tr><td>N05B</td><td>Anxiolytics</td><td/><td/><td>29.42%</td><td/></tr><tr><td>D04A</td><td colspan="3">Antipruritics, incl. antihistamines, anesthetics, etc.</td><td>27.62%</td><td/></tr><tr><td>C07A*</td><td colspan="2">Beta blocking agents</td><td/><td>27.08%</td><td/></tr><tr><td>L01X</td><td colspan="2">Other antineoplastic agents</td><td/><td>24.72%</td><td/></tr><tr><td>R05C</td><td colspan="3">Expectorants, excl. combinations with cough sup-</td><td>20.43%</td><td/></tr><tr><td/><td>pressants</td><td/><td/><td/><td/></tr><tr><td>N03A*</td><td colspan="2">Antiepileptics</td><td/><td>20.00%</td><td/></tr></table>

Table structure as in case 1.

Table content: <table><tr><td/><td>Top predictions</td><td>Prob.</td><td/><td>True labels</td><td>Prob.</td></tr><tr><td>N04B</td><td>Dopaminergic agents</td><td>97.66%</td><td>C10A</td><td>Lipid modifying agents, plain</td><td>13.90%</td></tr><tr><td>N03A</td><td>Antiepileptics</td><td>34.01%</td><td>C09A</td><td>Ace inhibitors, plain</td><td>9.21%</td></tr><tr><td>N02B</td><td>Other analgesics and antipyretics</td><td>32.81%</td><td>C01E</td><td>Other cardiac preparations</td><td>5.56%</td></tr><tr><td>N06A</td><td>Antidepressants</td><td>26.10%</td><td>C02C</td><td>Antiadrenergic agents, peripherally acting</td><td>0.72%</td></tr><tr><td>N02A</td><td>Opioids</td><td>20.33%</td><td>G03B</td><td>Androgens</td><td>0.32%</td></tr><tr><td/><td/><td/><td>A14A</td><td>Anabolic steroids</td><td>0.08%</td></tr></table>

Peer review (iclr2017):
In light of the detailed author responses and further updates to the manuscript, I am raising my score to an 8 and reiterating my support for this paper. I think it will be among the strongest non-traditional applied deep learning work at ICLR and will receive a great deal of interest and attention from attendees. ----- This paper describes modern deep learning approach to the problem of predicting the medications taken by a patient during a period of time based solely upon the sequence of ICD-9 codes assigned to the patient during that same time period. This problem is formulated as a multilabel sequence classification (in contrast to language modeling, which is multiclass classification). They propose to use standard LSTM and GRU architectures with embedding layers to handle the sparse categorical inputs, similar to that described in related work by Choi, et al. In experiments using a cohort of ~610K patient records, they find that RNN models outperform strong baselines including an MLP and a random forest, as well as a common sense baseline. The differences in performance between the recurrent models and the MLP appear to be large enough to be significant, given the size of the test set. Strengths: - Very important problem. As the authors point out, two the value propositions of EHRs -- which have been widely adopted throughout the US due to a combination of legislation and billions of dollars in incentives from the federal government -- included more accurate records and fewer medication mistakes. These two benefits have largely failed to materialize. This seems like a major opportunity for data mining and machine learning. - Paper is well-written with lucid introduction and motivation, thorough discussion of related work, clear description of experiments and metrics, and interesting qualitative analysis of results. - Empirical results are solid with a strong win for RNNs over convincing baselines. This is in contrast to some recent related papers, including Lipton & Kale et al, ICLR 2016, where the gap between the RNN and MLP was relatively small, and Choi et al, MLHC 2016, which omitted many obvious baselines. - Discussion is thorough and thoughtful. The authors are right about the kidney code embedding results: this is a very promising result. Weaknesses: - The authors make several unintuitive decisions related to data preprocessing and experimental design, foremost among them the choice NOT to use full patient sequences but instead only truncated patient sequences that each ends at randomly chosen time point. This does not necessarily invalidate their results, but it is somewhat unnatural and the explanation is difficult to follow, reducing the paper*s potential impact. It is also reduces the RNN*s potential advantage. - The chosen metrics seem appropriate, but non-experts may have trouble interpreting the absolute and relative performances (beyond the superficial, e.g., RNN score 0.01 more than NN!). The authors should invest some space in explaining (1) what level of performance -- for each metric -- would be necessary for the model to be useful in a real clinical setting and (2) whether the gaps between the various models are *significant* (even in an informal sense). - The paper proposes nothing novel in terms of methods, which is a serious weakness for a methods conference like ICLR. I think it is strong enough empirically (and sufficiently interesting in application) to warrant acceptance regardless, but there may be things the authors can do to make it more competitive. For example, one potential hypothesis is that higher capacity models are more prone to overfitting noisy targets. Is there some way to investigate this, perhaps by looking at the kinds of errors each model makes? I have a final comment: as a piece of clinical work, the paper has a huge weakness: the lack of ground truth labels for missing medications. Models are both trained and tested on data with noisy labels. For training, the authors are right that this shouldn*t be a huge problem, provided the label noise is random (even class conditional isn*t too big of a problem). For testing, though, this seems like it could skew metrics. Further, the assumption that the label noise is not systemic seems very unlikely given that these data are recorded by human clinicians. The cases shown in Appendix C lend some credence to this assertion: for Case 1, 7/26 actual medications received probabilities < 0.5. My hunch is that clinical reviewers would view the paper with great skepticism. The authors will need to get creative about evaluation -- or invest a lot of time/money in labeling data -- to really prove that this works. For what it is worth, I hope that this paper is accepted as I think it will be of great interest to the ICLR community. However, I am borderline about whether I*d be willing to fight for its acceptance. If the authors can address the reviewers* critiques -- and in particular, dive into the question of overfitting the imperfect labels and provide some insights -- I might be willing to raise my score and lobby for acceptance.