Title: GRAM: GRAPH-BASED ATTENTION MODEL FOR HEALTHCARE REPRESENTATION LEARNING

Abstract: Deep learning methods exhibit promising performance for predictive modeling in healthcare, but two important challenges remain:
• Data insufficiency: Often in healthcare predictive modeling, the sample size is insufficient for deep learning methods to achieve satisfactory results. • Interpretation: The representations learned by deep learning models should align with medical knowledge. To address these challenges, we propose a GRaph-based Attention Model, GRAM that supplements electronic health records (EHR) with hierarchical information inherent to medical ontologies. Based on the data volume and the ontology structure, GRAM represents a medical concept as a combination of its ancestors in the ontology via an attention mechanism. We compared predictive performance (i.e. accuracy, data needs, interpretability) of GRAM to various methods including the recurrent neural network (RNN) in two sequential diagnoses prediction tasks and one heart failure prediction task. Compared to the basic RNN, GRAM achieved 10% higher accuracy for predicting diseases rarely observed in the training data and 3% improved area under the ROC curve for predicting heart failure using an order of magnitude less training data. Additionally, unlike other methods, the medical concept representations learned by GRAM are well aligned with the medical ontology. Finally, GRAM exhibits intuitive attention behaviors by adaptively generalizing to higher level concepts when facing data insufficiency at the lower level concepts. 𝛼 ",$ × 𝐞 " 𝛼 (,$ × 𝐞 ( 𝛼 ),$ × 𝐞 ) 𝛼 $,$ × 𝐞 $ + i-th column of G . . .

1. INTRODUCTION:

The rapid growth in volume and diversity of health care data from electronic health records (EHR) and other sources is motivating the use of predictive modeling to improve care for individual patients. In particular, novel applications are emerging that use deep learning methods such as word embedding (Choi et al., 2016c; e) , recurrent neural networks (RNN) (Che et al., 2016; Choi et al., 2016a; b; Lipton et al., 2016) , convolutional neural networks (CNN) (Nguyen et al., 2016) or stacked denoising autoencoders (SDA) (Che et al., 2015; Miotto et al., 2016) , demonstrating significant performance enhancement for diverse prediction tasks. Deep learning models appear to perform significantly better than logistic regression or multilayer perceptron (MLP) models that depend, to some degree, on expert feature construction (Lipton et al., 2015; Razavian et al., 2016) .

