Estimating counterfactual treatment outcomes over time through adversarially balanced representations
Authors: Ioana Bica
Ahmed M Alaa
James Jordon
Mihaela van der Schaar
Venue:  (2020)
Keywords: treatment effects over time
causal inference
counterfactual estimation

Abstract: Identifying when to give treatments to patients and how to select among multiple treatments over time are important medical problems with a few existing solutions. In this paper, we introduce the Counterfactual Recurrent Network (CRN), a novel sequence-to-sequence model that leverages the increasingly available patient observational data to estimate treatment effects over time and answer such medical questions. To handle the bias from time-varying confounders, covariates affecting the treatment assignment policy in the observational data, CRN uses domain adversarial training to build balancing representations of the patient history. At each timestep, CRN constructs a treatment invariant representation which removes the association between patient history and treatment assignments and thus can be reliably used for making counterfactual predictions. On a simulated model of tumour growth, with varying degree of time-dependent confounding, we show how our model achieves lower error in estimating counterfactuals and in choosing the correct treatment and timing of treatment than current state-of-the-art methods.

Published as a conference paper at ICLR 2020
ESTIMATING COUNTERFACTUAL TREATMENT OUTCOMES OVER TIME THROUGH ADVERSARIALLY BALANCED REPRESENTATIONS
Ioana Bica
Ahmed M. Alaa
Department of Engineering Science University of Oxford, Oxford, UK The Alan Turing Institute, London, UK ioana. bi caleng ox ac.uk
Department of Electrical Engineering
University California, Los Angeles, USA
eduu
James Jordon
Mihaela van der Schaar
Department of Engineering Science University of Oxford, Oxford, UK james. jerdon@wolfsee ox ac
University of Cambridge, Cambridge, UK University of California, Los Angeles, USA The Alan Turing Institute, London, UK mv 4 720 cam ac.uk
ABSTRACT
Identifying when to give treatments to patients and how to select among multiple treatments over time are important medical problems with few existing solutions. In this paper, we introduce the Counterfactual Recurrent Network (CRN), a novel sequence-to-sern -esuueecc model that leverages the increasingly available patient ob- servational data to estimate treatment effects over time and answer such medical questions. To handle the bias from varying confounders, covariates affect- ing the treatment assignment policy in the observational data, CRN uses domain adversarial training to build balancing representations each timestep, CRN constructs treatment invariant representation which removes the association between patient history and treatment assignments be reliably used for making counterfactual predictions. On simulated model of tumour growth, with varying degree of time-dependent confounding, we show how model achieves lower error in estimating counterfactuals and in choosing the correct treatment and timing of treatment than current state-of-the-arr methods.
the patient history. At
thus can
INTRODUCTION
As clinical decision-makers are often faced with the problem of choosing between treatment alterna- tives for patients, reliably estimating their effects paramount. While clinical trials represent the gold standard for causal inference, they are expensive, have a few patients and narrow inclusion criteria (Booth Tannock (2014). Leveraging the increasingly available observational data about patients, such as electronic health records, represents more viable alternative for estimating treatment effects.
large number of methods have been proposed for performing causal inference using observational data in the static setting (Johansson et al. 2016 Shalit et al. 2017 Alaa & van der Schaar 2017 Li & Fu 2017 Yoon [ Alaa & van der Schaar 2018 Yao 2018) and only a few methods address the longitudinal setting (Xu et al. 2016 Roy et al. 2016 Soleimani et al.) 201
Schulam & Saria 2017 Lim et al. 2018). However, estimating the effects treatments over time poses unique opportunities such understanding how diseases evolve under different treatment plans, how individual patients respond to medication over time, but also which are optimal timings for assigning treatments, thus providing new tools to improve clinical decision support systems.
The biggest challenge when estimating the effects of time- demendeenddnt treatments from observational data involves correctly handling the time-dependent confounders: patient covariates that affected by past treatments which then influence future treatments and outcomes (Platt et al., 2009). For instance, consider that treatment A is given when certain patient covariate (e.g. white blood cell

---

Published as a conference paper at ICLR 2020
Chemotherapy
Radiodiengy
Patient history
Coumeriictaal ouscemes
Patient history
Counerbactual
Pattent histrry IT
Countrricttall outsemes
Tumour volume
Tumne volamo
T volumee
Rost
Beat oakoumec
Beet cutzom
Poi
Produtions
Prolitiee
time
Time
Current
time
Time
time
Currest
Time
(a) Decide treatment plan
(b) Decide optimal time treatment
(c) Decide when stop treatment
Figure 1: Applicability of CRN in cancer treatment planning. We illustrate 3 patients with different covariate and treatment histories H, For current time I, CRN can predict counterfactual trajectories (the coloured dashed branches) for planned treatments in the future. Through the counterfactual predictions, we can decide which treatment plan results the best patient outcome (in this case, the lowest tumour volume). This way, CRN can be used to perform all of the following: choose optimal treatments (a), find timing when treatment is most effective (b) decide when to stop treatment (c).
count) has been outside of normal range values for several consecutive timesteps. Suppose also that this patient covariate was itself affected by the past administration of treatment B. If these patients more likely die, without adjusting for the time-dependent confounding (e.g. the changes in the white blood cell count over time), we will incorrectly conclude that treatment is harmful to patients. Moreover, estimating the effect of different sequence of treatments on the patient outcome would require not only adjusting for the bias at the current step (in treatment A), but also for the bias introduced by the previous application of treatment
Existing methods for causal inference in the static setting cannot be applied in this longitudinal setting since they are designed to handle the cross-sectional set-up, where the treatment and outcome depend only static value of the patient covariates. If we consider again the above example, these methods would not be able to model how the changes in patient covariates over time affect the assignment of treatments and they would also not be able to estimate the effect of a sequence of treatments on the patient outcome (e.g. sequential application of treatment followed by treatment B). Different models that can handle these temporal dependencies in the observational data and varying-length patient histories are needed for estimating treatment effects over time.
dependenn confounders present in observational data because doctors follow policies: the history of the patients' covariates and the patients' response to past treatments are used to decide future treatments (Mansournia et al., 2012). The direct use of supervised learning methods will be biased by the treatment policies present in the observational data and will not be able to correctly estimate counterfactuals for different treatment assignment policies.
Standard methods for adjusting for time-varying confounding and estimating the effects of time- varying exposures are based on ideas from epidemiology. The most widely used among these Marginal Structural Models (MSMs) (Robins et al., 20000 Mansournia et al. (0122) which use the inverse probability of treatment weighting (IPTW) adjust for the time-dependent confounding bias. Through IPTW, MSMs create a ppcuuo-population where the probability of treatment does not depend on the time varying confounders. However, MSMs are not robust to model misspecification in computing the IPTWs. MSMs can also give high-variance estimates due to extreme weights; computing the IPTW involves dividing by probability of assigning a treatment conditional on patient history which can be numerically unstable if the probability is small.
We introduce the Counterfactual Recurrent Network (CRN), a novel sequence-to-sequeee archi- tecture for estimating treatment effects over time. CRN leverages recent advances in representation learning (Bengio et al. 2012) and domain adversarial training (Ganin et al., 2016) to overcome the follows.
problems of existing methods for causal inference
time. Our main contributions
Treatment invariant representations over time. CRN constructs treatment invariant represen- tations each timestep in order to break the association between patient history and treatment assignment and thus removes the bias from time- tine-dependent confounders. For this, CRN uses domain adversarial training (Ganin al. 2016 Li et al. 2018. Sehag et al. 2019) to trade-off between build-
2

---

Published as a conference paper at ICLR 2020
ing this balancing representation and predicting patient outcomes. We show that these representations remove the bias from Line-varying confounders and can be reliably used for estimating counterfactual outcomes. This represents the first work that introduces ideas from domain adaptation to the area of estimating treatment effects over time. In addition, by building balancing representations, we propose novel way of removing the bias introduced by time-varying confounders.
Counterfactual estimation of future outcomes. To estimate counterfactual outcomes for treatment plans (and not just single treatments), we integrate the domain adversarial training procedure as part of a sequence-to-sequeee invariant representations of the patient history that are used initialize the decoder. The decoder network estimates outcomes under an intended sequence of future treatments, while also updating the balanced representation. By performing counterfactual estimation of future treatment outcomes, CRN can be used to answer critical medical questions such as deciding when to give treatments to patients, when to start and stop treatment regimes, and also how to select from multiple treatments over time. the applicability of our method in choosing optimal cancer treatments.
architecture. CRN consists of an encoder network which builds treatment
We illustrate in Figure I
In our experiments, we evaluate CRN in a realistic set-up using a model of tumour growth (Geng et al. 2017). We show that CRN achieves better performance in predicting counterfactual outcomes, but also in choosing the right treatment and timing of treatment than current state-of-the-an methods.
2 RELATED WORK
We focus methods for estimating treatment effects over time and for building balancing represen- tations for causal inference. A more in-depth review of related work is in Appendix A
Treatment effects over time. Standard methods for estimating the effects of time-varying exposures were first developed in the epidemiology literature and include the g- g-computation formula, Structural
Nested Models and Marginal Structural Models (MSMs) (Robins 1986 1994 RRobins al. 2000 Robins & Hernán 2008). Originally, these methods have used predictors performing logistic/line linear regression which makes them unsuitable for handling complex dime-dependencies (Hernán al.) Mansourna et al. 2012, Mortimer et al. 2005). To address these limitations, methods that use Bayesian non-par non-parametries or recurrent neural networks as part of these frameworks have been proposed. (Xu et al. 2016 Roy et al.) 2016) Lim et al. 2018).
To begin with, Xu et al. (2016) use Gaussian processes to model discrete patient outcomes as generalized mixed-effects model and uses the g-computation method to handle time-varying confounders. [Soleimani et al.] (2017) extend the approach Xu et al. (2016) to the continuous time-setting and model treatment responses using linear lime-invariant dynamical systems. |Roy et al. (0116) use Dirichlet and Gaussian processes to model the observational data and estimate the IPTW in Marginal Structural Models. Schulam & Saria (017)) build upon work from Lok et al. 2008); Arjas & Parner (2004) and use marked point processes and Gaussian processes to learn causal effects in continuous time data. These Bayesian non- non-parametric methods make strong assumptions about model structure and consequently cannot handle well heterogeneous treatment effects arising from baseline variables (Soleimani et al.
& Saria 2017) and multiple treatment outcomes
(Xu et al. 2016 Schulam & Saria 2017).
The work most related to ours the one Lim et al. (2018) which improves on the standard MSMs by using recurrent neural networks to estimate the inverse probability of treatment weights (s..Ws). Lim et al. (2088) introduces Recurrent Marginal Structural Networks (RMSNs) which also use sequence-tee posennee-t--eqqeee seqqueecc deep learning architecture to forecast treatment responses in a similar fashion to our model. However, RMSNs require training additional RNNs to estimate the propensity weights and does not overcome the fundamental problems with IPTWs, such as the high-variance of the weights. Conversely, CRN takes advantage of the recent advances in machine learning, in particular, representation learning to propose novel way of handling time-varying confounders.
Balancing representations for treatment effect estimation. Balancing the distribution of control and treated groups has been used for counterfactual estimation in the static setting. The methods proposed in the static setting for balancing representations are based on using discrepancy measures in the representation space between treated and untreated patients, which do not generalize to multiple treatments (Johansoon et al. 2016;; (halit et al. 2017] [i & Ful 2017)) Yao et al.] 2018). Moreover, due to the sequential assignment of treatments in the longitudinal setting, and due to the change of
3

