Title: ESTIMATING COUNTERFACTUAL TREATMENT OUTCOMES OVER TIME THROUGH ADVERSARIALLY BALANCED REPRESENTATIONS

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.

1. INTRODUCTION:

As clinical decision-makers are often faced with the problem of choosing between treatment alternatives for patients, reliably estimating their effects is 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 a more viable alternative for estimating treatment effects.

A 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 et al., 2018; Alaa & van der Schaar, 2018; Yao et al., 2018) and only a few methods address the longitudinal setting (Xu et al., 2016; Roy et al., 2016; Soleimani et al., 2017; Schulam & Saria, 2017; Lim et al., 2018) . However, estimating the effects of treatments over time poses unique opportunities such as 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-dependent treatments from observational data involves correctly handling the time-dependent confounders: patient covariates that are affected by past treatments which then influence future treatments and outcomes (Platt et al., 2009) . For instance, consider that treatment A is given when a certain patient covariate (e.g. white blood cell We illustrate 3 patients with different covariate and treatment histories Ht . For a current time t, 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 in 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 are more likely to 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 A is harmful to patients. Moreover, estimating the effect of a 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 B.

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 on a 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 A 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.

Time-dependent confounders are 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 timevarying exposures are based on ideas from epidemiology. The most widely used among these are Marginal Structural Models (MSMs) (Robins et al., 2000; Mansournia et al., 2012) which use the inverse probability of treatment weighting (IPTW) to adjust for the time-dependent confounding bias. Through IPTW, MSMs create a pseudo-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-sequence architecture 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 problems of existing methods for causal inference over time. Our main contributions are as follows.

Treatment invariant representations over time. CRN constructs treatment invariant representations at each timestep in order to break the association between patient history and treatment assignment and thus removes the bias from time-dependent confounders. For this, CRN uses domain adversarial training (Ganin et al., 2016; Li et al., 2018; Sebag et al., 2019) to trade-off between build-ing this balancing representation and predicting patient outcomes. We show that these representations remove the bias from time-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 a 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-sequence architecture. CRN consists of an encoder network which builds treatment invariant representations of the patient history that are used to 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. We illustrate in Figure 1 the applicability of our method in choosing optimal cancer treatments.

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-art methods.

2. RELATED WORK:

We focus on methods for estimating treatment effects over time and for building balancing representations 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-computation formula, Structural Nested Models and Marginal Structural Models (MSMs) (Robins, 1986; 1994; Robins et al., 2000; Robins & Hernán, 2008) . Originally, these methods have used predictors performing logistic/linear regression which makes them unsuitable for handling complex time-dependencies (Hernán et al., 2001; Mansournia et al., 2012; Mortimer et al., 2005) . To address these limitations, methods that use Bayesian non-parametrics 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 a generalized mixed-effects model and uses the g-computation method to handle time-varying confounders. Soleimani et al. (2017) extend the approach in Xu et al. (2016) to the continuous time-setting and model treatment responses using linear time-invariant dynamical systems. Roy et al. (2016) use Dirichlet and Gaussian processes to model the observational data and estimate the IPTW in Marginal Structural Models. Schulam & Saria (2017) 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-parametric methods make strong assumptions about model structure and consequently cannot handle well heterogeneous treatment effects arising from baseline variables (Soleimani et al., 2017; Schulam & Saria, 2017) and multiple treatment outcomes (Xu et al., 2016; Schulam & Saria, 2017) .

The work most related to ours is the one of Lim et al. (2018) which improves on the standard MSMs by using recurrent neural networks to estimate the inverse probability of treatment weights (IPTWs). Lim et al. (2018) introduces Recurrent Marginal Structural Networks (RMSNs) which also use a sequence-to-sequence 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 a 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 (Johansson et al., 2016; Shalit et al., 2017; Li & Fu, 2017; Yao et al., 2018) . Moreover, due to the sequential assignment of treatments in the longitudinal setting, and due to the change of patient covariates over time according to previous treatments, the methods for the static setting are not directly applicable to the time-varying setting (Hernán et al., 2000; Mansournia et al., 2012) .

3. PROBLEM FORMULATION:

Consider an observational dataset D = {x

(i) t , a (i) t , y (i) t+1 } T (i)

t=1 ∪{v (i) } N i=1 consisting of information about N independent patients. For each patient (i), we observe time-dependent covariates

X (i) t ∈ X t , treatment received A (i) t ∈ {A 1 , . . . A K } = A and outcomes Y (i) t+1 ∈ Y t+1 for T (i) discrete timesteps.

The patient can also have baseline covariates V (i) ∈ V such as gender and genetic information. Note that the outcome Y (i) t+1 will be part of the observed covariates X (i) t+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 extended by (Robins & Hernán, 2008) to account for time-varying treatments. Let Y[ā] be the potential outcomes, either factual or counterfactual, for each possible course of treatment ā. Let Ht = ( Xt , Āt-1 , V) represent the history of the patient covariates Xt = (X 1 , . . . , X t ), treatment assignments Āt = (A 1 , . . . , A t ) and static features V. We want to estimate:

EQUATION

where ā(t, t + τ -1) = [a t , . . . a t+τ -1 ] represents a possible sequence of treatments from timestep t just until before the potential outcome Y t+τ is observed. We make the standard assumptions (Robins et al., 2000; Lim et al., 2018) needed to identify the treatment effects: consistency, positivity and no hidden confounders (sequential strong ignorability). See Appendix B for more more details.

4. COUNTERFACTUAL RECURRENT NETWORK:

The observational data can be used to train a supervised learning model to forecast: E(Y t+τ | Ā(t, t + τ -1) = ā(t, t + τ -1), Ht ). However, without adjusting for the bias introduced by timevarying confounders, this model cannot be reliably used for making causal predictions (Robins et al., 2000; Robins & Hernán, 2008; Schulam & Saria, 2017) . The Counterfactual Recurrent Network (CRN) removes this bias through domain adversarial training and estimates the counterfactual outcomes E(Y t+τ [ā(t, t + τ -1)]| Ht ), for any intended future treatment assignment ā(t, t + τ -1).

Balancing representations. The history Ht = ( Xt , Āt-1 , V) of the patient contains the timevarying confounders Xt which bias the treatment assignment A t ∈ {A 1 , . . . A K } in the observational dataset. Inverse probability of treatment weighting, as performed by MSMs, creates a pseudopopulation where the probability of treatment A t does not depend on the time-varying confounders (Robins et al., 2000) . In this paper, we propose instead building a representation of the history Ht that is not predictive of the treatment A t . This way, we remove the association between history, containing the time-varying confounders Xt , and current treatment A t . Robins (1999) shows that in this case, the estimation of counterfactual treatment outcomes is unbiased. See Appendix C for details and for an example of a causal graph with time-dependent confounders.

Let Φ be the representation function that maps the patient history Ht to a representation space R.

To obtain unbiased treatment effects, Φ needs to construct treatment invariant representations such that P (Φ( Ht )

| A t = A 1 ) = • • • = P (Φ( Ht ) | A t = A K ).

To achieve this and to estimate counterfactual outcomes under a 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 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, rather than the use of sequence-to-sequence models, which have already been applied to forecast treatment responses (Lim et al., 2018) . Figure 2 illustrates our model architecture.

Encoder. The encoder network uses an RNN, with LSTM unit (Hochreiter & Schmidhuber, 1997) , to process the history of treatments Āt-1 , covariates Xt and baseline features V to build a treatment invariant representation Φ( Ht ), but also to predict one-step-ahead outcomes Y t+1 . To achieve this, the encoder network aims to maximize the loss of the treatment classifier G a and minimize the loss of the outcome predictor network G y . This way, the balanced representation Φ( Ht ) is not predictive of the assigned treatment A t , but is discriminative enough to estimate the outcome Y t+1 . 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. During training, the decoder uses as input the outcomes from the observational data (Y t+1 , . . . Y t+τ -1 ), the static patient features V and the intended sequence of treatments ā(t, t + τ -1). The decoder is trained in a similar way to the encoder to update the balanced representation and to estimate the outcomes. During testing, we do not have access to ground-truth outcomes; thus, the outcomes predicted by the decoder ( Ŷt+1 , . . . Ŷt+τ-1 ) 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 Φ( Ht ) is built by applying a fully connected layer, with Exponential Linear Unit (ELU) activation to the output of the LSTM. The treatment classifier G a and the predictor network G y consist of a hidden layer each, also with ELU activation. The output layer of G a uses softmax activation, while the output layer of G y 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 E 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 et al., 2016) .

