Year: 2021

Abstract:
In many domains, including healthcare, biology, and climate science, time series are irregularly sampled with variable time between successive observations and different subsets of variables (sensors) are observed at different time points, even after alignment to start events. These data create multiple challenges for prevailing models that assume fully observed and fixed-length feature representations. To address these challenges, it is essential to understand the relationships between sensors and how they evolve over time. Here, we introduce RAINDROP, a graphguided network for learning representations of irregularly sampled multivariate time series. RAINDROP represents every sample as a graph, where nodes indicate sensors and edges represent dependencies between them. RAINDROP models dependencies between sensors using neural message passing and temporal selfattention. It considers both inter-sensor relationships shared across samples and those unique to each sample that can vary with time, and it adaptively estimates misaligned observations based on nearby observations. We use RAINDROP to classify time series and interpret temporal dynamics of three healthcare and human activity datasets. RAINDROP outperforms state-of-the-art methods by up to 11.4% (absolute points in F1 score), including methods that deal with irregular sampling using fixed discretization and set functions, and even in challenging leave-sensorout settings and setups that require generalizing to new patient groups.

heading: 1 INTRODUCTION
text: Multivariate time series are prevalent in a variety of domains including healthcare, space science, cybersecurity, biology, and finance (Ravuri et al., 2021; Sousa et al., 2020; Sezer et al., 2020; Fawaz et al., 2019; Abanda et al., 2019; Tang et al., 2018). Practical issues often exist in collecting sensor measurements that lead to various types of irregularities caused by missing observations, such as cost saving, sensor failures, external forces in physical scenarios, medical interventions, to name a few (Choi et al., 2020). While temporal machine learning models usually assume fully observable and fixed-size input data, irregularly sampled time series raise considerable challenges. For example, the observations of multiple sensors are not well-aligned; the time intervals among adjacent observations are different across sensors; and different samples have different numbers of observations for different subsets of sensors recorded at different time points.
Prior methods for dealing with irregularly sampled time series involve filling in missing values, using interpolation, kernel methods, and probabilistic approaches (Schafer & Graham, 2002). The missingness of observations carry informative power (Little & Rubin, 2014) and thus imputation of missing data is not always beneficial (Agniel et al., 2018). While modern techniques involve recurrent neural network architectures (e.g., RNN, LSTM, GRU) (Hochreiter & Schmidhuber, 1997; Cho et al., 2014) and transformers (Vaswani et al., 2017), they are restricted to regular sampling or assume aligned measurements across modalities. For misaligned measurements, existing methods tend to rely on a two-stage process of imputation to get a regularly-sampled version of a dataset, and then performing a downstream task such as classification. This decoupled approach might not fully exploit informative missingness patterns that could be essential for the downstream task thus achieving suboptimal performance (Wells et al., 2013; Li & Marlin, 2016). Thus, several approaches circumvent imputation and directly model irregularly sampled time series data (Che et al., 2018; Horn et al., 2020; Shukla & Marlin, 2021; 2018). However, few explicitly consider relational structure to
address the characteristics of irregularly sampled multivariate time series and none leverages graph neural networks.
Present work. To address the characteristics of irregularly sampled time series data, we propose to model temporal dynamics of sensor dependencies and how those relationships evolve over time. Previous studies (Wu et al., 2021; Li et al., 2020a; Zhang et al., 2019) have noted that the inter-sensor correlations bring rich information in modeling time series. Our intuitive assumption is that the observed sensors can indicate how the unobserved sensors currently behave, which can further improve the representation learning of irregular multivariate time series. We develop RAINDROP1, a graph-guided network that leverages relational structure to embed and classify irregularly sampled multivariate time series. RAINDROP takes samples as input, each sample containing multiple sensors and each sensor consisting of irregularly recorded observations (e.g., in clinical data, an individual patient’s state of health, recorded at irregular time intervals with different subsets of sensors observed at different times). RAINDROP model is inspired by the idea of raindrops falling into a pool at
sequential but non-uniform time intervals and thereby creating ripple effects that propagate across the pool. Mathematically, in RAINDROP, observations (i.e., raindrops) hit a sensor graph (i.e., pool) asynchronously and at irregular time intervals, and each observation is processed by passing messages to neighboring sensors (i.e., causes a ripple effect in the pool) taking into account the learned sensor dependencies (Figure 1). As such, RAINDROP can handle misaligned observations, varying time gaps, arbitrary numbers of observations, and produce embeddings via a novel hierarchical attention.
RAINDROP is the first to explicitly model sensor dependencies in learning representations of irregularly sampled time series. We represent dependencies with a separate sensor graph for every sample, wherein nodes indicate sensors and edges denote relationships between them. In addition to capturing distinct sensor dependencies within each sample, RAINDROP i) takes advantage of similarities between different samples via sharing of parameters in the calculation of attention weights, and ii) considers importance of successive observations via temporal attention. RAINDROP adaptively estimates missing observations based on recorded information and the learned graph structure. We compare RAINDROP to five state-of-the-art methods across three datasets and four experimental settings, including a setup where a subset of sensors in the test set have malfunctioned (i.e., have no readouts at all). Experiments show that RAINDROP outperforms baselines on all datasets with an average AUROC improvement of 3.5% in absolute points on classification tasks. Further, RAINDROP achieves a considerable margin (9.3% absolute points in accuracy on activity recognition) when a subset of sensors malfunction.