---

Published as a conference paper at ICLR 2020
patient covariates over time according to previous treatments, the methods for the static setting not directly applicable to the varring setting (Hernán et al. |000 Mansournia et al. 2012).
3 PROBLEM FORMULATION
Consider an observational dataset D = { [x] (1) a 1 17th ((v(")) consisting of information
about . independent patients. For each patient (i), we observe time-dependent covariates x!") E Xx treatment received A," E {AI,... Ax} } = A and outcomes Y The patient can also have baseline covariates V E V such as gender and genetic information. Note that the outcome y(f)
++1 V4+1 for T discrete timesteps.
++1 will be part of the observed covariates 2+1 For simplicity, the patient superscript (i) will be omitted unless explicitly needed.
We adopt the potential outcomes framework proposed by (Neyman 1923 Rubin 1978) and ex- tended by (Robins & Hernán 2008) account for time-varying treatments. Let Y a be the potential outcomes, either factual or counterfactual, for each possible course of treatment a. Let H. = (X,, A,- 111 V) represent the history of the patient covariates X, (X,,...,X, treatment assignments A, (Ai,..., Ar) and static features V. We want to estimate:
B B+Y++r (a + H)),
(I)
where a(t, a(t,t+T- 1) = |ay.................... a represents a possible sequence of treatments from timestep t just until before the potential outcome Y, ,+ is observed. We make the standard assumptions (Robins)
et al. 2000 Lim et al. (018) needed to identify the treatment effects: consistency, positivity and no hidden confounders (sequential strong ignorability). See Appendix B for more more details.
COUNTERFACTUAL RECURRENT NETWORK
The observational data can be used to train supervised learning model to forecast: (Y 1+-7 Att, + 1) a(t, t + 1), H,). However, without adjusting for the bias introduced by time- varying confounders, this model cannot be reliably used for making causal predictions (Robins et al.. Robins Robinss & Hernán 2008: [Schulam & Saria 2017). The Counterfactual Recurrent Network (CRN) removes this bias through domain adversarial training and estimates the counterfactual outcomes (YY E(Y++(((aa))))) (a(,((( H)), for any intended future treatment assignment ál,, 1).
Balancing representations. The history H. (X, A,-1, V) of the patient contains the time- varying confounders X, which bias the treatment assignment A, (A1.... AA) in the observational dataset. Inverse probability of treatment weighting, as performed by MSMs, creates a pseudo- population where the probability of treatment A, does not depend on the time-varying confounders (Robins et al. 2000). In this paper, we propose instead building a representation that is not predictive of the treatment A,. This way, we remove the association between history, containing the time-varying confounders X,, and current treatment A,, Robins (199) shows that in this case, the estimation of counterfactual treatment outcomes is unbiased. See Appendix Cor details and for an example of causal graph with time-dependent confounders.
the history H,
Let & be the representation function that maps the patient history H. a representation space R. To obtain unbiased treatment effects, & needs to construct treatment invariant representations such that (H)(H. A, A1)
P o (H.) A, Ax). To achieve this and to estimate counterfactual outcomes under planned sequence of treatments, we integrate the domain adversarial training framework proposed by Ganin et al. (2016) and extended by [Sebag et al.) (2019) to the multi-domain learning setting, into a sequence-to-sequence to-seqeence architecture. In our case, the different treatments at each timestep are considered the different domains. Note that the novelty here comes from the use of domain adversarial training to handle the bias from the time-dependent confounders, models, which have already been applied to forecast
rather than the use of peruuccce-to-scqueee treatment responses (Lim et al. 2018). Figure 2 illustrates model architecture.
Encoder. The encoder network uses an RNN, with LSTM unit (Hochreiter Schmidhuberi 1997). to process the history of treatments A,- covariates X, and baseline features V to build a treatment invariant representation & (H,). ). but also to predict mee-tee-ahead outcomes Y,41- To achieve this, the encoder network aims maximize the loss of the treatment classifier G, and minimize the loss
4

---

Published as a conference paper at ICLR 2020
Gredias sevenad laye
A
Y=41
A--1
V Ye42
A148
Yus
A
G
G
Gy
G
Gy
G
G
2
GBL
CBL
c H
O (++1)
d H
+ H
Patienl bistoy H,
RNN
RNN
RNN
(H.)
RNN
RNN
RNN
V X
V
X
V
Yoek
a.41 V You
BEiT-
Encoder building the balancing representation
Decoder for sequence prediction treaiments effects
Figure 2: CRN architecture. Encoder builds representation (H,,) ) that maximizes loss of treatment classifier G. and minimizes loss of outcome predictor Gy. ((H)) is used initialize the decoder, which continues to update it to predict counterfactual outcomes of a sequence of future treatments.
of the outcome predictor network Gy- This way, the balanced representation + (H,)) )) is not predictive of the assigned treatment A., but is discriminative enough estimate the outcome Y To train this model using gradient descent, we use the Gradient Reversal Layer (Ganin et al. 2016).
Decoder. The decoder network uses the balanced representation computed by the encoder to initialize the state of an RNN that predicts the counterfactual outcomes for a sequence of future treatments. Dur- ing training, the decoder uses input the outcomes from the observational data (Y++1.... Y Y++r-1). the static patient features V and the intended sequence of treatments ä(t, L +++ 7 - 1). The decoder is trained in similar way the encoder update the balanced representation and estimate the outcomes. During testing, we do not have access to ground-trul -rrth outcomes; thus, the outcomes predicted by the decoder (Y (Y....... 11) are auto-regressively used instead as inputs. By running the decoder with different treatment settings, and by auto-regressively feeding back the outcomes, we can determine when to start and end different treatments, which is the optimal time to give the treatment and which treatments to give over time to obtain the best patient outcomes.
The representation ((H,, ( is built by applying a fully connected layer, with Exponential Linear Unit (ELU) activation to the output of the LSTM. The treatment classifier G. and the predictor network Gy consist of a hidden layer each, also with ELU activation. The output layer of G. uses softmax activation, while the output layer of Gy uses linear activation for continuous predictions. For categorical outcomes, softmax activation can be used. We follow an approach similar to Lim et al. (2018) and we split the encoder and decoder training into separate steps. See Appendix for details.
The encoder and decoder networks use variational dropout (Gal & Ghahramani 2016) such that the CRN can also give uncertainty intervals for the treatment outcomes. This is particularity important in the estimation of treatment effects, since the model predictions should only be used when they have high confidence. Our model can also be modified to allow for irregular samplings of observations by using a PhasedLSTM (Neil al. 2016).
5 ADVERSARIALLY BALANCED REPRESENTATION OVER TIME
At each timestep I, let the K different possible treatments A, E (A ... Ax) represent our domains. As described in Section 4 to remove the bias from deeendependent confounders, we build represen-
tation of history H. that invariant across treatments: (()(H.) A1) = P @ (H)(H.) AA))
This requirement can be enforced by minimizing the distance in the distribution of (H,) between any two pairs of treatments. Kifer et al. (2004); David al. (2007), propose measuring the disparity between distributions based on their separability by discriminatively--raiin Let the symmetric hypothesis class H consist of the set of symmetric multiclass classifiers, such as neural network architectures. The H-divergencee between all pairs of two distributions is defined in terms of the capacity of the hypothesis class H. to discriminate between examples from the multiple distributions. Empirically, minimizing the H divieeegeence involves building representation where examples from the multiple domains are as indistinguishable as possible (Ben- David et al., 2207 let al. 2018 Sebag et al. 2019). Ganin et al. (016) use this idea propose an adversarial framework
classifier.
5

---

Published as a conference paper at ICLR 2020
for domain adaptation involving building representation which achieves maximum enror domain classifier and minimum eπror on an outcome predictor. Similarly, in our case, we use domain adversarial training to build representation of the patient history ((H,) that is both invariant the treatment given at timestep I, A, and that achieves low error in estimating the outcome Y (41-
Let ( (a)) be the treatment classifier with param- eters 8, and let G ( (H,) fa) be the output correspond- ing treatment Aj. Let G (H,,: 0yy) be the predictor network with parameters Oy, The representation func- tion parameterized by the parameters 0, in the RNN:
Forward propagation
Backpropsgatioon
A
L
L
1+11
& (He; 0,). Figure 3) shows the adversarial training proce- dure used.
A
Y
(+1
For timestep t and patient (i), let L (Or, , 0a) be the treat- ment (domain) loss and let clyy (Or, (Or, 02) the outcome loss, defined as follows:
DC.
acg
00
G
G
00,
Ceatiestroveraa
GRI
A
S
(H,)
(,,, (0,, 0a)
E j=1
I (af0 =a,) log log(
OC,
0C,
Baknnol
80,
70.
(2)
h
h:
RNN
(0,,,) (G,
(G,(
(3)
A
If the outcome is binary, the cross-entropy loss be used instead for Lt.g To build treatment invariant repre- sentations and to also estimate patient outcomes, aim
Figure 3: Training procedure for build- ing balancing representation.
to maximize treatment loss and minimize outcome loss.
Thus, the overall loss L.y at timestep is given by:
c (0, 0.)
where the hyperparameter À controls this trade-off between domain discrimination and outcome prediction. We use the standard procedure for training domain adversarial networks from Ganin (et al. (2016) and we start off with an initial value for À and use an exponentially increasing schedule during training. To train the model using backpropagation, (Ganin et al. 2016). For more details about the training procedure, see Appendix
the Gradient Reversal Layer (GRL)
By using the objective (,,( (0,, 0y, 0a), we reach the saddle point (ō., 0,, 0o) that achieves the equilibrium between domain di discrimination and outcome estimation.
(0r: Q,) = arg min c.0 (0r, 0y, 9,,9,
0, = arg max c.0h (Ō,, Oa).
(5)
6n
The result stated in Theorem proves that the treatment (domain) loss part of our objective (from equation 2 aims to remove the time-dependent confounding bias.
Theorem 1. Let E .1,2,...}. For each j = 1,.., K, let Pj denote the distribution of He conditional on A, = Aj and let ph denote the distribution of (H,)) conditional on A =A. Let
G denote the output of Ga corresponding treatment Aj. Then the minimax game defined by
K EAA,P, ([
K
min max G.
subject to
(o),) =1
(6)
==1
==1
has global minimum which is atlained if and only if Pf representations are invariant across all treatments.
p
pt i.e. when the learned
Proof. This result is restatement of the one in Li al. 2018 For details, see the Appendix D n

---

Published as a conference paper at ICLR 2020
A good representation allows us obtain low enror in estimating counterfactuals for all treatments, while at the same time to minimize the H-diiergencee between induced marginal distributions of all the domains. We use algorithm that directly minimizes combination of the H divvieegee and the empirical training margin.
6 EXPERIMENTS
In real datasets, counterfactual outcomes and the degree of time-dependent confounding are not known (Schulam Saria 2017 Lim et al. 2018). To validate the CRN" we evaluate on Pharmacokinetic- Sacikiinttii-Paaaaacoooyye Pharmacoddnnaii model of tumour growth (Geng al. 2017), which uses state-
of-the-art bio-mathematical model to simulate the combined effects of chemotherapy and radiotherapy
in lung cancer patients. The same model was used by Lim et al. (2018) to evaluate RMSNs.
Model of tumour growth The volume of tumour days after diagnosis is modelled as follows:
K
V(l 1)
plog
B,C(L)
(a,, + +,)(1))
+ Et )(())
V(t)
Chemotherapy
Radiotherapy
Noise
Tumar growth
where K, Ps Be, 3r are sampled as described in Geng et al. (2017). To incorporate heterogeneity in patient responses, the prior means for 3. and ax, adjusted to create patient subgroups, which are used as baseline features. The chemotherapy concentration C'(t) and radiotherapy dose d(t)
are modelled as described in Appendix (F) Time-varying confounding is introduced by modelling chemotherapy and radiotherapy assignment as Bernoulli random variables, with probabilities pe and Dmux
pr- depending on the tumour diameter: (t) = a "o (D(t) 5z)) and (t = a (22 (D(t)
6,)) where D(t) is the average diameter over the last 15 days, D = 13cm, af-) is the sigmoid and 6, = 6, = D mmxx/ 2. The amount of time-dependent confounding is controlled through Ye. 7ri the higher 7. is, the more important the history is in assigning treatments. At each timestep, there are four treatment options: no treatment, chemotherapy. radiotherapy, combined chemotherapy and radiotherapy. For details about data simulation, see Appendix [F]
Benchmarks We used the following benchmarks for performance comparison: Marginal Structural Models (MSMs) (Robins et al. 2000,, which use logistic regression for estimating the IPTWs and linear regression for prediction (see Appendix for details). We also compare against the Recurrent al. (2018), which is the current mfateoof-the-ar model in estimating treatment responses. RMSNs use RNNs to estimate the IPTWss and the patient outcomes (details in Appendix H.. To show that standard supervised learning models do not handle the time- varying confounders
Marginal Structural Networks (RMSNs)[Lim
compare against an RNN and a linear regression model, which receive as imput treatments and covariates to predict the outcome (see Appendix [ for details). Our model architecture follows the description in Sections and 5 with full training details and hyperparameter optimization in Appendix[] [ To show the importance of adversarial training, we also benchmark against CRN (A 0) model with the same architecture, but with À 0, i.e model architecture without adversarial training.
6.1 EVALUATE MODELS ON COUNTERFACTUAL PREDICTIONS
Previous methods focused on evaluating the error only for factual outcomes (observed patient outcomes) (Lim et al. 2018). However, build decision support systems, we need evaluate how well the models estimate the counterfactual outcomes, i.e patient outcomes under alternative treatment options. The parameters "e and 7/ control the treatment assigmment policy, i.e. the degree of time-dependent confounding present in the data. We evaluate the benchmarks under different degrees of time-dependent confounding by setting 7 = To = Tr. For each 7 we simulate 10000 patients for training, 1000 for validation (hyperparameter tuning) and 1000 for out-of-sample testing. For the patients in the test set, for each time t, we also simulate counterfactuals (1+1, represented by tumour volume (tt + 1), under all possible treatment options.
'The implementation of the model can be found https://bitbuukket wwtp://etttto//oooeooeeeeooo m1 for htt.ftteiiiaaapp/e//e////////see http:s///aaaacccoooooooooooooaan
/ttp///////////boooooooooo/oo//a /
and

---

Published as a conference paper at ICLR 2020
5%
5%
Lincar
RMN CRN
RMSN
RMEN CRN
Normalized RMSE
4%
No time dependent
3%
confounding
3%
2%
Large dime-dependent confounding
1%
17
9
10
Degree time- depeededddd confunding 7 (a) One-step ahead prediction
Degree dependend confunding
(b) Five-step ahead prediction
Figure 4: Results for prediction of patient counterfactuals.
Figure (a) shows the normalized root mean squared error (RMSE) for one-step ahead estimation of counterfactuals with varying degree of time-dependent confounding 7. The RMSE normalized by the maximum tumour volume: Vaaz = 11 150cm3. The linear and MSM models provide baseline for performance they achieve the highest RMSE. While the use of IPTW in MSMs helps when y increases, using linear modelling has severe limitations. When there is time-dependent confounding, the machine learning methods achieve similar performance, close to 0.6% RMSE, As the bias in the dataset increases, the harder it becomes for the RNN and the CRN (A 0) to generalize to estimate outcomes of treatments not matching the training policy. When y = 10, CRN improves by 48. 18.1% the same model architecture without domain adversarial training CRN (A 0).
Our proposed model achieves the lowest RMSE across all values of y. Compared to RMSNs, CRN improves by 17% when y >> 6. To highlight the gains of our method for smaller 7, Figure 4 (b) shows the RMSE for five-step ahead prediction (with counterfactuals generated as described in Section 6.2 and Appendix [). RMSNs also use decoder for sequence prediction. However, RMSNs require training additional RNNs to estimate the IPTW, which are used to weight each sample during the decoder training. For 7-step ahead prediction, IPTW involves multiplying weights which can result in high variance. The results in Figure 4 (b) show the problems with using IPTW to handle the denendennn confounding bias. See Appendix forr more results on multi-step ahead prediction.
Balancing representation: To evaluate whether the CRN has indeed learnt treatment invariant represenations, for ? 5, we illustrate in Figure the T-SNE embeddings of the balancing rep- resentations (HH) built by the CRN encoder for test patients. We color each point by the treat- ment A, E {no treatment, chemotherapy, radiotherapy, combined chemotherapy and radiotherapy received at timestep to highlight the invariance of ((H,) across the different treatments. In Figure b),, we show &b (H,) only for chemotherapy and radiotherapy for better understanding.
10
20
Radethenpy
-10
Cheme i Ride
tixe
Figure 5: TSNE embedding of the balancing representation (HH, learnt by the CRN encoder different timesteps t. Notice that (H,,) is not predictive of the treatment A, given at timestep t.
6.2 EVALUATE RECOMMENDING THE RIGHT TREATMENT AND TIMING OF TREATMENT
Evaluating the models just in terms of the RMSE on counterfactual estimation is also not enough for assessing their reliability when used as part of decision support systems. In this section we assess how well the models can select the correct treatment and timing of treatment for several forecasting horizons т. We generate test sets consisting of 1000 patients where for each horizon 7 and for each time in a patient's trajectory, there are options for giving chemotherapy one of f,.., and 7 options for giving radiotherapy at one of ff....t treatment is applied. These 27 treatment plans assessed in terms of the tumour volume outcome Y,+r- We select the treatment (chemotherapy or radiotherapy) that achieves lowest Y Y,tr. and within the correct treatment the timing with lowest Y e+.- We also compute the normalized RMSE for
1. At the rest of the future timesteps, no
8

---

Published as a conference paper at ICLR 2020
Table I: Results for recommending the correct treatment and timing of treatment.
%
5, 7r- =
5
Vo = 5, 7yr
Ve- = 0, Yr =
7
CRN
RMSN
MSM
CRN
RMSN
MSM
CRN
RMSN
MSM
Normalized RMSE
2.43% %
3.16%
6.75%
1.08%
1.35%
3.68%
1.54%
1.59%
3.23%
2.83%
3.95%
7.65%
1.21 %
1.81%
3.84%
1.81
2.25%
3.52%
3.18%
4.37%
7.95%
1.33%
2.13%
3.91%
2.03%
2.71%
3.63%
3.51 %
5.61%
8.19%
1.42%
2.41%
3.97%
2.23%
2.73%
3.71%
7 3.93%
6.21%
8.52%
1.53%
2.43%
4.04%
2.43%
2.88%
3.79%
Treatment Accuracy
3 83.1%
75.3%
73.9%
83.2
78.6%
77.1%
92.9%
87.3%
74.9%
82.5%
74.1%
68.5%
81.3%
77.7%
73.9%
85.7%
83.8%
74.1%
73.5%
72.7%
63.2%
78.3%
77.2%
72.3%
83.8%
82.1%
72.8%
69.4%
66.7%
62.7%
79.5%
76.3%
71.8%
78.6%
69.7%
64.5%
7 71.2%
68.8%
62.4%
72.7%
71.8%
71.6%
71.9%
69.3%
61.2%
Treatment
79.6%
78.1%
67.6%
80.5%
76.8%
77.5%
79.8%
75.7%
60.6%
Timing Accuracy
4 73.9%
70.3%
63.1%
79.0% %
77.2%
73.4%
75.4%
71.4%
58.2%
69.8%
68.6%
62.4%
78.3%
73.3%
63.6%
66.9%
31.3%
29.5%
6 66.9%
66.2%
62.6%
73.5%
72.1%
63.9%
65.8%
24.2%
15.5%
7
64.5%
63.6%
62.2%
70.6%
57.4%
44.2%
63.9%
25.6%
12.5%
predicting Y(... See Appendix L for more details about the test set. The models are evaluated for
settings % and ?/-
Table shows the results for this evaluation set-up. The treatment accuracy denotes the percentage of patients for which the correct treatment selected, while the treatment timing accuracy the percentage for which the correct timing was selected. Note that when %o = 0 and 7 5, RMSN 4. CRN performs similarly
and MSM select the wrong treatment timing for projection horizons among the different policies present in the observational data and achieve the lowest RMSE and highest accuracy in selecting the correct treatment and timing of treatment.
In Appendix M we also show the applicability of the CRN in more complex medical scenarios involving real data. We provide experimental results based on the Medical Information Mart for Intensive Care (MIMIC III) database (Johnson et al. 2016) consisting of electronic health records from patients the ICU.
CONCLUSION
Despite its wide applicability, the problem of causal inference for time-dependent treatments has been relatively less studied compared problem of causal inference in the static setting, Both new methods and theory are necessary to be able to harness the full potential of observational data for learning individualized effects of complex treatment scenarios. Further work in this direction is needed for proposing alternative methods for handling time-dependent confounders, for modelling combinations of treatments assigned treatments with associated dosage.
time or for estimating the individualized effects of time-dependent
In this paper, introduced the Counterfactual Recurrent Network (CRN), a model that estimates individualized effects of treatments over time using a novel way of handling the bias from time- dependent confounders through adversarial training. Using model of tumour growth, we validated CRN in realistic medical scenarios and we showed improvements over existing state-of state-of-the-art arr methods. We also showed the applicability of the CRN a real dataset consiting of patient electronic health records. The counterfactual predictions of CRN have the potential to be used as part of clinical decision support systems to address relevant medical challenges involving selecting the best treatments for patients over time, identify optimal treatment timings but also when the treatment is no longer needed. In future work, we will aim to build better balancing representations and provide theoretical guarantees for the expected error the counterfactuals.

---

Published as a conference paper at ICLR 2020
ACKNOWLEDGMENTS
We would like to thank the reviewers for their valuable feedback. The research presented in this paper was supported by The Alan Turing Institute, under the EPSRC grant EP/N510129/1 and by the US Office Naval Research (ONR).
REFERENCES
Alberto Abadie and Guido W Imbens. Matching on the estimated propensity score. Econometrica, 84(2):781-807 2016.
Ahmed Alaa and Mihaela van der Schaar. Limits of estimating heterogeneous treatment effects: Guidelines for practical algorithm design. In International Conference on Machine Learning, pp. 129-138, 2018.
Ahmed M Alaa and Mihaela van der Schaar. Bayesian inference of individualized treatment effects using multi-task gaussian processes. In Advances in Neural Information Processing Systems, pp- 3424-3432, 2017.
Muhammad Ali, Humaira Naureen, Muhammad Haseeb Tariq, Muhammad Junaid Farrukh, Abubakar Usman, Shahana Khattak, and Hina Ahsan. Rational use of antibiotics in an intensive care unit: retrospective study of the impact on clinical outcomes and mortality rate. Infection and Drug Rexistance, 12:493, 2019.
Elja Arjas and Jan Parner. Causal reasoning from longitudinal data. Scandinavian Journal of Statistics, 31(2):171-177 2004.
Onur Atan, William R Zame, and Mihaela van der Schaar. Learning optimal policies from observa- tional data. International Conference on Machine Learning CausalML workshop. 2018.
Peter C Austin. An introduction propensity score methods for reducing the effects of confounding in observational studies. Multivariate behavioral research, 46(3):399 2011.
Helmut Bartsch, Heike Dally, Odilia Popanda, Angela Risch, and Peter Schmezer. Genetic risk profiles for cancer susceptibility and therapy response. In Cancer Prevention, pp. 19-36. Springer, 2007.
Shai nen-David, John Blitzer, Koby Crammer, and Fernando Percira. Analysis of representations for domain adaptation. In Advances in neural information processing systems, pp. 137-144, 2007.
Y Bengio, A Courville, and P Vincent. Representation learning: a review and new perspectives. arxiv. org. 2012.
CM Booth and IF Tannock. Randomised controlled trials and population- observational research: partners in the evolution of medical evidence. British journal of cancer, 1 110(3):551, 2014.
Liesbet De Bus, Bram Gadeyne, Johan Steen, Jerina Boelens, Geert Claeys, Dominique Benoit, Jan De Waele, Johan Decruyenaere, and Pieter Depuydt. A complete and multifaceted overview of antibiotic use and infection diagnosis in the intensive care unit: results from prospective four-year registration. Critical Care, 22(1):241, 2018.
Shayan Doroudi, Philip Thomas. and Emma Brunskill. Importance sampling for fair policy selection. Grantee Submission, 2017.
Yarin Gal and Zoubin Ghahramani. A theoretically grounded application of dropout in recurrent neural networks. Advances in neural information processing systems, pp. 1019-1027, 2016.
Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario Marchand, and Victor Lempitsky. Domain-adversarial training of neural networks. The Journal of Machine Learning Research, (711):096-2033 -030, 2016.
Changran Geng, Harald Paganetti, and Clemens Grassberger. Prediction of treatment response for combined cne radiation therapy for non-small cell lung cancer patients using a bio- mathematical model. Scientific reports, 131:23 2017.
10

---

Published as a conference paper at ICLR 2020
Zhaohan Guo, Philip S Thomas, and Emma Brunskill. Using options and covariance testing for long horizon ofllcy policy evaluation. In Advances in Neural Information Processing Systems, pp. 2492-2501, 2017.
Assaf Hallak, François Schnitzler, Timothy Mann, and Shie Mannor. Off-policy model-based learning under unknown factored dynamics. In International Conference on Machine Learning, pp. 711-719, 2015.
Miguel A Hernán, Babette Brumback, and James M Robins. Marginal structural models to estimate the joint causal effect of nonrandomized treatments. Journal the American Statistical Association, 96(454):440-444 448, 2001.
Miguel Ángel Heránn, Babette Brumback, and James M Robins. Marginal structural models to estimate the causal effect of zidovudine on the survival of hiv-positive men. Epidemiology, pp. 561-570, 2000.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8): 1735-1780, 1997.
William Hoiles and Mihaela Van Der Schaar. A non non- con-parmetric learning method for confidently estimating patient's clinical state and dynamics. In Advances in Neural Information Processing Systems, pp. 2020-2028, 2016.
Chanelle Howe, Stephen R Cole, Shruti H Mehta, and Gregory D Kirk. Estimating the effects of multiple time-varying exposures using joint marginal structural models: alcohol consumption, injection drug use, and hiv acquisition. Epidemiology (Cambridge, Mass.), 23:4):574, 2012.
Kosuke Imai and Mare Ratkovic. Covariate balancing propensity score. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 76(1):243-22) 2014.
Kosuke Imai and David A Van Dyk. Causal inference with general treatment regimes: Generalizing the propensity score. Journal of the American Statistical Association, 99(467):854-866, 2004.
Nan Jiang and Lihong Li. Doubly robust off-policy value evaluation for reinforcement learning. arXiv preprint 1rXiv:1511.0372 2015.
Fredrik Johansson, Uri Shalit, and David Sontag. Learning representations for counterfactual inference. In International conference on machine learning, pp. 3020-3029, 2016.
Alistair EW Johnson, Tom J Pollard, Lu Shen, Lehman Li-wei, 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:160035, 2016.
Daniel Kifer, Shai Ben-David, and Johannes Gehrke. Detecting change in data streams. In Proceedings of the Thirtieth international conference on Very large data bases-Volume 30, pp. 180-191. VLDB Endowment, 2004.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Sheng Li and Yun Fu. Matching balanced nonlinear representations for treatment effects estimation. In Advances in Neural Information Processing Systems, pp. 929-939, 2017.
Ya Li, Xinmei Tian, Mingming Gong, Yajing Liu, Tongliang Liu, Kun Zhang, and Dacheng Tao. Deep domain generalization via conditional invariant adversarial networks. In Proceedings of the European Conference on Computer Vision (ECCV), pp. 624-639, 2018.
Bryan Lim, Ahmed Alaa, and Mihaela van der Schaar. Forecasting treatment responses over time using recurrent marginal structural networks. In Advances in Neural Information Processing Systems, pp. 7493-7503, 2018.
Judith Lok et al. Statistical modeling of causal effects in continuous time. The Annals of Statistics, 3663):144-11577, 2008.
11

---

Published as a conference paper at ICLR 2020
Mohammad Ali Mansournia, Goodarz Danaei, Mohammad Hossein Forouzanfar, Mahmood Mah- moodi, Mohsen Jamali, Nasrin Mansournia, and Kazem Mohammad. Effect of physical activity functional performance and knee pain in patients with ostcoarthritis: analysis with marginal structural models. Epidemiology, pp. 631-640, 2012.
Mohammad Ali Mansournia, Mahyar Etminan, Goodarz Danaei, Jay S Kaufman, and Gary Collins. Handling time varying confounding in observational research. bmj, 359:j4587, 2017.
Kathleen M Mortimer, Romain Neugebauer, Mark Van Der Laan, and Ira B Tager. An application of model-fitting procedures for marginal structural models. American Journal of Epidemiology. 162 (4):382-388, 2005.
Daniel Neil, Michael Pfeiffer, and Shih-Chii Liu. Phased Istm: Accelerating recurrent network training for long or event-based sequences. In Advances in Neural Information Processing Systems, pp. 3882-3890, 2016.
Jersey Neyman. Sur les applications de la théorie des probabilités aux experiences agricoles: Essai des principes. Roczniki Nauk Rolniezych, 10:1-51, 1923.
Cosmin Päduraru, Doina Precup, Joelle Pineau, and Gheorghe Comânici. An empirical analysis of off-policy learning in discrete mdps. In European Workshop Reinforcement Learning, pp. 89-102, 2013.
Judea Pearl et al. Causal inference in statistics: An overview. Statistics surveys, 3:96-146, 2009.
Robert W Platt, Enrique Schisterman, and Stephen R Cole. mide-fediiied confounding. American journal of epidemiology, 170(6):6876694 2009.
Doina Precup. Eligibility traces for off-policy policy evaluation. Computer Science Department Faculty Publication Series, pp. 80, 2000.
James Robins. A approach to causal inference in mortality studies with sustained exposure crriiodapplicattion to control of the healthy worker survivor effect. Mathematical modelling, (9-12):1993-1512, 1986.
James Robins. Correcting for non-compliance in randomized trials using structural nested mean models. Communications in Statistics-Thee Theory and methods, 23(8):2379-241) 1994.
James M Robins. Association, causation, and marginal structural models. Synthese, 121(1):151-17,, 1999.
James M Robins and Miguel Heráán. Estimation of the causal effects of time-varying exposures. In Longitudinal data analysis, pp. 547-593. Chapman and Hall/CRC, 2008
James M Robins, Miguel Angel Hernan, and Babette Brumback. Marginal structural models and causal inference epidemiology, 2000.
Jason Roy, Kirsten J Lum, and Michael J Daniels. A bayesian nonparametric approach to marginal structural models for point treatments and a continuous or survival outcome. Biostatistics, 18(1): 32-47, 2016.
Donald B Rubin. Bayesian inference for causal effects: The role of randomization. The Annals of statistics, pp. 34-58, 1978.
Enrique Schisterman, Stephen R Cole, and Robert W Platt. Overadjustment bias and unnecessary adjustment in epidemiologic studies. Epidemiology (Cambridge, Mass.), 20(4):488, 2009.
Peter Schulam and Suchi Saria. Reliable decision support using counterfactual models. In Advances in Neural Information Processing Systems, pp. 1697-1708, 2017.
Alice Schoenauer Sebag, Louise Heinrich, Mare Schoenauer, Michèle Sebag, Lani Wu, and Steven Altschuler. Multi-domain adversarial learning. In Seveentht annual International Confer- on Learning Representations, 2019.
12

---

Published as a conference paper at ICLR 2020
Uri Shalit, Fredrik D Johansson, and David Sontag. Estimating individual treatment effect: general- ization bounds and algorithms. In Proceedings of the 34th International Conference on Machine Learning-
Khamingg-Vlluee Volume 70, pp. 3076-3085. JMLR. org, 2017.
Hossein Soleimani, Adarsh Subbaswamy, and Suchi Saria. Treatment-response models for coun- terfactual reasoning with contimoous-time, continuous-valued interventions. arXiv preprint 02000, 2017.
Adith Swaminathan and Thorsten Joachims. Batch learning from logged bandit feedback through counterfactual risk minimization. Journal of Machine Learning Research, 16(1):1731-1755, 2015a.
Adith Swaminathan and Thorsten Joachims. The self-normalized estimator for counterfactual learning. In advances in neural information processing systems, pp. 3231 9231-3239, 2015b.
Philip Thomas, Georgios Theocharous, and Mohammad Ghavamzadeh. High-confidence off-policy evaluation. In Twenty-Ninth AAAI Conference on Artificial Intelligence, 2015.
U Waheed, P Williams, S Brett, G Baldock, and N Soni. White cell count and intensive care unit outcome. Anaesthexia, 58(2):180-18) 2003.
Yanbo Xu, Yanxun Xu, and Suchi Saria. A bayesian nonparametric approach for estimating individu- alized geatmenntresponne curves. In Machine Learning for Healthcare Conference, pp. 282-300, 2016.
Liuyi Yao, Sheng Li, Yaliang Li, Mengdi Huai, Jing Gao, and Aidong Zhang. Representation learning for treatment effect estimation from observational data. In Advances in Neural Information Processing Systems, pp.
2018.
Jinsung Yoon, James Jordon, and Mihaela van der Schaar. Ganite: Estimation of individualized treatment effects using generative adversarial nets. International Conference on Learning Repre- sentations (ICLR), 2018.
13

---

Published as a conference paper at ICLR 2020
APPENDIX
A EXTENDED RELATED WORK
Causal inference in the static setting: large number of methods have been proposed learn treatment effects from observational data in the static setting. In this case, it is needed to adjust for the selection bias; bias caused by the fact that, in the observational dataset, the treatment assignments depend on the patient features. Several ways of handling the selection bias involve using propensity matching (Austin 2011] [mai & Ratkovic 2014 Abaaie & Imbens| (016), building representations where treated and un-treated populations had similar distributions (Johansson al. 2016 Shalit et al. 2017; & Fu 2017, Yao et al.
2018) or performing propensity-awar hyperparameter tuning (Alaa
& van der Schaar 2017 2018). However, these methods for the static setting cannot be extended directly to time- varying treatments (Hernán et al. 2000 Schisterman et al. 2009).
Learning optimal policies: A related problem to ours involves learning the optimal treatment policies from logged data (Swaminathan & Joachims 2015a b Atan et al., 2018). That is, learning the treatment option that would give the best reward. Note the difference to the causal inference setting considered in this paper, where the aim is to learn the counterfactual patient outcomes under all possible treatment options. Learning all of the counterfactual outcomes is harder problem and can also be used for finding the optimal treatment.
A method for learning optimal policies, proposed by Atan al. 2018) uses domain adversarial training to build a representation that is invariant to the following two domains: observational data and simulated randomized clinical trial data, where the treatments have equal probabilities, Atan et al. (2018) only considers the static setting and aims to choose the optimal treatment instead of estimating all of the counterfactual outcomes. In our paper the aim is to eliminate the bias from the time-dependent confounders and reliably estimate all of the potential outcomes; thus, at each timestep we build a representation that is invariant to the treatment.
Off-policy evaluation in reinforcement learning: In reinforcement learning. a similar problem to ours is off-policy evaluation, which uses retrospective observational data, also known as logged bandit feedback (Hoiles & Van Der Schaar 2016 Pāduraru et al.) (2013 Doroudi al.] (2017). In this case, the retrospective observational data consists of sequences of states, actions and rewards which were generated by an agent operating under unknown policy. The off-policy evaluation methods aim to use this data to estimate the expected reward of target policy. These methods use algorithms based on importance sampling (Precup) 20000 Thomas al. 20155] Guo et al. (011)) action-value function approximation (model based) (Hallak et al. 2015) or doubly robust combination of both approaches (Jiang & Li, 2015). Nevertheless, these methods focus on obtaining average rewards of policies, while in our case the aim is to estimate individualized patient outcomes for future treatments.
B ASSUMPTIONS
The standard assumptions needed for identifying the treatment effects are (Robins & Hernán 2008 Lim et al. 2018 Schulam & Saria 2017):
Assumption 1: Consistency. If A, a, for a given patient, then the potential outcome for treatment a, is the same as the observed (factual) outcome: Y Y++1 [a]] Y(+1.
Assumption 2: Positivity (Overlap) (Imai & Van Dyk 2004
If P(A, = 1; X, = x,) * 0
then P(A, a, A a X, x2) for all a,
Assumption 3: Sequential strong ignorability. Y+ti IL A. A, Va, E Vt.
Assumption 2 means that, for each timestep. each treatment has non-zero probability of being assigned. Assumption 3 means that there are no hidden confounders, that is, all of covariates affecting both the treatment assignment and the outcomes are present in the the observational dataset. Note that while assumption 3 is standard across all methods for estimating treatment effects, it is not testable in practice (Robins et al. 2000, Pearl et al. 2009)
14

---

Published as a conference paper at ICLR 2020
C TIME-DEPENDENT CONFOUNDING
Figure 6 fillustrates the causal graphs for time-varying exposures with 2-steps (Robins et al. 2000). In Figure (a), the covariate is a time-dependent confounder because affects the treatment assignments and at the same time, its value is changed by past treatments (Mansournia et al. 2017). as illustrated by the red arrows. Thus, the treatment probabilities at each time / depend on the history of covariate X and past treatments. Note that Uo and U1 are hidden variables which only affect the covariates, i.c. they do not have arrows into the treatments. Thus, the no hidden confounders assumption (Assumption 3) is satisfied.
Figure (a) and (b) illustrate the two cases when there is no bias from time-dependent confounding. In Figure 6 (a) the treatment probabilities are independent, while in Figure they depend past treatments.
X
2
Uo
(a) Bias from time-dependent confounders.
Xo
Ao
Xo
A
2
Uo
Uo
U1
(a) Option for removing bias from time-dependent confounders.
(a) Option 2 for removing bias from time-dependent confounders.
Figure 6: Causal graphs for 2-step time-varying exposures (Robins al. 2000). Xo, Xy are patient covariates, Ag. A are treatments, Un. U, are unobserved variable and Y2 is the outcome.
Marginal Structural Models Robins et 2000). To remove the association between time- dependent confounders and time- varying treatments, Marginal Structural Models propose using inverse probability of treatment weighting (IPTW). Without loss of generality, consider the use of MSMs with univariate treatments, baseline variables and outcomes. The outcome after timesteps is parametrized as follows: E GY++1 ai.... a,, V] = g(a)..... a. V: where g(-) is usually linear function with parameters 6. To remove the bias from the dinendepnneent confounders present the observational dataset, in the regression model g(.) MSMs weights each patients using either stabilized weights:
f(Ar A [= f(A: X ALL V)
SW())
II
(8)
or unstabilized weights:
W(t) = II
(9)
f(A, A
1=1
where f(.) represents the conditional probability mass function for discrete treatments.
15

---

Published as a conference paper at ICLR 2020
Inverse probability of treatment weighting (IPTW) creates peeuulppopulation where each member consists of themselves and -1 (or SW 1) copies added though weighting. this pseudo- population, Robins Robins (1999) shows that X, does not predict treatment A,, thus removing the bias from time-dependent confounders.
When using unstabilized weights W, the causal graph in the pseudo-populane is the in Figure 6 (a) where P(A, X, A, V) = P(A,). On the other hand, when using stabilized weights SW, causal graph in the pouullppottttion is the in Figure 6 (b) where P(A, Xx; A1-1, V) P(Ar A.-
Counterfactual Recurrent Networks. Instead of using IPTW, we proposed building representation of X, A, V that is not predictive of treatment A.. At timestep t. we have l different possible treatments A {Ai.... Ak) We build representation the history and covariates and treatments that has the same distribution across the different possible treatments: P( At-1, V) A, A1)
(X,,X,, A V) A, = Ax). By breaking the association between past exposure and current treatments A,, we satisfy the causal graph in Figure (a) and thus we remove the bias from time-dependent confounders.
16

---

Published as a conference paper at ICLR 2020
D PROOF OF THEOREM
We first prove the following proposition.
Proposition For fixed i, let x' ((h)). Then the optimal prediction probabilities of G. are
given by
(a()) (x') Σ (x))
G (x')
(10)
Proof. For fixed 4, the optimal prediction probabilities are given by
K (=((=)=
G arg (og( (x')de'
subject to
=1
(1I)
Maximising the value function pointwise and applying Lagrange multiplies, we get
K
G arg max Ga
log log(G
(x))
Σ
-2(0)-1)
(12)
Setting the derivative (w.r.t. Gy (x')) to and solving for G% (a') we get
Po (x')
(13)
where À can now be solved for using the constraint to be X = - (x'). This gives the result.
Proof. (of Theorem By substituting the expression from Proposition into the minimax game defined in Eq. [6] the objective for & becomes
Ph (ar)
min
w' p/
pd
log
>K p.o
(xx)
(14)
We then that
K
===1 p.t (at((a) )
p pt((') (x²)
K Bxr
p Pa (x')
(p(())) (xx) )
Ex'n Pp
log
+ K log K =
log
+ log K
5K
5K
===1 P.
(15)
K Eg'rg
(e
log
(16)
1 EK P.O (zl))
K
Σ KL (f (z')
1 K
Pf (x')
(17)
K JSD(P ..., P
(18)
where K KL is the Kullback-Leine divergence and JSD( ) is the
Jensen-
Shannon Divergence eLi al. 2018). Since K log K is constant and the non-negative and 0 if and only all distributions are equal, we have that pt = :..= P
JSD is
17

---

Published as a conference paper at ICLR 2020
E TRAINING PROCEDURE FOR CRN
Let D a. -1t1lt= U (v())) / be an observational dataset consisting of information (=1 about N independent patients that we use to train CRN. The encoder and decoder networks part of CRN are trained into two separate steps.
To begin with, the encoder is trained to built treatment invariant representations of the patient history and to perform one-step ahead prediction. After the encoder is optimized, we use compute the balancing representation br for each timestep in the trajectory of patient (i). To train the decoder, modify the training dataset as follows. For each patient (i), split their trajectory into shorter sequences of the Tmax timesteps of the form:
{ br (y Bip 1444 INE1
1mmaz Uv(i),
(19)
for 1... for training the decoder. The different sequences obtained for all patents are randomly grouped into minibatches and used for training.
T Tmax- Thus, each patients contributes with T(1) Tmax examples in the dataset
The pseudocode in Algorithm T shows the training procedure used for the encoder and decoder networks part of CRN. The model was implemented in TensorFlow and trained on an NVIDIA Tesla K80 GPU. The Adam optimizer (Kingma & Ba 2014) was used for training and both the encoder and the decoder are trained for 100 epochs.
18

---

Published as a conference paper at ICLR 2020
Algorithm Pseudo-code for training CRN
Input: Training data: D (x
a .itUEE U Uv())
( ( T
i=1
(1) Encoder optimization: parameters Pg B, E Learning rate: for p 1,..., max epochs do
X 1 + exp(-10-2
for Batch B {x4 as 18tt1li=0 Uv()] [8|
{
(0) (1) ()) 17())
in epoch do (6)
Compute [ (frrr BE,A) = [8] EieB Et=1 t.a (fg,r, OF.0)
Compute C (F FF = [B
EieB E1=1 C.. F0 F..))
OC (.,.) Xp
D. (E.
Bg f Fgr - 1
80
00g.
ac
R. (0E...0E..)))
Dg + 0E
ac 1..
pg + EE
00
end for end for
(2) Compute the encoder balanced representation and use it to initialize the decoder hidden state.
for i = l.... NN do for t 1, .,.... ( do br. = encoder(x)) a v(0) 0rr)
end for end for
(3) Split dataset in sequences of Tmax timesteps:
{(D u (yll, a vieet1 14411 tmm Uv(i)
7 Tmaax
N
1=1
(4) Optimize decoder: parameters 0p 0p,a, OD. Learning rate: for p ... max epochs do
A + exp(-10) pp)
(1) (8) Tmaa U[r(0)] [8]
for Batch B br U (y a Y14e+1 11+4+1110 1=o
in epoch do
i=1
Compute C. OD,a) T3T Eies Etm (()) (0 OD.))
Compute C (0 0D.) [ EieB Σ [ (0D.r OD.))
ocg (PD..,. Ap
ac "t..,(D.... (D.a)
OD, -- 0p, - 1
80p.r
d0p.r
8C (DD...D.s))
0D.y +- 0D. |
D.
3C (...).
0n.e OD.a
dUD.a
end for end for
Output: Trained CRN encoder (parameters gg AE BE.,) and trained CRN decoder (parameters Orr 0 D, OD.y
19

---

Published as a conference paper at ICLR 2020
F PHARMACOKKNETIC-HHARAACCCMMM
MODEL OF TUMOUR GROWTH
To evaluate the CRN on counterfactual estimation, we need access to the data generation mechanism to build test that consists of patient outcomes under all possible treatment options. For this purpose, we use the state-of-the-ar pharmacokineti--pharnaaaa parre (PK-PD) model of tumour growth proposed by Geng al. (017)) and also used by Lim al. (2018) for evaluating RMSMs. The PK-PD model characterizes patients suffering from non soall cell lung cancer and models the evolution of their tumour under the combined effects of chemotherapy and radiotherapy. In addition, the model includes different distributions of tumour sizes based on the cancer stage at diagnosis.
Model of tumour growth The volume of tumour days after diagnosis modelled as follows:
V +
(1 plog(
K
3,0(()
(c,,d()) + B,)(1)²)
V(()
V(t)
(20)
Chemotherapy
Radiotherapy
Noise
Tumor growth
where the parameters K. p. Se. ar, B, are sampled from the prior distributions described in (Geng et al. 2017) and et ~ (V 0.013) is a noise term that accounts for randomness in the tumour growth.
To incorporate heterogeneity among patient responses, due to, for instance, gender or genetic factors Bartsch al. (2007), the prior means for 3e and O, are adjusted to create three patient subgroups g {1, 2, 3) as described in Lim et al. 2018)). This way, we incorporate in the model of tumour growth specific characteristics that affect the patient's individualized response to treatments. Thus, the prior mean M&a. of Be and the prior mean Na. of a, are augmented as follows.
Me. (1) = {
1. [1.1 if ((() 33
1. 1.1 its =1
Map (1)
(21)
p otherwise
Pa. otherwise
where s4g, and Ma, are the mean parameters from Geng al. (2017) and 1lg, (1) and Pla, (i) are the parameters used in the data simulation. The patient subgroup S()) E (1,2,3) is used as baseline features.
The chemotherapy drug concentration follows an exponential decay with half life of 1 day:
C(t) Č(t) + C(t 1)/2,
(22)
where C(t) 5.0mg//' of Vinblastine if chemotherapy is given at time t. d(t) 2.0Gy fractions of radiotherapy if the radiotherapy treatment applied at timestep t.
Time-varying confounding is introduced by modelling chemotherapy and radiotherapy assignment as Bernoulli random variables, with probabilities Pa and p, depending on the tumour diameter:
(t) = a ( To (D(1) ( ( (3)
( ( 7/- (D(I) (D())-8,)) 6.)). (Dias.
(t)
(23)
where D(t) is the average tumour diameter over the last 15 days, Dinx = 13em is the maximum tumour diameter and a al.) is the sigmoid activation function. The parameters & and 6, set to 5e = 8, = D such that there is 0.5 probability of receiving treatment when tumour is half of its maximum size. To: Tr control the amount of time-dependent confounding; the higher 7 is, the more important the history of tumour diameter is in assigning treatments. Thus, at each timestep, there are four treatment options options: no treatment (A1), chemotherpy (A2), radiotherapy (Ag). combined chemotherapy and radiotherapy (A,).
Since the work most relevant to ours the one of Lim et al. (018 we used the same data simulation and same settings for y % = 'ir as in their case. When Y O, there is no deeedeeendent confounding and the treatments are randomly assigned. By increasing Y we increase the influence of
the volume size history (encoded in D(t) on the treatment probability. For example, assume D(t)= 3Dmax From equation (7), the probability of chemotherapy in this case is (t) oi (D(1)
Dmaz )) = (0.25ry)) where a--) (-) is the sigmoid function. When 7 = 1, (t( = 0.56, when y = 5, be increased further to increase the bias. However, the values used in the experiments evaluate the model on a wide range of settings
(t) = 0.77 and when y 10, (t) 0.92 in this example.
for the time-dependent confounding bias.
20

---

Published as a conference paper at ICLR 2020
G MARGINAL STRUCTURAL MODELS
Marginal Structural Models (Robins et al. 2000 Hernán et al. 2001) have been widely used in epidemiology and part of follow up studies. In our case, we would like to estimate the effects of sequence of treatments in the future given the current patient history:
B(Y++T At,, 1) t 1), H.) = g(5, a H.)
(24)
where g a generic function and a(l,1+T-1) = .., represents a possible sequence of treatments from timestep just until before the potential outcome Y is observed. After removing the bias form time-dependent confounders, (Y A(t, t + T 1) = a(t, t + T 1), H+)
E(Y++ 1+r (
1)].
Note that for implementing MSMs, we encode the treatments at timestep & the model of tumour growth as A. [Ar.c At,al indicate the binary application of chemotherapy and radiotherapy. In order to remove the time-dependent confounding bias and estimate future outcomes, we use the stabilized weights of MSMs weight each patient in the dataset:
1+7
f(An An-1)
L+T
[N S(Ank An-1) IIEE(Ca) f(An, An-1, X,, V)²
SW(,T)
II
(25)
n f(A, An
n==
where f()) represents the conditional probability mass function for discrete treatments.
We adopt the implementation in (Hermán et al. 200T Howe et al. 2012 Lim et al. 2018 for MSMs and use logistic regression for estimating the propensity weights as follows:
f(Ar, A(-1) o
Σ Wk Σ Ac.)))
(26)
f(Ackk H₂) =
Σ óe!
Agk)) + +
(27)
KE (c,d)
where wi. o. and w. are regression coefficients, &o E {c, indicates the chemotherapy or radiother- apy treatments and (-)) is the sigmoid function.
For predicting the outcome, the following regression model is used, where each individual patient is weighted by its propensity score:
th-
g(T,
H.)
3A( (E(d)))
Σ
Σ
+1x,,
12x,1
(28)
where B, and 1, are regression coefficients.
MSMs do not require hyperparameter tuning so we use the patients from both the train and validation sets for training.
H RECURRENT MARGINAL STRUCTURAL NETWORKS
MSMs are very sensitive to model missspcciiication in computing the propensity weights and estimating the outcomes. Recurrent Marginal Structural Models (RMSNs) (Lim et al. 2018) overcome this problem by using recurrent neural networks to estimate the propensity scores and to build the outcome model. RNNs
robust to changes in the treatment assignment policy. RMSNs were implemented as descried in Lim et al. 22
For implementing RMSNs, we also encode the treatments timestep in the model of tumour growth as A, = [Arer Ar,a) to indicate the binary application of chemotherapy and radiotherapy. The propensity weights are estimated using recurrent neural networks as follows:
fAA,k A,-1) = RNNsw. (A+-1) f(Ackk
X,, A,-1) = RNNsw, ( X,
(29)
We used the publicly available implementation from
https://g whw...thbbooooo..o rmsn
2018
21

---

Published as a conference paper at ICLR 2020
For predicting one-step-ahed outcome, R-MSNs encoder network: g(L, at, t), H,) = RRN, ş (a,, A, V),
(30)
where in the loss function, each patient weighted by their stabilized IPTW.
For estimating the treatment responses for sequence of treatments in the future, RMSNs use decoder network:
g(T, (t,a(t,t+7-1))
RNND 1: A V).
(31)
See Lim et al. (0088) for more details about the R-MSNs model architecture and training procedure of the propensity weights, encoder and decoder networks. Tables 2 and show the hyperparameter search ranges used optimize this model for evaluation in our paper. The hyperparameterss were
selected in the same way as proposed by Lim et al. (2018), based on the error on the factual outcomes in the validation dataset. All of the models are trained using Adam optimizer for 100 epochs.
Table 2: Hyperparameter search range for propensity networks and encoder (same as in Lim et al. ())18)). C is the size of the input.
Hyperparameter
Search range
Iterations of Hyperparameter Search
Learning rate
0.01, 0.005, 0.001 64, 128, 256 0.5C, 1C, 2C, 3C, 4C 0.1, 0.2, 0.3, 0.4, 0.5 0.5, 1.0,
Minibatch size
RNN state size
Dropout rate
Max Gradient Norm
Table 3: Hyperparameter search range for decoder (same as in Lim et al. (0088) ).
C the input size.
Hyperparameter
Search range
Iterations of Hy plperparraeeerr Search Learning rate Minihatch size
20
0.01, 0.001, 0.0001
256, 512, 1024
RNN state size
IC, 2C, 4C, 8C, 16C
Dropout Rate
0.1, 0.2, 0.3, 0.4, 0.5
Max Gradient Norm
0.5, 1.0, 2.0. 4.0
I BASELINE RNN AND LINEAR MODEL
For the baseline linear model, we fit the same regression model used for Marginal Structural Networks, but without using the IPTW. The baseline RNN uses an LSTM unit and, at each timestep, receives as input the current treatment, the patient covariates and the patient static features to perform one-step- ahead prediction. To have model of similar capacity to the CRN (similar number of parameters), we add a fully connected layer on top of the output of the LSTM unit in order to obtain the outcomes. Table 4 shows the hyperparameter search range used to optimize this model. The hyperparameters were selecting according to the error on the factual outcomes in the validation set. We train the baseline RNN using the Adam optimizer for 100 epochs.
22

---

Published as a conference paper at ICLR 2020
Table 4: Hyperparameter search range for baseline RNN model. C is the size of the input.
Hyperparameter
Search range
Iterations of Hyperparameler Search
50
Learning rate Minibatch size
0.01, 0,001, 0,0001
64, 128, 256 0.5C, IC, 2C, 4C 0,5C, 1C, 2C, 3C, 4C
RNN hidden units
hidden units
RNN dropout probability
0.1, 0.2, 0.3, 0.4,
J HYPERPARAMETER OPTIMIZATION FOR CRN
As described in Appendix C, the dataset for training the decoder are used by splitting the sequences of the patients in the training set. This creates a larger dataset for training (where each patient (i) contributes T')) Tmax times to the dataset) which requires a different hyperparameter search range. Moreover, the balancing representations computed by the encoder are used to initialize the state of the RNN for the decoder. Thus, the decoder RNN size equal to the size of the balancing representation size of the encoder. Table shows the hyperparameter search ranges for the encoder and decoder networks CRN. We selected hyperparameters hased on the error of the model on the factual outcomes in the validation dataset. All models trained for 100 epochs.
In addition, Tables and 7 illustrate the optimal hyperparameters chosen.
Table 5: Hyperparameter search range for CRN encoder. C is the size of the input and R is the size of the balancing representation.
Hyperparameter
Search range encoder
Search range decoder
Iterations of Hyperparameter Search Learning rate Minibatch size
50
30
0.01, 0.001, 0.0001
0.01, 0.001, 0.0001
64, 128, 256
256, 512, 1024
RNN hidden units
0.5C, 2C, 3C, 4C
Balancing representation size of encoder 0.5C, IC, 2C, 3C, 0.5R, 1R, 2R, 3R, 4R 0.1, 0.2, 0.3, 0.4, 0.5
Balancing representation size
0.5C. IC, 2C, 3C,
FC hidden units
0.5R, IR, 2R, 3R, 4R
RNN dropout probability
0.1, 0.2, 0.3, 0.4, 0.5
23

---

Published as a conference paper at ICLR 2020
Table Optimal hyperparameters for the CRN encoder when different degrees of time-dependent confounding are applied in the model of tumour growth. The parameters Ye and "/- measure the degree of dime-dependent confounding applied. When % and " are to the same value, we denote this with 7..
7
7
7
2
7/
3
7/4
7/4
5
Learning rate Minibatch size
0.001
0.1
0.001
0.01
0.01
0.001
64
128
128
RNN hidden units Balancing representation size
12
18
18
24
24
18
18
12
18
6
12
FC hidden units RNN dropout probability
18
18
36
54
24
48
0.1
=
7
8
7/+ 9
% 10
Learning rate Minibatch RNN hidden units Balancing representation size
0.001
0.001
0.01
0.001
0.01
64
128
128
128
24
18
12
24
24
12
18
24
18
12
FC hidden units RNN dropout probability
48
72
12
36
12
0.1
0.2
0.1
0.1
0.1
= %r
5,
0
Learning rate Minibatch size
0.01
0.001
128
64
RNN hidden units Balancing representation size
18
FC hidden units
36
96
RNN dropout probability
0.1
0.1
Table 7: Optimal hyperparameters for the CRN decoder when different degrees of time-dependent confounding are applied in the model of tumour growth. The parameters Ye and "rr measure the degree of time-dependent confounding applied. When % and y set to the same value, denote this with 7.
7/4
7/4
7.
7
5
Learning Minihatch size RNN hidden units Balancing representation size
0.001
0.001
0.001
0.001
0.001
1024
1024
512
1024
1024
18
12
18
6
12
18
18
6
18
FC hidden units RNN dropout probability
18
36
18
72
0.1
0.2
0.3
0.1
0, 7rr
%e
5, 7r
Learning rate Minibatch size
0.01
0.001
512
1024
RNN hidden units
18
24
Balancing representation size
18
12
FC hidden units
36
24
RNN dropout probability
0.1
0.03
24

---

Published as a conference paper at ICLR 2020
K FULL RESULTS FOR COUNTERFACTUAL PREDICTION
K.1 MULTI-STEP AHEAD PREDICTION OF COUNTERFACTUALS
Figure shows the normalized RMSE for multiple step-ahead prediction of counterfactuals. The RMSE is normalized by the maximum tumour volume: Vmaz 1150cm The counterfactuals in this case are generated as described in Section 6.3 and Appendix I. We notice that performance gains of CRN compared to RMSN increase with the number of future timesteps for which the counterfactuals are estimated.
5%
RMSU
RMSS
RMSN
CRV
CRN
4%
Nomaal lzed RMSE
Noma lzze RMSE
Norallzdd
39
2%
Degree time-dependend confunding
Degree of time-dependend confunding
Degree of time-dependend corfunding
(a) Two-step ahead prediction
(b) Three-step ahead prediction
(c) Four-step ahead prediction
8%
RMSN
RMSN
CRN
6%
CRN
4%
Normaliee RMSSE
Normllzed
0%
2%
Degres dine-dependend confunding 7
Dagis of time-eeeerdnnd confunding
(d) Five-step ahead prediction
(e) Six-step ahead prediction
Figure Results for prediction of patient counterfactuals for multiple steps ahead.
25

---

Published as a conference paper at ICLR 2020
K.2 DETAILED RESULTS FOR THE COUNTERFACTUAL PREDICTIONS
Tables 8 and 9 show detailed results for the counterfactual predictions.
Table 8: Normalized RMSE for one-step-ahead prediction of counterfactuals. The parameter measures the degree of time-dependent confounding applied.
2
Y 3
Y =
Y = 5
Linear (no IPTW)
0.99%
1.08%
1.36%
1.68%
2.11%
2.77%
MSM
0.99%
1.08%
1.34%
1.63%
2.02%
2.61%
RNN
0.70%
0.70%
0.84%
1.05%
1.24%
1.69%
CRN (0=
0.66%
0.77%
0.92%
0.95%
1.24%
1.54%
RMSN
0.60%
0.61%
0.72%
0.81%
0.94%
1.23%
CRN
0.56 %
0.57% %
0.62%
0.67 %
0.87%
1.20%
7
YV
7 8
Y
7 10
Linear (no IPTW)
3.55%
4.15%
4.80%
5.09%
5.22%
MSM
3.30%
3.79%
4.30%
4.47%
4.47%
RNN
2.03%
2.52%
2.88%
3.79%
4.01%
CRN (A=0)
1.98%
2.42%
2.73%
3.17%
3.57%
RMSN
1.70%
2.18%
2.37%
2.77%
2.83%
CRN
1.48% %
1.56 %
2.05%
2.36%
2.41%
Table 9: Normalized RMSE for 7-step-ahead prediction of counterfactuals. The parameter y measures the degree of time-dependent confounding applied.
==1
2
3
Y
5
2
RMSN
0.90%
1.15%
1.53%
2.14%
2.91%
CRN
0.84%
0.96%
1.21
1.46%
2.45%
3
RMSN
0.97%
1.36%
1.87%
2.44%
3.47%
CRN
0.86% %
0.96%
1.47% %
1.51 %
2.84%
RMSN
1.24%
1.79%
2.60%
3.33%
3.88%
CRN
0.91%
1.08
1.74%
1.76%
2.82%
RMSN 1.51%
2.1%%
3.06%
4.07%
4.58%
CRN
0.85%
1.10%
1.73%
2.00%
3.43%
6
RMSN
2.10%
2.89%
3.06%
4.16%
6.32%
CRN
1.16%
1.52%
2.29 %
2.66 %
4.91
26

---

Published as a conference paper at ICLR 2020
L TEST SET GENERATION FOR EVALUATING TIMING OF TREATMENT
In order to evaluate how well the models select the correct treatment and timing of treatment we simulate counterfactual outcomes as follows. We generate 1000 test samples using the model of tumour growth described in Section 6. Let H, be the current history of the patient and let T be future time horizon. For each timestep in the future, we have 4 treatment options at: no treatment (Ap), chemotherapy (A1), radiotherapy (A2), chemotherapy and radiotherapy. (Ag).
Using the model of tumour growth where the outcome u+r is given by the volume of the tumour, we generate the following 27 counterfactuals:
Chemotherapy application
Yetr
A1, Ap,.. a Ap, H,
(32)
Yotr
a. An, = -A,,.../477-1 An, H,
a1+7-1
(33)
Yorr
a. = ab,1 = A.......aau = A1, H
(34)
Radiotherapy application
Yotr
a a241 = Ap.... a = An, H
(35)
Y r+r
a, Ao, ac+1 Ag....
ac+r-1 . Ag, Hr
(36)
Y YeT
Ao, Ao,ac+1 ac+1
An,.. a A2, H,
(37)
We perform this for each patient in the test set and at each time in the history. For instance, for patient with timesteps the model of tumour growth and for time horizon T = we generate 2 50 300 counterfactuals.
Using the true generated couterfactual data, we select the treatment that has the lowest Y among the 7 options generated for each treatment. Then, select the time of applying treatment (among t, 1, -++-11111 1) that resulted in the lowest +- For each model, we generate the counterfactuals under the same treatment plans and patient histories. Then, perform the selection of treatment and timing of treatment in the same way and we compare these with the true ones. Note that in order to account for numerical instability (two outcomes Y having very similar values), we consider two outcomes the same if they are within = 0.001 of each other.
27

---

Published as a conference paper at ICLR 2020
M RESULTS ON FACTUAL PREDICTION ON MIMIC III
Using the Medical Information Mart for Intensive Care (MIMIC III) (Johnson et al.) 2016) database consisting of electronic health records from patients in the ICU, we also show how the CRN can be used on a real medical dataset. From MIMIC III we extracted the patients on antibiotics, with trajectories up 30 timesteps. thus obtaining a dataset with 3487 patients. For each patient, we extracted 25 patient covariates including lab tests and vital signs measured over time, as well as static patient features such as age and gender.
We used a binary treatment at each timestep indicating whether the patient was administered antibi- oties or not. Note that for the longitudinal covariates we used aggregate value for each day since the ICU admission. The reason for this is because antibiotic treatment is decided daily for the patient. We split the dataset into 2826/33 2/16/31//// /313/348 patients for training, validation and testing respectively. We performed hyperparameter optimization on the validation patient set, using the search ranges in Table S and we again selected hyperparameters based on the error on the factual outcomes.
We estimate the individualized effect of antibiotics assigned over time on the patient's 's white blood cell count. A high white blood cell count associated with severe illness and poor outcome for ICU patients (Waheed et al. 2003). Antibiotic administration in the ICU aims reduce the white blood cell count. However, the effectiveness of the antibiotics treatment in reducing the white blood cell count highly dependent on the time they are administered with respect to the history of the patient covariates. In this context we again have time-dependent confounders: the patient features change over time and affected by the previous administration of antibiotics. Moreover, the history of the patient features also determines antibiotics administration and affects future patient outcomes (De Bus et al. 2018 Ali et al. 2019).
In Table 10 we report the root mean squared error for factual prediction of the patients' white blood cell count for multiple prediction horizons т. Note that for this dataset we do not have access to counterfactual data, which is why we report error on factual predictions.
Table 10: RMSE for r-step-ahead prediction of factual outcomes on MIMIC III.
T
T = 2
T=
T
RMSN
2.84
3.87
4.46
4.79
CRN
2.68
3.54
4.07
4.67
We notice that CRN also achieves better performance than RMSN estimating factual outcomes in a real-world dataset containing electronic health records. In this context, where couterfactual data is not available, domain expert knowledge is required to validate the model's counterfactual predictions under other antibiotic treatment alternatives. This further medical validation is outside the scope of this paper.
28


Reviews:
[{"review_id": "rylO65p7Kr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper adapts domain adversarial training to construct treatment invariant representation for adjusting time-varying confounding in counterfactual treatment outcome predictions over time. \n\nThe proposed method extends the recent static settings of balancing representation in counterfactual inference to longitudinal settings, and it also overcomes the problem of high-variance of IPTW weighting in MSN based methods such as the current state-of-art method RMSN.\n\nThe paper is very well written. Method is novel to me. Experiments are sufficient. However since I am not very familiar with this area, so there can be things I miss in the evaluation of this paper.\n\nThere is one problem with the illustrations in Figure 1: predictions on the potential outcomes before the first treatment was given should have stayed on the same path, and then depart when different treatment was initiated (at different time)", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I did not assess the derivations or theory.", "review_assessment:_checking_correctness_of_experiments": "I did not assess the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}, {"review_id": "S1lBUXOaKr", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This work addresses the problem of causal inference in time-dependent treatment regimes. To address the problem, the authors propose an extension of the balancing representation for causal inference framework that seeks render the current treatment independent from a representation of the history of treatment and confounders. This is sensibly actualized within an RNN. The authors provide empirical results that demonstrate the proposed method performing very well in comparison to prior art. \n\nThis paper is well written, proposes a sensible solution to a difficult problem, and performs well empirically. It would be nice to see theory that connects the imbalance to the expected error of the causal estimand, and a notion of hyperparameter tuning (the authors say they search over hyperparameters but do not specify how the parameters are ultimately selected. However, I don't feel that either of those points outweigh the benefit of the paper. ", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "H1e1Ni-3FB", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper introduces Counterfactual Recurrent Network (CRN) that is able to estimate the effects of various treatments from longitudinal data. The claim is that the model can decide (i) treatment plan; (ii) optimal time of treatment; and (iii) when to stop treatment. The proposed method attempts to learn time-invariant representations that are not predictive of the next treatment by borrowing ideas from Ganin, et al. (2016)’s work on for domain adversarial training. In fact, this paper is an extension of (Atan et al., 2018) to be applicable for longitudinal data. \n\nThis paper should be rejected due to the following arguments:\n\t- Page 1, par. 4, line 4-7: The example given in sentence “For instance, if ...” is true regardless of time. Therefore, it does not describe a situation where a time-variant confounder might make causal inference hard. In fact, nowhere in the paper there is a solid description of how time-variant confounders would break the current approaches for causal inference. \n\n\t- Page 3, par. 4, last sentence: I don’t understand why sequential treatments and change of covariates through time can stop us from using the conventional methods for learning balanced representations. Why not consider each time point as an instance (thus representing each patient with multiple instances) and learn the representation from the collection of all these instances?\nAll in all, I don’t understand why the authors assume that the world is non-Markovian and that they have to consider the entire \\bar{H}_t. Especially when in the their Experiments section the model used for tumor growth that generates the synthetic data is Markov, and consequently, does not check if the proposed method works in a non-Markovian world.\nIn summary, if the goal is to solve the non-Markov case, the presented experiments don’t show that the proposed method is addressing that; and if the  goal is to solve the Markov case: (i) why insist on framing the problem as non-Markov; and (ii) (Atan et al., 2018) already have solved this (since we established that there are no time-dependent confounders and as the authors state on page 4, par. 3, line -4 to -3 “the novelty here comes from the use of domain adversarial training to eliminate bias from the time-dependent confounders”).\n\n\t- Page 6, Theorem 1 and 2 lines above it: This theorem does not prove that the confounding bias is completely removed because the objective function in Eq. (4) being optimized also includes a loss term for outcome prediction. In fact, IF there is a confounding bias (that [partially] determines both treatment and outcome), it would be wrong to remove it from the learned representation.\n\nThings to improve the paper that did not impact the score:\n\t- Page 2, par. 1, line -2: Use of adverb “Moreover” is wrong. Extreme weights are due to division by small Pr(t|x). The consequence of this numerical instability is high-variance estimates.\n\t- Page 2, par. 4, 1st sentence: incomprehensible\n\t- Page 2, par. 4, line -2: “but also” should come after a “not only”.\n\t- Page 3, par. 4, line -1: incorrect usage of citation; both author name(s) and publication year must be in parentheses.\n\t- Page 4, par. 5, last line: sentence “For testing, the … “ needs elaboration.\n\t- Page 6, par. 1, line -3: “start of” should be “start off”.\n\nReferences:\n\t- Atan, O., Zame, W. R., & Van Der Schaar, M. (2018). Counterfactual Policy Optimization Using Domain-Adversarial Neural Networks. In ICML CausalML workshop.\n\t- Ganin, Y., Ustinova, E., Ajakan, H., Germain, P., Larochelle, H., Laviolette, F., ... & Lempitsky, V. (2016). Domain-adversarial training of neural networks. The Journal of Machine Learning Research, 17(1), 2096-2030.\n\n\n********UPDATE after reading the rebuttal********\nThe authors have addressed my major concerns in their rebuttal and therefore I have increased my rating form “reject” to “weak accept”. \n\n", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers uniformly suggest acceptance. Please look carefully at reviewer comments and address in the camera-ready. Great work!
Decision: Accept (Spotlight)