5. ADVERSARIALLY BALANCED REPRESENTATION OVER TIME:

At each timestep t, let the K different possible treatments A t ∈ {A 1 , . . . A K } represent our domains.

As described in Section 4, to remove the bias from time-dependent confounders, we build a representation of history Ht that is invariant across treatments:

P (Φ( Ht ) | A 1 ) = • • • = P (Φ( Ht ) | A K ).

This requirement can be enforced by minimizing the distance in the distribution of Φ( Ht ) between any two pairs of treatments. Let G a (Φ( Ht ); θ a ) be the treatment classifier with parameters θ a and let G j a (Φ( Ht ); θ a ) be the output corresponding to treatment A j . Let G y (Φ( Ht ); θ y ) be the predictor network with parameters θ y . The representation function Φ is parameterized by the parameters θ r in the RNN: Φ( Ht ; θ r ). Figure 3 shows the adversarial training procedure used.

For timestep t and patient (i), let L (i) t,a (θ r , θ a ) be the treatment (domain) loss and let L (i) t,y (θ r , θ y ) the outcome loss, defined as follows:

L (i) t,a (θ r , θ a ) = - K j=1 I {a (i) t =aj } log(G j a (Φ( Ht ; θ r ); θ a ))

(2)

EQUATION

)

If the outcome is binary, the cross-entropy loss can be used instead for L t,y . To build treatment invariant representations and to also estimate patient outcomes, we aim to maximize treatment loss and minimize outcome loss.