heading: 2 RELATED WORK
text: Learning with irregularly sampled multivariate time series. Irregular time series indicate that the time intervals between adjacent observations are varying (Zerveas et al., 2021; Tipirneni & Reddy, 2021; Chen et al., 2020). In the multivariate case, irregularity means that observations are misaligned across different sensors. Further, because of a multitude of sampling frequencies and varying time intervals, the number of observations can vary considerably across samples (Fang & Wang, 2020; Kidger et al., 2020). Predominant downstream tasks for time series are classification (i.e., predicting a label for a given sample, e.g., Tan et al. (2020); Ma et al. (2020)) and forecasting (i.e., anticipating future observations based on historical observations, e.g., Wu et al. (2020a)). The above mentioned data characteristics create considerable challenges for machine learning models that expect wellaligned and fixed-size input (Shukla & Marlin, 2020). An intuitive way to deal with irregular time series is to impute missing values and process them as a regular time series dataset (Mikalsen et al., 2021; Li & Marlin, 2020; Shan & Oliva, 2021). However, imputation methods can distort underlying
1Code and datasets are available at https://anonymous.4open.science/r/Raindrop.
distributions and introduce unwanted bias in the dataset. To this end, recent methods directly learn from irregularly sampled time series. For example, Che et al. (2018) propose a decay mechanism based on gated recurrent units (GRU-D), incorporating binary masking indicators and time intervals to capture long-term temporal dependencies. SeFT (Horn et al., 2020) re-formulates irregularly sampled time series into a set of observations and learn representations through set functions that are insensitive to alignment. Chen et al. (2018) model continuous-time hidden dynamics by latent ordinary differential equations (Latent-ODE), combined with an RNN for temporal representations. mTAND (Shukla & Marlin, 2021) leverages multi-time attention mechanism to learn temporal similarity from non-uniformly collected measurements and produces continuous-time embeddings. IP-Net (Shukla & Marlin, 2018) and DGM2 (Wu et al., 2021) adopt similar imputation paradigms to interpolate the irregular time series against a set of reference points using a kernel-based method. The learned inter-sensor relations are static without considering sample-specific and time-specific characteristics. In contrast with these methods, RAINDROP leverages dynamic graphs to address the characteristics of irregular time series and improve the quality of learned representations.
Learning with graphs and neural message passing. There has been a surge of interest in applying neural networks to graphs, leading to the development of graph embeddings (Zhou et al., 2020; Li et al., 2021), graph neural networks (Wu et al., 2020b), and message passing neural networks (Gilmer et al., 2017). Our approach is mainly related to methods that perform message passing along edges to update node representations through neural transformations (Riba et al., 2018; Nikolentzos et al., 2020; Galkin et al., 2020; Fey et al., 2020; Lin et al., 2018; Zhang et al., 2020). However, in contrast to message passing used to make predictions on graphs, we are interested in leveraging it to meet the challenges of irregularly sampled time series. In particular, we consider message passing where nodes are sensors, describing a particular sample (e.g., patient, Figure 1), and we design a message-passing network with learnable adjacency matrices. The key difference with the predominant use of message passing is that RAINDROP uses it to estimate edges (dependencies) between sensors rather than applying it on a fixed, apriori-given graph. To the best of our knowledge, no prior work utilized sensor dependencies for learning representations of irregularly sampled multivariate time series. Finally, while prior work used message passing for regular time series (Wang et al., 2020; Wu et al., 2020c; Kalinicheva et al., 2020), no modeling of irregularly sampled time series has been attempted.
3 RAINDROP
Let D = {(Si, yi) | i = 1, . . . , N} denote an irregular time series dataset with N labeled samples (Figure 2). Every sample Si is an irregular multivariate time series with a corresponding label yi ∈ {1, . . . , C}, indicating which of the C classes Si is associated with. Each sample contains M non-uniformly measured sensors that are denoted as u, v, etc. RAINDROP can also work on samples with a subset of sensors (see Sec. 4.1). Each sensor is given by a sequence of observations ordered by time. For sensor u in sample Si, we denote a single observation as a tuple (t, xti,u), meaning that sensor u was recorded with value xti,u ∈ R at timestamp t ∈ R+. We omit sample index i and sensor index u in timestamp t. Sensor observations are irregularly recorded, meaning that time intervals between successive observations can vary across sensors. For sensor u in sample Si, we use Ti,u to denote the set of timestamps that u, or at least one of u’s L-hop neighbors (L is the number of layers in RAINDROP’s message passing) is recorded. We use || and T to denote concatenation and transpose, respectively. We omit layer index l ∈ {1, . . . , L} for simplicity when text is clear. Problem (Representation learning for irregularly sampled multivariate time series). A datasetD of irregularly sampled multivariate time series is given, where each sample Si has multiple sensors and each sensor has a variable number of observations. RAINDROP learns a function f : Si → zi that maps Si to a fixed-length representation zi suitable for downstream task of interest, such as
classification. Using learned zi, RAINDROP can predict label ŷi ∈ {1, . . . , C} for Si.
RAINDROP learns informative embeddings for irregularly samples time series. The learned embeddings capture temporal patterns of irregular observations and explicitly consider varying dependencies between sensors. While we focus on time-series classification in this work, the proposed method can be easily extended to broader applications such as regression, clustering and generation tasks.

heading: 3.1 OVERVIEW OF RAINDROP
text: RAINDROP aims to learn a fixed-dimensional embedding zi for given Si, and then predict the associated label ŷi. RAINDROP learns sample embeddings in a hierarchical architecture composed of observation embedding, sensor embedding and sample embedding (Figure 2). We describe the procedure of RAINDROP in the context of at most one observation at a specific time t (one sensor has observation while all residual sensors do not have observation). If there are multiple observations at the same time, RAINDROP can effortlessly process them all in parallel.
We first construct a graph for each sample where nodes represent sensors and edges indicate relations between sensors (Sec. 3.2). Each sample Si corresponds to a graph Gi, where we use ei,uv to denote the weight of directed edge going from sensor u to sensor v. The graph structures are automatically learned while considering sample-wise and time-wise specificity.
The key idea of RAINDROP is to borrow information from u’s neighbors based on estimated relationships between u and other sensors. This is achieved via message passing carried out on Si’s dependency graph and initiated at node u in the graph. When an observation (t, xti,u) is recorded for sample Si at time t, RAINDROP first embeds the observation at active sensor u (i.e., sensor whose value was recorded) and then propagates messages (i.e., the observation embeddings) from u to neighboring sensors along edges in sensor dependency graph Gi. As a result, recording the value of u can affect u’s embedding as well as embeddings of other sensors that related to u (Sec. 3.3). Finally, RAINDROP generates sensor embeddings by aggregating all observation embeddings for each sensor (across all timestamps) using temporal attention weights (Sec. 3.4). At last, RAINDROP embeds every sample Si based on the embeddings of all sensors (Sec. 3.5) and feeds the learned sample embedding into downstream classifier.