Training deep learning models typically requires large amounts of data that often cannot be met by a single health system or provider organization. Sub-optimal model performance can be particularly challenging when the focus of interest is predicting onset of a specific disease (e.g. heart failure) or related events such as accelerated disease progression. For example, using Doctor AI (Choi et al., 2016a) , we discovered that RNN alone was ineffective to predict the onset of diseases such as cerebral degenerations (e.g. Leukodystrophy, Cerebral lipidoses) or developmental disorders (e.g. autistic disorder, Heller's syndrome), partly because their rare occurrence in the training data provided little learning opportunity to the flexible models like RNN.

The data requirement of deep learning models comes from having to assess exponential number of combinations of input features. This can be alleviated by exploiting medical ontologies that encodes hierarchical clinical constructs and relationships among medical concepts. Fortunately, there are many well-organized ontologies in healthcare such as the International Classification of Diseases (ICD), Clinical Classifications Software (CCS) (Stearns et al., 2001) or Systematized Nomenclature of Medicine-Clinical Terms (SNOMED-CT) (Project et al., 2010) . Nodes (i.e. medical concepts) close to one another in medical ontologies are likely to be associated with similar patients, allowing us to transfer knowledge among them. Therefore, proper use of medical ontologies will be helpful when we lack enough data for the nodes in the ontology to train deep learning models.

In this work, we propose GRAM, a method that infuses information from medical ontologies into deep learning models via neural attention. Considering the frequency of a medical concept in the EHR data and its ancestors in the ontology, GRAM decides the representation of the medical concept by adaptively combining its ancestors via attention mechanism. This will not only support deep learning models to learn robust representations without large amount of data, but also learn interpretable representations that align well with the knowledge from the ontology. The attention mechanism is trained in an end-to-end fashion with the neural network model that predicts the onset of disease(s). We also propose an effective initialization technique in addition to the ontological knowledge to better guide the representation learning process. We compared predictive performance (i.e. accuracy, data needs, interpretability) of GRAM to various models including the recurrent neural network (RNN) in two sequential diagnoses prediction tasks and one heart failure (HF) prediction task. We demonstrate that GRAM is up to 10% more accurate than the basic RNN for predicting diseases less observed in the training data. After discussing GRAM's scalability, we visualize the representations learned from various models where GRAM provides more intuitive representations by grouping similar medical concepts close to one another. Finally, we show GRAM's attention mechanism can be interpreted to understand how it assigns the right amount of attention to the ancestors of each medical concept by considering the data availability and the ontology structure.

2. METHODOLOGY:

We first define the notations describing EHR data and medical ontologies, followed by a description of GRAM (Section 2.2), the end-to-end training of the attention generation and predictive modeling (Section 2.3), and the efficient initialization scheme (Section 2.4).

2.1. BASIC NOTATION:

We denote the set of entire medical codes from the EHR as c 1 , c 2 , . . . , c |C| ∈ C with the vocabulary size |C|. The clinical record of each patient can be viewed as a sequence of visits V 1 , . . . , V T where each visit contains a subset of medical codes V t ⊆ C. V t can be represented as a binary vector x t ∈ {0, 1} |C| where the i-th element is 1 only if V t contains the code c i . To avoid clutter, all algorithms will be presented for a single patient.

We assume that a given medical ontology G typically expresses the hierarchy of various medical concepts in the form of a parent-child relationship, where the medical codes C form the leaf nodes. Ontology G is represented as a directed acyclic graph (DAG) whose nodes form a set

D = C + C . C = {c |C|+1 , c |C|+2 , . . . , c |C|+|C | }

defines the set of all non-leaf nodes (i.e. ancestors of the leaf nodes), where |C | represents the number of all non-leaf nodes. We use knowledge DAG to refer to G. A parent in the knowledge DAG G represents a related but more general concept over its children. Therefore, G provides a multi-resolution view of medical concepts with different degrees of specificity. While some ontologies are exclusively expressed as parent-child hierarchies (e.g. ICD-9, CCS), others are not. For example, in some instances SNOMED-CT also links medical concepts to causal or treatment relationships, but the majority relationships in SNOMED-CT are still parent-child. Therefore, we focus on the parent-child relationships in this work.

2.2. KNOWLEDGE DAG AND THE ATTENTION MECHANISM:

GRAM leverages the parent-child relationship of G to learn robust representations when data volume is constrained. GRAM balances the use of ontology information in relation to data volume in determining the level of specificity for a medical concept. When a medical concept is less observed in the data, more weight is given to its ancestors as they can be learned more accurately and offer general (coarse-grained) information about their children. The process of resorting to the parent concepts can be automated via the attention mechanism and the end-to-end training as described in Figure 1 .

In the knowledge DAG, each node c i is assigned a basic embedding vector e i ∈ R m , where m represents the dimensionality. Then e 1 , . . . , e |C| are the basic embeddings of the codes c 1 , . . . , c |C| while e |C|+1 , . . . , e |C|+|C | represent the basic embeddings of the internal nodes c |C|+1 , . . . , c |C|+|C | . The initialization of these basic embeddings is described in Section 2.4. We formulate a leaf node's Figure 1 : The illustration of GRAM. Leaf nodes (solid circles) represents a medical concept in the EHR, while the non-leaf nodes (dotted circles) represent more general concepts. The final representation g i of the leaf concept c i is computed by combining the basic embeddings e i of c i and e g , e c and e a of its ancestors c g , c c and c a via an attention mechanism. The final representations form the embedding matrix G for all leaf concepts. After that, we use G to embed patient visit vector x t to a visit representation v t , which is then fed to a neural network model to make the final prediction ŷt . final representation as a convex combination of the basic embeddings of itself and its ancestors:

EQUATION

where g i ∈ R m denotes the final representation of the code c i , A(i) the indices of the code c i and c i 's ancestors, e j the basic embedding of the code c j and α ij ∈ R the attention weight on the embedding e j when calculating g i . The attention weight α ij in Eq. ( 1) is calculated by the following Softmax function,

EQUATION

f (e i , e j ) is a scalar value representing the compatibility between the basic embeddings of e i and e k . We compute f (e i , e j ) via the following feed-forward network with a single hidden layer (MLP),

f (e i , e j ) = u a tanh(W a e i e j + b a )

where W a ∈ R l×2m is the weight matrix for the concatenation of e i and e j , b ∈ R l the bias vector, and u a ∈ R l the weight vector for generating the scalar value. The constant l represents the dimension size of the hidden layer of f (•, •). Note that we always concatenate e i and e j in the child-ancestor order.

Remarks: The example in Figure 1 is derived based on a single path from c i to c a . However, the same mechanism can be applicable to multiple paths as well. For example, code c k has two paths to the root c a , containing five ancestors in total. Another scenario is where the EHR data contain both leaf codes and some ancestor codes. We can move those ancestors present in EHR data from the set C to C and apply the same process as Eq. ( 1) to obtain the final representations for them.

2.3. END-TO-END TRAINING WITH A PREDICTIVE MODEL:

We train the attention mechanism together with a predictive model such that the attention mechanism improves the predictive performance. Once the final representations g 1 , g 2 , . . . , g |C| of all medical codes are obtained, we can convert visit V t to a visit representation v t by using the embedding matrix G ∈ R m×|C| where g i is its i-th column as in Figure 1 . We continue the mathematical formulation under the assumption that we are using the RNN to perform sequential diagnoses prediction (Choi et al., 2016a; b) with the objective of predicting the disease codes of the next visit V t+1 given the visit records up to the current timestep V 1 , V 2 , . . . , V t , which can be expressed as follows, where x t ∈ R |C| denotes the t-th visit; v t ∈ R m the t-th visit representation; h t ∈ R r the RNN's hidden layer at t-th time step (i.e. t-th visit); W ∈ R |C|×r and b ∈ R |C| the weight matrices and the bias vector of the Softmax function; r denotes the dimension size of the hidden layer. We use "RNN" to denote any recurrent neural network variants that can cope with the vanishing gradient problem (Bengio et al., 1994) , such as LSTM (Hochreiter & Schmidhuber, 1997) , GRU (Cho et al., 2014), and IRNN (Le et al., 2015) , with any varying numbers of hidden layers. The prediction loss for all time steps is calculated using the cross entropy as follows,

y t = x t+1 = Softmax(Wh t + b), where h 1 , h 2 , . . . , h t = RNN(v 1 , v 2 , . . . , v t ), where (4) v 1 , v 2 , . . . , v t = tanh(G[x 1 , x 2 , . . . , x t ])

L(x 1 , x 2 . . . , x T ) = -1 T -1 T -1 t=1 y t log( y t ) + (1 -y t ) log(1 -y t )

where we sum the cross entropy errors from all dimensions of y t , T denotes the length of the visit sequence. Note that the above loss is defined for a single patient. But we can take the average of the individual loss for multiple patients.

2.4. INITIALIZING BASIC EMBEDDINGS:

The attention generation mechanism in Section 2.2 requires basic embeddings e i of each node in the knowledge DAG. The basic embeddings of ancestors, however, pose a difficulty because they are often not observed in the data.To better initialize them, we use co-occurrence information to learn the basic embeddings of medical codes and their ancestors. Co-occurrence has proven to be an important source of information when learning representations of words or medical concepts (Mikolov et al., 2013; Choi et al., 2016c; e) . To train the basic embeddings, we employ GloVe (Pennington et al., 2014) , which uses the global co-occurrence matrix of words to learn their representations. In our case, the co-occurrence matrix of the codes and the ancestors was generated by counting the co-occurrences within each visit V t , where we augment each visit with the ancestors of the codes in the visit. Details of training the basic embeddings are described in the Appendix A. Note that, with or without the initialization, the basic embeddings e i 's of both leaf nodes (i.e. medical codes) and non-leaf nodes (i.e. ancestors) are fine-tuned when training our model, since the error signal flows from the output y t to the final representations g i 's which are convex combinations of e i 's.

3. EXPERIMENTS:

We conduct three experiments to determine if GRAM offered superior prediction performance when facing data insufficiency. We first describe the experimental setup followed by results comparing predictive performance of GRAM with various baseline models. After discussing GRAM's scalability, we qualitatively evaluate the interpretability of the resulting representation. The source code of GRAM is publicly available at https://github.com/mp2893/gram.

3.1. EXPERIMENT SETUP:

Prediction tasks and source of data: We conduct two sequential diagnoses prediction tasks, which aim at predicting all diagnosis categories in the next visit, and one heart failure (HF) prediction task, which is a binary prediction task for predicting a future HF onset where the prediction is made only once at the last visit x T . Two sequential diagnoses predictions are respectively conducted using 1) Sutter Palo Alto Medical Foundation (PAMF) dataset, which consists of 18-years longitudinal medical records of 258K patients between age 50 and 90. This will determine GRAM's performance for general adult population with long visit records. 2) MIMIC-III dataset (Johnson et al., 2016; Goldberger et al., 2000) , which is a publicly available dataset consisting of medical records of 7.5K intensive care unit (ICU) patients over 11 years. This will determine GRAM's performance for high-risk patients with very short visit records. We utilize all the patients with at least 2 visits. We prepared the true labels y t by grouping the ICD9 codes into 283 groups using CCS single-level diagnosis grouper 1 . This is to improve the training speed and predictive performance for easier analysis, while preserving sufficient granularity for each diagnosis. Each diagnosis code's varying frequency in the training data can be viewed as different degrees of data insufficiency. We calculate Accuracy@k for each of CCS single-level diagnosis codes such that, given a visit V t , we get 1 if the target diagnosis is in the top k guesses and 0 otherwise. We conduct HF prediction on Sutter heart failure (HF) cohort, which is a subset of Sutter PAMF data for a heart failure onset prediction study with 3.4K HF cases and 27K controls chosen by a set of criteria (see Appendix B). This will determine GRAM's performance for a different prediction task where we predict the onset of one specific condition. We randomly downsample the training data to create different degrees of data insufficiency. We use area under the ROC curve (AUC) to measure the performance. A summary of the datasets are provided in Table 1 .We used CCS multi-level diagnoses hierarchy2 as our knowledge DAG G. We also tested the ICD9 code hierarchy3 , but the performance was similar to using CCS multi-level hierarchy. For all three tasks, we randomly divide the dataset into the training, validation and test set by .75:.10:.15 ratio, and use the validation set to tune the hyper-parameters. Further details regarding the hyper-parameter tuning are provided in Appendix C. The test set performance is reported in the paper.

