Stabilized Neural Prediction of Potential Outcomes in Continuous Time
Authors: Konstantin Hess
Stefan Feuerriegel
Venue: ICLR 2025 Poster (2025)
Keywords: causal inference
potential outcomes
continuous time
treatment effects
inverse propensity weighting
confounding
medicine
Primary area: causal reasoning
TL;DR: We predict potential outcomes with a novel stabilized inverse propensity weighting network in order to adjust for time-varying confounders in continuous time.

Abstract: Patient trajectories from electronic health records are widely used to estimate conditional average potential outcomes (CAPOs) of treatments over time, which then allows to personalize care. Yet, existing neural methods for this purpose have a key limitation: while some adjust for time-varying confounding, these methods assume that the time series are recorded in discrete time. In other words, they are constrained to settings where measurements and treatments are conducted at fixed time steps, even though this is unrealistic in medical practice. In this work, we aim to estimate CAPOs in continuous time. The latter is of direct practical relevance because it allows for modeling patient trajectories where measurements and treatments take place at arbitrary, irregular timestamps. We thus propose a new method called stabilized continuous time inverse propensity network (SCIP-Net). For this, we further derive stabilized inverse propensity weights for robust estimation of the CAPOs. To the best of our knowledge, our SCIP-Net is the first neural method that performs proper adjustments for time-varying confounding in continuous time.

Published as a conference paper at ICLR 2025
STABILIZED NEURAL PREDICTION OF POTENTIAL OUTCOMES IN CONTINUOUS TIME
Konstantin Hess & Stefan Feuerriegel
Munich Center for Machine Learning LMU Munich (k. hess, feuerriege mum de
ABSTRACT
Patient trajectories from electronic health records are widely used estimate con- ditional average potential outcomes (CAPOs) of treatments over time, which then allows personalize care. Yet, existing neural methods for this purpose have a key limitation: while some adjust for time-varying confounding, these methods assume that the time series are recorded discrete time. other words, they are constrained to settings where measurements and treatments are conducted at fixed time steps, even though this is unrcalistic in medical practice. this work, we aim to estimate CAPOs in continuous time. The latter is of direct practical rele- vance because it allows for modeling patient trajectories where measurements and treatments take place at arbitrary, irregular timestamps. We thus propose new method called stabilized continuous time inverse propensity network (SCIP-Net). For this, we further derive stabilized inverse propensity weights for robust esti- mation of the CAPOs. To the best of our knowledge, our SCIP-Net the first neural method that performs proper adjustments for time-varying confounding continuous time.
INTRODUCTION
Estimating conditional average potential outcomes (CAPOs) of treatments is crucial personalize treatment decisions in medicine (Feuerriegel et al. 2024). Such CAPOs are increasingly predicted based on patient data from electronic health records (Bica et al. 202T}. This thus requires methods that can model the time dimension in patient trajectories and, therefore, estimate CAPOs over time.
Existing neural methods for estimating CAPOs over time primarily model the patient trajectory in discrete time (e.g.. [Bia et al.) 2020 Li al. 2021 Lim et al. 2018 Melnychuk et al. 2022). As such, these methods make wnrealistic assumptions that both health measurements and treatments oc- a fixed, regular schedule (such as, e.g., daily hourly). However, both health measurements and treatments take place at arbitrary, irregular timestamps based on patient needs. For example, patients in critical state may be subject to closer monitoring, that measurements recorded more frequently (Allam et al. 2021).
To account for arbitrary, irregular timestamps of both health measurements and treatments, methods are needed that correctly model the patient trajectory in continuous time (Lok] 2008 Røysland 2011 Rytgaard al. 2022). However, neural methods that operate in continuous time scarce (see Sec. 2. Crucially, existing neural methods have a key limitation in that they fail to properly account for time-varying confounding (e.g.. [Seedat al., (002)). This means that, for a sequence of future treatments, the corresponding confounders lie also in the future, are thus unobserved, and therefore need to be adjusted for. Yet, existing neural methods rely only on heuristies such as balancing, which targets improper estimand and thus leads to estimates that biased. To the best of our knowledge, there no neural model that estimates CAPOs in continuous time while properly adjusting for time-narying confounding.
In this paper, we aim to estimate CAPOs for sequences of treatments in continuous time while prop- erly adjusting for time-varying confounding. However, this is non-trivial challenge, this requires a method that can perform adjustments at arbitrary timestamps. While there are methods to adjust for varying confounding in discrete time, similar methods for continuous time are still lacking.

---

Published as a conference paper at ICLR 2025
Correct timestamps?
Adjustment time-varying confounding?
Existing
X
x
CRN Biga el al. 2020 CT Molnyclokk at 2022
Neural mothods
Limn
Lr 2001
Neural methods continuous time
x
TE-CDE Seedat al. 2022
SCIP-Net (ours)
Table I:
Comparison of key neural methods for estimating CAPOs over time. Our SCIP-Net the first method to perform proper adjustments for time-varying confounding in continuous time.
Therefore, we first derive a tractable expression for inverse propensity weighting (IPW) in continu- ous time. However, a direct application of IPW may suffer from severe overlap violations and thus lead to extreme weights. As a remedy, we further derive stabilized IPW in continuous time. We then use our stabilized IPW to propose novel method, which we call stabilized continuous time inverse propensity network (SCIP-Net). Unlike existing methods, ours is the first neural method to estimate CAPOs in continuous time while properly adjusting for time-varying confounding.
We make the following contributions D (1) We introduce SCIP-Net, novel neural method for esti- mating conditional average potential outcomes in continuous time. (2) We derive a tractable version of IPW in continuous time, which provides the theoretical foundation of our paper for proper adjust- ments for time-varying confounding. Further, we propose stabilized IPW in continuous time, which then in our SCIP-Net. (3) We demonstrate through extensive experiments that our SCIP-Net outperforms existing neural methods.
2 RELATED WORK
Table presents an overview of key neural methods for estimating CAPOs over time. An extended related work in Supp. A
Average vs. individualized estimation: Estimating average potential outcomes over time is a well-
studied problem in classical statisties (e.g.. Bang & Robins 2005 Lok, 2008 Robins 1986 1999
Robins& & Hernán 2009] Røysland 2011: Rytgaard et al. 2022; 2023 van der Laan & Gruber
2012). However, these methods are population-level approaches and thus do not make individu- alized estimates at the patient level. Put simply, the observed history of an individual patient is ignored. Therefore, they are not suitable for personalized medicine. In contrast, our work (and the following overview) focuses on potential outcome estimation conditional on the observed patient history, which thus allows us to make individual-level estimates for personalized medicine.
① Neural methods in discrete time: Some neural methods for estimating CAPOs over time impose discrete time model on the data (e.g.,Bica et al. 2020 Li et al. 202T: Lim et al. 2018 Melnychuk] al. 2022). As such, these methods operate under the assumption of both fixed observation and treatment schedules, yet which unrealistic in clinical settings. Instead, patient health is typically monitored at arbitrary, irregular timestamps, and the timing of treatments also takes place at arbi- trary, irregular timestamps, which may directly depend on the health condition of a patient. Hence, methods discrete time rely on data model that is not flexible enough account for arbitrary, irregular monitoring and treatment times, because of which their suitability in medical practice is limited.
Only few neural methods have been developed for es- timating CAPOs continuous time. Yet, existing methods have key limitations. One stream of methods (Hess et al. 2024b; Vanderschueren et al. 2023) ignores time-varying confounding and is thus not applicable to our setting.
2 Neural methods in continuous time:
To the best of our knowledge, there is only one neural method that works in continuous time and that applicable our setting: TE-CDE (Seedat et al., (002). This method tries handle time-varying confounding through balancing. However, balancing fact, balancing was originally proposed for variance reduction (Johansson) al. 2016) and may even increase bias (Melnychuk et al.) 2024). Therefore, TE-CDE suffers from an infinite data bias that comes from the fact that is does not properly adjust for time-varying confounding.
heuristic approach to adjust for time-
varying confounding:
Cooee available at httpaa/:///ggttbb.ccoo/////////t
2

---

Published as a conference paper at ICLR 2025
Research gap: To the best of our knowledge, none of the above neural methods performs proper adjustments for vareingying confounding in continuous time. As a remedy, we propose SCIP-Net, which is the first neural method that estimates CAPOs in continuous time while properly adjusting for time-varying confounding.
PROBLEM FORMULATION
Let [0, 7] be the time window. the following, we assume every stochastic process Vi2 (tt) defined on [0, to
Notation:
bed Intevvntion
be cådlág, and we let V,- = lim V, denote the left time limit. Further, we write V {V.)agr and Vr, {Va)a2r.
Treatmat
Seatmert etechs
Setup (see Fig. 1]: We consider outcomes Y, E Rd,, discrete treatments A, E (0, 1)d1, and covariates X₁ E RRd,, where we assume that X, contains Y, for all < t. Without loss of generality, we assume that static covariates are included in X,, At time t, we let (Ae) and (X,e denote the ob- servational treatment propensity and the covariate distribution, respectively. Both measurement times and treatment times typically dynamic and do not follow fixed schedules. Rather, both are recorded at arbitrary, irregular timestamps.
ouzome
Figure Setup. Shown are treatment and covariate trajecto- ries in continuous time. OB- servational treatment assignments are confounded by covariates and, hence, estimating CAPOs requires adjustments.
To formalize the above in continuous time, we need to be able to model arbitrary timestamps, which increases the complexity compared to the discrete time setting considerably. Following Ryyggard et al. (2022; 2023), we let the counting processes N₂ and N₁ govern the times at which covariates are measured and which treatments may be administered, respectively. For both E {a, 2). let - Tx;) } denote the set of jumping times of the process Nē with T%; S т. Further, we let Aÿ denote the cumulative intensity of Nj, and let (t) denote the corresponding intensity function. Further, we use the short-hand notation dt = [t, + df) as in (Gill & Johansen 1990. Rytgaard et al. 2022). Then, have that
Xt( dt (f) Ez, [N' (d)) Po (N (dt) 1).
Observational likelihood: We write the full observed history up to time as
H. {Ayy, Xry, Ny, NE,NE TỤ T E T.
Then, following (Rytgaard al. 2022), we can write the observational likelihood of the data Hr via
(Fo(H))) (Xo) π
(dAA(s Hx-))o,,(X,, H)- (=(ss) (1-dAī(s | 2- NE(D,T]
HHa)))
(dAg(s) H,,) H)-)) Nd(ss) ( dAg(s H.-)) I- H
where π is the geometric product integral. Intuitively, the geometric product integral π can be thought of as the infinitesimal limit of the discrete product operator II. Importantly, the geometric product integral Equation 3 is the natural way to describe joint likelihoods continuous time. For more details, we provide a brief overview of product integration in Supp. B
Objective: We are interested in estimating the response of the outcome variable Y, when intervening on the treatment sequence starting time given an observed history H;, h, For this, adopt the potential outcomes framework (Neyman 1923 Rubin 1978). That is, we seek to estimate the CAPO
E Lo,,, [.,] H. Tr.
under interventions on both the treatment propensity 4, = a,,t and the treatment frequency Nº = n= given the history H, h,_. Here, the interventions are hard interventions, which standard in the literature (e.g., Bica et al. 2020 Lim al.) 2018 Melnychuk et al.) [2022] Seedat
et al. 2022). That is, for 2 I, we are interested in deterministic interventions of the form
A, (A, H h 1 +A.=a,.)
No dsy(s H, h ) = I Le..)_. (8) ds, (5)
3