heading: 3.2 CONSTRUCTING SENSOR DEPENDENCY GRAPHS
text: We build a directed weighted graph Gi = {V, Ei} for every sample Si and refer to it as the sensor dependency graph for Si. Nodes V represent sensors and edges Ei describe dependencies between sensors in sample Si that RAINDROP infers. As we show in experiments, RAINDROP can be directly used with samples that only contain a subset of sensors in V . We denote edge from u to v as a triplet (u, ei,uv, v), where ei,uv ∈ [0, 1] represents the strength of relationship between sensors u and v in sample Si. Edge (u, ei,uv, v) describes the relationship between u and v: when u receives an observation, it will send a neural message to v following edge ei,uv. If ei,uv = 0, there is no exchange of neural information between u and v, indicating that the two sensors are unrelated. We assume that the importance of u to v is different than the importance of v to u, and so we treat sensor dependency graphs as directed, i.e., ei,uv 6= ei,vu. All graphs are initialized as fully-connected graphs (i.e., ei,uv = 1 for any u, v and Si) and edge weights ei,uv are updated following Eq. 3 during model training. If available, it is easy to integrate additional domain knowledge into graph initialization.

heading: 3.3 GENERATING EMBEDDINGS OF INDIVIDUAL OBSERVATIONS
text: Let sensor u be activated at time t ∈ Ti,u and receive observation xti,u and let u be connected to v through edge (u, ei,uv, v). We next describe how to produce observation embeddings hti,u ∈ Rdh and hti,v ∈ Rdh for sensors u and v, respectively (Figure 3a). We omit layer index l and note that the proposed strategy applies to any number of layers.
Embedding observation of an active sensor. We treat u as an active sensor that has received observation xti,u. For sufficient expressive power (Veličković et al., 2018), we map observation x t i,u to a high-dimensional space using a nonlinear transformation: hti,u = σ(x t i,uRu). Here, we use sensorspecific transformations because values recorded at different sensors can follow different distributions, which is achieved by trainable weight vectors Ru depending on what sensor is activated (Li et al., 2020b). Alternative functions (such as a multilayer perceptron) can be considered for transforming
i,uv from previous layer and the learned inter-sensor
attention weights in all time steps. We explicitly show layer index l as multiple layers are involved.
xti,u into h t i,u. As h t i,u represents information brought on by measuring x t i,u, we regard h t i,u as u’s observation embedding at t. The sensor-specific weight vectorsRu are shared across all samples.
Message passing in sensor dependency graphs. For sensors that are not active at timestamp t but are neighbors of the active sensor u in the sensor dependency graph Gi, RAINDROP uses relationships between u and those sensors to estimate observation embeddings for them. We proceed by describing how RAINDROP generates observation embedding hti,v for sensor v assuming v is a neighbor of u in Gi. Given hti,u and edge (u, ei,uv, v), we first calculate inter-sensor attention weight αti,uv ∈ [0, 1], representing how important u is to v via the following equation:
αti,uv = σ(h t i,uD[rv||pti]T ), (1)
where rv ∈ Rdr is a trainable weight vector that is specific to the sensor receiving the message (i.e., hti,u). Vector rv allows the model to learn distinct attention weights for different edges going out from the same sensor u. Further, pti ∈ Rdt is the time representation obtained by converting a 1-dimensional timestamp t into a multi-dimensional vector pti by passing t through a series of trigonometric functions (Horn et al., 2020). See Appendix A.1 for details. RAINDROP uses pti to calculate attention weights that are sensitive to time. Finally,D is a trainable weight matrix mapping hti,u from dh dimensions to (dr+dt) dimensions. Taken this together, we can estimate the embedding hti,v for u’s neighbor v as follows:
hti,v = σ(h t i,uwuw T v α t i,uvei,uv), (2)
wherewu,wv ∈ Rdh are trainable weight vectors shared across all samples. Thewu is specific to active sensor u and wv is specific to neighboring sensor v. In the above equation, ei,uv denotes edge weight shared across all timestamps. The above message passing describes the processing of a single observation at a single timestamp. In case multiple sensors are active at time t and connected with v, we normalize αti,uv (with softmax function) across active sensors and aggregate messages at v.
Overall, RAINDROP produces observation embedding hti,v for sensor v through its relational connection with u, even though there is no direct measurement of v at time t. These message passing operations are performed to adaptively and dynamically estimate missing observations in the embedding space based on recorded information and learned graph structure.
Updating sensor dependency graphs. We describe the update of edge weights and prune of graph structures in the situation that stacks multiple RAINDROP layers (Figure 3). Here we explicitly show layer index l because multiple layers are involved in the computation. As no prior knowledge is assumed, we initialize the graph as all sensors connected with each other. However, the fully connected edges may bridge sensors that should be independent, which will introduce spurious correlations and prevent the model from paying attention to the truly important connections. Addressing this issue, RAINDROP automatically updates edge weights and prunes out less important edges. Based on the aggregated temporal influence driven by the inter-sensor attention weights α(l),ti,uv , we update edge weights e(l)i,uv in each layer l ∈ {1, . . . , L} by:
e (l) i,uv =
e (l−1) i,uv |Ti,u| ∑ t∈Ti,u α (l),t i,uv, (3)
where Ti,u denotes the set of all timestamps where there is message passes from u to v. In particular, we set e(0)i,uv = 1 in the initialization of graph structures. We use L = 2 in all our experiments. In every layer, we order the estimated values e(l)i,uv for all edges in sample Si and prune bottom K% edges with smallest edge weights (Yang et al., 2021). Pruned edges will not re-appear in later layers.