Implementation details: We implemented GRAM with Theano 0.8.2 (Team, 2016). For training models, we used Adadelta (Zeiler, 2012) with a mini-batch of 100 patients, on a machine equipped with Intel Xeon E5-2640, 256GB RAM, four Nvidia Titan X's and CUDA 7.5.

Models for comparison are the following. The first two GRAM+ and GRAM are the proposed methods and the rest are baselines. Hyper-parameter tuning is configured so that the number of parameters for the baselines would be comparable to GRAM's. Further details are provided in Appendix C.

• GRAM: Input sequence x 1 , . . . , x T is first transformed by the embedding matrix G, then fed to the GRU with a single hidden layer, which in turn makes the prediction, as described by Eq. ( 4).

The basic embeddings e i 's are randomly initialized.

• GRAM+: We use the same setup as GRAM, but the basic embeddings e i 's are initialized according to Section 2.4. • RandomDAG: We use the same setup as GRAM, but each leaf concept has five randomly assigned ancestors from the CCS multi-level hierarchy to test the effect of correct domain knowledge. • RNN: Input x t is transformed by an embedding matrix W emb ∈ R k×|C| , then fed to the GRU with a single hidden layer. The embedding size k is a hyper-parameter. W emb is randomly initialized and trained together with the GRU. • RNN+: We use the same setup as RNN, but we initialize the embedding matrix W emb with GloVe vectors trained only with the co-occurrence of leaf concepts. This is to compare GRAM with a similar weight initialization technique. • SimpleRollUp: We use the same setup as RNN. But for input x t , we replace all diagnosis codes with their direct parent codes in the CCS multi-level hierarchy, giving us 578, 526 and 517 input codes respectively for Sutter data, MIMIC-III and Sutter HF cohort. This is to compare the performance of GRAM with a common grouping technique. • RollUpRare: We use the same setup as RNN, but we replace any diagnosis code whose frequency is less than a certain threshold in the dataset with its direct parent. We set the threshold to 100 for Sutter data and Sutter HF cohort, and 10 for MIMIC-III, giving us 4,408, 935 and 1,538 input codes respectively for Sutter data, MIMIC-III and Sutter HF cohort. This is an intuitive way of dealing with infrequent medical codes.