Thus, the overall loss

L (i)

t,y at timestep t is given by:

EQUATION

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, we use the Gradient Reversal Layer (GRL) (Ganin et al., 2016) . For more details about the training procedure, see Appendix E.

By using the objective L (i) t (θ r , θ y , θ a ), we reach the saddle point ( θr , θy , θa ) that achieves the equilibrium between domain discrimination and outcome estimation.

( θr , θy ) = arg min θr,θy

EQUATION

The result stated in Theorem 1 proves that the treatment (domain) loss part of our objective (from equation 2) aims to remove the time-dependent confounding bias.

Theorem 1. Let t ∈ {1, 2, . . . }. For each j = 1, ..., K, let P j denote the distribution of Ht conditional on A t = A j and let P Φ j denote the distribution of Φ( Ht ) conditional on A t = A j . Let G j a denote the output of G a corresponding to treatment A j . Then the minimax game defined by

min Φ max Ga K j=1 E Ht∼Pj log(G j a (Φ( Ht ); θ a )) subject to K j=1 G j a (Φ( Ht )) = 1 (6)

has a global minimum which is attained if and only if P Φ 1 = P Φ 2 = ... = P Φ K , i.e. when the learned representations are invariant across all treatments.

Proof. This result is a restatement of the one in Li et al. (2018) . For details, see the Appendix D.

A good representation allows us to obtain a low error in estimating counterfactuals for all treatments, while at the same time to minimize the H-divergence between induced marginal distributions of all the domains. We use an algorithm that directly minimizes a combination of the H-divergence 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 CRN1 , we evaluate it on a Pharmacokinetic-Pharmacodynamic model of tumour growth (Geng et al., 2017) , which uses a stateof-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 t days after diagnosis is modelled as follows:

EQUATION

where K, ρ, β c , α r , β r , e t are sampled as described in Geng et al. (2017) . To incorporate heterogeneity in patient responses, the prior means for β c and α r are 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 p c and p r depending on the tumour diameter: p c (t) = σ γc Dmax ( D(t)δ c ) and p r (t) = σ γr Dmax ( D(t)δ r ) where D(t) is the average diameter over the last 15 days, D max = 13cm, σ(•) is the sigmoid and δ c = δ r = D max /2. The amount of time-dependent confounding is controlled through γ c , γ r ; the higher γ 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 G for details). We also compare against the Recurrent Marginal Structural Networks (RMSNs) Lim et al. (2018) , which is the current state-of-the-art model in estimating treatment responses. RMSNs use RNNs to estimate the IPTWs and the patient outcomes (details in Appendix H). To show that standard supervised learning models do not handle the timevarying confounders we compare against an RNN and a linear regression model, which receive as input treatments and covariates to predict the outcome (see Appendix I for details). Our model architecture follows the description in Sections 4 and 5, with full training details and hyperparameter optimization in Appendix J. To show the importance of adversarial training, we also benchmark against CRN (λ = 0) a model with the same architecture, but with λ = 0, i.e our 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, to build decision support systems, we need to evaluate how well the models estimate the counterfactual outcomes, i.e patient outcomes under alternative treatment options. The parameters γ c and γ r control the treatment assignment 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 γ = γ c = γ r . For each γ we simulate a 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 Y t+1 , represented by tumour volume V (t + 1), under all possible treatment options. Balancing representation: To evaluate whether the CRN has indeed learnt treatment invariant represenations, for γ = 5, we illustrate in Figure 5 the T-SNE embeddings of the balancing representations Φ( Ht ) built by the CRN encoder for test patients. We color each point by the treatment A t ∈ {no treatment, chemotherapy, radiotherapy, combined chemotherapy and radiotherapy} received at timestep t to highlight the invariance of Φ( Ht ) across the different treatments. In Figure 5 (b), we show Φ( Ht ) only for chemotherapy and radiotherapy for better understanding.

-100 0

x-tsne Notice that Φ( Ht ) is not predictive of the treatment A t 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 τ and for each time t in a patient's trajectory, there are τ options for giving chemotherapy at one of t, . . . predicting Y t+τ . See Appendix L for more details about the test set. The models are evaluated for 3 settings of γ c and γ r .

Table 1 shows the results for this evaluation set-up. The treatment accuracy denotes the percentage of patients for which the correct treatment was selected, while the treatment timing accuracy is the percentage for which the correct timing was selected. Note that when γ c = 0 and γ r = 5, RMSN and MSM select the wrong treatment timing for projection horizons τ > 4. CRN performs similarly 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 in the ICU.

7. CONCLUSION:

Despite its wide applicability, the problem of causal inference for time-dependent treatments has been relatively less studied compared to 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 over time or for estimating the individualized effects of time-dependent treatments with associated dosage.

In this paper, we 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 timedependent confounders through adversarial training. Using a model of tumour growth, we validated CRN in realistic medical scenarios and we showed improvements over existing state-of-the-art 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 to provide theoretical guarantees for the expected error on the counterfactuals.

Inverse probability of treatment weighting (IPTW) creates a pseudo-population where each member consists of themselves and W -1 (or SW -1) copies added though weighting. In this pseudopopulation, Robins Robins (1999) shows that Xt does not predict treatment A t , thus removing the bias from time-dependent confounders.

When using unstabilized weights W , the causal graph in the pseudo-population is the one in Figure 6 (a) where P (A t | Xt , Āt-1 , V ) = P (A t ). On the other hand, when using stabilized weights SW , causal graph in the pseudo-population is the one in Figure 6 (b) where

P (A t | Xt , Āt-1 , V ) = P (A t | Āt-1 ).

Counterfactual Recurrent Networks. Instead of using IPTW, we proposed building a representation of Xt , Āt-1 , V that is not predictive of treatment A t . At timestep t, we have k different possible treatments A t ∈ {A 1 , . . . A K }. We build a representation of the history and covariates and treatments that has the same distribution across the different possible treatments:

P (Φ( Xt , Āt-1 , V) | A t = A 1 ) = • • • = P (Φ( Xt , Āt-1 , V) | A t = A K ).

By breaking the association between past exposure and current treatments A t , we satisfy the causal graph in Figure 6 (a) and thus we remove the bias from time-dependent confounders.

D PROOF OF THEOREM 1:

We first prove the following proposition. Proposition 1. For fixed Φ, let x = Φ( ht ). Then the optimal prediction probabilities of G a are given by

EQUATION

)

Proof. For fixed Φ, the optimal prediction probabilities are given by

EQUATION

)

Maximising the value function pointwise and applying Lagrange multiplies, we get

EQUATION

)

Setting the derivative (w.r.t. G j a * (x )) to 0 and solving for G j a * (x ) we get

EQUATION

)

where λ can now be solved for using the constraint to be λ = -

K i=1 P Φ i (x )

. This gives the result.

Proof. (of Theorem 1) By substituting the expression from Proposition 1 into the minimax game defined in Eq. 6, the objective for Φ becomes

EQUATION

)

We then note that

EQUATION

)

where KL(•||•) is the Kullback-Leibler divergence and JSD(•, ..., •) is the multi-distribution Jensen-Shannon Divergence (Li et al., 2018) . Since K log K is a constant and the multi-distribution JSD is non-negative and 0 if and only if all distributions are equal, we have that P Φ 1 = ... = P Φ K .

Algorithm 1 Pseudo-code for training CRN

Input: Training data: D = {x (i) t , a (i) t , y (i) t+1 } T (i) t=1 ∪ v (i) N i=1

(1) Encoder optimization: parameters θ E,r , θ E,a , θ E,y . Learning rate: µ for p = 1, . . . , max epochs do

EQUATION

t , y

(i) t+1 } T (i) t=0 ∪ v (i) |B| i=1 in epoch do Compute L B E,a (θ E,r , θ E,a ) = 1 |B| i∈B T (i) t=1 L (i) t,a (θ E,r , θ E,a ) Compute L B E,y (θ E,r , θ E,y ) = 1 |B| i∈B T (i) t=1 L (i) t,y (θ E,r , θ E,y ) θ E,r ← θ E,r -µ ∂L B E,y (θ E,r ,θ E,y ) ∂θ E,r -λ p ∂L B E,a (θ E,r ,θ E,a ) ∂θ E,r θ E,y ← θ E,y -µ ∂L B E,y (θ E,r ,θ E,y ) ∂θ E,y θ E,a ← θ E,a -µ ∂L B E,a (θ E,r ,θ E,a ) ∂θ E,a

end for end for

(2) Compute the encoder balanced representation and use it to initialize the decoder hidden state. for i = 1, ..., N do for t = 1, . . . , T (i) do br i) ; θ E,r ) end for end for

EQUATION

(3) Split dataset in sequences of τ max timesteps:

br (i) l ∪ {y (i) l+t , a (i) l+t , y (i) l+t+1 } τmax t=1 ∪ v (i) T (i) -τmax l=1 N i=1

(4) Optimize decoder: parameters θ D,r , θ D,a , θ D,y Learning rate: µ for p = 1, . . . , max epochs do