heading: 3.4 GENERATING SENSOR EMBEDDINGS
text: Next we describe how to aggregate observation embeddings into sensor embeddings, taking sensor v as an example (Figure 3b). Previous step outlined in Sec. 3.3 outputs observation embeddings for every timestamp when either v or v’s neighbor in Si’s sensor dependency graph are observed. The following step aggregates those embeddings into a sensor embedding zi,v via temporal attention.
We concatenate observation embedding hti,v with time representation p t i to include information of timestamp. Then, we stack the concatenated embeddings [hti,v||pti] for all t ∈ Ti,v into a matrixHi,v . The Ti,v = {t1, t2, . . . , tT } includes all timestamps when a readout is observed in v (we can directly generate hti,v) or v’s neighbor (we can generate h t i,v through message passing). We calculate the temporal attention weight βti,v using self-attention (Hu et al., 2020; Yun et al., 2019; Vaswani et al., 2017). The βti,v represents the importance of observation embedding to the whole sensor embedding. The βti,v is the corresponding element of vector βi,v that is calculated through:
βi,v = softmax
( Qi,vK
T i,v√
dk s
) , (4)
where query matrix Qi,v = Hi,vWQ and key matrix Ki,v = Hi,vWK are two matrices linearly mapped fromHi,v parameterized byWQ andWK , respectively. The √ dk is a scaling factor where dk is the dimension after linear mapping. Our calculation of βi,v is different from the typical selfattention by having mapping weight vector s. The standard dot-product self-attention generates an attention matrix with dimension of T × T (where T = |Ti,v| can vary across samples) that has an attention weight for each pair of observation embeddings. In our case, we only need a single attention vector where each element denotes the attention we should pay to an observation embedding when generating the sensor embedding. Thus, we modify the self-attention model to fit our case and use a trainable s ∈ RT×1 to map the results of dot-product from matrix (RT×T ) to T -dimensional vector βi,v (RT×1) through matrix product (Appendix A.2).
We calculate sensor embedding zi,v through: zi,v = ∑ t∈Ti,v (βti,v[h t i,v||pti]W ), (5)
where weight matrixW is a linear projector shared by all sensors and samples. It’s worth to mention that all attention weights (such as αti,uv and βi,v) can be multi-head. In this work, we describe the model using a single head for brevity.
Using attentional aggregation, RAINDROP can learn a fixed-length sensor embedding for arbitrary number of observations. Meanwhile, RAINDROP is capable of focusing on the most informative
observation embeddings. We process all observation embeddings as a whole instead of sequentially, which allows parallel computation for faster training and also mitigates the performance drop caused by modeling long dependencies sequentially. In the case of sensors with very large number of observations, we can reduce the length of time series by subsampling or splitting a long series into multiple short series.

heading: 3.5 GENERATING SAMPLE EMBEDDINGS
text: Finally, we aggregate sensor embeddings zi,v (Eq. 5) to obtain an embedding zi ∈ Rdz for sample Si, as follows: zi = [zi,1||zi,2|| · · · ||zi,M ] across all sensors v = 1, 2, . . . ,M . Our preliminary experiments show that concatenation outperforms other popular aggregation functions such as averaging (Errica et al., 2021) and squeeze-excitation readout function (Kim et al., 2021; Hu et al., 2018). While any of those aggregation functions can be considered, we used concatenation throughout all experiments in this manuscript (Appendix A.3). Taken together Sec. 3.2-3.5, given an input of sample Si, RAINDROP returns a sample embedding zi that can in turn be fed into a downstream task.

heading: 3.6 IMPLEMENTATION AND PRACTICAL CONSIDERATIONS
text: Loss function. The RAINDROP’s loss function is formulated as: L = LCE + λLr, where Lr = 1 M2 ∑ u,v∈V ∑ i,j∈V ||ei,uv − ej,uv||2/(N − 1)2, where LCE is a cross-entropy classification loss, and Lr is a regularizer that encourages the model to learn similar sensor dependency graphs for similar samples. The Lr measures averaged Euclidean distance of edge weights across all samples pairs, in all sensor pairs (including self-connections). The λ is a user-defined coefficient. Practically, as N can be large, we calculate Lr only for samples in each batch. Downstream tasks. If a sample has auxiliary attributes (e.g., a patient’s demographics) that do not change over time, we can project the attribute vector to a da-dimensional vector ai with a fullyconnected layer and concatenate it with the sample embedding, getting [zi||ai]. At last, we feed [zi||ai] (or only zi if ai is not available) into a neural classifier ϕ : Rdz+da → {1, . . . , C}. In our experiments, ϕ is a 2-layer fully-connected network with C neurons at the output layer returning prediction ŷi = ϕ([zi||ai]) for sample Si. Sensor dependencies. While modeling sensor dependencies, we involve observation embedding (hti,u, Eq. 1) of each sample in the calculation of attention weights. Similarly, to model time-wise specificity in graph structures, we consider time information (pti, Eq. 1) when measuring α t i,uv. RAINDROP can capture similar graph structures across samples from three aspects (Appendix A.4): (1) the initial graphs are the same in all samples; (2) the parameters in message passing (Ru; wu, wv, Eq. 2), inter-sensor attention weights calculation (D, Eq. 1), and temporal attention weights calculation (s, Eq. 4; W , Eq. 5) are shared by all samples; (3) we encourage the model to learn similar graph structures by adding a penalty to disparity of structures (Lr). Scalability. RAINDROP is efficient because embeddings can be learned in parallel. In particular, processing of observation embeddings is independent across timestamps. Similarly, sensor embeddings can be processed independently across different sensors (Figure 3). While the complexity of temporal self-attention calculation grows quadratically with the number of observations, it can be practically implemented using highly-optimized matrix multiplication.

heading: 4 EXPERIMENTS
text: Datasets. Below we briefly overview healthcare and human activity datasets. (1) P19 (Reyna et al., 2020) includes 38,803 patients that are monitored by 34 sensors. Each patient is associated with a binary label representing the occurrence of sepsis. (2) P12 (Goldberger et al., 2000) records temporal measurements of 36 sensors of 11,988 patients in the first 48-hour stay in ICU. The samples are labeled based on hospitalization length. (3) PAM (Reiss & Stricker, 2012) contains 5,333 segments from 8 activities of daily living that are measured by 17 sensors. Details are in Appendix A.5.
Baselines. We compare RAINDROP with five state-of-the-art time series classification methods: Transformer (Vaswani et al., 2017), Trans-mean, GRU-D (Che et al., 2018), SeFT (Horn et al., 2020), and mTAND (Shukla & Marlin, 2021). The Trans-mean is an imputation method combining
transformer architecture with commonly used average interpolation (i.e., missing values are replaced by average observations in each sensor). Further details are in Section 2. The mTAND (Shukla & Marlin, 2021) method has been shown to outperform numerous recurrent models including RNNImpute (Che et al., 2018), RNN-Simple, Phased-LSTM (Neil et al., 2016), and IP-Nets (Shukla & Marlin, 2018), along with ordinary differential equations (ODE)-based models such as LATENT-ODE and ODE-RNN (Chen et al., 2018). For this reason, we compare with mTAND and do not report comparison with those techniques in this paper. Details on hyperparameter selection of RAINDROP and baselines are in Appendix A.6. Chosen evaluation metrics are presented in Appendix A.7.