3.2. PREDICTION PERFORMANCE AND SCALABILITY:

Tables 2a and 2b show the sequential diagnoses prediction performance on Sutter data and MIMIC-III. Both figures show that GRAM+ outperforms other models when predicting labels with significant data insufficiency (i.e. less observed in the training data).The performance gain is greater for MIMIC-III, where GRAM+ outperforms the basic RNN by 10% in the 20th-40th percentile range. This seems to come from the fact that MIMIC patients on average have significantly shorter visit history than Sutter patients, with much more codes received per visit. Such short sequences make it difficult for the RNN to learn and predict diagnoses sequence. The performance difference between GRAM+ and GRAM suggests that our proposed initialization scheme of the basic embeddings e i is important for sequential diagnosis prediction.

Table 2c shows the HF prediction performance on Sutter HF cohort. GRAM and GRAM+ consistently outperforms other baselines (except RNN+) by 3∼4% AUC, and RNN+ by maximum 1.8% AUC. These differences are quite significant given that the AUC is already in the mid-80s, a high value for HF prediction, cf. (Choi et al., 2016d) . Note that, for GRAM+ and RNN+, we used the downsampled training data to initialize the basic embeddings e i 's and the embedding matrix W emb with GloVe, respectively. The result shows that the initialization scheme of the basic embeddings in GRAM+ gives limited improvement over GRAM. This stems from the different natures of the two prediction tasks.

While the goal of HF prediction is to predict a binary label for the entire visit sequence, the goal of sequential diagnosis prediction is to predict the co-occurring diagnosis codes at every visit. Therefore the co-occurrence information infused by the initialized embedding scheme is more beneficial to sequential diagnosis prediction. Additionally, this benefit is associated with the natures of the two prediction tasks than the datasets used for the prediction tasks. Because the initialized embedding shows different degrees of improvement as shown by Tables 2a and 2c , when Sutter HF cohort is a subset of Sutter PAMF, thus having similar characteristics. Additional prediction results when varying the k of Accuracy@k are discussed in the Appendix D.

Overall, GRAM showed superior predictive performance under data insufficiency in three different experiments, demonstrating its general applicability in predictive healthcare modeling. Now we briefly discuss the scalability of GRAM by comparing its training time to RNN's. Table 3 shows the number of seconds taken for the two models to train for a single epoch for each predictive modeling task. GRAM+ and RNN+ showed the same behavior as GRAM and RNN. GRAM takes approximately 50% more time to train for a single epoch for all prediction tasks. This stems from calculating attention weights and the final representations g i for all medical codes. GRAM also generally takes about 50% more epochs to reach to the model with the lowest validation loss. This is due to optimizing an extra MLP model that generates the attention weights. Overall, use of GRAM adds a manageable amount of overhead in training time to the plain RNN.

3.3. QUALITATIVE EVALUATION OF INTERPRETABLE REPRESENTATIONS:

To qualitatively assess the interpretability of the learned representations of the medical codes, we plot on a 2-D space using t-SNE (Maaten & Hinton, 2008) the final representations g i of 2,000 (Figure 2a ). The colors represent the highest disease categories and the text annotations represent the detailed disease categories in CCS multi-level hierarchy. For comparison, we also show the t-SNE plots on the strongest results from RNN+ (Figure 2b ), and GloVe (Figure 2c ), the same embedding technique in initializing the basic embeddings e i . Figures 2b and 2c confirm that interpretable representations cannot simply be learned only by co-occurrence or supervised prediction without medical knowledge. GRAM+ learns disease representations that are significantly more consistent with the given knowledge DAG G. Therefore the neural network predictive model that accepts g i is using accurate representations that lead to higher predictive performance. Additional scatterplots of other models are provided in Appendix E for comparison. An interactive visualization tool can be accessed at http://www.sunlab.org/research/gram-graph-based-attention-model/.

3.4. ANALYSIS OF THE ATTENTION BEHAVIOR:

Next we show that GRAM's attention can be interpreted to understand how it considers data availability and knowledge DAG's structure when performing a prediction task. Using Eq. ( 1), we can calculate the attention weights of individual disease. Figure 3 shows the attention behaviors of four representative diseases when performing HF prediction on Sutter HF cohort.