λ p = 2 1 + exp(-10 • p) -1 for Batch B = br (i) l ∪ {y (i) l+t , a (i) l+t , y (i) l+t+1 } τmax t=0 ∪ {v (i) } |B| i=1 in epoch do Compute L B D,a (θ D,r , θ D,a ) = 1 |B| i∈B τmax t=1 L (i) t,a (θ D,r , θ D,a ) Compute L B D,y (θ D,r , θ D,y ) = 1 |B| i∈B τmax t=1 L (i) t,y (θ D,r , θ D,y ) θ D,r ← θ D,r -µ ∂L B D,y (θ D,r ,θ D,y ) ∂θ D,r -λ p ∂L B D,a (θ D,r ,θ D,a ) ∂θ D,r θ D,y ← θ D,y -µ ∂L B D,y (θ D,r ,θ D,y ) ∂θ D,y θ D,a ← θ D,a -µ ∂L B D,a (θ D,r ,θ D,a ) ∂θ D,a

end for end for Output: Trained CRN encoder (parameters θ E,r , θ E,a , θ E,y ) and trained CRN decoder (parameters θ D,r , θ D,a , θ D,y . )

F PHARMACOKINETIC-PHARMACODYNAMIC MODEL OF TUMOUR GROWTH:

To evaluate the CRN on counterfactual estimation, we need access to the data generation mechanism to build a test set that consists of patient outcomes under all possible treatment options. For this purpose, we use the state-of-the-art pharmacokinetic-pharmacodynamic (PK-PD) model of tumour growth proposed by Geng et al. (2017) and also used by Lim et al. (2018) for evaluating RMSMs. The PK-PD model characterizes patients suffering from non-small 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.

The volume of tumour t days after diagnosis is modelled as follows:

V (t + 1) = 1 + ρlog( K V (t) )

Tumor growth

EQUATION

where the parameters K, ρ, β c , α r , β r are sampled from the prior distributions described in (Geng et al., 2017 ) and e t ∼ N (0, 0.01 2 ) 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 et al. (2007) , the prior means for β c and α r are adjusted to create three patient subgroups S (i) ∈ {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 µ βc of β c and the prior mean µ αr of α r are augmented as follows.

EQUATION

where µ βc and µ αr are the mean parameters from Geng et al. (2017) and µ βc (i) and µ αr (i) are the parameters used in the data simulation. The patient subgroup S (i) ∈ {1, 2, 3} is used as baseline features.

The chemotherapy drug concentration follows an exponential decay with half life of 1 day:

EQUATION

where C(t) = 5.0mg/m 3 of Vinblastine if chemotherapy is given at time t. d(t) = 2.0Gy fractions of radiotherapy if the radiotherapy treatment is applied at timestep t.

Time-varying confounding is introduced by modelling chemotherapy and radiotherapy assignment as Bernoulli random variables, with probabilities p c and p r depending on the tumour diameter:

EQUATION

where D(t) is the average tumour diameter over the last 15 days, D max = 13cm is the maximum tumour diameter and σ(•) is the sigmoid activation function. The parameters δ c and δ r are set to δ c = δ r = D max /2 such that there is 0.5 probability of receiving treatment when tumour is half of its maximum size. γ c , γ r control the amount of time-dependent confounding; the higher γ 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 (A 1 ), chemotherpy (A 2 ), radiotherapy (A 3 ), combined chemotherapy and radiotherapy (A 4 ).

Since the work most relevant to ours is the one of Lim et al. (2018) we used the same data simulation and same settings for γ = γ c = γ r as in their case. When γ = 0, there is no time-dependent confounding and the treatments are randomly assigned. By increasing γ we increase the influence of the volume size history (encoded in D(t)) on the treatment probability. For example, assume D(t) = 3Dmax 4

. From equation ( 7), the probability of chemotherapy in this case is p c (t) = σ( γc Dmax ( D(t) -Dmax 2 )) = σ(0.25γ c ), where σ(•) is the sigmoid function. When γ = 1, p c (t) = 0.56, when γ = 5, p c (t) = 0.77 and when γ = 10, p c (t) = 0.92 in this example. γ can be increased further to increase the bias. However, the values used in the experiments evaluate the model on a wide range of settings for the time-dependent confounding bias.

G MARGINAL STRUCTURAL MODELS:

Marginal Structural Models (Robins et al., 2000; Hernán et al., 2001) have been widely used in epidemiology and as part of follow up studies. In our case, we would like to estimate the effects of a sequence of treatments in the future given the current patient history:

EQUATION

)

where g is a generic function and ā(t, t + τ -1) = [a t , . . . 

a

EQUATION

)

where f (•) represents the conditional probability mass function for discrete treatments.

We adopt the implementation in (Hernán et al., 2001; Howe et al., 2012; Lim et al., 2018) for MSMs and use logistic regression for estimating the propensity weights as follows:

EQUATION

EQUATION

where ω , φ and w are regression coefficients, k ∈ {c, d} indicates the chemotherapy or radiotherapy 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:

g(τ, a(t, t + τ -1), Ht ) = k∈{c,d}

EQUATION

where β and l 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 mis-specification 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 are more robust to changes in the treatment assignment policy. RMSNs were implemented as descried in Lim et al. (2018) 2 .