heading: 4.1 RESULTS ACROSS A MULTITUDE OF EVALUATION SETTINGS
text: Setting 1: Classic time series classification. (1) Setup.t .et . We randomly split the dataset into training (80%), validation (10%), and test (10%) set. The indices of these splits are fixed across all methods. (2) Results.lt .s se lt . As shown in Table 1, RAINDROP obtains the best performance across three benchmark datasets, implying our superiority in standard time series classification scenario. In particular, in binary classification (P19 and P12), RAINDROP outperforms the strongest baselines by 5.3% in AUROC and 4.8% in AUPRC on average. In a much more challenging task of 8-class classification, RAINDROP also outperforms existing approaches by 5.7% in accuracy and 5.5% in F1 score on PAM. More experiments on exploration and comparison are shown in Appendix A.9-A.10. Comparison with extensive baselines is in Appendix A.11.

heading: P19 P12 PAM
text: Setting 2: Leave-fixed-sensors-out. (1) Setup.t .et . We evaluate whether RAINDROP can learn dependencies among sensors, thus achieve good performance even some sensors are completely missing. This setting is practically meaningful, such as facing sensor failure. or some sensors are unavailable in specific scenes. Our intuition is that RAINDROP can compensate for the missing information from nearby observations by exploiting relational dependencies. In this setting, we select a proportion of sensors, and set all their observations as zero in validation and test set (training samples are not changed). We mask out the most informative sensors which are selected through information gain (Appendix A.8). The selected sensors are fixed across samples and models. (2) Results.lt .s se lt . We report results taking PAM as an example. In Table 2 (left columns), we observe that RAINDROP achieves the highest values in 18 out of 20 settings while the missing data ratio ranges from 10% to 50%. The missing data ratio refers to the proportion of sensors we mask out. With the increase of missing data ratio, our model has greater relative gains. RAINDROP outperforms competitive baselines by up to 24.9% in accuracy, 50.3% in precision, 29.3% in recall, and 42.8% in F1 score.
Setting 3: Leave-random-sensors-out. (1) Setup.t .et . Setting 3 is similar to Setting 2 except that the missing sensors in this setting are randomly selected instead of fixed. In each test sample, we select a subset of sensors and regard them as missing through replacing all of their observations with zeros. (2) Results.lt .s se lt . We provide results of PAM in Table 2 (right columns). Similar to Setting 2, we observe that RAINDROP achieves better performance than baselines in 16 out of 20 settings while Trans-mean and GRU-D are competitive methods.
Setting 4: Group-wise time series classification. (1) Setup.t .et . To understand whether RAINDROP can adaptively adjust its structure and generalize well to other groups of samples which were not observed while training the model. In this setting we split the data into two groups, based on a specific static attribute. The first split attribute is age, where we classify people into young (< 65 years) and old (≥ 65 years) groups. We also split patients into male and female by gender attribute. Given the split attribute, we use one group as a train set and randomly split the other group into equally sized validation and test set. (2) Results.lt .s se lt . Taking P19 as an example, we present the classification results when the training and testing samples are from different groups. As shown in Table 3, RAINDROP achieves the best results over all of the four given cross-group scenarios (Appendix A.13). For instance, RAINDROP claims large margins (with 4.8% in AUROC and 13.1% in AUPRC absolute improvement) over the second best model while training on males and testing on female patients.

heading: 4.2 ABLATION STUDY AND VISUALIZATION OF LEARNED GRAPH STRUCTURE
text: Ablation study. Taking PAM in the most common task (Setting 1) as an example, we conduct an ablation study to evaluate the necessity of three fundamental steps of RAINDROP: inter-sensor dependency (further decomposed to key weights including ei,uv , rv , pti, andα t i,uv), temporal attention, and sensor-level concatenation. In Appendix A.14 (Table 9), we observe that all of the components and the regularization term Lr are important to improve model performance. Visualizing learned sensor dependency graphs. We investigate whether samples with the same label have similar graph structure by visualizing the learned inter-sensor dependencies (P19; Setting 1). As shown in Figure 4, we can observe distinguishable patterns between the graph structures of negative and positive samples, indicating RAINDROP can adaptively learn relationships that are sensitive to the classification task. The differential analysis can also provide inspirations to domain experts such as intensivists. More details in Appendix A.15.

heading: 5 CONCLUSION
text: We develop RAINDROP, a graph-guided network for learning with irregularly sampled multivariate time series. RAINDROP learns a distinct sensor dependency graph for every sample that captures time-varying dependencies between sensors in the sample. The ability to leverage the relational structure gives RAINDROP unique ability to naturally handle misaligned observations, non-uniform time intervals between successive observations, and sensors with a varying number of recorded observations. Our novel findings have implications for using message passing as a powerful approach to take advantage of relational structure in multivariate time series.

heading: REPRODUCIBILITY STATEMENT
text: We ensure the reproducibility of our work by clearly presenting the model and provide publicly accessible code and data. For all datasets used in this work, we share the downloadable links to the raw sources, and well-processed and ready-to-run datasets, with the research community through link https://anonymous.4open.science/r/Raindrop. We specify all the training details (e.g., preprocessing, data splits, hyperparameters, sensor selection) in the main text and Appendix. Python implementation of RAINDROP and all baseline methods is available at the aforementioned link. Detailed description of all data and scripts, configurations, along with examples of usage, will also be provided.