Other pneumothorax (ICD9 512.89) in Figure 3a is rarely observed in the data and has only five siblings. In this case, most information is derived from the highest ancestor. Temporomandibular joint disorders & articular disc disorder (ICD9 524.63) in Figure 3b is rarely observed but has 139 siblings. In this case, its parent receives a stronger attention because it aggregates sufficient samples from all of its children to learn a more accurate representation. Note that the disease itself also receives a stronger attention to facilitate easier distinction from its large number of siblings. In each figure, the leaf node represents the disease and upper nodes are its ancestors. The size of the node shows the amount of attention it receives, which is also shown by the bar charts. The number in the parenthesis next to the disease is its frequency in the training data. We exclude the root of the knowledge DAG G from all figures as it did not play a significant role.

Unspecified essential hypertension (ICD9 401.9) in Figure 3c is very frequently observed but has only two siblings. In this case, GRAM assigns a very strong attention to the leaf, which is logical because the more you observe a disease, the stronger your confidence becomes. Need for prophylactic vaccination and inoculation against influenza (ICD9 V04.81) in Figure 3d is quite frequently observed and also has 103 siblings. The attention behavior in this case is quite similar to the case with fewer siblings (Figure 3b ) with a slight attention shift towards the leaf concept as more observations lead to higher confidence.

4. RELATED WORK:

We introduce recent studies related to GRAM that learn the representations of graphs and discuss their relationship with GRAM. Several studies focused on learning the representations of graph vertices by using the neighbor information. DeepWalk (Perozzi et al., 2014 ) and node2vec (Grover & Leskovec, 2016) use random walk while LINE (Tang et al., 2015) uses breadth-first search to find the neighbors of a vertex and learn its representation based on the neighbor information. Graph convolutional approaches (Yang et al., 2016; Kipf & Welling, 2016 ) also focus on learning the vertex representations to mainly perform vertex classification. These works focus on solving the graph data problems whereas GRAM focuses on solving EHR data problems using the knowledge DAG as supplementary information.

Several researchers tried to model the knowledge DAG such as WordNet (Miller, 1995) or Freebase (Bollacker et al., 2008) where two entities are connected with various types of relation, forming a set of triples. They aim to project entities and relations (Bordes et al., 2013; Socher et al., 2013; Wang et al., 2014; Lin et al., 2015) to the latent space based on the triples or additional information such as hierarchy of entities (Xie et al., 2016) . These works demonstrated tasks such as link prediction, triple classification or entity classification using the learned representations. More recently, Li et al. (2016) learned the representations of words and Wikipedia categories by utilizing the hierarchy of Wikipedia categories. GRAM is fundamentally different from the above studies in that it aims to design intuitive attention mechanism on the knowledge DAG as a knowledge prior to cope with data insufficiency and learn medically interpretable representations to make accurate predictions.

A classical approach for incorporating side information in the predictive models is to use graph Laplacian regularization (Weinberger et al., 2006; Che et al., 2015) . However, using this approach is not straightforward as it relies on the appropriate definition of distance on graphs which is often unavailable.

5. CONCLUSION:

Data insufficiency, either due to less common diseases or small datasets, is one of the key hurdles in healthcare analytics, especially when we apply deep neural networks models. To overcome this hurdle, we leverage the knowledge DAG, which provides a multi-resolution view of medical concepts.

We propose GRAM, a graph-based attention model using both a knowledge DAG and EHR to learn an accurate and interpretable representations for medical concepts. GRAM chooses a weighted average of ancestors of a medical concept and train the entire process with a predictive model in an end-to-end fashion. We conducted three predictive modeling experiments on real EHR datasets and showed significant improvement in the prediction performance, especially on low-frequency diseases and small datasets. Analysis of the attention behavior provided intuitive insight of GRAM.

Figure 4 : Creating the co-occurrence matrix together with the ancestors. Here we exclude the root node, which will be just a single row (column). We first create an augmented dataset by adding the ancestors of the code to the dataset. Then, we count the co-occurrence of the codes. Performing GloVe on this matrix produces the embedding vectors e i .

Leaf medical codes First ancestor Second ancestor Third ancestor Fourth ancestor

A GENERATING GLOVE EMBEDDINGS:

We learn the basic embeddings e i 's of medical codes and their ancestors using GloVe (Pennington et al., 2014) , which uses global co-occurrence matrix of words to learn their representations. We generate the co-occurrence matrix of the codes and the ancestors by counting the co-occurrence within each visit V t . However, since visits only contain the leaf codes c ∈ C, we augment each visit with the ancestors of the codes in each visit, then count the co-occurrence of codes and ancestors altogether.

We describe the details of the algorithm with an example. We borrow the parent-child relationships from the knowledge DAG of Figure 1 . Given a visit V t ,

EQUATION

we augment it with the ancestors of all the codes to obtain the augmented visit V t ,

V t = {c d , c b , c a , c i , c g , c c , c a , c k , c j , c f , c c , c b , c a } (6)

where the added ancestors are underlined. Note that a single ancestor can appear multiple times in V t . In fact, the higher the ancestor is in the knowledge DAG, the more times it is likely to appear in V t . We count the co-occurrence of two codes in V t as follows,

EQUATION