For implementing RMSNs, we also encode the treatments at timestep t in the model of tumour growth as A t = [A t,c , A t,d ] to indicate the binary application of chemotherapy and radiotherapy. The propensity weights are estimated using recurrent neural networks as follows: (i) τ max 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 is equal to the size of the balancing representation size of the encoder. Table 5 shows the hyperparameter search ranges for the encoder and decoder networks in CRN. We selected hyperparameters based on the error of the model on the factual outcomes in the validation dataset. All models are trained for 100 epochs.

f (A t,k | Āt-1 ) = RNN SWn ( Āt-1 ) f (A t,k | Xt , Āt-1 ) = RNN SW d ( Āt-1 , Xt , V) (29)

In addition, Tables 6 and 7 illustrate the optimal hyperparameters chosen. Figure 7 shows the normalized RMSE for multiple step-ahead prediction of counterfactuals. The RMSE is normalized by the maximum tumour volume: V max = 1150cm 3 . 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. 

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. 

γ = 0 γ = 1 γ = 2 γ = 3 γ = 4 γ =

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 Ht be the current history of the patient and let τ be a future time horizon. For each timestep in the future, we have 4 treatment options at: no treatment (A 0 ), chemotherapy (A 1 ), radiotherapy (A 2 ), chemotherapy and radiotherapy. (A 3 ).

Using the model of tumour growth where the outcome Y t+τ is given by the volume of the tumour, we generate the following 2τ counterfactuals: 

We perform this for each patient in the test set and at each time t in the history. For instance, for a patient with 50 timesteps in the model of tumour growth and for time horizon τ = 3, we generate 2 • 3 • 50 = 300 counterfactuals.

Using the true generated couterfactual data, we select the treatment that has the lowest Y t+τ among the τ options generated for each treatment. Then, we select the time of applying treatment (among t, t+1, . . . t+τ -1) that resulted in the lowest Y t+τ . For each model, we generate the counterfactuals under the same treatment plans and patient histories. Then, we 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 t+τ having very similar values), we consider two outcomes the same if they are within = 0.001 of each other.

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 to 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 antibiotics 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/313/348 patients for training, validation and testing respectively. We performed hyperparameter optimization on the validation patient set, using the search ranges in Table 5 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 white blood cell count. A high white blood cell count is associated with severe illness and poor outcome for ICU patients (Waheed et al., 2003) . Antibiotic administration in the ICU aims to reduce the white blood cell count. However, the effectiveness of the antibiotics treatment in reducing the white blood cell count is 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 are 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. We notice that CRN also achieves better performance than RMSN in 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.

:

The implementation of the model can be found at https://bitbucket.org/mvdschaar/ mlforhealthlabpub/src/master/alg/counterfactual_recurrent_network/ and at https://github.com/ioanabica/Counterfactual-Recurrent-Network.

We used the publicly available implementation from https://github.com/sjblim/rmsn_nips_ 2018.

Figures and Tables:

Figure1: Applicability of CRN in cancer treatment planning. We illustrate 3 patients with different covariate and treatment histories Ht . For a current time t, 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 in 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).

Figure 2: CRN architecture. Encoder builds representation Φ( Ht ) that maximizes loss of treatment classifier G a and minimizes loss of outcome predictor G y . Φ( Ht ) is used to initialize the decoder, which continues to update it to predict counterfactual outcomes of a sequence of future treatments.

Figure 3: Training procedure for building balancing representation.

Figure 4: Results for prediction of patient counterfactuals.

Figure4(a) shows the normalized root mean squared error (RMSE) for one-step ahead estimation of counterfactuals with varying degree of time-dependent confounding γ. The RMSE is normalized by the maximum tumour volume: V max = 1150cm 3 . The linear and MSM models provide a baseline for performance as they achieve the highest RMSE. While the use of IPTW in MSMs helps when γ increases, using linear modelling has severe limitations. When there is no 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 (λ = 0) to generalize to estimate outcomes of treatments not matching the training policy. When γ = 10, CRN improves by 48.1% on the same model architecture without domain adversarial training CRN (λ = 0).Our proposed model achieves the lowest RMSE across all values of γ. Compared to RMSNs, CRN improves by ∼ 17% when γ > 6. To highlight the gains of our method even for smaller γ, Figure4(b) shows the RMSE for five-step ahead prediction (with counterfactuals generated as described in Section 6.2 and Appendix L). RMSNs also use a 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 τ -step ahead prediction, IPTW involves multiplying τ weights which can result in high variance. The results in Figure4(b) show the problems with using IPTW to handle the time-dependent confounding bias. See Appendix K for more results on multi-step ahead prediction.

Figure 5: TSNE embedding of the balancing representation Φ( Ht ) learnt by the CRN encoder at different timesteps t. Notice that Φ( Ht ) is not predictive of the treatment A t given at timestep t.

