Title: CLINICAL RISK: WAVELET RECONSTRUCTION NETWORKS FOR MARKED POINT PROCESSES

Abstract: Timestamped sequences of events, pervasive in domains with data logs, e.g, health records, are often modeled as point processes with rate functions over time. Leading classical methods for risk scores such as Cox and Hawkes processes use such data but make strong assumptions about the shape and form of multivariate influences, resulting in time-to-event distributions irreflective of many real world processes. Recent methods in point processes and recurrent neural networks capably model rate functions but may be complex and difficult to interrogate. Our work develops a high-performing, interrogable model. We introduce wavelet reconstruction networks, a multivariate point process with a sparse wavelet reconstruction kernel to model rate functions from marked, timestamped data. We show they achieve improved performance and interrogability over baselines in forecasting complications and scheduled care visits in patients with diabetes.

1. INTRODUCTION:

Clinical risk scores are commonly used analytic devices in health care. There are risk scores for predicting strep throat from sore throats (Centor et al., 1981) , mortality from vital signs (Gardner-Thorpe et al., 2006) , heart attacks from routine clinic visits (D'Agostino et al., 2008) , and many more. Policy is implemented around these risk scores, from rates of reimbursement to physician compensation (Asch et al., 2015) . When used for early warning, risk scores have been associated with reduced mortality (Seymour et al., 2017) .

Underlying these approaches is the formulation of risk over time given some set of features. For example, in Cox models, a study time t 0 = 0 is defined and the risk model is λ(t; t , x), where x are defined by features timestamped before time t = t 0 , i.e., a time invariant model. In Hawkes processes, the risk model is also λ(t; t , x) but where x contains all history up to time t = t, i.e., a nowcasting model. Cox and Hawkes models can be limited by their assumptions, which often are inappropriate for the health care setting, and include, for example, the assumptions of proportional hazards and summation over kernel activations. We formulate a point process model to address these limitations and develop multi-forecasting, the forecasting task across the two dimensions of time: t and t.

To motivate our specific formulation, consider the limitations of the Hawkes process in health care. First, the Hawkes process encodes an additive relationship of change in rate from recurring precursors, i.e., burstiness, whereas in health care, the repeated measurement of an event beyond the first, say of glucose, might be irrelevant. Second, clinical event timing may be routine, scheduled, or emergent, which suggests that kernel learning will improve model performance because changes in the rate may be time-dependent and not immediate. Third, clinical event processes are marked, with marks that could be categorical, real, or null values: e.g., bacterial culture: staphylococcus aureus, glucose: 200, and ketoacidosis: NA Our model addresses each of these limitations. To address the first limitation, summation, we adopt a reduction layer where we allow for reductions other than "sum" of kernel contributions from recurring events. To address the second and third limitations, non-specific timing and lack of marks, we propose a kernel learning method over one dimension (time) and two dimensions (time-value) using wavelet reconstructions. The motivation for wavelets is illustrated in Figure 1 , where a discrete wavelet reconstruction encodes the relationship of time-delayed events identified through cross-correlation. While cross-correlations capture all relative timings, many may be spurious or coincidental. The wavelet representation of relative timing instead is learned through likelihood optimization. To capture the effect of the value distribution of marks on events, two-dimensional wavelet reconstructions are value indexed, producing a one-dimensional reconstruction that is passed forward to the reduction layer. We show that we can encode the wavelets, the relative-to-absolute mapping, and the reduction step on a computation graph to conduct learning. We apply our model first in simulations of heart attacks and scheduled hemoglobin A1c checks. We then evaluate the model in a real cohort to forecast complications and adherence to clinical participation in patients with diabetes. Empirically, we show that our model has improved performance over baselines in single forecasts and explore its performance in multi-forecasts in terms of prediction and representation.

Point process:

Our results provide benefit for diabetes risk assessment. Whereas the central focus of the diabetes adherence literature is on adherence to medication and therapeutic regimen (García-Pérez et al., 2013; Edelman and Polonsky, 2017) , our method provides a lens to investigate adherence to continued diabetes care participation, which is important because, for one, participation is associated with better diabetes control (Schectman et al., 2008) . Additionally, forecasts of complications and adherence to participation, at different times in the care process and at different forecast distances, are central concerns in clinical decision-making. For example, many diabetes medications require regular clinical monitoring, and the risk of complications and noncompliance to clinical participation both affect regimen choice. By focusing on these outcomes, our forecasting tools provide relevant information for the clinical decision-making process.

Related work. Both neural networks and Hawkes process variants are used for rate modeling in health care: a few recent examples include Choi et al. (2015) , Du et al. (2016 ), Alaa et al. (2017 ), and Bao et al. (2017) . To effectively model the hazard two key properties are used: (1) the ability of Hawkes processes to capture relative timing of interdependent events, and (2) the flexible functional forms of neural networks that are able to capture relative timings, albeit somewhat opaquely. The closest work is likely that of Bao et al. (2017) , where the authors adopt dyadic influence functions. However in that work marks are not used and the dyads selected are a subset of the Haar wavelet basis. Outside of health applications, related literature includes Hawkes kernel learning, e.g. Zhou et al. (2013) ; Linderman and Adams (2014) ; Lee et al. (2016) . Our method follows these approaches, instead using a wavelet representation across value and time to capture long-range dependencies. Compared to Mei and Eisner (2017) that generalizes the Hawkes process with neural networks, our method (1) allows events to have marks, (2) enables forecasting and multi-forecasting (not just nowcasting), and (3) can be seen to have a generalized linear form.

Our method of mapping relative timing hazard components onto absolute time possesses the advantages of the Hawkes approach and adopts a simple but well-performing neural network architecture. To achieve the mapping, a stepwise hazard approximation is made, as is done in Jing and Smola (2017) and Weiss (2017) , however, instead of LSTMs and forests that have challenging interpretations, our method remains interpretable for small data sets, where exploration with visualization similar to that of Caruana et al. (2015) can be performed. Our method uses wavelets to represent event contributions, and several survival analysis approaches have also adopted them in univariate models, e.g., in Antoniadis et al. (1994) and Brillinger (1997) . Outside of survival analysis and point processes, wavelet-inspired neural networks have seen success, with Wave-net using wavelets to classify time series (Bakshi and Stephanopoulos, 1993) , and Wavenet adopting a multi-layer hidden neural architecture to connect distant time steps (Van Den Oord et al., 2016) .

Contributions. The contributions of this work are as follows: our work generalizes multivariate Hawkes processes to allow for non-additive event rate relationships. Like other works, e.g. Linderman and Adams (2014) , our work learns the kernel function that relates multivariate event histories to the rate. However, in our work we use wavelets as the kernel, akin to a multivariate development from Brillinger (1997) . We leverage the scaling property of wavelets to formulate a regularization that balances spatiotemporal generalizability with deterministic or near-deterministic event timing. Unlike many sequence models, e.g. Hochreiter and Schmidhuber (1997) , which are affected strongly by choice of time step, our work adopts an absolute and relative time frame, and therefore the granularity of the absolute time domain need not be determined a priori. Additionally, unlike some point process formulations, our work models marks that are 1-d (event times) and 2-d (event times and their category or real value). Our work is explicit about types of forecasting tasks, and the methods are adapted for several purposes: forecast performance, interrogability, and representation. We show that our method performs as well as or better than comparison algorithms in predicting complications and forecasting adherence.

2. BACKGROUND:

Let E be the set of events with target event y ∈ E the event we want to forecast. Associated with each event e is a value v ∈ V . An example consists of a sequence of (time, event, value) tuples and a period of interest for forecasting. For the n-th example, n ∈ {1 . . . N }, define T n as the number of tuples. Then the sequence can be written as (t in , e in , v in ) for i ∈ {1 . . . T n }, with the period of interest denoted as τ ny .

Let λ y (t) be the rate functions of interest, dropping the subscript n for ease of notation. The multivariate Hawkes process can then be written as follows:

λ y (t) = λ 0 (t) + |E| e=1 β e T i=1 g e (t -t i )1(t i < t, e i = e)

where λ 0 (t) is a baseline population rate function, g e (•) is a kernel function for event e relating its effect on the rate of y, β e are event-specific parameters, and 1(•) is the indicator function. Typically g e (•) is an event-specific exponential decay function with a learnable decay parameter. Self-exciting processes are defined by g y (•) > 0, bursty processes by g e (•) > 0, and inhibitory processes by g e (•) < 0. A few recent variations include Linderman and Adams (2014) where g e is a Bayesian graph kernel and Xu et al. (2017) where g e is an infectivity function and triggering kernel product.

Given λ ny (t), the log likelihood of the data is:

EQUATION

The form of the Hawkes process is limiting, however, because (1) the effect of g e (•) decays over time, (2) the effect over g e (•) is additive, (3) the value associated with each event is not considered, and (4) the time restriction in the indicator function implies nowcasting (1(t i < t)) not forecasting (1(t i < t -c) for some c > 0). Making modifications to achieve these characteristics is desirable to effectively model many real-world processes. For example, a patient with new-onset diabetes schedules an appointment with a typical gap interval of 3 months, and the presence of 1 or 10 elevated readings may not affect the timing of the scheduled appointment. The former suggests the utility of kernel learning, and the latter suggests summation over g e is not the appropriate reduction.

The proposed method addresses these concerns. In particular, we adopt discrete wavelet reconstructions, which both allows kernel learning and the use of marks. Additionally, by representing the point process as a neural network, we are able to (1) use maximization alongside summation in a reduction layer (the formulation enables specification of any number of reductions), and (2) conduct time-dependent censoring. We formalize the model below.

3. WAVELET RECONSTRUCTION NETWORKS:

We now define wavelet reconstruction networks (WRNs) shown in Figure 2 . We specify the form of the rate function, define our kernel function, and impose restrictions on the kernel function for forecasting and multi-forecasting. indicates the reduction over T functions over the interval τ y . We propose the following rate function:

h j=2 (w,vi) w j=1 w j=2 ɸ j=1 (w) ɸ j=2 (w,∘) t Reduce (Tmax → R) v ɸ q, C h j=1 (w) t -ti E E E, C, Tmax Permute pool / linear (𝛽) Reduce (E,R → 1) + 𝜆0 X = {vi, ei, ti} NLL Y = {ei = y, ti}

λ y (t) = λ 0 (t) + |E| e=1 R,J r,j β erj r τy i={1,...,T } h ej (t; t i , τ d )1(e i = e) .

(2)

Note that the kernel function h ej is on absolute time, whereas g e is defined on relative time in the previous section. This is done to specify the translation of the wavelet reconstructions on discrete relative-time intervals onto discrete absolute-time intervals, which requires additional treatment to prevent causal leakage. We define the kernel function on relative time below, followed by the definition of h ej using the translation function.

Discrete wavelet reconstruction kernel. Recall that the discrete wavelet transform (DWT) is an invertible transform of a signal between a time space and a time-frequency space used in multiresolution analysis and signal compression (Mallat, 1989 ). Here we use the inverse DWT to encode the parameters in the time-frequency space to reconstruct the signal as the kernel function as follows.

For each event, we use one-to two-dimensional wavelets, with j = 1 referring to time reconstructions, and j = 2 referring to time and event value reconstructions. We use discrete wavelets of size (α) and (α, |S|), with α time intervals on the interval [0, max n (max tiny t iny -min tine t ine )] on the time dimension, and with s ∈ S disjoint value intervals in [min(v), max(v)] on the value dimension. For point events, we use the notation s ∈ S where |S| = 1, and categorical events are treated as separate point events. Note that j = 1 reconstructions are the temporal analogues of missingness indicators.

Let parameters w ej be the wavelet coefficient tensors for event e and wavelet reconstruction dimension j, and let g ejs be the kernel function on relative time for event e and interval s, with g ejs (t -t i ) = 0 for t -t i < 0. Define the set of wavelet reconstruction functionals by: Φ = {φ ej : (w ej , v e ) → g ejs }.

Conceptually, given event value v e , the wavelet reconstruction functional φ ej reconstructs the signal from w ej and indexes the value dimension with v e , producing function g ejs , a function with inputs of relative time.

Relative-to absolute-time transformations. To relate the absolute-time kernel h ej with the relativetime discrete kernel g ejs , we define the following causally-protective translation function. Let τ d denote disjoint caglad intervals that comprise τ y the target interval in absolute time, and let τ e be the relative-time intervals of the wavelet reconstruction. We denote lower and upper endpoints with • and • . Then, for event time t i , the absolute wavelet reconstruction intervals have endpoints τ ie = τ e + t i and τ ie = τ e + t i . The transformation q is given by:

h ej (t; t i , τ d ) = q g ejs (t-t i ), τ d = τie∧τ d =∅, τie ≥ τ d (max( τ d , τ ie ) -τ ie )g ejs (t -t i ) τ d -τ d ,

where ∧ denotes interval intersection. The second condition for inclusion in summation, τ ie ≥ τ d , prevents causal leakage by ignoring intervals τ ie that affect an interval τ d that both precedes and intersects it. One advantage of our relative-time specification is that the granularity over absolute time can be adjusted with small effect on the hazard. By comparison, RNNs and their analogues need to be retrained if there are changes to the time-step specification. Figure 3 (left) illustrates this with an absolute-time hazard function of a single trajectory using different absolute time steps and applying censorship at different times, resulting in absolute-time hazards similar but not identical. Sharp discontinuities near target event times can result in relatively large likelihood differences, and the ability to choose the absolute time granularity τ d post-training facilitates hazard recovery as needed. Figure 3 (left) also acts as a causality leakage check by demonstrating that the hazard is unaffected by the presence or censorship of future events.

Forecasting. Thus far we have specified h ej as a nowcasting kernel because g ejs is zero and noncontributory when t-t i < 0 for all e and s. For forecasting, we incorporate time-dependent censoring with functional C, the Hadamard censor, and hyperparameter c the censoring distance, to prevent recent events, i.e., risk modifiers, from affecting the rate. Let C(c) t; t i equal 1 if t -t i > c and 0 otherwise, and let ψ ejc (w ej , v e ) = ψ ej (c)(w ej , v e ) = C(c) • φ ej (w ej , v e ) where • is the Hadamard product. We can specify the forecasting rate function analogous to Equation 2 as follows:

EQUATION

For multi-forecasting, we choose a vector of desired relative forecast times {c} and maximize the average log-likelihood over all c. This may be distinguished from training separate forecasting models because the parameters of the model are tied. Figure 3 (right) contrasts the forecasting tasks.

Compared to single forecasts along any given line, multi-forecasting captures more of the valid forecast region (shaded gray) that is relevant for clinical decision making.

Learning. The parameters of the model are Θ = {w ej , β erj }. Because the system may be overdetermined, we add regularization terms. The first is γ β e rj ||β erj || 1 akin to the LASSO (elastic-net regularization is equally straightforward). The second is the regularizer γ w ej ||u(w ej )|| 1 akin to sparse shrinkage on the wavelet tensor with a choice for u.

We define u(w ej ) = k∈{1,...,j} 2 l k /2 • w ej , where l k is the wavelet scale parameter of the k-th dimension. The idea is that regularization on wavelets for point events corresponds to smoothing a function of Dirac deltas over time, and we want the log loss effect of a Dirac delta (an element of the first term in Equation 1) to be in proportion to the activation of the unnormalized wavelet basis function so that the data drive the choice of smoothness. To do so, the regularization must be in proportion to k∈{1,...,j} 2 l k /2 . An example is the orthonormal two-level Haar wavelet, where the orthonormal transformation matrix is written as the Hadamard product of expanded, exponentiated scale parameters and unnormalized basis functions, where 1 is a column vector of ones:

1 • 1 1 1 1 1 1 -1 -1 1 -1 0 0 0 0 1 -1

. NA q q q q q Time 0 1 2 3 4 5 7 markStart Troponin ACS markEnd NA q q q q q Time 0 1 2 3 4 markStart Troponin ACS markEnd q qqqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Time 0 5 10 15 20 Improving prediction. The formulation in Equation 3 can be seen as a generalized linear model of censored, add-/max-reductions of wavelet reconstructions. While the generalized linear form lends itself to interpretation, we consider whether non-linearities will further improve predictive performance. Unlike images, where local pixels are spatially related, any given ordering of clinical event types may not capture valuable event relationships. Therefore we introduce permute-and-pool layers (WRN-PPL) that randomly permute event ordering within time step, randomly select sign, perform max-pool, and project linearly to the next layer. In place of the double summation in Equation 3, we apply a random sign ({-1, 1}) Hadamard tensor Z and pass the result to P parallel permutation layers with max pools of size min(2 p , |E||J||R|) for p = 0 to P -1. The outputs of the max pool are then linearly combined and output to the next layer.

4. EXPERIMENTS:

We conduct tests in two simulations and on one real health care data set. For direct comparison, we evaluate performance against methods in the nowcasting framework using negative log likelihood and ranking measures. Then we explore WRN and WRNPPL in forecasting and multi-forecasting.

Setup. We divide the data into a train, tune, and held out test set. Model development is performed on train and tune sets with parameters determined by early stopping. Models are then evaluated on the held out test set. We use the Goodman-Kruskal γ statistic as a measure of concordance among non-tied pairs: (agree -disagree -prior)/(agree + disagree + prior), where the prior penalizes algorithms providing identical predictions. Agreement occurs when the predicted average rate and the empirical average rate are both greater in one example compared to the other, ties when either is the same, and disagreement otherwise. A crude interpretation of γ is the difference (1 -γ)/2 gives the width of a band corresponding to random guessing of correct ordering, and correctly predicting all other pairs. Further details of parameter setting, e.g. preprocessing, bin widths, optimizer settings, are in the Appendix. The code is written in PyTorch 0.4.1 and will be released upon publication.

Comparison methods. We compare wavelet reconstruction networks (WRNs) with homogeneous Poisson processes, time-invariant and nowcasting Fourier basis functions, multivariate Hawkes processes, and two long short-term memory (LSTM) networks. Briefly, the Fourier methods are given by f (t|t 0 , x) = k l w kl sin((2πl/τ )(t -t 0 )) + v kl cos((2πl/τ )(t -t 0 ))x k , where x corresponds to features, i.e., event by value interval occurrences, at or before time 0 (t 0 = 0, time-invariant) and at t 0 = t (nowcasting), k indexes the event by value interval by time step, and l indexes the basis function component. L2 regularization are applied to w kl and v kl . Then the rate function is defined as λ(•) = w 0 + f (•) 2 . For nowcasting, the Fourier method is given features from 16 previous time steps, and since in nowcasting t -t 0 = 0, the formula reduces to a generalized linear model of v kl x k terms. The multivariate Hawkes process we use includes a kernel with event-specific exponential decay parameter γ e > 0: i.e, g e (t -t i ) = e -γe(t-ti) . We use a learnable constant baseline rate λ 0 . We learn β e without constraint, rather than β e ≥ 0 or β e ≤ 0 of Hawkes and inhibitory processes respectively. We apply a positivity constraint to ensure the rate is non-negative.

The first LSTM method is a variant of the multi-task healthcare LSTM from Lipton et al. (2016) where the preprocessing involves zero-or last-value carry forward-imputation, mean-reducing, and adding missing indicators. Because our task is nowcasting not multi-label classification, we modify the loss function accordingly. The second LSTM is a WRN preprocessing LSTM analogue. The LSTM includes a linear-embedded input (i × h), two LSTM hidden layers (h × h), and output to a rectified linear layer (h × 1) where h is the hidden unit width. For each model, the output is a hazard per time step λ ny , and the loss is the point process log likelihood in Equation 1. Simulations. The first simulation is of heart attack diagnoses denoted by acute coronary syndrome (ACS). In this simulation, it is the elevation in value of troponin, a heart enzyme measurement, outside the normal range (less than 0.01 ng/mL) that indicates ACS will occur in the next time unit uniformly at random. Figure 4 (left, middle) illustrates two trajectories; note both the mark and the timing are important for ACS determination. The second simulation is of diabetes care: patients with diabetes undergo semi-regular appointments, e.g., annual eye and foot exams, quarterly hemoglobin A1c measurements, and pre-and post-prandial glucose measurements. These patients are often non-adherent with worsening adherence as a function of increasing time from adverse events. Figure 4 (right) illustrates the timings of an example trajectory.

Diabetes visits. We partnered with a regional health care system to investigate the risk of adverse outcomes of diabetes and adherence to the care those patients received. From the regional cohort followed from 2010 to 2017, we selected those at risk of diabetes as defined by an outpatient measurement of hemoglobin A1c or glucose, or a diagnosis of hyperglycemia. Among those, we excluded any individuals without at least two clinic encounters more than six months apart. We additionally applied a censor date at the time of the last clinical event before a 30-month gap in care, where there is uncertainty that the patient is lost to follow-up or is receiving care outside of network.

Application of the inclusion and exclusion criteria resulted in 798,818 timestamped events in a study population of 4,732 individuals each representing a single example. We divided the population into thirds: {train, tune, test} sets. We focused on two outcomes: (1) hemoglobin A1c measurements, as a proxy for scheduled diabetes care, and (2) a combined outcome of {ketoacidosis, neuropathy, retinopa-thy} as defined by ICD 9 and ICD 10 codes. Features included were extracted with string matching on event descriptions of events documented as putative risk factors in clinical guidelines from the ADA, AHA, and UpToDate, and included events from demographics, medications, encounters, laboratory, diagnosis, and procedures tables. The extraction resulted in 575 features. Hemoglobin A1c was measured at least once in 820 individuals (21%), and an adverse event occurred at least once in 137 individuals (3%). Additional details are given in the Appendix.

5. RESULTS:

Nowcasting. method outperformed the WRN-PPL method at the KNR task, which could be due to the relatively low rate of target events (0.025 per year) where complexity may lead to overfitting; however the difference in negative log likelihood is small. The WRN and WRN-PPL γ statistics showed effective risk stratification of individuals with diabetes. A γ of 0.98 corresponds to a band of 0.01 under the crude interpretation and suggests strong overall risk stratification.

At the other end, the time-invariant method performed marginally better than homogeneous Poisson. The nowcast method performance shows the utility of incorporating time-varying information but substantially underperformed compared to the other methods. The Hawkes process also lacked performance with comparable NLL only in the A1c* data set. LSTM 1 had good ranking characteristics on the simulations but mediocre predictive performance with notably poor performance on the simulation A1c*, where the timing of events, not their value distribution, determines the ground truth rates. LSTM 2 mostly outperformed LSTM 1 , suggesting the usefulness of WRN preprocessing, but mostly underperformed against WRN-PPL.

Figure 5 shows the hemoglobin A1c predicted hazards profile for random test set patients using WRN-PPL. The step function represents predicted hazard over time, points indicate true event times and dotted lines show the difference from the baseline (homogeneous Poisson) rate. The WRN-PPL algorithm makes predictions that anticipate appointments where hemoglobin A1c will be measured in quasi-periodic fashion (right). Similarly, Figure 6 (left) illustrates the ability to model the rates of complications. Medical guidelines do not specify scheduling for regular follow-up of the adverse events, and this is congruent with the lack of periodicity in the KNR hazard predictions.

Forecasting, multi-forecasting, and interrogability. As one would expect, a trade-off occurs between early prediction and predictive performance. The effect of WRN-PPL forecast distance c on KNR prediction is shown in Figure 6 (right). Notably, the 3-month censored WRN-PPL has approximately the same performance as the nowcasting LSTM 2 . Similarly, effects of single-model, multiple-c prediction are shown in Figure 7 (left), illustrating the WRN-PPL improvement over WRN for nowcasting (up to c < 1) but not for c ≥ 1. The coefficient profile as a function of c in multiple-c prediction is also shown (middle left), demonstrating that relative-time attributions, which are commonly used in association statements in health literature, appear to depend on censor time c. 

6. DISCUSSION:

The performance of WRN-PPL in Table 1 and Figure 5 illustrates the utility of our model, in particular in identifying the near-periodicity of recurring events. For example, the rate prediction for the individual denoted in green in Figure 5 (right) suggests that individual may have skipped, missed, or rescheduled 5 to 6 appointments over the last decade. The peaks reach hazards of approximately 3, indicative of a mixture of belief and uncertainty-belief that in those months the event should occur at a rate above three per year, and uncertainty about the occurrence of the appointment.

For multi-forecast learning, a comparison of the results in Table 1 and Figure 7 (left) demonstrates the value of model expressivity. In particular, Table 1 shows that single forecasting outperforms multi-forecasting at c = 0 in Figure 7 (right). However, Figure 7 (middle right and right) illustrates that multi-forecasting improves the learned wavelet representation. These findings suggest that the layering between the wavelet reconstruction (WRN: {reduction layer, linear}, and WRN-PPL: {reduction, permute and pool, linear}) and the hazard output is not adequately expressive to map the true wavelet reconstruction to the true hazard. We argue the solution is not in simplification nor abandonment of the multi-forecast setting, but in leveraging the multi-forecast setting to facilitate recovery of the wavelet reconstruction by using an even more expressive mapping.

Conclusion.:

Wavelet reconstruction networks is a forecasting method tailored for health care settings. Its advantages include multi-resolution representation of relative time dependencies in 1 and 2 dimensions and enables time step specification at test time. The combinations of relative time functions result in a model with powerful predictive capacity, flexibility, and interrogability. We demonstrated improved performance of our system over competing methods in a analysis of diabetes and showed the ability to capture quasi-periodic events that could be used to measure adherence and forecast risk of complications. APPENDIX C LOSS-BY-EPOCH CURVES Epoch curves were used for early stopping-the parameters used on the held out test set were those for which the tune set negative log likelihood was minimized, only for the best performing hyperparameter setting (identified through tune set performance). 

:

-1 2 -1 2 -1/2 2 -1/2

Figures and Tables:

Figure 1: Illustrative 1-d cross-correlation motivating the discrete wavelet reconstruction kernel for relative time dependencies.

Figure 2: Wavelet reconstruction network architecture

Figure 3: Left: three absolute time hazards (green, blue, and purple) for one trajectory with different steps sizes and censor times. Right: forecasting tasks. The vertical distance from the line t = t is censor distance c. Common loss definitions are along solid lines.

Figure 4: Examples: acute coronary syndrome-angina (left), ACS (middle); diabetes simulation (right).

Figure7(middle right and right) shows the wavelet reconstruction for the effect of troponin level and timing on rate. Both reconstructions demonstrate recovery that acute coronary syndrome is diagnosed within the next time unit after a troponin greater than 0.01 ng/mL. The multiple-c reconstruction on the right more accurately reflects the uniform distribution hazard, namely, increasing hazard if the event has not yet occurred.

Figure 6: Left: ketoacidosis, retinopathy, or polyneuropathy diagnoses and rate predictions for thirty random test set patients. Right: combined outcome (KNR) negative log likelihood as a function of censor distance (in years).

Figure 7: Best viewed in color. Left: negative log likelihood as a function of forecast censoring distance c for multi-forecasts. The permute and pool layer expresses greater expressivity to model the hazard than WRN. Middle left: coefficient profile as a function of forecast censoring distance c for WRN. The reduction layer comprises additions and maximums of 1-d and row-indexed 2-d reconstructions. Middle right: wavelet reconstruction of troponin contribution to ACS hazard. A preceding troponin above 0.01 indicates increased rate of ACS occurrence within the next hour. Right: WRN-PPL reconstruction image for multi-forecasting at c = {0, 0.25, 0.5, 0.75, 1}.

Figure 8: Hemoglobin A1c hazard prediction with nowcasting

Figure 10: Loss by epoch curves for real data, target: {ketoacidosis, polyneuropathy, retinopathy} used for early stopping (left) and target: hemoglobin A1c (right).

Negative log likelihood on the held out test set. Asterisk (*) denotes simulation; KNR: {ketoacidosis, neuropathy, retinopathy}; H. Poission: Poisson; Time-invariant: Fourier prediction at time t = 0; Nowcast: 16-step history as nowcast features. LSTM 1 : Lipton et al. (2016), best of zero-imputation or lastvalue carry forward, with missingness indicators; LSTM 2 : LSTM with WRN preprocessing; WRN: wavelet reconstruction network; PPL: with permute pool layer. Best performer in bold.

Table content: <table><tr><td>Dataset</td><td/><td/><td/><td colspan="2">Method (NLL)</td><td/><td/><td/></tr><tr><td/><td colspan="8">H. Poisson Time-invarant Nowcast Hawkes LSTM 1 LSTM 2 WRN WRN-PPL</td></tr><tr><td>ACS*</td><td>0.44</td><td>0.43</td><td>0.36</td><td>0.39</td><td>0.21</td><td>0.13</td><td>0.23</td><td>0.15</td></tr><tr><td>A1c*</td><td>18.54</td><td>19.20</td><td>13.56</td><td>3.87</td><td>11.80</td><td>4.10</td><td>3.93</td><td>3.78</td></tr><tr><td>A1c</td><td>2.86</td><td>2.76</td><td>2.52</td><td>1.67</td><td>1.15</td><td>1.29</td><td>1.23</td><td>1.13</td></tr><tr><td>KNR</td><td>0.75</td><td>0.71</td><td>0.58</td><td>0.31</td><td>0.46</td><td>0.35</td><td>0.24</td><td>0.26</td></tr></table>

Goodman-Kruskal γ, a measure of concordance, on held out test set. Asterisk (*) denotes simulation.

Table content: <table><tr><td>Dataset</td><td/><td/><td colspan="3">Method (Goodman-Kruskal γ)</td><td/><td/><td/></tr><tr><td/><td colspan="8">H. Poisson Time-invarant Nowcast Hawkes LSTM 1 LSTM 2 WRN WRN-PPL</td></tr><tr><td>ACS*</td><td>-1.00</td><td>-0.81</td><td>0.08</td><td>-0.97</td><td>0.98</td><td>0.91</td><td>0.80</td><td>0.85</td></tr><tr><td>A1c*</td><td>-1.00</td><td>-0.02</td><td>0.81</td><td>0.64</td><td>0.84</td><td>0.85</td><td>0.78</td><td>0.77</td></tr><tr><td>A1c</td><td>-1.00</td><td>0.23</td><td>0.71</td><td>0.79</td><td>0.83</td><td>0.87</td><td>0.93</td><td>0.93</td></tr><tr><td>KNR</td><td>-1.00</td><td>0.25</td><td>0.81</td><td>0.91</td><td>0.84</td><td>0.91</td><td>0.95</td><td>0.98</td></tr></table>

Table 1 reports the negative log likelihoods for the experiments on the held out test sets.Overall, the proposed wavelet reconstruction network WRN-PPL outperformed the other algorithms. The WRN-PPL method excelled particularly in tasks with many target occurrences (A1c* and A1c experiments) and performed near to the best in rare occurrence data (ACS* and KNR). The WRN Figure5: Left: hazards and hemoglobin A1c events for A1C* simulation. Right: hazards and hemoglobin A1c events for five random, test set patients in real cohort.

Table content: <table/>

Peer review (iclr2019):
 [PROS] [originality] The paper proposed to construct the intensity function of a point process using wavelet, in order to improve its expressiveness, e.g. allowing non-additivity. The authors did extensive experiments to investigate their model performance compared to many appropriate baselines, on both synthetic and real-world datasets. [CONS] [clarity] The major drawback of this submission is its clarity. The paper is vague at various important points in both method and experiments, thus leaving their correctness and soundness undetermined. In the method section, the authors did not specify a well-defined and self-consistent notation system. This makes the paper really hard to understand. For example, one may be easily confused with things like: 1) How q, q(g(t), t_i), g(t), g_{es}, g_{d} are connected and distinguished? 2) The function g_{es} maps from t to R, then is t a space or a variable? 3) The state s seems crucial in the function g_{es}, but why it is only mentioned one time in the paper? How it is defined and how it is used? 4) How the Hadamard product is applied to two matrices of different sizes? 5) j=1 is time dimension and j=2 is time and value dimension, then why j=1 is not part of j=2? Time is needed in both cases and it seems natural that the associated parameters are shared. 6) Figure-3 has t_i in figure, t’ in caption, but the text in main paper mentions t_0 for Figure-3. How are they related? Are they actually the same? The most confusing part is the censoring distance c. Its introduction around eqn-(2) suggests that c > 0 and the censoring section clearly mentions that. But c is also used to denote the forecast distance, which is clearly < 0 according to Figure-6 and Figure-7. What’s worse, there is also a term called forecast censoring distance. What are the relationships among these terms? If they are all the same, then is the c actually a model parameter or an evaluation control knob? Such things are very important to clarify. Moreover, the paper did not clearly explain how the model is trained in each case, especially for (multi-)forecasting. In details: 1) What is the training objective?. 2) What is the optimization method for this objective? 3) How is it implemented and would the code be released? It is good that the experimental section lists many appropriate baseline models and multiple evaluation metrics, but it is not clear how they are used. For example: 1) Fourier methods and Hawkes process do not deal with the value v, then how are they fairly compared to the proposed model which takes v into account as in eqn-(2)? 2) How is the Goodman-Kruskal gamma exactly computed? On all the instances of the held-out set? What is exactly the rank in this case? 3) The authors also leave out the positiveness constraints of a Hawkes process to incorporate inhibitory interactions, but how the positivity of the intensity function is ensured in this case? [quality and significance] The method is well-motivated and innovative. But details of the model and experiments are very unclear, so its overall soundness is hard to judge. For example: 1) The authors claim that, compared to neural models, their model has the advantage of the interpretability (for small datasets), but they also have neural components in their model. So why their model is more interpretable than others (e.g. Mei and Eisner 2017 as they cited) is not clear to me. 2) It is not clear why the interpretability is associated with the size of the dataset (quote `remains interpretable for small data sets’). What’s worse, the interpretability seems the only advantage of the model over other neural models (please correct me if I am wrong). If this edge could not scale up to large datasets, then does it mean on large datasets, a neural model should always be preferred over this model, because they are supposed to be more expressive?