---

Published as a conference paper at ICLR 2025
where 0x,s RR" + is step-wise constant function with jumping points {+h
Estimating the CAPO for treatment sequence is notoriously challenging due to the fundamental problem of causal inference (Imbens & Rubin 2015). This means, only factual outcomes Y, in are observed in the data, but not the potential outcomes when intervening on the treatment. In the fol- lowing, first define the interventional distribution for objective and then ensure identifiability.
Interventional distribution: We now define the interventional distribution for our objective. For this, rewrite Equation by reweighting Y, under the observational distribution dlF For this, follow Rytgaard et al. 2023) and first split the likelihood into two separate parts as
π =E(0,7]
dlo (+)) = dlp Qo (H+) = (X))(Xo)
(H.)
(6)
where
dGo,, (dAg(s) To,) H,-)) (1ss)) dAgss H₂-)) (ds)
the treatment part the intervene and where
(ds) dQo.: (dA7 HX,) H,-)) ( (1-dA Ha-)) 1-
(8)
remains unchanged. Then, we can write the interventional distribution as
(H.) π
(H,, dl' Q..G. (Hlr) = (o,o(Xoo
xE(0,7)
where
(H.) = (dA; H.-) o H1 {) N*
dA-s HI.))) 1-
H)
(10)
Identifiability: To ensure identifiability, we need to make the following assumptions (Lok 2008
Robins & Hernán 2009: Rytgaard al. 2022) that are standard in the literature for estimat- ing CAPOs over time (e.g.. SSedat et al. 2022).. (i) Consistency: Given an intervention on the treatment propensity and the frequency [e. p.,,,, the observed outcome Y, coincides with the potential outcome Y, [2.. under this intervention. (ii) Positivity: Given any history H,_, the Radon-Nikodym derivative dGG., /.G.e/dGo, exists. (iii) Unconfoumdednees:: Given any his- tory H._, the potential outcome is independent of the treatment assignment probability, that is, Y- Y- e: n (A,, No) H...
Proposition 1. Under assumptions (i)-iii), we can estimate the CAPO from observational data from data sampled under dP,)) via inverse propensity weighting, that is,
π W. (A. (a.,;, n²,) H. = Th- (11)
E Y- [V n
H:- = To
E Y-
where the inverse propensity weights for s 2 are defined as
dG.,s (Ha)
W. = w. (H.)) = dGo,s (H.)
(12)
Proof. See Supp.
Proposition important for the rest of our paper: it tells us that we can estimate CAPOs in continuous time from data sampled under dPp. For this, we need to quantify the change in measure from the observational distribution dlPo to the interventional distribution dP.,, which is given by Equation I2
However, the above formulation is based on a product integral J. which is not computationally tractable. + We later derive novel, tractable formulation where we rewrite Equation using the product operator I[ (Sec. 4.1 Further, inverse propensity weights W, may lead to extreme weights and, hence, unstable performance. This is known issue because settings over time prone to low overlap (Frauen et al. 20255 [im et al.) 2018). + We later derive novel stabilized weights that are tailored
continuous time setting (Sec.
4

---

Published as a conference paper at ICLR 2025
4 SCIP-NET
In this section, we introduce our SCIP-Net. It is designed to perform proper adjustments for time- varying confounding in continuous time.
Objective: Our objective is to find the optimal parameter 2 Φ of neural network ms via
x arg min Ee.
(V-le... 20, mo (A., No. ) He- = he-
(13)
-arg min Ero (y=+ ma (A, Ar,
I W. (A,, e. Hr- = he- (14)
Note that the above objective makes use of inverse propensity weights W.. However, these weights suffer from two drawbacks: they are (1) intractable as they rely on product integrals Л[. and (2) they may lead to unstable performance. As remedy, we (1) derive tractable expression for this product integral (Sec. [.11. and we further (2) introduce stabilized weights (Sec. 4.2). Finally, we present our neural architecture (Sec. 4.3) and how perform inference (Sec.
4.1
REWRITING THE OBJECTIVE FOR COMPUTATIONAL TRACTABILITY
We now derive a tractable expression to compute our unstabilized weights W, in Equation 12 Proposition 2. Let t2 = t for notational convenience. The unstabilized weights in Equation 12 scatssyy
((A,, = (E,,, Hr- Hh π w. (IA,. Ni) (0., Hi. hr-)= We, 22
he-).
==1
(15)
where
exp SSE 122,,) 28(s HT cl
We,
(16)
A0(12 Her To.t. (an.e Hee
Proof. See Supp.
Proposition 2 has important implications for tractability: Given history H, h4- and the future sequence treatments (A,, No) (a.,t: .. ) the product integral π of the unstabilized inverse propensity weights W, reduces to finite product I[. This product includes both treatment propen- sities and treatment intensities. In Sec. 4.3 we show how to learn these quantities from data.
Importantly, the unstabilized weights in Equation 16 already sufficient to adjust for time-varying confounding, However, as they may lead to extreme weights, now propose stabilized weights.
4.2 STABILIZED WEIGHTS
In the following, we first define our stabilized weights (Def. DT. Then, we show that the optimal parameters in Equation 14 are the same, regardless of whether the original, unstabilized inverse propensity weights from above used or our stabilized weights (Proposition Finally, present a tractable expression to compute the stabilized weights (Proposition 4
Definition 1. For t. let the scaling factor E, be given by the ratio of the marginal transition probabilities of treatment, that is,
E. (A No) = dGo,- No)
(17)
dG a (Ag, Na)
We define the stabilized weights W. as
W. =EWW..
(18)
Throughout, refer to the weights of neural nets parameters to make the distinction to inverse propen- sity weights clear.
5

---

Published as a conference paper at ICLR 2025
The idea of stabilized weights (e.g., Lim et al. 2018) is that the marginal transition probabilities E., on average over the population, downscale the inverse propensity weights. Importantly, the scaling factors are not conditioned on the individual history and, therefore, do not change the objective.
To formalize this, we make use of the fact that the optimal parameters in Equation 14 are invariant to multiplicative scaling of the optimization problem with respect to constant scaling factors. We summarize this in the following proposition.
Proposition 3.
The optimal parameters & in Equation 14 can equivalently be obtained by
é = arg min B2. (y- mo . Ni. H H.-))
W (A,
(a. (1. 2i,,, H hi-
(19)
Proof. See Supp. D
Proposition guarantees that we can substitute the original, unstabilized weights W. with the stabi- lized version W As the scaling factors Eg downscale W., reduce the risk of receiving extreme inverse propensity weights, and, thus, obtain a more stable objective.
The results from Proposition still rely product integrals Л[ However, in Proposition 2 now derive an equivalent, tractable expression that relies on the product operator II instead.
Proposition Let t²20 = t for notational convenience. The scaling factor E, from Equation then satisfies
π 82t
E.
(A,, (a.,:: (A,_, N; (ae- nT-))
(20)
II j=1
Beg ((d,. (a (A-, (a._,
(21)
where
A Aee, -, No k0,ee, (are Are No
(22)
exp S A8(s Ao_ No ds
Proof. See Supp. D
Together, Propositions and yield a (1) tractable objective function that relies on (2) stabilized inverse propensity weights. As show in the following Sec. 4.3] we can estimate the stabilized weights from data. Thereby, we present our SCIP-Net, which adjusts for time-varying confounding continuous time.
4.3 NEURAL ARCHITECTURE
Overview: We now introduce the neural architecture of our SCIP-Net, which consists four com- ponents (see Fig.
2): The S stabilization network learns estimator of the scaling factors (.(-)
from Equation 17 The (w) weight network learns an estimator of the unstabilized inverse propensity
weights wx(-) from Equation 12 Combining both, we have an estimator for the stabilized weights
w. (-) from Equation 18 The E encoder learns representation of the observed history, which is then passed to the decoder. Finally, the D decoder takes the learned representations and the stabi-
lized inverse propensity weights as input to estimate the CAPOs Y, [a.,, n H,- = h,-
and D use neural controlled differential equations (CDEs) backbones. Neural CDEs have several benefits for our SCIP-Net. First, neural CDEs process data in continuous time. Second, neural CDEs update their hidden states data becomes available over time. We provide brief introduction in Supp.
Backbones: All components S W (Kidger al. 2020 Morril et 2021)

---

Published as a conference paper at ICLR 2025
Stability network
A,, A
Meuri GIDE
A(t.J)
Decoder
f2
*(Auo)
Xry N
CDE
Weight network
W
Wir
aCE
{(trg)
f
CDE
W
Mg
#(Aag)
P
Legend
CC
Cross entropy
Encoder
BCE
Binery GE
M
Y
MSE
Meen squared emrrr
WWSE:
Welghted MSE
Figure 2: Neural architecture of our SCIP-Net.
In the following, we denote the training samples by {he, 1--1 and the test samples by ((hg,- Ba 2e )= Reassuringly, we emphasize that hy - are realizations from the observational distribution dlPo, whereas (2,, is the interventional treatment sequence.
The stability network learns an estimator Ea(-) for the scaling factors E.g from Proposition consists of a linear input layer vy. neural vector field fo. and two linear output layers 1 and s.s which estimate the treatment intensity and propensity, respectively.
S Stability network:
Training: The stability network receives treatment decisions and treatment times Ng) = (āirr no We distinguish two cases. For 0 < t < the latent representation zi, evolves as
= 2000 1 fo (21.. 8) dia,..-. nf,.-
(23)
where dlai.s-,n,, n9,,- denotes Riemann -Stilltees integration with respect to the control path and <0 = v (ar,0, n) Case (1): The latent representation is then passed to the intensity layer M in order to estimate the probability whether a treatment decision is made at time L. For this, our SCIP-Net optimizes the binary cross entropy (BCE) loss
c 5.1 (0) = BCE
(1 S.1 (=..). dn
(24)
Thereby, the intensity layer p learns an estimator the treatment intensity function via
Ag(t A,_,N; 1 S.1 (Z2))
(25)
Case (2): If treatment decision is made time T = f the latent representation is additionally passed through the propensity output layer i to estimate which treatment is administered. For this, our SCIP-Net minimizes the cross entropy (CE) loss
5 C (0)
CE 1
(2 ( ), (.a..r.))
(26)
Hence, the propensity layer learns an estimator of the propensity score via
(ot Ng u S.P (2,).
(27)
Scaling factor: After training, the stability network again receives the training samples (Ar, No) = (āi.., n. Following Proposition it then computes
π ((rT,,s
(28)
At Aeg,- a..f, Nga No.,- ni,1e -0,,__1 )) No.r9, (arg, Arg,- = 0r.09 No. nleg, II
exp Joe Ag(s As- ai,o-, No n...-) ds

---

Published as a conference paper at ICLR 2025
where
arbitrary quadrature scheme to compute the integral in the denominator.
W Weight network: The weight network learns an estimator w.)() for the unstabilized weights It also consists linear input layer vo a neural vector field fø and
We from Proposition 2
two linear output layers M48 propensity, respectively.
W.P which are trained to estimate the treatment intensity and W.J and p
Training: The weight network receives samples H, = hirr- After 14,0 is transformed into =W via W. the latent state zi,t of the weight network evolves &
3,0
2,,0 + ( [ ()
d
(29)
where d[hs,e-] denotes Riemann- Stinlttess integration w.r.t. the control path. Case (1): As for the stability network, the weight network estimates the probability of treatment decision at time t through the intensity layer Ha
W.J via
C.W. BCE WiJ (= (
(30)
and, hence, learns estimator of the treatment intensity function via
Àg H1-) - 1 T.J (Z.").
(31)
Case (2): If treatment decision is made at time Ty = 12g the latent state is also passed through the propensity layer H and jointly trained via
C
W.P = CE
W.P J, CE
12 (
(32)
such that our SCIP-Net learns an estimator of the propensity score as
(a He-) Hb
W.P
)z
(33)
Inverse propensity weight: After training, the weight network again receives the training samples H, = hi.r and estimates the unstabilized inverse propensity weights according to Proposition via
exp SaEtnn 1155)) \i [s| H. hhsss ds Heg,- = hic, -) fo,.e (ar,r,, His,- = hvg
π u. (hi..-) II
(34)
1=1 Ag(t
The encoder computes a latent representation of the history, which is then passed to the decoder. It consists of a linear input layer E g. a neural vector field f and an output layer u
E Encoder:
Training: The encoder receives samples H, h.. First, it transforms h4o0 into 2 via ug. The latent state 2,, then evolves according to
= 2.. + [ f (= sh
E 2
(35)
where d[hr,s-] denotes Riemann- -temannnStiiltjes integration w.r.t. the control path. At jumping times ni,,- we pass the latent state 2ı, to the encoder output layer 2 and minimize the mean squared error (MSE) loss for outcomes at the jumping times via
c = MSE He (2 ar,)), v..e)
(36)
D Decoder: The decoder receives input: (i) the encoded history of the encoder and (ii) the future sequence of treatments. Then, it outputs an estimate of the CAPOs by adjusting for time-varying confounding. has linear input layer vD neural vector field 1 and output layer ug
8

---

Published as a conference paper at ICLR 2025
Training: During training, the decoder receives the final latent representation 2 of the encoder as well the observed treatments (A. At) (ei 21.)) It then transforms the encoder representa- tion through z5. = v (zf,)) and computes
2 = + S. s)
(37)
where diai,.-, n,,- denotes Riemann Riemann-Stieltjes integration w.r.t. the control path. At time т, we pass zfPr through the linear output layer ng. Importantly, the decoder is trained by minimizing the MSE loss weighted by the stabilized weights, i.c.,
C
MSE ue (2 ( ar.)).1..))
(38)
where
W We, (he,r,)) Ext,, (ar,2,, Tiey,) ) (bs,02,)
(39)
using the stability network (.) and the weight network (.)))), respectively. By Proposition 3 thereby target the optimal model parameters which, unlike existing methods, explicitly adjust for time-varying confounding in continuous time.
4.4 INFERENCE
order to estimate CAPOs for an observed history H,- hu,-- and a future sequence of treat- ments NT) ( ;) first encode the history via
= = =k,0 +
1 (= 8) d[hx,s-)-
(40)
This latent representation then passed to the decoder along with (a.,: n.,)) in order to compute the final representation
2,, - J 1 () d(on,,,- n
(41)
Finally, the output layer ng of the decoder is used to estimate the CAPO at time т, given the history H1- hg.e-, via
n8 H, hg.c-] M (= a..)).
(42)
NUMERICAL EXPERIMENTS
Baselines: We now demonstrate the performance of our SCIP-Net against key neural baselines for estimating CAPOs over time (see Table [ Importantly. our choice of basclines and datasets is consistent with prior literature (c.g.. Bica al. 2020 Lim et al. 2018 Melnychuk Seedat et al. (022). Further, we report the performance of the CIP-Ner Net ablation, where we directly train the decoder with the unstabilized weights. This allows us to understand the performance gain of our stabilized weights. Note that the CIP-Net ablation is still new method (as no other neural method performs proper adjustments for vire-ng confounding in continuous time).
al., 2022
Datasets: We use (i) synthetic dataset based on tumor growth model (Geng et al. 2017 and a (ii) semi-synthetic dataset based on the MIMIC-II dataset (Johnson et al. 2016). For both datasets, the outcomes are simulated, so that we have access the ground-truth potential outcomes, which allows for comparing the performance in terms of root mean squared error (RMSE). We report the mean + the standard deviation five runs with different seeds. We perform rigorous hyperparameter tuning for all baselines to ensure a fair comparison (see Supp. G
Tumor growth data: Tumor growth models are widely used for estimating CAPOs model for the evolution lung cancer Y, under radio therapy treatment Af and chemo therapy treat- ment AC (details in Supp. E. Here, we follow Vanderschueren et al. (2023), where observation
time. is

---

Published as a conference paper at ICLR 2025
times at arbitrary, irregular timestamps. However, our setup has two differences: (i) We not primarily interested in informative sampling times. Instead, we consider a scenario with observation times completely at random. Hence, observation times follow Hawkes process with constant inten- sity. (ii) Further, we do not consider an RCT setting. Instead, treatment assignments are confounded via AE; AT Ber((n/ Dmax Dis (Y maxx /2)x where ? controls the confounding strength, and Dmxx and D15 are the maximum and average tumor diameter over the last 15 days, respectively. For evaluation, we then generate ground truth potential outcomes under hard interventions that would be unobserved in real-world data. For this, we randomly sample treatment sequences, irrespective of the history, as is done in (Melnychuk et al., 2022). Below, we vary the prediction horizons (up to three days ahead) and the confounding strength. We provide details in Supp. E.
CEN
CRN
CRN
IMEN
BMRN
RMRN
Araug 1000
Gáta
Cobca
Caka
TEC36
TRCDE
CIF-Ne (pem)
Confeandine srenth
Confeandinz sreneh PM
Conferndiaa swc.alaa PV)
(a) One day ahead prediction
(b) Two days ahead prediction
(c) Three days ahead prediction
Figure 3: Performance for the tumor growth model. We compare different forecast horizons and different confound strengths. Shown is the average RMSE of the potential outcomes under hard interventions over five seeds. Our proposed SCIP-Net performs best, followed by our CIP-Net.
The results are in Fig. 3 We make the following observations: (1) Our proposed SCIP-Net per- forms best. The performance gains become especially obvious for increasing levels of vire-ng confounding because ours is the first method to perform proper adjustments in continuous time. (2) Our proposed SCIP-Net performs more robust than the CIP-Net ablation, which demonstrates the effectiveness of our stabilized weights over the unstabilized weights. The stabilized weights help estimating more accurately especially for larger prediction horizons and strong confounding. (3) Nevertheless, our CIP-Net, has a competitive performance. (4) The only baseline designed for continuous time is TE-CDE (shown in green), which we outperform by a large margin. (5) All other neural baselines are instead designed for discrete time (shown in gray) and are outperformed clearly.
Our experiments are based on the MIMIC-III extract from Wang al. (2020). Here, we use real-world covariates at irregular measurement timestamps. and then simulate treat- ments and outcomes, respectively. This is done analogous to (Melnychuk et al., (022). such that we have access to the ground-truth potential outcomes. For the outcome variable, we additionally apply a random observation mask in order to mimic observations at arbitrary, irregular timestamps. We provide more details in Supp. D
MIMIC-III data:
Prediction wandow
CRN
RMSN:
G-Na
TE-CDE
CIP-Net
SCIP-Net
Rel. improvement
1,022= 0.003
1.000
0.065
0.074
0.060
0.915
0,676 0.041
0.272
0.274
0.385
0,784
0.765 0.117
0.034
19.1%
hours
0.232
0.355
1.300
1,300
0.242
1,201
Table 2: Performance for MIMIC-II. Reported is the average RMSE of the potential outcomes under hard interventions over five seeds. We highlight the relative improvement of our SCIP-Net over existing baselines. Again, our proposed SCIP-Net performs best.
Table 2 shows the results for different prediction windows. (1) Our SCIP-Net has the lowest error among all methods and performs thus again best. (2) Our CIP-Net ablation, which is a new method in itself, again has highly competitive performance. (3) Yet, the ablation shows that using stabi- lized weights as in SCIP-Net (as opposed to the unstabilized weights as in CIP-Net) makes a large contribution to the overall performance.
To the best of our knowledge, SCIP-Net is the first neural method for estimating con- ditional average potential outcomes through proper adjustments for time-varying confounding in continuous time. For this, we first derive a tractable expression for inverse propensity weighting in continuous time. Then, we propose stabilized weights in continuous time to stabilize the training objective. Our experiments show that our SCIP-Net has clear benefits over existing baselines when observation times and treatment times take place at arbitrary, irregular timestamps
Conclusion:
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
This work has been supported by the German Federal Ministry of Education and Research (Grant: 011S24082).
REFERENCES
Ahmed Allam, Stefan Feuerriegel, Michael Rebhan, and Michael Krauthammer. Analyzing patient trajectories with artificial intelligence. Journal of Medical Internet Research, 23(12)::29812, 2021.
Heejung Bang and James M. Robins. Doubly robust estimation in missing data and causal inference models. Biometrics, 644)9662,, 2005.
Agamirza E. Bashirov, Emine Misrh, Yücel Tandoğdu, and Ali ÖzyapicI. On modeling with multi- plicative differential equations. Applied Mathematics Journal of Chinese Universities, 26:425 - 438, 2011.
Samuel L. Battalio, David E. Conroy, Walter Dempsey, Peng Liao, Marianne Menictas, Susan Mur- phy, Inbal Nahum-Shani, Tianchen Qian, Santosh Kumar, and Bonnie Spring. Sense2Stop: A micro-randomized trial using wearable sensors optimize just-inttime-ddapiie stress manage- ment intervention for smoking relapse prevention. Contemporary Clinical Trials, 109:106534, 2021.
Inci Baytas, Cao Xiao, Xi Zhang, Fei Wang, Anil Jain, and Jiayu Zhou. Patient subtyping via time-aware Istm networks. In KDD, 2017.
Ioana Bica, Ahmed M. Alaa, James Jordon, and Mihaela van der Schaar. Estimating counterfactual treatment outcomes over time through adversarially balanced representations. In ICLR, 2020.
Ioana Bica, Ahmed M. Alaa, Craig Lambert, and Mihaela van der Schaar. From real-world patient data to individualized treatment effects using machine learning: Current and future methods to address underlying challenges. Clinical Pharmacology and Therapeutics, 109(():87-100, 2021.
Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. GRU-ODE-Bayes: Continuous
modeling of sporadically-obsevee time series. In NeurlPS, 2019.
Zhengping Che, Sanjay Purushotham, KyungHyun Cho, David Sontag, and Yan Liu. Recurrent neural networks for multivariate time series with missing values. Scientific Reports, 8(1):6085, 2018.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary dif- ferential equations. In NeurIPS, 2018
Yuqi Chen, Kan Ren, Yansen Wang, Yuchen Fang, Weiwei Sun, and Dongsheng Li. Contiformer: Continuous-time transformer for irregular time series modeling. NeurlPS, volume 36, 2023.
Alicia Curth and Mihaela van der Schaar. Nonparametric estimation of heterogeneous treatment effeetts: From theory to learning algorithms. In AISTATS, 2021.
Nan Du, Hanjun Dai, Rakshit Trivedi, Utkarsh Upadhyay, Manuel Rodrez-Rodriguez and Le Song. Recurrent marked temporal point processes: Embedding event history to vector. In KDD, 2016.
Stefan Feuerriegel, Dennis Frauen, Valentyn Melnychuk, Jonas Schweisthal, Konstantin Hess, Ali- cia Curth, Stefan Bauer, Niki Kilbertus, Isaac S. Kohane, and Mihaela van der Schaar. Causal machine learning for predicting treatment outcomes. Nature Medicine, 30:958 2024.
Dennis Frauen, Tobias Hatt, Valentyn Melnychuk, and Stefan Feuerriegel. Estimating average causal effects from patient trajectories. In AAAI. 2023.
Dennis Frauen, Konstantin Hess, and Stefan Feuerriegel. Model-agnostic meta-learnen for estimat- ing heterogeneous treatment effects over time. In ICI.R, 2025.
11

---

Published as a conference paper at ICLR 2025
Crispin W. Gardiner. Handbook of stochastic methods for physics, chemistry and the natural sci- ences. Springer-Verie Verlag, Berlin; Heidelberg, 1985. ISBN 3540156070.
Changran Geng, Harald Paganetti, and Clemens Grassberger. Prediction of treatment response for combined chemo- and radiation therapy for non-small cell lung cancer patients using bio- mathematical model. Scientific Reports, 7(1):13542, 2017.
Richard D. Gill and Soren Johansen. A survey product-integration with a view toward application in survival analysis. Annals of Statistics, 18(4), 1990.
Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34(1):014004, 2018.
Tobias Hatt and Stefan Feuerriegel. Sequential deconfounding for causal inference with unobserved confounders. arXiv preprint, uxXvv2104.093233 2021a.
Tobias Hatt and Stefan Feuerriegel. Estimating average treatment effects via orthogonal regulariza- tion. In C/KM, 2021b.
Konstantin Hess, Dennis Frauen, Valentyn Melnychuk, and Stefan Feuerriegel. G-transformer for estimating conditional average potential outcomes over time. arXiv preprint, arXiv:2405.21012, 2024a.
Konstantin Hess, Valentyn Melnychuk, Dennis Frauen, and Stefan Feuerriegel. Bayesian neural controlled differential equations for treatment effect estimation. In ICLR, 2024b.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Computation, 9(8): 1735-1780, 1997.
Çaglar HIzh, ST John, Anne Juuti, Tuure Saarinen, Kirsi Pietiläinen, and Pekka Marttinen. Causal modeling of policy interventions from sequences of treatments and outcomes. In ICML, 2023.
Guido W. Imbens and Donald B. Rubin. Causal inference for statistics, social, and biomedical sci- ences: An introduction. Cambridge University Press, Cambridge, 2015. ISBN 9781 9781139025751.
Fredrik D. Johansson, Uri Shalit, and David Sonntag. Learning representations for counterfactual inference. In ICML, 2016.
Alistair E. W. Johnson, Tom J. Pollard, Lu Shen, Li-wei H. Lehman, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G. Mark. MIMIC-III, freely accessible critical care database. Scientific Data, 3(1):160035, 2016.
Nathan Kallus, Xiaojie Mao, and Angela Zhou, Interval estimation of individual causal effects under unobserved confounding. In AISTATS, 2019.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled diffcrential equa- tions for irregular time series. In NeurIPS, 2020.
Patrick Kidger, James Foster, Xuechen Li, Harald Oherhauser, and Terry Lyons. Neural SDEs as Deminte-Dimensionall GANs. In ICML, 2021.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
Milan Kuzmanovic, Tobias Hatt, and Stefan Feuerriegel. Deconfounding temporal autoencoder:
Estimating treatment effects time using noisy proxies. In ML4H, 2021.
Rui Li, Stephanie Hu, Mingyu Lu, Yuria Utsumi,, Prithwish Chakraborty, Daby M. Sow, Piyush Madan, Jun Li, Mohamed Ghalwash, Zach Shahn, and Li-wei Lehman. G-Net: A recurrent net- work approach to G-computation for counterfactual prediction under a dynamic treatment regime. In ML4H, 2021.
Bryan Lim, Ahmed M. Alaa, and Mihacla van der Schaar. Forecasting treatment responses over time using recurrent marginal structural networks. In NeurlPS, 2018.
12

---

Published as a conference paper at ICLR 2025
Judith J. Lok. Statistical modeling of causal effects in continuous time. Annals of Statistics, 36(3), 2008.
Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong. Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. In ICML, 2018.
Yuchen Ma, Valentyn Melnychuk, Jonas Schweisthal, and Stefan Feuerriegel DiffPO: A causal diffusion model for learning distributions of potential outcomes. In NewrIPS, 2024.
Farokh Marvasti. Nonuniform sampling. Springer US, 2001, 2001. ISBN 9781461512295.
Valentyn Melnychuk, Dennis Frauen, and Stefan Feuerriegel. Causal transformer for estimating counterfactual outcomes. In ICML, 2022.
Valentyn Melnychuk, Dennis Frauen, and Stefan Feuerriegel. Bounds on representation-innuued confounding bias for treatment effect estimation. In ICLR, 2024.
James Morrill, Patrick Kidger, Lingyi Yang, and Teny Lyons. Neural controlled differential equa- tions for online prediction tasks. arXiv preprint, 2106. 11028, 2021.
Elizabeth Murray, Eric B. Hekler, Gerhard Andersson, Linda M. Collins, Aiden Doherty, Chris Hol- lis, Daniel E. Rivera, Robert West, and Jeremy C. Wyatt. Evaluating Digital Health Interventions: Key Questions and Approaches. American Journal of Preventive Medicine, 51(5):843-851, 2016.
Jerzy Neyman. the application of probability theory to agricultural experiments. Annals of Agriculural Sciences, 10:1-51, 1923.
Yilmazcan Özyurt, Mathias Kraus, Tobias Hatt, and Stefan Feuerriegel. AtDMM: An attentive deep
Markov model for risk scoring in intensive care units. In KDD. 2021.
James M. Robins. A new approach to causal inference in mortality studies with sustained exposure period: Application to control of the healthy worker survivor effect. Mathematical Modelling, 7: 1393-1512, 1986.
James M. Robins. Robust estimation in sequentially ignorable missing data and causal inference models. Proceedings of the American Statistical Association Bayesian Statistical Science, pp- 6-10, 1999.
James M. Robins and Miguel A. Hernán. Estimation of the causal effects of time-varying exposures. Chapman & Hall/CRC handbooks of modern statistical methods. CRC Press, Boca Raton, 2009.
Kjetil Røysland. A martingale approach to continuous-time marginal structural models. Bernoulli, 17(3):895 915, 2011.
Donald B. Rubin. Inference and missing data. Biometrika, 63(3):581, 1976.
Donald Rubin. Bayesian inference for causal effects: The role of randomization. Annals of Statistics, 6(1):34-58, 1978.
Helene C. Rytgaard, Thomas A. Gerds, and Mark J. der Laan. Continuous-time targeted min- imum lass-based estimation of ciarventionnsseciife mean outcomes. The Annals of Statistics, 2022.
Helene C. Rytgaard, Frank Eriksson, and Mark J van der Laan. Estimation of time-specific interven- tion effects continuously distributed time-to-event outcomes by targeted maximum likelihood estimation. Biometrics, 79(4):3038-300 -0499 2023.
Maresa Schröder, Dennis Frauen, Jonas Schweisthal, Konstantin Hess, Valentyn Melnychuk, and Stefan Feuerriegel. Conformal prediction for causal effects of continuous treatments. arXiv preprint, arXiv:2407.03094, 2024.
Peter Schulam and Suchi Saria. Reliable decision support using counterfactual models. NeurlPS, 2017.
13

---

Published as a conference paper at ICLR 2025
Nabeel Seedat, Fergus Imrie, Alexis Bellot, Zhaozhi Qian, and Mihaela der Schaar. Continuous- time modeling of counterfactual outcomes using neural controlled differential equations. In ICML, 2022.
Uri Shalit, Fredrik D. Johansson, and David Sontag. Estimating individual treatment effect: Gener- alization bounds and algorithms. In ICML, 2017.
Oleksandr Shchur, Ali Caner Türkmen, Tim Januschowski, and Stephan Günnemann. Neural tem- poral point processes: A review. IJCAI, 2021.
Yi Shirakawa, Toru; Li, Yulun Wu, Sky Qiu, Yuxuan Li. Mingduo Zhao, Hiroyasu Iso, and Mark van der Laan. Longitudinal targeted minimum loss-based estimation with temporal-difference heterogeneous transformer. In JCML, 2024.
Hossein Soleimani, Adarsh Subbaswamy, and Suchi Saria. Tresoment-response models for counter- factual reasoning with continuous-time, continuous-valuee interventions. In UAI, 2017.
Mark J. van der Laan and Susan Gruber. Targeted minimum loss based estimation of causal effects of multiple time point interventions. The International Journal of Biostatistics, 8(1):Article 9, 2012.
Toon Vanderschueren, Alicia Curth, Wouter Verbeke, and Mihaela van der Schaar. Accounting for informative sampling when learning forecast treatment outcomes over time. In ICML, 2023.
Shirly Wang, Matthew B.A. McDermott, Geeticka Chauhan, Marzyeh Ghassemi, Michael C. Hughes, and Tristan Naumann. MIMIC-extract: A data extraction, preprocessing, and repre- sentation pipeline for MIMIC-III CHIL, 2020.
Yanbo Xu, Yanxun Xu, and Suchi Saria. non-parametric bayesian approach for estimating treatment-response curves from sparse time series, In ML4H, 2016.
14

---

Published as a conference paper at ICLR 2025
A EXTENDED RELATED WORK
Irregular time series: Outside of the causal literature, stochastic processes in continuous time and irregular time series are heavily studied area of research, in particular in terms of stochastic processes [Gardiner 1985). Directly related to the continuous time literature are non-uniformly spaced observations and their impact on inference (Marvasti (001). Irregular observation times can also be seen as a special case of missing data (Rubin 1976). More recently, neural network architectures handling irregularity involve explicit incorporation of time gaps into recurrent units [Baytas et al. 2017). or otherwise handling missing values in recurrent neural networks (Che et al., 2018 Brouwer et al. 2019). Further, transformers have been adapted for irregular time series (Chen) al. 2023) and neural point processes Du et al. (2016); Shchur et al. (2021). Finally, and directly related to our work, neural controlled differential equations (Kidger et al. 2020 2021 Morrill et al, 2021) can directly handle irregular time series by design. However, these works focus traditional estimation and not causal inference tasks. In other words, the above methods solve a different task and would be biased in our setting.
Conditional average potential outcomes in the static setting: large body of research focuses
on methods for estimating conditional average potential outcomes in the static setting (e.g.. g..... Curth & van der Schaar [021]] Hatt & Feuerriegel, (2221b] Johansson et al. [016)] Kallus et al. 2019;
| Ma
al. 2024; Schröder et al. 2024] (hhllit al. 2017). However, they cannot adequately address the complexities of time-varying data, which are crucial in healthcare applications where patient
conditions evolve over time (e.g., EHRs (Allam al. 2021: Bica al., 2021). wearable devices
(Battalio et al. 2021: Murray et al. 2016)).
Nonparametrie methods in continuous time: Some nonparametric methods have been proposed (Hatt& & Feuerriegel| 2021aa] Hizh et al. 2023) Schulam & Saria 2017, /oolemani et al.] 2017; Xu al. 2016), but they suffer from scalability issues. Further, some struggle with static and high- dimensional covariates, complex outcome distributions, or impose additional identifiability assump- tions. As result, we focus on neural methods for their flexibility and scalability.
Additional research directions: Hess et al. (2024b) propose neural method for uncertainty quan- tification of conditional average potential outcomes continuous time. Further, Vanderschueren
et al. (023) develop framework to account for bias due informative observation times. Both approaches are orthogonal to our work and do not focus on time-varying confounding. Further, additional research directions in discrete time consider nonparametric learners (Frauen et al.) 2025) and robust estimation of the G-computation formula (Hess et al.] (202a)). Other works focus on learning treatment effects from noisy proxies (Kuzmanovic al. [002)). Finally, there are efforts to leverage neural methods for estimating average potential outcomes (Frauen et al. [2023] Shirakawa al., 2024]. However, these methods are not applicable our setting, as the do not operate in continuous time and do not focus on estimating CAPOs.
15

---

Published as a conference paper at ICLR 2025
B PRODUCT INTEGRATION
Product integration is useful tool for describing stochastic processes and their joint likelihoods in contimuous time (e.g., survival analysis; see Gill Johansen 1990. Rytgaard et al. 2022 2023).
In the following, we provide brief introduction to product integration. Importantly, the type that we refer to is also known as the geometric integral or the multiplicative integral (Bashirov et al, 2011
The product integral π. can intuitively be thought of as the infinitesimal limit of the product operator II, or, equivalently, as the multiplicative version of the (remann-integras f. That is, the standard Riemann-integral is defined as
/ d(z) = lim E(t,)A. la.0]
(43)
Ax
where Ar = T++1 x2 and F E (I2, x411) for a partition U, 1441] of [a, b].
The product-integral is analogously defined as
U (a.b]
[())
lim Ax
[Z:)
(44)
For the proof's in Supp. D we make use of the identity
π [a,b]
[ax)d exp
(
log (xx dx))
(45)
[a.M]
which holds since
π f(x)dr= = [a.b]
lim
II f())
(46)
Ar
lim Az-+
cxp
log
f(z,)
(47)
lim exp Ax+0
Σ log (Tx) (m.))a))
(48)
= =exp
lim Ar
Σ log f(t) )00)
(49)
(ʃ
= =exp
log (xx dx)
(50)
(a.0)
Similar to the additivity of the standard integral f, the product integral ][ is multiplicative. That is, we have that for a <
π f(
π f(x)dzz Troont
(51)
[0,0]
(0.0] [0,0]
16

---

Published as a conference paper at ICLR 2025
which can be shown with the above exp-log identity
π
(f
f()
exp
log (xx) dr)
(52)
[a.c]
[a,0]
(f
= exp
log f(z) dr +
log f(x) dz)
(53)
[a.0]
[b.c]
(f
exp ( f
exp
log (x) dr)
log f(x) dx)
(54)
[a.b]
[b.c]
(a))de dz π f(x)dr
(55)
[0,0]
[b,c]
17

---

Published as a conference paper at ICLR 2025
C NEURAL DIFFERENTIAL EQUATIONS
We provide a brief summary of neural ordinary differential equations and neural controlled differ- ential equations, similar that in (Hess al. 2024b).
Neural ODEs: Neural ordinary differential equations (ODEs) (Chen et al. 2018 Haber & Ruthotto 2018 Lu al. 2018 integrate neural networks with ordinary differential equations. In neural ODE, the neural network f-(-) defines the vector field of the initial value problem
Z = L// 80( Zo X.
(56)
Thereby, neural ODE captures the continuous evolution of hidden states Ze over time scale. Thereby, learns continuous flow of transformations, where the input = Zo is passed through an ODE solver to obtain the output Ý = Z, (possibly, after another output transformation).
While neural ODEs are generally suitable for describing continuous time evolution, they have an important limitation: all data need to captured the initial value. Hence, they are not capable of updating their hidden states new data becomes available time, is the case for, e.g., electronic health records.
Neural controlled differential equations (CDEs) Kidger al.) 2020 overcome the above limitation. Put simply, one can think of them as a continuous-time counterpart to recurrent neural networks. Given path of data X, TRd,, t [0, r]. neural CDE consists of an embedding network (.).), readout network (-)-), and a neural vector field fa. Then, the neural CDE is defined as
Neural CDEs:
( [6(2,,)) EE (0, with Zo ==>((o)
= (7+)) Z:
(57)
where Ze Rd and tZ, Rd.xdd, The integral here is Chemann-Stieltie integral, where fo (Z,, 8) d [[]] corresponds to matrix multiplication. In this context, the neural differential equation controlled by the process [X.]. Importantly, one rewrite it (under some regularity conditions)
as
Zi z- f ds
dx. ds. tE (0,
(58)
Calculating the time derivative requires C -e-representation of the data X, for all E [), As result, irregularly sampled observations ((lo, Xo), (41, X.)...., (In, Xn)) must be interpolated over time, producing a continuous representation X,. Here, we can use any interpolation scheme such as linear interpolation (Morrill et al. 2021) in TE-CDE (Seedat et al. (202). Then, the neural CDE essentially reduces to neural ODE for optimization.
The key difference to neural ODEs is, however, that the neural vector field is controlled by quentially incoming data and, therefore, updates the hidden states as data becomes available over time. This is clear advantage over ODEs, which require capturing all data in the initial value and because of which neural ODES are not suitable.
18

---

Published as a conference paper at ICLR 2025
D PROOFS OF PROPOSITIONS
Proposition 1. Under assumptions (i)-(ii), we can estimate the conditional average potential out- come from observational data (i.e..from data sampled under ol via inverse propensity weighting, that is,
π W (A. (e., He- he-) (59)
E Y 1. D.,
Hi- hi-
E Y
where the inverse propensity weights for s 2 are defined as
W₂ = w. (H.) = dG.,,HH,
(60)
oG (H₂)
Proof. The proof follows similar ideas as in (Rytgaard et al. 2022
Y- la., r.
Hi_ Trc]]
(61)
/ la., e π (h.) (Go,a(ha)
(62)
f vr vela,, a
(h.)
(63)
62e
f Mr (1..)
π
dQu.a (A) = (..)) dGna(hs (A) = (g.,er mi,,))
(64)
(1)
y1 dQo (h, Ni) (. dG. (A) (a..:
(65)
dG.,, (h, Ni) = (g, .: ,.,,,)
Ur
dGo, (ha (A, @,, 1(,,,))
x
dQo.e (A,, No) = dGo.o NA) = (a., 1.
(66)
EE Y-
(He (A, (0. (o.. (e, flr hr-
(67)
= B Y=
Wa (Hz) N2) = (a.,,, n. Hr- = he-]
(68)
F.
W A..N5 (a.,: He- Al- hr-
(69)
22
where Equation 66 uses the multiplicativity of the product integral.
19

---

Published as a conference paper at ICLR 2025
Proposition 2. Let t2,0 t for notational convenience. The unstabilized weights in Equation 72 satisfy
π w. (A4). = (I,,, H. Ho_ = -h--)= We, ((A,. = (L, Hy- = Tn-))
(70)
where
exp SeE 18(s Ha_ ds
We:, 18(12)) Hig Toss (an,1 Heg
(71)
Proof.
The inverse propensity weight in Equation 12 is by definition given by
π
dG., (Ht) (H,))HH.) 22
π
dA II,-)T,, (A. H (dsss)
dA (1.-)
1-N
dAg(s H2-) H
- dAg(s Hs-)
(72)
We simplify the first part via
π
dA H
Hs-) N2(ds)
(73)
dAg(s H ) Hg-))
#2
S (A, (H,-) ds (ds)
π
(1 H (
(74)
(8(s H-- To, (A, H
#2
(5 11. (1.))To, (As H4s) No (das)
exp
log
(75)
Xg(s H HH,-)roo ) H1s-)
s21
(T Hrg )mo Ty (Ary Hrg-)) )
(I) Hig- H TT, (Ary Hrg-))
exp
log
(76)
TEETSNTE
XT(T)
H ( (Ary Hrg-
(77)
ATTT
Hry Mo (Arz Hry-)
TTT
Now, since we condition on N D² we know that the jumping times are given by
T TP ((2g)/=
(78)
Hence, we have that
XT
Hy (Ar; (Hrg-) ((A,. = (i,,, Hr -hr-)
(79)
(TT HT (Arę HTg-)
ET
X (t
Ale, [ )T,, (ante, (
Rr,e): Mee -oog-)Joore (ar,t (y ((dr, = (a,, (.,.).1...). Hir ho-)
(80)
II
(t -ee To, (ar.e, ((a--)) (AA,, = (a.,::27.)). Hi = hr-))
==1
(81)
20

---

Published as a conference paper at ICLR 2025
Further, can simplify the second part via
1 (s H.-) (de)
π
(82)
dAg (1-))
Pa(N" Hx-))
(da)
(83)
(os Ha-)
#2f
K
P. P. 1x 1)) = H(r-))
1- ([6,tlt)))))))
lim
II
(84)
K -000
(PoNa(lle, =0 Hip
A=1
exp
(s Ho- dss
exp
No(s H
(85)
Jazt
a22
exp
( H
Ha- ds
(86)
x2t
#22
where U [tx, tx+1) is a disjoint partition of [t, r] with limg-sso max (tx+11 txx [x)==0 = Fi- nally, again by conditioning No nh,r, have that
[N n2.)) = - A A0 (s) (Nr, 22,1 ) = =0, (87)
(s H
which leaves
(A Hs-) (dN=(do)
π #2!
((4,, ((d,, He- = h
(88)
(A
(1,-)
exp
(s
HH-)dd ds
(14, (L 1i,)) H Tre.).
(89)
Combining Equation 81 with Equation 89 the unstabilized weights are then given by
π w. (AA,.
Hi. he-)
(90)
82!
π .2t
dG.,a (As))
|
(A4,. (a.,.) :,,,); Fr. = hr-))
(91)
dGo,
exp
(a Ha_) ds
x II
(88(1
Hee, tro.e, (amte Hey
-1 11 ((A, = (e.,): Hi_ = hr-))
1=1
(92)
exp Secice (e Ag(s Hg-) ds
((4,..7)) = (i,,), Hr = hu-)
(93)
Xi(te H -To,e,, (a..1., Hos
1=1
II Wes (A) (a.,1. H Hi hv-).
(94)
1=1
21

---

Published as a conference paper at ICLR 2025
Definition 1. For t, let the scaling factor Eg be given by the ratio of the marginal transition probabilities of treatment, that is,
E Eo Noy = (As,
(95)
dG., As, No)
We define the stabilized weights Ws as
E,W
(96)
Proposition 3.
The optimal parameters & in Equation 14 can equivalently be obtained by
arg min B (y.-mo(A,, (o N?, H--))
W. NT) (00,,) Fr. hi-
(97)
#2
Proof. The scaling factor is constant conditionally on the future and past treatment propensity and frequency, i.c.,
E, (A,, (a.,t: n2,)) Ho_ hi-
(98)
=E. (A, Not 1), (Ac-, NN_ = (az_, AT-)
(99)
Eo (1 ] [ 1)
(100)
Ecnnt..,
(101)
where we use
[āa,s) [c-
UU 15rsi
U
and)
(Uar) and r<t
[nC n
(U nnr 15rSt
(0(Ung) rst
m
(102)
for the concatenation of interventional and observational treatments at time t. Hence, with E. E IRt can use linearity of the expectation and multiplicativity of the product integral, such that
arg min Ero ((r, (0% (Y5-mm((4 (A No, H. ))²
π
Ws) N3) (a.,. n.,,), Hi. he- a2
(103)
arg min Ero ((x, - my (4,, No, H)-
I E.W.) (A, = (1, Flo = hr-) (104)
arg min {T Eallar..: ae-1, E
(105)
π
x EFo (8= (A, No, fl. 11(-))
We
(A, = (a.,, 1,,,,, Fy. For
(106)
arg min Bro (8- mo(A, Nr. ))²
W. (A,. (a..: Fr. Tu-]
(107)
62
(108)
22

---

Published as a conference paper at ICLR 2025
Proposition 4. Let t2,0 then satisfies
t for notational convenience. The scaling factor Eg from Equation 77
π $21
NA) (a.,; (AL,, Ni_) (āc- nT-))
(109)
1 j=1
Eet ((4 (en, (A,, = (Q)=)())))
(110)
where
1g(t
Ar, - Ni. ko,e%, (ar:, , No.
(111)
exp exp]eetee Ssetes ,, 28(s) Ax- No. N ds
Proof. For the proof, we follow the steps as in the proof of Proposition By definition, E, satisfies
π E π +2
dGo. (A,, No)
(112)
(Ao,
+2
= -π
dAS(s A No. DTO. (Ao
Aoo, No_)
(²(a))
dA
A-_ No
(dx
(113)
ds Aa_,, No )o (As As_, No.
(A A._, No.
First, we simplify
π
dAg(s A.-, No.. J
Nx_ (dsds)
(114)
dA A,_, No. (T, Ag_, No
22
π
Aa A._, - TA,
- (de)
(115)
(2 Ao_, Ny ) d.
To.
Ao_, N
22
A A._, No. (A,
exp
log
(2 Â._, Ni. (T.
No_,No_ N (ds)
(116)
16(TT Ang - ( (Ara Ang- Nio
exp
log
(T Arg No. IT (Arg
Arg Nip
(117)
TEETPNTE
XB(TT
Arg- Nig. To (Arg Ang Nig
II TEETENT
(118)
(TT [Arg-. N ,T (Arg
AArg- Noj.
Conditionally on N , t. the jumping times are fixed, i.e.
T Te = (12.)))-a
(119)
Therefore, it follows that
X6(T)) Arg- Nig - " ry (Arg Arg Ang--. NNgg)
((A,, = (i,) No_) = (ā.-, n2-
TgeTs XS Ary- Ng.. - Jo (Ary Arg-, N
(120)
(1 A N2 JTT (an, Aee,-, No. )
((A,. = (a.,: (A,,, No) = (ā._, T(p_)
5=1 (t Ace Ng.
JT (an.e Aee -, No.
(121)
1 AL | A - -)oo.e (ar,e
((4, = 1e,,) (Ae_, N;-) = (a.-,
(Ae,,
(122)
23

---

Published as a conference paper at ICLR 2025
For the second part, we also follow the steps as in Proposition 2 that is,
π
ds A._, No-
1-N7(da)
(123)
ds As-; No_
1 Po(N (ds) = 1 A._, N.-) (do)
(124)
1 P. (No = 1 N..,
x2
Po(N ([te (x+1)) )) = AA,_,,,,, )
lim
(125)
K
P.(No (++1)) =0 | N1x-)
exp
(s As-, No- )d]] /
exp
/ \o Ao_, N2-
(126)
exp
AC(s| | A.-, No_ ds AB(s | Ao_, N2_ N
(127)
.21
- where U [1x: (+11 is a disjoint partition of [t, r] with kim max (1x+1 (1++1-4)) 0. Fi- nally, again by conditioning on No- n2 have that
(s A._, N.- (Ni ni..) = / 1 112,)321 (Nz n..1) =0, (128)
s2t
which leaves
π #21
1 - dA Ao_, No-) (ds) ((4,. = (e,,: (A,, (ai, (129)
I ds Ac_, No-)
[-
1b | Ac_, No ds (44, = (a D. (Ap_, No_) = (āz_, m2 m()) (130)
exx
Finally, we combine equation 122 with equation 130 Hence, the scaling factors satisfy
π E. (A = (e., Di
= (āc_, mi-))
(131)
π dGp, dG (A,,
dG.,. (A,, N2) (A,, = m1 (Ae_, Nz. = (a-,nnt))
(132)
exp
Xo(s| No ds}
II At Are, Ae No., -moor,, (anes, No (44, / = (.p Nº-) = (oe_ 138
(133)
X6 (t2.))) Are -, No. Iuarre . (aner,, N ((A,, Ni) = (a.): (Ac_, = (ā.., m7.
exp SopSsete sepe Nb(s Ao-, No_ )
(134)
Beg, (AA.. (Q., ni,) (Ae_, Bab,,No NE) = (ār., nÿ nz-)).
(135)
1=1
24

---

Published as a conference paper at ICLR 2025
E DATA GENERATION
E. TUMOR GROWTH DATA
The tumor data used in Sec. was simulated based on the lung cancer model proposed by Geng ct al. (2017), which has been previously used in several works (Lim et al. 2018 Bica et al. 2020
Li al. 2021 Melnychuk et al. 2022 Seedat et al. 2022 Vanderschueren al. 2023).
Specifically, we adopt the simulation framework introduced by Vanderschueren et al. (2023). which includes irregularly spaced observations. However, different their work, we explicitly add con- founding bias to the treatment assignment (that is, both treatment times and treatment choice).
The tumor volume is the outcome variable. evolves over time according to the ordinary differential equation
K V,
dy;- = + log
Cr
-a,,d + Br + 6t
de,
(136)
Chemotherpyy
Radiotherapy
Noise
Tumor growth
where P is the tumor growth rate, K represents the carrying capacity, and Cee. Or, and B, control the effects of chemotherapy and radiotherapy, respectively. The term e, introduces randomness into
the dynamics. The parameters were drawn following the distributions with details provided in Table The variables and de represent chemotherapy and radiotherapy
in [Geng al. (2017).
treatments, respectively, and follow previous works (Lim et al. 2018 Bica al. 2020 Seedat al.
2022). Time t is measured in days.
Variable
Parameter
Distribution
Value (μ, o²)
Growth parameter
Normal
(7.00 x 10-3, 7.23 x 10-3) 30
7.23
Tumor growth
Carrying capacity
K
Constant
Radiotherapy
Radio cell kill
Or
Normal
(0.0398, 0.168) Set to 8, 10 X
Radio cell kill
Br
Chemotherapy
Chemo cell kill
Normal
(0.028, 7.00 10
Noise
Normal
(0, 0.01²)
Table 3: Parameter details for the synthetic data generating process.
The radiation dosage de and chemotherapy drug concentration e, are applied with probabilities The treatments are administered according the following, history dependent treatment probabilities
AC, AC Ber (a (0(- -
7 Dmax
-Dmax/
(137)
where y controls the confounding strength, Dmux is the maximum tumor volume, D15 the average tumor diameter of the last 15 time steps, and & controls the confounding strength. For test data, want evaluate the potential outcomes under hard interventions. Hence, uniformly sample random treatment sequence and apply it to the outcome, irrespective of the history, as is done in (Melnychuk et 2022).
Importantly, this treatment assignment process is only used for training and validation. For testing, we randomly sample hard interventions as described in Sec. 3 Thereby, we can directly investigate how all baselines perform under time-varying confounding.
We add an observation process that randomly masks away observations of the outcome variable. Hence, at some days, the tumor diameter remains unobserved. Irregular observations times are the domain that continuous time methods are tailored for.
For this, our setup is consistent with Vanderschueren et al. cess as a Hawkes process with intensity (1) given by
2023 who define the observation pro-
Dr
1)
(t) sigmoid
(138)
Dret
2
25

---

Published as a conference paper at ICLR 2025
where determines the informativeness of the sampling, Duet = 13 represents the reference tumor diameter, and D, is the average tumor diameter over the past 15 days. In this work, however, main focus is not informative sampling, which is orthogonal research direction. Therefore, we opted for setting the informativeness parameter to w = Thereby, we are in the setting that is known sampling completely random.
Following Kidger et al. (2020) we added a multivariate counting variable that counts the number of observations up each day, respectively. We normalized this counting variable with the maximum time scale 7.
Finally, consistent with (Lim et al.
2018 Bica al. 2020 Seedat et al. 2022
Vanderschueren
et al., 2023), we introduced patient heterogeneity by modeling distinct subgroups. Each subgroup differs in their average treatment response, characterized by the mean of the normal distributions. Specifically, for subgroup A, we increased the mean of a, by 10%, and for subgroup B, we increased the mean of ag by 10%.
The observed time window for training, validation, and testing is set to = 30 days. We generate 1000 observations for training, validation, and testing, respectively.
26

---

Published as a conference paper at ICLR 2025
E.2 MIMIC- III DATA
For our semi-synthetic experiments in Sec. 5. we upon the MIMIC-extract dataset (Wang et al. 2020), which is based on the MIMIC-III database (Johnson al. 2016) and widely used in research (e.g., Ozyurt et al., 2021). Importantly, measurements in this dataset have irregular timestamps for different covariates. Therefore, we can directly use this missingness without artificially introduces any masking process for covariates.
time-varying covariates (i.e., vital signs) alongside the static covariates gen- der, ethnicity, and age. As we are interested in conditional average potential outcomes, we need to introduce a synthetic data outcome generation process. For this, we simulate two-dimensional out- come variable for training and validation purposes and generate interventional outcomes for testing. As defined Sec. 3 we add past observed outcomes to list of covariates. Hence, have have dg 14-dimensional covariate space. The outcome-generatin process follows (Melnychuk et al. 2022
In our setup,
Simulating untreated outcomes: We first simulate two untreated outcomes Y22 for j = 1,2 as fol- lows:
Y = B-splinee + og + of t
(139)
where ag, og. and a² are weight parameters. Here, B-spline(t) is drawn from mixture of three cubic splines, and ()-()) is a random function approximated using random Fourier features from a Gaussian process.
Simulating treatment assignments: We simulate da = 3 synthetic treatments A for = 1,2,3 according to:
A Ber(p)), = (N Th (Xi) +b))
(140)
where 7t and 7 are parameters that control the influence of past treatments and covariates on treatment assignment. Y.A. represents summary of previously treated outcomes, b' is a bias term, and f--(-) is another random function sampled using a random Fourier features approximation of Gaussian process. For test data, we want to evaluate the potential outcomes under hard interventions. Hence, uniformly sample a random treatment sequence and apply it the outcome, irrespective of the history, as is done in (Melnychuk et al. (022)).
Applying treatments to outcomes: Finally, for training and validation, treatments are applied to untreated outcomes Yt using:
min/=1.....d de 1 14 (wl - 8)2
Y = Vt
(141)
im Empa
where wi defines the duration of the treatment effect window, and ghd determines the maximum effect of treatment A' on outcome Y2. Importantly, we do not follow this treatment assignment mechanism for testing. Instead, as we are interested in estimating conditional average potential
outcomes, we randomly assign hard interventions as in Sec. 3
Masking the outcome: Finally, we add an observation mask to the outcome variable Y?. For this, we randomly mask away the outcome variable with observation probability p = 0.15
In our experiments Sec. we used 1000 samples for training, validation and testing, respectively. For testing, we simulate 50 different intervention sequences per patient. The time window was set between 30 < < 50.
27

---

Published as a conference paper at ICLR 2025
F ADDITIONAL RESULTS
F.1 VARIATION OF OBSERVATION INTENSITIES
In the following, evaluate the stability of our SCIP-Net for different sampling intensities. For this, we use the tumor growth model as in Section 5 In our main study, observation times follow intensity process as informed by prior literature (Vanderschueren et al. 2023
a
Here, the observation probabilities (t) are given by
De -
(1) sigmoid
(142)
where determines the informativeness of the sampling, Diet 13 cm represents the reference tumor diameter, and D, is the average tumor diameter over the past 15 days. More details are provided in Supplement E.
In our main study in Section S we focused on the sampling completely at random setting, where we set the informativeness parameter = O. In the following, increase the informativeness param- eter up to w = 0.5. Tables and [ show the performance of our SCIP-Net against the baselines. We find that our SCIP-Net has very robust performance and consistently outperforms the baselines.
Informalivenes Prediction window observation days
G-Net
RMSNs
TE&CDE
SCIP-Net
Melnychuk
2022
1.09
0.95
1.29
7,85
4.33
1.64
2.35
1.45
3.20
1.72
6.10
6.65
0.86
0.88
0.68
7,95
5.63
1.44
1.43
1.19
7.75 1.66
7.60
1,00
0.30
0.51
1.16
5.47
1.40
1.26
1:47
4.69
3.05
0..73
7.64::00.00
8:65
5.48
1.64
1.72
7.76
4.74
0.77
7,06
5.97
1.50
1.57
1.53
1,22
5.97
0.883
1.008
8,96
1,52
1,10
1.54
1.09
Table 4: Informative sampling: Performance for the tumor growth model with irregular sampling times and confounding strength 7 8. We vary the informative sampling parameter w. Our SCIP-Net has robust performance and consistently outperforms the baselines.
Informativenss.o Prediction window days
G-Net D0T
RMSNs
CRN
TE-CDE
SCIP-Nit
observation times
Li
Meinychuk
2022
2018
Dica 2000
Seedat
8.76--00 0.10
9.55 0.75
9.52 0.05
12.15++0.71
10.02:+
5.13 0.44
1.29
1.16
$.95110.91
8.09.10.00
5,30 0.08
0.38
5.61
6.74
6,75
4.47 0.39
9:10 0.32
8.71 0.73
10.08 0.00
11.11 1.06
11.15+0.05
4.99 0.74
7.45 1.08
1 0.54
1.35 1.04
1.011.111
8.42.1..22
4.84
6.17 0.05
5.55
0.71
5.01
7.14
6,95
4.53
0.90
5.00 0.32
0.999
0.85
11.30 1.300.00.00
10.92.10.11
5.37
1.47
7.32 1.11
0.32
7.69 1: 1.0.5
1.02
8.55110005
5.03
0.50
0.58
0.45
6.994
8.81 0.43
8.65 0.00
0.72
11.28 0.000
10.01
6.33
2.31
2
7.25 1.14
1:11::00 0.41
1.89 1: 1.05
6.00
8.74 0.64
6.04
2.14
0.99
5,80
1.79
7.25 1.791.0.74
2.30
50 0.32
5.42 1.00
8.74 1.02
11.06 0.39
11.50:+:100 1.04
4.50 0.38
1.08
0.82
0.89
0.97
4.79
0.58
0.44
0.64
11101000 0.78
6,91
0.42
1.43
1.25
11,08
4.80
1.11
0.77
0.49
1 1.00
8.60
4.04
5.60 0.05
5.57
0.38
5.89
2.221.0..8 0.80
7.07 107.1.0.0
4.72 0.73
Table 5: Informative sampling: Performance for the tumor growth model with irregular sampling times and confounding strength 7 6. As in Table we vary the informative sampling parameter w. Our SCIP-Net again demonstrates superior performance and outperforms all baselines.
28

---

Published as a conference paper at ICLR 2025
F.2
ABLATION STUDY: SDIP-NET FOR DISCRETE TIME
strength days
Molryyau ctal,
1.401 1.11
2.33
1.75 0..7
2:4.1:0.2
1.53 = 0.12
2.37
2:11 2: 0.002
0.60
2778-4 D:
3.40 0.550
2.22:6...0
1.97 2 0.47
2 0.57
255 4: 0.64
2:24
0.66
2.91 6.5.5
2.79 00.47 0.47
5.00 21 0.36
5.6 0.55
3.57 0.32
2.50 6.6.0
1.37
2,80 22 0.42
4.14 0.02
1.00 0.25
7:
0.56
0:777
3.00 1-0.3
4:15+00:5
3.54 0.52
5,084--000
2:20
7-
9.61 0.54
1.99
6,81 1.07
5.68 1.18
5.06
L.30
5.51 0.67
0,76 L.17
4.31 22 0.98
4.70 1.08
4,70
L.00
3.70-
4
5.10 6.1-4-0077
3.52 0:771
4.72.10.57
3,88 0.50
6.54-1 0.81
2,88 0.84
4.41 1.18
4.7% 1.36
7.00 =1 1,50
7.17 1.47
5,M1+4
3.81 1.0.6
5:35 0.36
1.35 1.000
1.00
2:74 0.67
4.72 d: 0.68
5,07 1,30
8.32 1.22
7.47- 3.30
7.02 1.16
3.54 0..12
4.69 0.93
0.35 2 0.30
6:25 0.3.1.00 0.62
5.94 0.79
Table 6: Ablation study: Performance for the tumor growth model with regular sampling times. Our SDIP-Net ablation has comparable performance to the state-of-the-art baselines for estimating CAPOs in discrete time.
We perform am ablation study, where we use our stabilized inverse propensity weights in the discrete- time setting. For this, observation mask. Instead, all timestamps are observed.
the tumor growth data in Section [ However, do not apply an
LSTM (Hochreiter & Schmidhuber 1997) the neural backhone in order to demonstrate that our approach is also applicable to other neural backbones. Our stabilized discrete time inverse propensity network (SDIP-Net) has comparable performance to the baselines. Table 6 shows the results, which thus confirm the effectiveness of our approach. Nevertheless, we empha- size that approach tailored for the continuous time setting, and not for the more unrealistic discrete-time time setting.
Here,
29

---

Published as a conference paper at ICLR 2025
G HYPERPARAMETER TUNING
In order to ensure a fair comparison of all methods, we close follow hyperparameter tuning as in [Melnychuk et al. 2022) and (Hess al. 2024a). In particular, we performed random grid search. Below, we report the tuning grid for each method. Importantly, all methods are only tuned on factual data. For optimization, we use Adam (Kingma & Ba 2013). Both TE-CDE (Seedat al. 2022) and our SCIP-Ne Net used a simple Euler quadrature and linear interpolation for the Neural CDE control path (Morrill et al. 202T). For the neural CDE of both methods and the integrated intensity in our SCIP-Net, we did not tune the grid size of the solver. Further, we emphasize that different quadrature schemes may impact training time. Fortunately, we did not encounter large differences in performance between schemes of different orders and, hence, opted for Euler integration.
Runtime: All methods were trained on 1x NVIDIA A100-PCIE-4OGB. On average, training our SCIP-Net took 29.4 minutes on tumor growth data and approximately 1.7 hours on MIMIC-III data, which is comparable to the baselines.
Coepern:
64. 136 256
Bmodit
(BC
apocie O4.)
56
Chaane de sic
bullin
Maatio apocia (4.1
50
64. 126
11,
Pegerad
Charne sec
Scrvet
[A51 8:311 BOKI
laley
nec
512 3D
LS533 degend
spocka
L5533 apc
256
Sound
Kilger cull 2 botc
Exoddt
Soal Carvel
Merthe spode (N1)
50
Neaal
sppooct andit (6-)
Montic spece (n.)
Norwa Bodge (102 oi L1)
A4 136 236
8:31 BOK1
Moai (1.8 ddlio
Eardil
Mar ual
14
L.3
3D+
Moad cop
) we let dyza = dy + de + d, be the overall input size. Further, d, is the hidden representation size of our SCIP-Net, and corresponds to the balanced representation size of TE-CDE (Seedat et al., 2022), CRN (Bica et al. 2020), and CT (Melnychuk et al. 2022). and the LSTM output size of G-Net (Li 2021).
Table 7: Following MMenycchk et al., 2022),
30


Reviews:
[{"review_id": "604bU5Wrnn", "reviewer": "Reviewer_MNJ4", "summary": "Models that can predict patient outcomes over time from EHR and insurance claims trajectories are increasingly important. This problem is technically non-trivial since patient events take place at irregular times and often have burst-like patterns around a major event. The key here is to model outcomes in continuous time and account for confounding that varies in time. The paper considers methods for adjusting time-varying confounding in continuous time, and proposes a new framework that permits doing so i.e. for irregular measurements. The presentation is heavy in formalism, but quite well-written. However, despite the formalism, the method by itself is fairly simple, at least in principle. The main problem as mentioned earlier is to predict patient outcome with interventions on both the treatment propensity and the treatment frequency given the history of the patient. The authors first write a separable interventional distribution (involving a geometric product integral) and clarify how the process also ensures identifiability. The general formulation for patient outcomes involves inverse propensity scoring, which is non-trivial to work with for various reasons, including computational tractability and instability. The paper first proposes a simple method to make the objective tractable and then formulate a scaling factor to weigh inverse propensity weights to improve stability. The architecture description looks quite complicated at first, but under the formalism is quite straightforward (although I do have concerns about how easy it is to train). The architecture is an encoder-decoder type network that involves learning the unstabilized inverse propensity measures with another network to obtain scaling weights, all of which are implemented using neural controlled differential equations. The experiments are rather brief, but in my opinion, sufficient to illustrate the efficacy of the method. \n\nOverall, I thought the paper was innovative and I enjoyed reading it.", "questions": "See the above filed under \"weakness\": I am curious about how easy is it to train the proposed model compared to the competition. I understand that the approach intends to improve stability, but the actual training exercise could possibly be quite involved given multiple moving parts.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "-- The paper is generally well-written and rigorous. It is also well-motivated and provides a solution of a problem of quite some importance not just in healthcare (the focus of this paper), but also applicable to other areas such as banking (where customer interactions are modeled in a similar manner irregularly, where interactions can also happen in bursts i.e. such as around a fraudulent transaction).\n-- The contribution, to the best of my knowledge is novel and interesting.", "weaknesses": "-- The paper obviously needs a lot of formalism, which is fine, but on the surface it seems like the network construction is essentially quite simple. But still, it is not clear how easy it is to train. It would be great to include more details on training details and difficulties, if any, in the appendix. On looking at the paper, one might think that the approach is complicated.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "crvnez2EsG", "reviewer": "Reviewer_utuj", "summary": "This work proposed the inverse propensity weight adjustment and stabilized inverse propensity weight adjustment to the training of continuous time modeling of counterfactual outcomes using neural controlled differential equations. The weights and stabilized weights are learned by separate neural networks, and used to adjust the target of outcome model, yielding a weighted MSE loss. The approach was evaluated on synthetic dataset based on tumor growth model and semi-synthetic dataset based on MIMIC-III, on showed noticeable improvement to compared prior work.", "questions": "1. Corresponding to Weaknesses 2., the authors should provide more robust argument or experimental evidence that balancing is biased and \"not proper\" compared to IPW adjustment, and change the column name \"adjustment for time-varying covariates\" in Table 1 to \"propensity weight adjustment for time-varying covariates\".\n\n2. Other than balancing and IPW adjustment, the more recent temporal-difference learning approach [3] adjust for time-varying covariates by constructing targets with weighted average of neural pesudo-outcomes. Such method should be discussed in related work. \n\n[3] Shirakawa, Toru, et al. \"Longitudinal Targeted Minimum Loss-based Estimation with Temporal-Difference Heterogeneous Transformer.\" International Conference on Machine Learning. PMLR, 2024.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "1. To the best of my knowledge, the use of inverse propensity weight to adjust the neural estimation of conditional average potential outcome in continuous time has not been proposed, so there is some non-trivial novel contribution in this work.\n\n2. The paper is well-written and easy to follow. It is quite notation-heavy but I'm not sure if there is much room for improvement.\n\n3. Good experimental results are shown on MIMIC-III and one-day ahead prediction of the tumor growth model.", "weaknesses": "1. Corresponding to Strengths 1., although the use of inverse propensity weight to adjust the neural estimation of conditional average potential outcome in continuous time has not been proposed, I do not see it as a significant contribution. It is a combination of TE-CDE and the line of work by Lok, Roysland and Rytgaard, i.e. A+B style contribution. Therefore, I think a score in the acceptance range is warranted, but not high acceptance.\n\n2. There are quite a few statements in this paper that I think are either inaccurate or in lack of evidence. In their statement of contributions (L100-107), they claim that balancing is a heuristic approach to adjust for time-varying confounding and was originally proposed for variance reduction, however, it is in fact first proposed in [1] for counterfactual inference. [2] was cited to state that balancing may increase bias, however, bias adjustment in [2] is for marginal parameter (averaging over covariates) and I don't think it is related to the conditional estimation in the scope of this work. Overall, I don't see any real argument or evidence that balancing is more biased than the proposed method and hence, I definitely don't think it is appropriate to say that it does not \"properly\" adjust for time-varying confounding. Besides, the parameter of interest is the conditional average potential outcome instead of potential outcome, and it should be stated correctly in Proposition 1 and the title.\n\n[1] Johansson, Fredrik, Uri Shalit, and David Sontag. \"Learning representations for counterfactual inference.\" International conference on machine learning. PMLR, 2016.\n\n[2] Valentyn Melnychuk, Dennis Frauen, and Stefan Feuerriegel. Normalizing flows for interventional density estimation. In ICML, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZeGxJYo1iL", "reviewer": "Reviewer_bgav", "summary": "This study proposed SCIP-Net, a model designed to predict potential outcomes in continuous time. SCIP-Net adopts controlled differential equations (CDEs) as its architecture and primarily relies on inverse propensity weighting with stabilizing weights to estimate potential outcomes.", "questions": "Figure 3 and Table 2 have low readability due to the small fonts and figure legends. The authors might consider improving the readability of them.\n\nIn Figure 3, the performance of TE-CDE is very different from the original publication since TE-CDE was outperformed by RMSN and CRN. It is different from the experiment on MIMIC data. Can authors provide reasonable guesses?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper's contributions are practically significant, as detailed in the introduction. \n\nThe authors provide a thorough explanation of the problem formulation and model architecture. \n\nThe paper is well-structured with (1) visual emphasis on key points and (2) a coherent organization of sections.", "weaknesses": "The numerical experiment section needs further details, particularly about the experiment setup (e.g., how the data are randomly sampled). I see that some of the details were referenced in the previous study (Vanderschueren et al.) and supplementary (perhaps most likely due to the page limit). However, considering different sampling strategies from Vanderschueren et al. and novelties compared to the previous studies that aimed to POs in discrete time, I believe it will be worthwhile to put more details about the experiment setup.\n\nI believe the length of the forecast horizon (i.e., prediction window) affects the irregularity of the timestamps. Looking at the results when the forecast horizon is small (i.e., one-day prediction in tumor growth and 1 hour in MIMIC), SCIP and CIP did not show meaningful differences, which may be due to not significant irregularity of timestamps since the prediction horizon is small. Therefore, it will be very informative how irregular the timestamps (e.g., distributions of the length of the timestamps) are to analyze the performance and benefits of the proposed method. \n\nThe ablation studies that compared CIP-Net and SCIP-Net are appropriate to show the benefits of using stabilized weights. Related to the comment above, it will be also informative if authors can show the performance of the proposed methods and baselines on regularly sampled timestamps.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1S6yNDFYgy", "reviewer": "Reviewer_mKFR", "summary": "This paper studies the modeling of EHR as continuous time series. In particular, it focuses on modeling the time-varying confounding as well. This paper proposes a stabilized target to fit and a neural network architecture, and validated the proposed method on two synthetic datasets.", "questions": "1. How sensitive is SCIP-Net’s performance to the specific confounding strength? If confounding is minimal (and the data abundant), would SCIP-Net still provide benefits over simpler models? I saw that Figure 3 is related to this, but it's hard for me to make association between $\\gamma$ and more realistic dataset. \n\n2. Can the authors elaborate how the baselines were chosen?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The problem studied in this paper is important. Being able to perform \"what if\" analysis is critical in medical (and other) times series modeling tasks.", "weaknesses": "1. Experiment: Fundamentally I found it a bit difficult to conclude that the proposed method is effective. There are a lot of literature on similar topics (they may use the term \"irregular time series\" instead of \"contiuous\"). Although they do not model the confounding variables, in reality it is difficult to  know if such additional modeling helps or hurts. The paper now uses simulated data to show the point, but simulated data is fundamentally not reliable. One possible way to show that modeling confounding actually helps might be tweaking the distribution by oversampling certain interventions. If the proposed modeling actually helps, I'd imagine it not noly improves metrics measured on the original marginal distribution, but should also be better on the re-sampled distribution. \n\n\ta. To be very clear, if the goal is to predict the effect of drug A on symptom B, and the proposed method actually models 1. when A will be applied, given EHR history x 2. result of A on B, then I'd expected that it might not always beat something that learns the marginal distribution B|x, but if we focus on only B|A=1 or B|A=0 (using real, unsimulated data), the covariate shift should have little to no impact on the proposed method but greatly hurts other baselines. \n\tb. Please also include more literature in the \"irregular time series\" modeling space.\n\n2. Presentation and writing bothers me a bit. While I appreciate the use of coloring to distinguish certain quantities, things like \"=>Takeaway\" in L509 makes the manuscript a bit too casual in my opinion. Minor issues on some of the expressions as well (such as \"unconditional of\" at L272). \n\t\ta. Also related to 1, the experiment section is very unclear right now. The paper only says \"RMSE\" is the evaluation metric, but \"RMSE of what\"? This is particular problematic given that the dataset itself is synthetic. I'd suggest moving some details from the Appendix to the main paper and but very clear about the setup of the experiment (and the motivation behind it).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)