Figure 7: Results for prediction of patient counterfactuals for multiple steps ahead.

Chemotherapy applicationY t+τ | a t = A 1 , a t+1 = A 0 , . . . a t+τ -1 = A 0 , Ht (32) Y t+τ | a t = A 0 , a t+1 = A 1 , . . . a t+τ -1 = A 0 , Ht (33) . . . Y t+τ | a t = A 0 , a t+1 = A 0 , . . . a t+τ -1 = A 1 , Ht (34) Radiotherapy application Y t+τ | a t = A 2 , a t+1 = A 0 , . . . a t+τ -1 = A 0 , Ht (35) Y t+τ | a t = A 0 , a t+1 = A 2 , . . . a t+τ -1 = A 0 , Ht (36) . . . Y t+τ | a t = A 0 , a t+1 = A 0 , . . . a t+τ -1 = A 2 , Ht

Kifer et al. (2004);Ben-David et al. (2007), propose measuring the disparity between distributions based on their separability by a discriminatively-trained classifier.

Table content: <table/>

Results for recommending the correct treatment and timing of treatment.

Table content: <table><tr><td/><td>γ c = 5, γ r = 5</td><td>γ c = 5, γ r = 0</td><td>γ c = 0, γ r = 5</td></tr><tr><td/><td>τ CRN RMSN MSM</td><td colspan="2">CRN RMSN MSM CRN RMSN MSM</td></tr><tr><td colspan="4">Normalized 3 2.43% 3.16% 6.75% 1.08% 1.35% 3.68% 1.54% 1.59% 3.23%</td></tr><tr><td>RMSE</td><td colspan="3">4 2.83% 3.95% 7.65% 1.21% 1.81% 3.84% 1.81% 2.25% 3.52%</td></tr><tr><td/><td colspan="3">5 3.18% 4.37% 7.95% 1.33% 2.13% 3.91% 2.03% 2.71% 3.63%</td></tr><tr><td/><td colspan="3">6 3.51% 5.61% 8.19% 1.42% 2.41% 3.97% 2.23% 2.73% 3.71%</td></tr><tr><td/><td colspan="3">7 3.93% 6.21% 8.52% 1.53% 2.43% 4.04% 2.43% 2.88% 3.79%</td></tr><tr><td colspan="4">Treatment 3 83.1% 75.3% 73.9% 83.2% 78.6% 77.1% 92.9% 87.3% 74.9%</td></tr><tr><td colspan="4">Accuracy 4 82.5% 74.1% 68.5% 81.3% 77.7% 73.9% 85.7% 83.8% 74.1%</td></tr><tr><td/><td colspan="3">5 73.5% 72.7% 63.2% 78.3% 77.2% 72.3% 83.8% 82.1% 72.8%</td></tr><tr><td/><td colspan="3">6 69.4% 66.7% 62.7% 79.5% 76.3% 71.8% 78.6% 69.7% 64.5%</td></tr><tr><td/><td colspan="3">7 71.2% 68.8% 62.4% 72.7% 71.8% 71.6% 71.9% 69.3% 61.2%</td></tr><tr><td colspan="4">Treatment 3 79.6% 78.1% 67.6% 80.5% 76.8% 77.5% 79.8% 75.7% 60.6%</td></tr><tr><td>Timing</td><td colspan="3">4 73.9% 70.3% 63.1% 79.0% 77.2% 73.4% 75.4% 71.4% 58.2%</td></tr><tr><td colspan="4">Accuracy 5 69.8% 68.6% 62.4% 78.3% 73.3% 63.6% 66.9% 31.3% 29.5%</td></tr><tr><td/><td colspan="3">6 66.9% 66.2% 62.6% 73.5% 72.1% 63.9% 65.8% 24.2% 15.5%</td></tr><tr><td/><td colspan="3">7 64.5% 63.6% 62.2% 70.6% 57.4% 44.2% 63.9% 25.6% 12.5%</td></tr></table>