where count(c i , V t ) is the number of times the code c i appears in the augmented visit V t . For example, the co-occurrence between the leaf code c i and the root c a is 3. However, the co-occurrence between the ancestor c c and the root c a is 6. Therefore our algorithm will naturally make the ancestor codes have higher co-occurrence with other codes compared to leaf medical codes. We repeat this calculation for all pairs of codes in all augmented visits of all patients to obtain the co-occurrence matrix depicted by Figure 4 . For training the embedding vectors using the co-occurrence matrix, we use the same procedure and hyper-parameter as described in Pennington et al. (2014) .

B HEART FAILURE COHORT CONSTRUCTION:

For the heart failure (HF) case patients, we select patients between 40 to 85 years of age at the time of HF diagnosis. HF diagnosis (HFDx) criteria are defined as: 1) Qualifying ICD-9 codes for HF appeared in the encounter records or medication orders. Qualifying ICD-9 codes are listed in Table 4 .

2) at least three clinical encounters with qualifying ICD-9 codes had to occur within 12 months of each other, where the date of HFDx was assigned to the earliest of the three dates. If the time span between the first and second appearances of the HF diagnosis code was greater than 12 months, the date of the second encounter was used as the first qualifying encounter. Up to ten eligible controls (in terms of sex, age, location) were selected for each case, yielding average 9 controls per case. Each control was also assigned an index date, which is the HFDx date of the matched case. Controls are selected such that they did not meet the HF diagnosis criteria prior to the HFDx date plus 182 days of Table 4 : Qualifying ICD-9 codes for heart failure their corresponding case. Control subjects were required to have their first office encounter within one year of the matching HF case patient's first office visit, and have at least one office encounter 30 days before or any time after the case's HFDx date to ensure similar duration of observations among cases and controls.

C HYPER-PARAMETER TUNING:

We define five hyper-parameters for GRAM:

• dimensionality m of the basic embedding e i : [100, 200, 300, 400, 500] • dimensionality r of the RNN hidden layer h t from Eq. ( 4): [100, 200, 300, 400, 500] • dimensionality l of W a and b a from Eq. ( 3): [100, 200, 300, 400, 500] • L 2 regularization coefficient for all weights except RNN weights: [0.1, 0.01, 0.001, 0.0001]

• dropout rate for the dropout on the RNN hidden layer: [0.0, 0.2, 0.4, 0.6, 0.8]

We performed 100 iterations of the random search by using the above ranges for each of the three prediction experiments. For sequential diagnoses prediction on Sutter data, we used 10% of the training data to tune the hyper-parameters to balance the time and search space. To match the baselines' number of parameters to GRAM's, we add 550 to the list of m's possible values. This will make the baseline's largest possible number of parameters comparable to the GRAM's largest possible number of parameters.

For SimpleRollUp and RollUpRare, the number of input codes is smaller than other models due to the grouping. Therefore, to match their largest possible number of parameters to GRAM's, we need to add much larger values to m. However, after preliminary experiments, as expected, setting m to too large a value degraded the performance due to overfitting. Since the number of input codes decreased due to the grouping, increasing the dimensionality of e i is not a logical thing to do. Therefore, for SimpleRollUp and RollUpRare, we use the same list of values for m as other baselines.

Table 5 describes the final hyper-parameter settings we used for all models for each prediction experiments.

D PREDICTION RESULTS USING DIFFERENT k'S IN ACCURACY@K:

We show Accuracy@k using k = 5, 10, 20, 30 for sequential diagnoses prediction on Sutter data (Tables 6a, 6b, 6c and 6d) and MIMIC-III (Tables 7a, 7b, 7c and 7d). We can see from the tables that GRAM+ consistently outperforms other models under 40th percentile range, except when k = 20, 30 for sequential diagnoses prediction on Sutter data where SimpleRollUp shows similar performance. We can also see that GRAM+ performs significantly better than other models for all k = 5, 10, 20, 30 when predicting infrequently observed diseases on MIMIC-III. As discussed in Section 3.2, this seems to come from the short visit sequences of MIMIC patients.

E T-SNE 2-D PLOTS OF VARIOUS MODELS:

For further comparison, we display t-SNE scatterplots of GRAM (Figure 5a ) RandomDAG (Figure 5b , RNN (Figure 5c ), and Skip-gram (Figure 5d ). GRAM, RandomDAG and RNN were trained for sequential diagnoses prediction on Sutter data, and Skip-gram (Mikolov et al., 2013) was trained on Sutter data as it is an unsupervised method. For Skip-gram, we used each visit V t as the context window. As we do not distinguish between the target concept and the neighbor concepts, we calculated the Skip-gram objective function using all possible pairs of codes within a single visit.

We can see from Figure 5a that the quality of the final representations g i of GRAM is quite similar to GRAM+ (Figure 2a ). Compared to other baselines, GRAM demonstrates significantly more structured representations that align well with the given knowledge DAG. It is interesting that Skip-gram shows the most structured representation among all baselines. We used GloVe to initialize the basic embeddings e i in this work because it uses global co-occurrence information and its training time is dependent only on the total number of unique concepts |C|. Skip-gram's training time, on the other hand, depends on both the number of patients and the number of visits each patient made, which makes the algorithm generally slower than GloVe. However, considering both Figures 2c and 5d , initializing e i 's with Skip-gram vectors might give us additional performance boost. 

:

https://www.hcup-us.ahrq.gov/toolssoftware/ccs/AppendixASingleDX.txt

https://www.hcup-us.ahrq.gov/toolssoftware/ccs/AppendixCMultiDX.txt

http://www.icd9data.com/2015/Volume1/default.htm

The scatterplots of models trained for sequential diagnoses prediction on MIMIC-III and HF prediction for Sutter HF cohort were similar but less structured due to smaller data size.

Figures and Tables:

Scatterplot of the final representations gi's of GRAM+ (b) Scatterplot of the trained embedding matrix W emb of RNN+ (c) Scatterplot of the disease representations trained by GloVe

Figure 2: t-SNE scatterplots of medical concepts trained by GRAM+, RNN+ and GloVe

Figure 3: GRAM's attention behavior during HF prediction for four representative diseases (each column). In each figure, the leaf node represents the disease and upper nodes are its ancestors. The size of the node shows the amount of attention it receives, which is also shown by the bar charts. The number in the parenthesis next to the disease is its frequency in the training data. We exclude the root of the knowledge DAG G from all figures as it did not play a significant role.

Figure 5: Scatterplot of medical concepts trained by various models. We used t-SNE to reduce the dimension to 2-D.

Basic statistics of Sutter PAMF, MIMIC-III and Sutter heart failure (HF) cohort. Note that for all datasets, we selected patients who made at least two hospital visits.

Table content: <table><tr><td>Dataset</td><td colspan="2">Sutter PAMF MIMIC-III</td><td>Sutter HF cohort</td></tr><tr><td># of patients</td><td>258,555  †</td><td>7,499  †</td><td>30,727  † (3,408 cases)</td></tr><tr><td># of visits</td><td>13,920,759</td><td>19,911</td><td>572,551</td></tr><tr><td>Avg. # of visits per patient</td><td>53.8</td><td>2.66</td><td>38.38</td></tr><tr><td># of unique ICD9 codes</td><td>10,437</td><td>4,893</td><td>5,689</td></tr><tr><td>Avg. # of codes per visit</td><td>1.98</td><td>13.1</td><td>2.06</td></tr><tr><td>Max # of codes per visit</td><td>54</td><td>39</td><td>29</td></tr><tr><td>†</td><td/><td/><td/></tr></table>

Performance of three prediction tasks. The x-axis of (a) and (b) represents the labels grouped by the percentile of their frequencies in the training data in non-decreasing order. For (c), we vary the size of the training data to train the models. (b) uses Accuracy@20 because MIMIC-III has a large average number of codes per visit (see Table1).

Table content: <table><tr><td>Model</td><td>0-20</td><td colspan="4">20-40 40-60 60-80 80-100</td><td colspan="2">Model</td><td>0-20</td><td colspan="3">20-40 40-60 60-80 80-100</td></tr><tr><td>GRAM+</td><td colspan="5">0.0150 0.3242 0.4325 0.4238 0.4903</td><td colspan="2">GRAM+</td><td colspan="4">0.0672 0.1787 0.2644 0.2490 0.6267</td></tr><tr><td>GRAM</td><td colspan="5">0.0042 0.2987 0.4224 0.4193 0.4895</td><td colspan="2">GRAM</td><td colspan="4">0.0556 0.1016 0.1935 0.2296 0.6363</td></tr><tr><td colspan="6">RandomDAG 0.0050 0.2700 0.4010 0.4059 0.4853</td><td colspan="6">RandomDAG 0.0329 0.0708 0.1346 0.1512 0.4494</td></tr><tr><td>RNN+</td><td colspan="5">0.0069 0.2742 0.4140 0.4212 0.4959</td><td colspan="2">RNN+</td><td colspan="4">0.0454 0.0843 0.2080 0.2494 0.6239</td></tr><tr><td>RNN</td><td colspan="5">0.0080 0.2691 0.4134 0.4227 0.4951</td><td colspan="2">RNN</td><td colspan="4">0.0454 0.0731 0.1804 0.2371 0.6243</td></tr><tr><td colspan="6">SimpleRollUp 0.0085 0.3078 0.4369 0.4330 0.4924</td><td colspan="6">SimpleRollUp 0.0578 0.1328 0.2455 0.2667 0.6387</td></tr><tr><td>RollUpRare</td><td colspan="5">0.0062 0.2768 0.4176 0.4226 0.4956</td><td colspan="2">RollUpRare</td><td colspan="4">0.0454 0.0653 0.1843 0.2364 0.6277</td></tr><tr><td colspan="5">(a) Accuracy@5 of sequential diagnoses</td><td/><td/><td colspan="5">(b) Accuracy@20 of sequential diagnoses</td></tr><tr><td colspan="3">prediction on Sutter data</td><td/><td/><td/><td/><td colspan="4">prediction on MIMIC-III</td></tr><tr><td/><td>Model</td><td>10%</td><td>20%</td><td>30%</td><td>40%</td><td>50%</td><td>60%</td><td>70%</td><td>80%</td><td>90%</td><td>100%</td></tr><tr><td/><td>GRAM+</td><td colspan="10">0.7970 0.8223 0.8307 0.8332 0.8389 0.8404 0.8452 0.8456 0.8447 0.8448</td></tr><tr><td/><td>GRAM</td><td colspan="10">0.7981 0.8217 0.8340 0.8332 0.8372 0.8377 0.8440 0.8431 0.8430 0.8447</td></tr><tr><td/><td colspan="11">RandomDAG 0.7644 0.7882 0.7986 0.8070 0.8143 0.8185 0.8274 0.8312 0.8254 0.8226</td></tr><tr><td/><td>RNN+</td><td colspan="10">0.7930 0.8117 0.8162 0.8215 0.8261 0.8333 0.8343 0.8353 0.8345 0.8335</td></tr><tr><td/><td>RNN</td><td colspan="10">0.7811 0.7942 0.8066 0.8111 0.8156 0.8207 0.8258 0.8278 0.8297 0.8314</td></tr><tr><td/><td colspan="11">SimpleRollUp 0.7799 0.8022 0.8108 0.8133 0.8177 0.8207 0.8223 0.8272 0.8269 0.8258</td></tr><tr><td/><td>RollUpRare</td><td colspan="10">0.7830 0.8067 0.8064 0.8119 0.8211 0.8202 0.8262 0.8296 0.8307 0.8291</td></tr><tr><td/><td/><td colspan="8">(c) AUC of HF onset prediction on Sutter HF cohort</td><td/></tr></table>

Scalablity result in per epoch training time in second (the number of epochs needed).

Table content: <table><tr><td>Model</td><td>Sequential diagnosis prediction (Sutter data)</td><td>Sequential diagnosis prediction (MIMIC-III)</td><td>HF prediction (Sutter HF cohort)</td></tr><tr><td>GRAM</td><td>525s (39 epochs)</td><td>2s (11 epochs)</td><td>12s (7 epochs)</td></tr><tr><td>RNN</td><td>352s (24 epochs)</td><td>1s (6 epochs)</td><td>8s (5 epochs)</td></tr></table>

Hyper-parameters used by the models in each predictive modeling experiments

Table content: <table><tr><td>Experiment</td><td>Model</td><td>m</td><td>r</td><td>l</td><td>L2</td><td>Dropout rate</td></tr><tr><td/><td>GRAM+</td><td colspan="4">500 500 100 0.0001</td><td>0.6</td></tr><tr><td/><td>GRAM</td><td colspan="4">500 500 100 0.0001</td><td>0.6</td></tr><tr><td>Disease progression modeling (Sutter data)</td><td>RandomDAG RNN+ RNN</td><td colspan="4">500 500 100 0.0001 550 500 0.0001 550 500 0.0001</td><td>0.6 0.6 0.6</td></tr><tr><td/><td colspan="3">SimpleRollUp 500 500</td><td/><td>0.0001</td><td>0.4</td></tr><tr><td/><td>RollUpRare</td><td colspan="2">500 500</td><td/><td>0.0001</td><td>0.2</td></tr><tr><td/><td>GRAM+</td><td colspan="4">400 400 100 0.0001</td><td>0.6</td></tr><tr><td/><td>GRAM</td><td colspan="3">400 400 100</td><td>0.001</td><td>0.6</td></tr><tr><td>Disease progression modeling (MIMIC-III)</td><td>RandomDAG RNN+ RNN</td><td colspan="3">400 400 100 550 400 550 400</td><td>0.001 0.001 0.001</td><td>0.6 0.8 0.8</td></tr><tr><td/><td colspan="3">SimpleRollUp 400 400</td><td/><td>0.001</td><td>0.6</td></tr><tr><td/><td>RollUpRare</td><td colspan="2">400 400</td><td/><td>0.0001</td><td>0.0</td></tr><tr><td/><td>GRAM+</td><td colspan="3">200 100 100</td><td>0.001</td><td>0.6</td></tr><tr><td/><td>GRAM</td><td colspan="3">200 100 100</td><td>0.001</td><td>0.6</td></tr><tr><td>HF prediction (Sutter HF cohort)</td><td>RandomDAG RNN+ RNN</td><td colspan="3">300 100 200 200 100 200 100</td><td>0.001 0.0001 0.001</td><td>0.6 0.6 0.6</td></tr><tr><td/><td colspan="3">SimpleRollUp 300 200</td><td/><td>0.001</td><td>0.4</td></tr><tr><td/><td>RollUpRare</td><td colspan="2">100 100</td><td/><td>0.001</td><td>0.6</td></tr></table>

Peer review (iclr2017):
I read the authors* response and maintain my rating. --- This paper introduces an approach for integrating a direct acyclic graph structure of the data into word / code embeddings, in order to leverage domain knowledge and thus help train an RNN with scarce data. It is applied to codes of medical visits. Each code is part of an ontology, which can be represented by a DAG, where codes correspond to leaf nodes, and where different codes may share common ancestors (non-leaf nodes) in the DAG. Instead of embedding merely the leaf nodes, one can also embed the non-leaf nodes, and the embeddings of the code and its ancestors can be combined using a convex sum. That convex sum can be seen as an attention mechanism over the representation. The attention weights depend on the embeddings and the weights of an MLP, meaning that the model can separate learning the code embeddings and the interaction between the codes. Embedding codes are pretrained using GloVe, then fine-tuned. The model is properly evaluated on two medical datasets, with several variations to isolate the contribution of the DAG (GRAM or GRAM+ vs. RNN or RandomDAG) and of pretraining the embeddings (RNN+ vs RNN, GRAM+ vs GRAM). Both are shown to help achieve the best performance and the evaluation methodology seems thorough. The paper is also well written, and the case for MLP attention instead of a plain dot product of embeddings was made by the authors. My only two comments would be: 1) Why is there a softmax in equation 4, given that the loss is multivariate cross-entropy (in the predicted visit, several codes could be equal to 1), not a a single-class cross-entropy? 2) What is the embedding dimension m?