heading: ETHICS STATEMENT
text: The ability of RAINDROP model to learn powerful information about sensors’ representations and dependencies creates new opportunities for applications, where time series data is predominant, e.g., in healthcare, biology, and finance. In all these fields, especially in patient care, our method should be used with caution. Although our model can gain useful insights from time series data, users must be aware of the restrictions in terms of machine-guided predictions. The use of our model in real-life clinical settings should always be verified by qualified physicians. As with all data-driven solutions, it is possible that our model makes biased predictions. In the case of biomedical data, biases can exist within the data itself, which can be, for example, caused by considering static attributes, such as age, weight, and gender. When the ratio of target labels is highly imbalanced, our model mitigates the possible bias by performing minority class upsampling for every processed batch.
All datasets in this paper are publicly available and are not associated with any privacy or security concerns. Further, all data are anonymized to guard against breaching patient’s protected health information (PHI). We have followed PhysioNet privacy policy and guidelines (https://archive. physionet.org/privacy.shtml) when experimenting with P12 and P19 datasets.

heading: A APPENDIX

heading: A.1 TIMESTAMP ENCODING
text: For a given time value t, we pass it to trigonometric functions with the frequency of 10,000 (Vaswani et al., 2017) and generate time representation pt ∈ Rξ (omit sample index i for brevity) through (Horn et al., 2020):
pt2k = sin( t
100002k/ξ ), pt2k+1 = cos(
t
100002k/ξ ), (6)
where ξ is the expected dimension. In this work, we set ξ = 16 in all experimental settings for all models. Please note, we encode the time value which is a continuous timestamp, instead of time position which is a discrete integer indicating the order of observation in time series.
A.2 INTUITION OF VARIABLE s
Here we report the intuition and the meaning for the designed transformation variable s in Eq. 4. The target of Eq. 4 is to learn the temporal attention weights vector βi,v for sensor v, taking advantage of self-attention mechanism. However, the standard self-attention generates an attention matrix but we ask for an attention vector. In the attention matrix, each row represents the importance of the observation embedding at a certain time to the observation embeddings at all other time steps.
Our intuition is to aggregate a row in the self-attention matrix into a scalar that denotes the importance of the observation embedding at a specific time step to the whole sensor embedding. We apply the weighted aggregation, parameterized by s, to every row in the attention matrix and concatenate the generated scalars into an attention vector. Next, we give a concrete example to specifically describe the meaning of s. Each row, j, of the self-attention matrix captures relationships of observation embedding htji,v to all observation embeddings {h tk i,v : k = 1, ..., T}. Then, using the learnable weight vector s, these correlations between observations are aggregated across time to obtain temporal importance weight βtji,v. The β tj i,v represents the importance of the corresponding observation to the whole sensor embedding.

heading: A.3 ADDITIONAL INFORMATION ON SAMPLE EMBEDDING
text: As we generate sample embedding by concatenating all sensor embeddings, the sample embedding could be relatively long when there is a large number of sensors. To alleviate this issue, on one hand, we can reduce the dimension of sample embeddings by adding a neural layer (such as a simple fully-connected layer) after the concatenation. On the other hand, when the number of sensors is super large, our model is flexible and can effortlessly switch the concatenation to other readout functions (such as averaging aggregation): this will naturally solve the problem of long vectors. We empirically show that concatenation works better than averaging in our case. We see a boost in the AUROC score by 0.6% using concatenation instead of averaging for generating sample embeddings(P19; Setting 1).

heading: A.4 ADDITIONAL INFORMATION ON SAMPLE SIMILARITIES
text: In this work, we assume all samples share some common characteristics to some extent. When modeling the similarities across samples, we do not consider the situation where the samples are similar within latent groups and different across groups.
Our study focuses on the question of irregularity rather than the question of distribution shifts in time series. To this end, in our experiments, we first rigorously benchmark Raindrop using a standard evaluating setup (Setting 1, which is classification of irregular time series). This is the only setup that most existing methods consider (e.g., Shukla & Marlin (2021); Che et al. (2018)) and we want to make sure our comparisons are fair. In order to provide a more rigorous assessment of Raindrop’s performance, we also consider more challenging setups in our experiments (i.e., Settings 2-4) when the dataset is evaluated in a non-standard manner and the split is informed by a select data attribute. Our results on Setting 1 are consistent with those on Settings 2-4. Results on harder Settings 2-4 show that Raindrop can perform comparably better than baselines. Results across these diverse settings increase our confidence that Raindrop is quite flexible and widely applicable.

heading: A.5 FURTHER DETAILS ON DATASETS
text: P19: PhysioNet Sepsis Early Prediction Challenge 2019. P19 dataset (Reyna et al., 2020) contains 38,803 patients and each patient is monitored by 34 irregularly sampled sensors including 8 vital signs and 26 laboratory values. The original dataset has 40,336 patients, we remove the samples with too short or too long time series, remaining 38,803 patients (the longest time series of the patient has more than one and less than 60 observations). Each patient is associated with a static vector indicating attributes: age, gender, time between hospital admission and ICU admission, ICU type, and ICU length of stay (days). Each patient has a binary label representing occurrence of sepsis within the next 6 hours. The dataset is highly imbalanced with only ∼4% positive samples. P12: PhysioNet Mortality Prediction Challenge 2012. P12 dataset (Goldberger et al., 2000) includes 11,988 patients (samples), after removing 12 inappropriate samples following (Horn et al., 2020). Each patient contains multivariate time series with 36 sensors (excluding weight), which are collected in the first 48-hour stay in ICU. Each sample has a static vector with 9 elements including age, gender, etc. Each patient is associated with a binary label indicating length of stay in ICU, where negative label means hospitalization is not longer than 3 days and positive label marks hospitalization is longer than 3 days. P12 is imbalanced with ∼93% positive samples. PAM: PAMAP2 Physical Activity Monitoring. PAM dataset (Reiss & Stricker, 2012) measures daily living activities of 9 subjects with 3 inertial measurement units. We modify it to suit our scenario of irregular time series classification. We excluded the ninth subject due to short length of sensor readouts. We segment the continuous signals into samples with the time window of 600 and the overlapping rate of 50%. PAM originally has 18 activities of daily life. We exclude the ones associated with less than 500 samples, remaining 8 activities. After modification, PAM dataset contains 5,333 segments (samples) of sensory signals. Each sample is measured by 17 sensors and contains 600 continuous observations with the sampling frequency 100 Hz. To make time series irregular, we randomly remove 60% of observations. To keep fair comparison, the removed observations are randomly selected but kept the same for all experimental settings and approaches. PAM is labelled by 8 classes where each class represents an activity of daily living. PAM does not include static attributes and the samples are approximately balanced across all 8 categories.
To feed given data into neural networks, we set the input as zero if no value was measured. In highly imbalanced datasets (P19 and P12) we perform batch minority class upsampling, which means that every processed batch has the same number of positive and negative class samples. The dataset statistics including sparse ratio are provided in Table 4.

heading: A.6 FURTHER DETAILS ON MODEL HYPERPARAMETERS
text: Baseline hyperparameters. The implementation of baselines follows the corresponding papers including SeFT (Horn et al., 2020), GRU-D (Che et al., 2018), and mTAND (Shukla & Marlin, 2021). We follow the settings of Transformer baseline in (Horn et al., 2020) while implementing Transformer in our work. For average imputation in Trans-mean, we replace the missing values by the global mean value of observations in the sensor (Shukla & Marlin, 2020). We use batch size of 128 and learning rate of 0.0001. Note that we upsample the minority class in each batch to make the batch balance (64 positive samples and 64 negative samples in each batch).
The chosen hyperparameters are the same across datasets (P19, P12, PAM), models (both baselines and RAINDROP), and experimental settings. Remarkably, we found that all the baselines make
dummy predictions (classify all testing samples as the majority label) on PAM in Setting 2-3 while RAINDROP makes reasonable predictions. For the comparison to make sense (i.e., the baselines can make meaningful predictions), we use learning rate of 0.001 for baselines on PAM. GRU-D has 49 layers while other models have 2 layers. We run all models for 20 epochs, store the parameters that obtain the highest AUROC in the validation set, and use it to make predictions for testing samples. We use the Adam algorithm for gradient-based optimization (Kingma & Ba, 2014).
RAINDROP hyperparameters. Next, we report the setting of unique hyperparameters in our RAINDROP. In the generation of observation embedding, we setRu as a 4-dimensional vector, thus the produced observation embedding has 4 dimensions. The dimensions of time representation pt and rv are both 16. The trainable weight matrixD has shape of 4× 32. The dimensions of wu and wv are the same as the number of sensors: 34 in P19, 36 in P12, and 17 in PAM. We set the number of RAINDROP layers L as 2 while the first layer prunes edges and the second layer does not. We set the proportion of edge pruning as 50% (K=50), which means we remove half of the existing edges that have the lowest weights. The dk is set to 20, while the shape of W is 20 × 20. All the activation functions, without specific clarification, are sigmoid functions. The da is set equal to the number of sensors. The first layer of ϕ has 128 neurons while the second layer has C neurons (i.e., 2 for P19 and P12; 8 for PAM). We set λ = 0.02 to adjust the scale of regularization term Lr. All the preprocessed datasets and implementation codes are made available online. More details can be found through our publicly accessible link.

heading: A.7 PERFORMANCE METRICS
text: Since P19 and P12 datasets are imbalanced, we use the Area Under a ROC Curve (AUROC) and Area Under Precision-Recall Curve (AUPRC) to measure performance. As the PAM dataset is nearly balanced, we also report accuracy, precision, recall and F1 score. We report mean and standard deviation values over 5 independent runs. Model parameters that achieve the best AUROC value on the validation set are used for test set.

heading: A.8 FURTHER DETAILS ON SETUP DETAILS FOR SETTING 2
text: In Setting 2, the selected missing sensors are fixed across different models and chosen in the following way. First, we calculate the importance score for each sensor and rank them in a descending order. The importance score is based on information gain, which we calculate with feeding the observations into a Random Forest classifier with 20 decision trees. In particular, we treat each sample as only having one sensor, then feed the single sensor into random forest classifier and record the AUROC. The higher AUROC indicates the sensor provides higher information gain. When we have sensors ranked by their AUROC values, we choose the first n sensors (the ones with highest AUROC values) and replace all observations in these sensors by zeros in all samples in validation and test set. The number of missing sensors is defined indirectly from the user with the sensors’ missing ratio which ranges from 0.1 to 0.5.

heading: A.9 ADDITIONAL INFORMATION ON MISSING PATTERN
text: This work propose RAINDROP which is a novel solution for irregularity in multivariate time series through inter-sensor dependencies. RAINDROP is not in conflict with other solutions (such as missing pattern and temporal decay) for irregularity. However, as the missing pattern is widely discussed in modelling incomplete time series (Che et al., 2018), we explore how to combine the advantages of relational structures and missing pattern. We adopt mask matrix as a proxy of missing pattern as in Che et al. (2018). Taking the architecture of RAINDROP, we concatenate the observation xti,u with a binary mask indicator bt[ i, u] as input. The indicator b t [ i, u] is set as 1 when there is an observation of sensor i at time t and set as 0 otherwise. All the experimental settings and hyperparameters are the same as in RAINDROP (P19; Setting 1). The experimental results show that taking advantage of missing pattern can slightly boost the AUROC by 1.2% and AUPRC by 0.9% in P19. This empirically shed the light for future research on integrating multiple characteristics in representation of irregularly time series.

heading: A.10 COMPARISON BETWEEN TEMPORAL ATTENTION AND LSTM
text: We conduct extensive experiments to compare the effectiveness of temporal attention and LSTM. To this end, we replace the temporal attention in sensor embedding generation (Eq 4-5) in RAINDROP by LSTM layer which processes all observation embeddings sequentially. We use zero padding to convert the irregular observations into fixed-length time series so the data can be fed into LSTM architecture. We regard the last output of LSTM as generated sensor embedding. The number of LSTM cells equal to the dimension of observation embedding. All the model structures are identical except in the part of temporal attention and LSTM. We keep all experimental settings (P19; Setting 1) and hyperparameter selections the same. The experimental results show that the temporal self-attention outperform LSTM by 1.8% (AUROC) and additionally saved 49% of the training time. One potential reason is that the self-attention mechanism avoids recursion and allows parallel computation and also reduces performance degradation caused by long-term dependencies (Ganesh et al., 2021; Vaswani et al., 2017).

heading: A.11 ADDITIONAL INFORMATION ON METHOD BENCHMARKING
text: We conduct extensive experiments to compare Raindrop with ODE-RNN (Chen et al., 2020), IPNet (Shukla & Marlin, 2018), DGM2-O (Wu et al., 2021), EvoNet (Hu et al., 2021), and MTGNN (Wu et al., 2020c). For the five newly added baselines, we follow the settings as provided in their public codes. For methods, which cannot deal with irregular data (i.e., EvoNet and MTGNN), we first impute the missing data using mean imputation and then feed data into the model. For forecasting models (i.e., MTGNN) which are strictly not comparable with the proposed classification model, we formulate the task as a single-step forecasting, concatenate the learned representations from all sensors and feed into a fully-connected layer (work as classifier) to make classification, and use cross-entropy to quantify the loss. The results are from the PAM dataset under experimental Setting 1. As shown in the following Table 5, the experimental results confirm that Raindrop outperforms the strongest baseline in accuracy by 4.7%. The added experiments increase our confidence in the effectiveness of RAINDROP.

heading: A.12 RESULTS FOR P19 (SETTINGS 2-3)
text: Here we report the experimental results for P19 in Setting 2 (Table 6) and Setting 3 (Table 7).

heading: A.13 ADDITIONAL INFORMATION ON GROUP-WISE CLASSIFICATION
text: In experimental Setting 4, we evaluate RAINDROP under a challenging scenario. Although RAINDROP is not designed to address domain adaptation explicitly, the results show that RAINDROP performs better than baselines when transferring from one group of samples to another. One reason for our good performance is that the learned inter-sensor weights and dependency graphs are sample-specific and their learning is based on the sample’s observations. Thus, the proposed RAINDROP has the power, to some extent, to adaptively learn the inter-sensor dependencies based on the test sample’s measurements. RAINDROP is not generalizing to new groups, but generalizing to new samples, which leads to a good performance even though our model is not designed for domain adaptation.
The reason is validated empirically. We remove the inter-sensor dependencies (set all sensors isolated in the dependency graph; set all αti,uv and e t i,uv as 0) in RAINDROP and evaluate the model in Setting 4. The experimental results show that the performance drops a lot when excluding dependency graphs and message passing in RAINDROP (Table 8). Without inter-sensor dependencies our model is on par with other baselines and does not outperform them by a large margin.

heading: A.14 FURTHER DETAILS ON ABLATION STUDY
text: We provide ablation study, taking PAM at Setting 1 as an example, in Table 9. In the setup of ‘W/o sensor level concatenation’, we take the average of all sensor embeddings (in stead of concatenating them together) to obtain sample embedding. Experimental results show that the full RAINDROP model achieves the best performance, indicating every component or designed structure is useful to the model. For example, we find that excluding inter-sensor attention weights αti,uv will cause a decrease of 3.9% in accuracy while excluding edge weights ei,uv (i.e., dependency graphs) will drop the accuracy by 7.1%.
A.15 VISUALIZATION OF INTER-SENSOR DEPENDENCY GRAPHS LEARNED BY RAINDROP
We visualize the learned inter-sensor dependencies (i.e., ei,uv before the averaging operation in Eq. 3) on P19 in early sepsis prediction. The visualizations are implemented with Cytoscape (Shannon et al., 2003). The data shown are for testing set of P19 including 3,881 samples (3708 negative and 173 positive). As RAINDROP learns the specific graph for each sample, we take average of all positive samples and visualize it in Figure 4b; and visualize the average of all negative samples in Figure 4b. As we take average, the edges with weights smaller than 0.1 (means they rarely appear in graphs) are ignored. The averaged edge weights range from 0.1 to 1. We initialize all sample graphs as complete graph that has 1,156 = 34× 34 edges, then prune out 50% of them in training phase, remaining 578 edges. The 34 nodes in figures denote 34 sensors measured in P19, as listed https://physionet.org/content/challenge-2019/1.0.0/. We list the sensor names here: 0: HR; 1: O2Sat; 2: Temp; 3: SBP; 4: MAP; 5: DBP; 6: Resp; 7: EtCO2; 8: BaseExcess; 9: HCO3; 10: FiO2; 11: pH; 12: PaCO2; 13: SaO2; 14: AST; 15: BUN; 16: Alkalinephos; 17: Calcium; 18: Chloride; 19: Creatinine; 20: Bilirubin_direct; 21: Glucose; 22: Lactate; 23: Magnesium; 24: Phosphate; 25: Potassium; 26: Bilirubin_total; 27: TroponinI; 28: Hct; 29: Hgb; 30: PTT; 31: WBC; 32: Fibrinogen; 33: Platelets.
We also visualize the differential inter-sensor connections between the learned dependency graphs from patients who are likely to have sepsis and the graphs from patients who are unlikely to suffer from sepsis. Based on the aggregated graph structures of positive and negative samples, we calculate
the divergence between two groups of patients and report the results in Figure 5. In detail, we sort edges by the absolute difference of edge weights across negative and positive samples. On top of the visualization of the 50 most distinctive edges, we can have a series of concrete insights. For example, the dependency between node 6 (Respiration rate) to node 25 (Potassium) is important to the early prediction of sepsis. Note these data-driven observations could be biased and still need confirmation and future analysis from healthcare professionals. The edges in both Figure 4 and Figure 5 are directed. The edge arrows might be difficult to recognize due to the small figure size. We will provide high-resolution figures to our public repository.
Furthermore, we statistically measure the similarities across samples within the same class and dissimilarities across samples from different classes. Specifically, for every sample, we calculate: 1) the average Euclidean distance between its dependency graph and the dependency graphs of all samples from the same class; 2) the average distance with all samples from the different classes. The P19 dataset has 38,803 samples including 1,623 positive samples and 37,180 negative samples. For a fair comparison, we randomly select 1,623 samples from the negative cohort, then mixed them with an equal number of positive samples to measure the averaged Euclidean distances intra- and inter-classes. We select the cohort for 5 independent times with replacement. We find that the distance ((8.6 ± 1.7) × 10−5) among dependency graphs of positive samples is smaller than the distance ((12.9 ± 3.1) × 10−5) across samples. The results show that the learned dependency graphs are similar within the same class and dissimilar across classes, which demonstrates RAINDROP can learn label-sensitive dependency graphs.

Review summaries:
This paper introduces a GNN based method for classifying irregular multivariate time series (MTS). The proposed method explicitly models the dependency between sensors, and propagates impacts of intermittent time series observations from sensors to related sensors for learning MTS embeddings. It jointly learns embeddings and dependency graphs through several attentive mechanisms. The paper provides a good perspective of applying GNN methods on MTS related tasks.

This manuscript describes a method that models multiple time series data (sensors) across different individuals (samples). A graph is used to model the dependency between the sensors, such that the data of one sensor could be used to infer another, which is the key point of this ms. How one sensor affects another is modelled as a message passing problem on the graph. At each time point, an observation embedding is generated, either from a sensor or its neighbours in the graph. Observation embeddings of a sensor is then aggregated into a fixed-length sensor embedding, using temporal attention. These sensor embeddings are concatenated into a sample embedding for an individual for downstream analysis such as classification.

The core idea of the paper is transforming multivariate time series data into a vector z and using z for classification. The transformation has two steps. Firstly, each time record of each sensor is mapped into a  space of 20 dimensions. The first 4 dimensions are used for feature extraction, and the last 16 dimensions are used for positional encoding. Secondly, self-attention is used to transform the embedding into a vector z and use it for classification.