t+τ -1 ] represents a possible sequence of treatments from timestep t just until before the potential outcome Y t+τ is observed. After removing the bias form time-dependent confounders,E(Y t+τ | Ā(t, t + τ -1) = ā(t, t + τ -1), Ht ) = E(Y t+τ [ā(t, t + τ -1)].Note that for implementing MSMs, we encode the treatments at timestep t in the model of tumour growth as A t = [A t,c , A t,d ] to 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 to weight each patient in the dataset:

Table content: <table/>

Hyperparameter search range for baseline RNN model. C is the size of the input. 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

Table content: <table><tr><td>Hyperparameter</td><td>Search range</td></tr><tr><td>Iterations of Hyperparameter Search</td><td>50</td></tr><tr><td>Learning rate</td><td>0.01, 0.001, 0.0001</td></tr><tr><td>Minibatch size</td><td>64, 128, 256</td></tr><tr><td>RNN hidden units</td><td>0.5C, 1C, 2C, 3C, 4C</td></tr><tr><td>FC hidden units</td><td>0.5C, 1C, 2C, 3C, 4C</td></tr><tr><td>RNN dropout probability</td><td>0.1, 0.2, 0.3, 0.4, 0.5</td></tr><tr><td colspan="2">J HYPERPARAMETER OPTIMIZATION FOR CRN</td></tr></table>

Hyperparameter search range for CRN encoder. C is the size of the input and R is the size of the balancing representation.

Table content: <table><tr><td>Hyperparameter</td><td>Search range encoder</td><td>Search range decoder</td></tr><tr><td>Iterations of Hyperparameter Search</td><td>50</td><td>30</td></tr><tr><td>Learning rate</td><td>0.01, 0.001, 0.0001</td><td>0.01, 0.001, 0.0001</td></tr><tr><td>Minibatch size</td><td>64, 128, 256</td><td>256, 512, 1024</td></tr><tr><td>RNN hidden units</td><td colspan="2">0.5C, 1C, 2C, 3C, 4C Balancing representation size of encoder</td></tr><tr><td>Balancing representation size</td><td>0.5C, 1C, 2C, 3C, 4C</td><td>0.5C, 1C, 2C, 3C, 4C</td></tr><tr><td>FC hidden units</td><td>0.5R, 1R, 2R, 3R, 4R</td><td>0.5R, 1R, 2R, 3R, 4R</td></tr><tr><td>RNN dropout probability</td><td>0.1, 0.2, 0.3, 0.4, 0.5</td><td>0.1, 0.2, 0.3, 0.4, 0.5</td></tr></table>

Normalized RMSE for τ -step-ahead prediction of counterfactuals. The parameter γ measures the degree of time-dependent confounding applied.

Table content: <table><tr><td/><td>γ = 1</td><td>γ = 2</td><td>γ = 3</td><td>γ = 4</td><td>γ = 5</td></tr><tr><td colspan="3">τ = 2 RMSN 0.90% 1.15%</td><td>1.53%</td><td>2.14%</td><td>2.91%</td></tr><tr><td>CRN</td><td colspan="2">0.84% 0.96%</td><td colspan="2">1.21% 1.46%</td><td>2.45%</td></tr><tr><td colspan="3">τ = 3 RMSN 0.97% 1.36%</td><td>1.87%</td><td>2.44%</td><td>3.47%</td></tr><tr><td>CRN</td><td colspan="2">0.86% 0.96%</td><td colspan="2">1.47% 1.51%</td><td>2.84%</td></tr><tr><td colspan="3">τ = 4 RMSN 1.24% 1.79%</td><td>2.60%</td><td>3.33%</td><td>3.88%</td></tr><tr><td>CRN</td><td colspan="2">0.91% 1.08%</td><td colspan="2">1.74% 1.76%</td><td>2.82%</td></tr><tr><td colspan="3">τ = 5 RMSN 1.51% 2.13%</td><td>3.06%</td><td>4.07%</td><td>4.58%</td></tr><tr><td>CRN</td><td colspan="2">0.85% 1.10%</td><td colspan="2">1.73% 2.00%</td><td>3.43%</td></tr><tr><td colspan="3">τ = 6 RMSN 2.10% 2.89%</td><td>3.06%</td><td>4.16%</td><td>6.32%</td></tr><tr><td>CRN</td><td colspan="4">1.16% 1.52% 2.29% 2.66%</td><td>4.91%</td></tr></table>

RMSE for τ -step-ahead prediction of factual outcomes on MIMIC III.

Table content: <table><tr><td/><td colspan="4">τ = 1 τ = 2 τ = 3 τ = 4</td></tr><tr><td>RMSN</td><td>2.84</td><td>3.87</td><td>4.46</td><td>4.79</td></tr><tr><td>CRN</td><td>2.68</td><td>3.54</td><td>4.07</td><td>4.67</td></tr></table>

Peer review (iclr2020):
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. This paper should be rejected due to the following arguments: - 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. - 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? All in all, I don’t understand why the authors assume that the world is non-Markovian and that they have to consider the entire _x0008_ar{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. In 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”). - 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. Things to improve the paper that did not impact the score: - 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. - Page 2, par. 4, 1st sentence: incomprehensible - Page 2, par. 4, line -2: “but also” should come after a “not only”. - Page 3, par. 4, line -1: incorrect usage of citation; both author name(s) and publication year must be in parentheses. - Page 4, par. 5, last line: sentence “For testing, the … “ needs elaboration. - Page 6, par. 1, line -3: “start of” should be “start off”. References: - Atan, O., Zame, W. R., & Van Der Schaar, M. (2018). Counterfactual Policy Optimization Using Domain-Adversarial Neural Networks. In ICML CausalML workshop. - 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. ********UPDATE after reading the rebuttal******** The authors have addressed my major concerns in their rebuttal and therefore I have increased my rating form “reject” to “weak accept”.