{'fraction_non_alphanumeric': 0.030544046935010737, 'fraction_numerical': 0.00855838146792875, 'mean_word_length': 4.746958637469587, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
G-Net: A Deep Learning Approach to G-computation for Counterfactual Outcome Prediction Under Dynamic Treatment Regimes

Counterfactual prediction is a fundamental task in decision-making. G-computation is a method for estimating expected counterfactual outcomes under dynamic time-varying treatment strategies. Existing G-computation implementations have mostly employed classical regression models with limited capacity to capture complex temporal and nonlinear dependence structures. This paper introduces G-Net, a novel sequential deep learning framework for G-computation that can handle complex time series data while imposing minimal modeling assumptions and provide estimates of individual or population-level time varying treatment effects. We evaluate alternative G-Net implementations using realistically complex temporal simulated data obtained from CVSim, a mechanistic model of the cardiovascular system.


Introduction
Counterfactual prediction is a fundamental task in decisionmaking. It entails the estimation of expected future trajectories of variables of interest under alternative courses of action (or treatment strategies) given observed history. Treatment strategies of interest are usually time varying (meaning they comprise decisions at multiple time points) and dynamic (meaning the treatment decision at each time point is a function of history up to that time point).
As an example, consider the problem of fluid administration in Intensive Care Units (ICUs) [1]. It is frequently necessary for physicians to adopt strategies that administer large volumes of fluids to increase blood pressure and promote blood perfusion through organs in septic patients. However, such strategies can lead to fluid overload, which can have serious adverse downstream effects such as pulmonary edema. Fluid administration strategies are time varying and dynamic because at each time point physicians decide the volume of fluid to administer based on observed patient history (e.g. blood pressure and volume of fluid already administered) up to that time point. To aid in the choice between alternative dynamic fluid administration strategies, it would be desirable to obtain counterfactual predictions of a patient's probability of developing fluid overload (and other outcomes of interest) were they to follow each alternative strategy going forward given their observed covariate history up to the current time.
Counterfactual prediction is an inherently causal task in that it must account for the causal effects of following different treatment strategies. When treatment strategies of interest are time-varying, so-called "g-methods" [2,3] are required to estimate their effects. G-methods include gcomputation [4][5][6], structural nested models [7,8], and marginal structural models [9,10]. Of these methods, gcomputation is best suited for estimating effects of general dynamic treatment strategies conditional on high dimensional patient histories [11]. G-computation works by estimating the conditional distribution of relevant covariates given covariate and treatment history at each time point, then producing Monte Carlo estimates of counterfactual outcomes by simulating forward patient trajectories under treatment strategies of interest. Regression model(s) for the covariates and outcomes at each time point conditional on observed history are a critical component of this method. While any regression models could in theory be input to the G-computation algorithm, most existing G-computation implementations have employed simple regression models with limited capacity to capture complex temporal and nonlinear dependence structures. In recent years, sequential deep learning methods such as Recurrent Neural Networks (RNNs) have achieved state of the art performance in predictive modeling of complex time series data while imposing minimal modeling assumptions. In this paper, we propose G-Net, a sequential deep learning framework for G-computation. G-Net admits the use of recurrent networks such as LSTMs to model covariates in a manner suitable for G-computation. The G-Net framework is flexible and allows for various configurations depending on the problem at hand. To the best of our knowledge, this is the first work to investigate a RNN based approach to G-computation.
Unfortunately, it is impossible to reliably evaluate counterfactual predictions on real data, since only the outcomes corresponding to treatment strategies that were actually followed can be observed. Consequently, to explore and evaluate various implementations of G-Net, we used simulated data in which counterfactual ground truth can be known. We used CVSim [12], a well established mechanistic model of the cardiovascular system, to estimate counterfactual simulated patient trajectories under various fluid and vasopressor administration strategies. These experiments provide a template for causal model evaluation using complex and physiologically realistic simulated longitudinal data.

Related Work
Several recent works have proposed a deep learning framework for counterfactual prediction from observational data, including [13][14][15]. However, these works have mostly focused on learning point exposure as opposed to timevarying treatment effects, which are the focus of this paper. G-computation for estimating time-varying treatment effects was first proposed by Robins [4]. Illustrative applications of the general approach are provided in [6,16], and summaries of g-computation (and other "g-methods" for estimating time-varying treatment effects) can be found in [2,3]. The g-computation algorithm takes arbitrary regression models as inputs. While most applications (e.g. [6,16]) have thus far employed classical generalized linear models, there is no conceptual barrier to using more complex machine learning regression models. RNNs, and in particular LSTMs, have achieved state of the art performance on a wide variety of time series regression tasks, including healthcare related tasks [17][18][19]. However, despite the popularity and success of RNNs for time series regression, we have not seen in the literature any "deep" implementation of g-computation.
Recently, Lim et al. [20] plugged RNN regression models into history adjusted marginal structural models (MSM) [21] to make counterfactual predictions. However, these MSMs can only make counterfactual predictions under static time-varying treatment strategies that do not depend on recent covariate history. For example, a history adjusted MSM [20] could estimate the probability of fluid overload given patient history under the (static) treatment strategy "give 1 liter fluid each hour for the next 3 hours", but it could not estimate the probability of fluid overload given patient history under the (dynamic) treatment strategy "each hour for the next 3 hours, if blood pressure is less than 65 then give 1 liter fluids, otherwise give 0 liters". History adjusted MSMs cannot estimate effects of time-varying treatment strategies that respond to changes in the patient's health history, but g-computation can. Further, g-computation is able to straightforwardly estimate the distribution of a counterfactual outcome under a timevarying treatment strategy. This is not straightforward to do with history adjusted MSMs.
Schulam and Saria [22] propose Counterfactual Gaussian Processes, an implementation of continuous time gcomputation, and like us apply their method to ICU data. They only consider static time-varying treatment strategies, though it appears that their method might straightforwardly be extended to handle dynamic strategies as well. An advantage of Gaussian processes is interpretable uncertainty quantification. However, GPs are intractable for large datasets since they have time complexity of O(N 3 ), where N is the number of observations [23]. Sparse GPs, which introduce M inducing points, have at least O(M 2 N ) [23] time complexity. RNNs are more scalable. Recurrent dropout can also be employed in RNN based implementations to produce uncertainty estimates that approximate posterior distributions from a Gaussian process [24].

G-computation for counterfactual prediction
Our goal is to predict patient outcomes under various future treatment strategies given observed patient histories. Let: • t ∈ {0, . . . , K} denote time, assumed discrete, with K being the end of followup; • A t denote the observed treatment action at time t; • Y t denote the observed value of the outcome at time t • L t denote a vector of covariates at time t that may influence treatment decisions or be associated with the outcome; •X t denote the history X 0 , . . . , X t and X t denote the future X t , . . . , X K for arbitrary time varying variable X.
At each time point, we assume the causal ordering (L t , A t , Y t ). Let H t ≡ (L t ,Ā t−1 ) denote patient history preceding treatment at time t. A dynamic treatment strategy g is a collection of functions {g 0 , . . . , g K }, one per time point, such that g t maps H t onto a treatment action at time t. A simple dynamic strategy for fluid volume might be g t (H t ) = .5 × 1{bp t < 65}, i.e. give .5 liters of fluid if mean arterial blood pressure is less than 65 at time t.
Let Y t (g) denote the counterfactual outcome that would be observed at time t had, possibly contrary to fact, treatment strategy g been followed from baseline [4]. Further, let Y t (Ā m−1 , g m ) with t ≥ m denote the counterfactual outcome that would be observed had the patient received their observed treatmentsĀ m−1 through time m − 1 then followed strategy g from time m onward.
In counterfactual point prediction, our goal is to estimate expected counterfactual patient outcome trajectories given observed patient history through time m for any m and any specified treatment strategy g. We might also be interested in estimating the counterfactual outcome distributions at future time points If we do not condition on anything in H m , then (1) is an expectation (and (2) a distribution) over the full population.
If we condition on a small subset of variables contained in patient history, then (1) is an expectation (and (2) a distribution) over a sub-population. If we condition on all elements of a patient history, then (1) is still technically only an expectation (and (2) a distribution) over a hypothetical sub-population with the exact patient history conditioned on, but in this case (1) and (2) practically amount to what is usually meant by personalized prediction.
Assumption 1 states that the observed outcome is equal to the counterfactual outcome corresponding to the observed treatment. Assumption 2 states that there is no unobserved confounding of treatment at any time and any future outcome. Assumption 2 would hold, for example, under the conditions depicted in Figure 1. Positivity states that the counterfactual treatment strategy of interest has some nonzero probability of actually being followed. Under the assumption that we specify certain predictive models correctly such that their predictions extrapolate to parts of a joint distribution that they were not trained on, positivity is not strictly necessary.
Under assumptions 1-3, for t = m we have simply that i.e. the conditional distribution of the counterfactual is simply the conditional distribution of the observed outcome given patient history and given that treatment follows the strategy of interest. For t > m, things are slightly more complex because we need to adjust for time-varying confounding. With X i:j = X i , . . . , X j for any random variable X, under Assumptions 1-3 the g-formula yields Figure 1: A causal DAG representing a data generating process in which Assumption 2 (sequential exchangeability) holds. This DAG represents a simple two time step process where the outcome is measured only after the final time step. However, its salient property is that all variables influencing treatment (i.e. with arrows directly into treatment) and associated with future outcomes are measured.
It is not generally possible to compute this integral in closed form, but it could be approximated through Monte-Carlo simulation. We repeat the recursive process shown in Algorithm 1 M times. (Here the outcome Y t is without loss of generality deemed to be a variable in the vector L t+1 .) At the end of this process, we have M simulated draws of the counterfactual outcome for each time t = {m, . . . , K}.
For each t, the empirical distribution of these draws constitutes a Monte-Carlo approximation of the counterfactual outcome distribution (2). The sample averages of the draws at each time t are an estimate of the conditional expectations (1) and can serve as point predictions for Y t (Ā m−1 , g m ) in a patient with history H m .
Key to the g-computation algorithm is the ability to simulate from joint conditional distributions p(L t |L t−1 ,Ā t−1 ) of the covariates given patient history at time t. Of course, in practice we do not have prior knowledge of these conditional distributions and need to estimate them from data. Most implementations use generalized linear regression models to estimate the conditional distributions of the covariates. Often, these models do not capture temporal dependencies present in the patient data. We propose the G-Net for this task.

G-Net Design
The proposed G-Net framework depicted in Figure 2 enables the use of sequential deep learning models to estimate conditional distributions p(L t |L t−1 ,Ā t−1 ) of covariates given history at each time point and perform the G- Algorithm 1: G-Computation (One simulation) Continue simulations through time K computation algorithm described in Algorithm 1 to simulate covariates under various treatment strategies. In this setting, without loss of generality, we set Y t as one of the co-variates in L for notational simplicity.
Let L 0 t , . . . , L p−1 t denote p components of the vector L t , where each component L j t could be multivariate. We impose an arbitrary ordering and estimate the conditional distributions of each L j t given all variables preceding it in this ordering. At simulation time, we exploit the basic probability identity . There are at least two reasons to allow for subdivision of the covariates. First, if covariates are of different types (e.g. continuous, categorical, count, etc.), it is difficult to simultaneously simulate from their joint distribution. Second, customizing models for each covariate component can potentially lead to better performance. Figure 2 illustrates this decomposition where at each time point the components are depicted via ordered (grey shaded) boxes that are responsible for the estimation of the various terms needed to compute the conditional distributions. One could set p = 1 and model all covariates simultaneously or at the other extreme set p to be the total number of variables.
The sequential model used in G-Net provides us with estimates of the conditional expectations , we proceed as follows. If L j t is multinomial, its conditional expectation defines its conditional density. If L j t has a continuous density, there are various approaches we might take to simulate from its conditional distribution. Without making parametric assumptions, we could simulate from where j t is a draw from the empirical distribution of the residuals L j t −L j t in a holdout set not used to fit the model parameters used to generateL j t as an estimate of . This method makes the simplifying assumption that the covariate error distribution does not depend on patient history. This is the approach we take in the experiments in this paper, and is depicted in the simulation noise nodes at the top of Figure 2. Alternatively, we might specify a parametric distribution for g. a Gaussian, and directly estimate its parameters by maximum likelihood.
To obtain the estimates for the conditional distribution of covariates, G-Net admits sequential modeling of patients' data using the group decomposition and arbitrary ordering outlined before. As shown in the yellow boxes in Figure 2, at each time t, a representation R t of the patient history can be computed as where Θ represents model parameters learned during training. In its simplest form, r t may just be an identity func-tion passing the covariates. In other configurations (e.g. using the selector in Figure 2), r t can be used to provide abstractions of histories using sequential learning architectures such as RNN. This formulation of r t allows for a great deal of flexibility in how information is shared across variables and time.
Estimates from each of the p covariate groups can then be obtained, as shown in Figure 2, by successive estimation of conditional expectations of covariates. Specifically, the conditional expectation of each L j t+1 , 0 ≤ j < p given the representation of patient history R t and the other variables from time t + 1 that precede it in the arbitrary predefined ordering is estimated by the functions f j t . The complete sequence can be given as follows: where, each of the f j t represents the specialized estimation function for group j and Λ j are the learnable parameters for the same. Depending on modeling choices, f j t can be sequential models specialized for group j or models focusing only on the time step t (e.g. linear models).
Given G-Net parameters, the distribution of the Monte Carlo simulations produced by G-computation Algorithm 1 constitute an estimate of uncertainty about a counterfactual prediction. But this estimate ignores uncertainty about the G-Net parameter estimates themselves. One way to incorporate such model uncertainty would be to fit a Bayesian model and, before each Monte Carlo trajectory simulation in G-computation, draw new network parameters from their posterior distribution. These Monte Carlo draws would be from the posterior predictive distribution of the counterfactual outcome. Bayesian deep learning can be prohibitively computationally intensive, but can be approximated through dropout [24]. First, we fit the RNN component of the G-Net using recurrent dropout, with dropout masks at each layer held constant across time as described in [24]. Let M denote a dropout mask constant across time for each layer of the RNN,Θ(M ) denote estimated RNN parameters with mask M applied, and p M denote the distribution from which M was sampled during training. Then, during g-computation, we add step 0 to the beginning of Algorithm 1 before each simulation: Step  Therefore, the Monte Carlo simulations obtained from g-computation incorporating dropout (i.e. adding step 0 as above) approximate draws from the posterior predictive distribution of the counterfactual outcome.
The parameters, Θ and Λ, are learned by optimizing a loss function forcing the G-Net to accurately estimate covariates L t at each time point t using standard gradient descent techniques. It is to be noted that it is necessary to use teacher-forcing (using observed values of L j t+t as arguments to f t in equation ) to obtain unbiased estimates of f t , as shown in Figure 2.

Simulation Experiments Using CVSim
To evaluate counterfactual predictions, it is necessary to use simulated data in which counterfactual ground truth for outcomes under alternative treatment strategies is known. To this end, we performed experiments on data generated by CVSim [25], a program that simulates the dynamics of the human cardiovascular system. Data Generation: We generated an 'observational' dataset D o under treatment regime g o and two 'counterfactual' datasets D c1 and D c2 under treatment regimes g c1 and g c2 . The data generating processes producing D o and D cj were the same except for the treatment assignment rules. For each j, g cj was identical to g o for the first m − 1 simulation time steps before diverging to a different treatment rule for time steps m to K as illustrated in Figure 3.
A CVSim 6-compartment circulatory model takes as inputs 28 variables that together govern a hemodynamic system. It then deterministically simulates forward in time a set of 25 output variables according to a collection of differential equations (parameterized by the input variables) modeling hemodynamics. Important variables in CVSim include arterial pressure (AP), central venous pressure (CVP), total blood volume (TBV), and total peripheral resistance (TPR). In real patients, physicians observe AP and CVP and seek to keep them above a clinically safe threshold. They do this by intervening on TBV (through fluid administration) and TPR (through vasopressors).
We defined simulated treatment interventions that were designed to mimic the impact of fluids and vasopressors. These simulated interventions alter the natural course of the simulation by increasing either TBV (in the case of the simulated fluids intervention) or TPR (in the case of the simulated vasopressor intervention). We generated patients by randomly initiating baseline inputs (which we hid from our G-Nets to make this a stochastic modeling problem) within plausible physiologic ranges, then using CVSim to simulate covariates forward, intervening according to the relevant treatment strategy at each timestep. Full details of the simulation process can be found in the Appendix.
Under (stochastic) observational treatment strategy g o , the probability of receiving a non-zero vasopressor or fluid dose at a given time increases as MAP and CVP decrease according to a logistic regression function. Given that a dose is non-zero, the exact amount is drawn from a normal distribution with mean inversely proportional to MAP and CVP. Since all drivers of treatment under g o are observed in our data, the sequential exchangeability assumption holds and g-computation may be validly applied. g c1 is similar to g o , except it is a deterministic treatment strategy and the functions linking treatment and dose to covariates have different coefficents. Under g c2 , treatment is always withheld. Again, details are in the Appendix.
Experiment Setup: We set ourselves the task of training a G-Net on D o and using it to predict the trajectories of patients in D cj for time steps m to K for each j. This setup is designed to evaluate the performance of a G-Net in a situation in which we observe data from past patients (D o ) who received usual care (g o ) for K timesteps and would like to predict how a new patient who has been observed for m timesteps would fare were they to follow a different treatment strategy of interest (g cj ) for timesteps m to K. This is a standard use case for counterfactual prediction. D cj provides ground truth for a collection of patients whose trajectories follow just the path we are interested in predicting. Thus, by aggregating predictive performance metrics across simulated patients in D cj we generate measures of the population level performance of our G-Net at the counterfactual prediction task for which it was intended.
As shown in Table 1, we explore two specific criteria: (a) using sequential vs. Identity functions for r t and (b) using sequential models of entire patient history vs. linear models  The best parameters for each model found from grid search are as follows. For M2 (LSTM ) the hidden dimension for the representational LSTM is 30, the hidden dimension for the categorical LSTM is 5, the hidden dimension for the continuous LSTM is 30, and the learning rate is .001. For M3 (LSTM ), the hidden dimension for the categorical LSTM is 10, the hidden dimension for the continuous LSTM is 75, and the learning rate is .005. For M4 (LSTM ), the hidden dimension for the representation LSTM is 30, the hidden dimension for the categorical LSTM is 5, the hidden dimension for the continuous LSTM is 30, and the learning rate is .001. For all models, the batch size used was 64. These parameters were the ones used to achieve the results presented in the experiments and results section.
Evaluation: We evaluate the accuracy (Mean Squared Error -MSE) and calibration of the counterfactual simulations generated by our G-Nets as follows. Say D cj comprises N c trajectories of random variable (L cj K ,Ā cj K ). Given observed history H cj If L t has dimension d, we compute the MSE of counterfactual predictions by a G-Net G in the dataset D c as We assess the calibration of a G-Net G as follows. Given lower and upper quantiles α low and α high , the calibration measures the frequency with which the actual counterfactual covariate L h,cj ti is between the α low and α high quantiles of the M simulations {L h,cj ti (H cj mi , G, k) : k ∈ 1 : M }. If this frequency is approximately α high − α low , then G is well calibrated.
Experiments/Results: Using CVSim, we generated a total of 10,000 trajectories in D o (N o = 10, 000), of which 80% were used for training, and the remaining 20% for validation. For testing, we generated 500 observations in the D cj datasets (N c = 500). We included a total of 18 output variables (including all variables influencing treatment as-signment under g o ) from CVSim to construct D o and D cj ; each trajectory is of length 64 time steps (d=18, K=64). In each D cj , the switching time point m from g o to g c is fixed at 34 for all trajectories (m = 34).
We fit the four models described in Table 1 to the training portion of D o (80%), and used the remaining portion as validation to tune our model hyperparameters. Next, given observed covariate history through 34 time steps and treatment history through 33 time steps of each trajectory in each D cj , we computed the MSE and calibration of the G-Nets' counterfactual predictions for time steps 34 to 64. Overall, for this experiment we found M3 (Identity function for r t and LSTM for the f i functions), performed best over both treatment strategies. Also, a comparative analysis of M1 vs M3 (both with Identity representation function for r t ) and M2 vs M4 (both with LSTM representation functions for r t ) showed that G-Net provide better estimates of conditional distributions by admitting sequential prediction models focusing on the entire patient history compared to prediction models focusing on a single time point.   Table 1) and actual population average trajectories under gc1 for select variables. G-Nets again perform better than the linear model implementation. This is in part because the counterfactual predictive density estimates in these experiments do not take into account uncertainty about model parameter estimates, which could be addressed with dropout as discussed in Section 4.
The G-Net can also be used to estimate population average counterfactual outcomes and treatment effects, quantities sometimes more relevant to policy decisions than individual level counterfactual predictions. Figure 6 displays G-Net (M3 from Table 1) estimates and true values of population average trajectories for select variables under g c1 . Figure 7 shows estimates and true values of the population average treatment effect of following g c1 as opposed to g c2 on select variables. We see that the G-Net does a good job of estimating these population level quantities of interest.

Discussion and Future Work
In the G-Net we have introduced a novel and flexible framework for counterfactual prediction under dynamic treatment strategies through G-computation using sequential deep learning models. We have illustrated several implementations of the G-Net in a realistically complex simulation setting where we had access to ground truth counterfactual outcomes. In particular, we considered alternative approaches to representation learning that either share representations of patient history across predictive tasks or keep them separate. In the particular implementation we considered, shared representations seemed to aid simple linear classifiers but harm LSTMs.
The G-Net framework's flexibility means that there are many other alternative implementations to explore. For example, we might consider alternative architectures for representing patient history, such as attention mechanisms or memory networks.
Another direction of future work is incorporation of prior causal knowledge. For example, if it is known which variables are confounders and which merely predictive of the outcome, we might include weights in the loss function emphasizing faithful representation and prediction of confounding variables compared to non-confounders.

A.1 CVSim
CVSim is a open-source cardiovascular simulator, aiming for education and research purpose developed by Thomas et al [12]. In this work, we focus on CVSim-6C which consists of 6 components, functioning as pulmonary and systemic veins, arteries, and micro-circulations respectively. CVSim-6C is regulated by arterial baroreflex system to simulate Sahs lumped hemodynamic model [12]. The aggregate model is capable of simulating pulsatile waveforms, cardiac output and venous return curves, and spontaneous beat-to-beat hemodynamic variability. In this work, we modified and built on CVSim by adding stochasticcomponents to it for the purposes of evaluating our counterfactual simulators. We call our stochastic simulation engine S-CVSim.

A.2 Inputs of S-CVSim
By varying hemodynamic parameters of CVSim, it can simulate cardiovascular system under various conditions. We, therefore, sample values of a subset of model parameters while initiating simulation at time 0 to obtain ideal distribution of trajectories. These model parameters are listed at Table 2. Note that this does not necessarily mean values of input covariates would not exceed or drop below such range at any time t, where t > 0..

A.3 Outputs of S-CVSim
At each time t, CVSim-6C generates hemodynamic data including vascular resistance and varieties of type of flow, pressure, and volume. In addition to original 25 hemodynamic outputs, we introduce 3 new outputs, including systolic blood pressure, diastolic blood pressure, and mean arterial pressure, bringing the total number of output variables to 28. For this work, we only predict a subset of output covariates highlighted in Table 3. A complete list of output are listed in Table 3.
• Systoblic Blood Pressure (SBP) is defined as the highest measured arterial blood pressure while heart contracting. • Diastolic Blood Pressure (DBP) is defined as the lowest measured arterial blood pressure while heart contracting.

A.4 Simulation Process
To obtain observational data for our intended purpose, we implemented and built two extra events, Disease, S t and Treatment, A t . Note that at each time T , both S t and A t could happen simultaneously but S t also happens before A t . With these stochastic events, A t and S t , D obs and D CF can be simulated and obtained on S-CVSim.

A.4.1 Data Generation, D o , D c
To simulate a patient trajectory under treatment strategy g, we obtain D o and D c with the following algorithm for each individual trajectory. D o , training set, consists of 10k trajectories based on g o . D c , test set, consists of 1k trajectories based on g c .
• Initialize input variables V 1 , ..., V N by drawing from independent uniform distributions using a predefined plausible physiological ranges for each variable as Table 2 suggested.
We introduce the concept S t to simulate hemodynamic instability of cardiovascular system such as sepsis and bleeding at each timestep. S t consists of two events, sepsis and blood loss. In module of S t , we denote that P (S t |L t ) = 0.05, and P (Sepsis|S t ) = P (BloodLoss|S t ) = 0.5; Note that blood loss and sepsis events are mutually exclusive in the simulation process at any time t.

A.4.3 Treatment Simulation, A t
We developed two treatment strategy g o , policy for observational regime and g c , policy for counterfactual policy, to simulate clinical treatment and validate our model. Under any given g, A any t is defined as g(L t ). Similar to disease simulation, A any t is either A 1 t , increasing quantity of total blood volume or A 2 t increasing quantity of arterial resistance. The probability of choosing fluids is equal to vasopressor but will not be administered at the same time. The dosage of A t depends on a subset of L t which indicates hemodynamic balance. More specific, since adequate blood pressure is an important clinical goal [26], we denote mean arterial pressure, MAP, of 65 mmHg and central venous pressure, CVP, of 10 mmHg as target goals. Therefore, We define ∆ map,t ≡ 65 − L map t and ∆ cvp,t ≡ 10 − L cvp t as proxies of how much dosage should be delivered. The fol-lowing section will discuss the difference of A t between g CF and g obs .

A.4.4 Observational Regime, g o
Under g obs , probability and dosage of treatment are denoted as the followings • Probability of treatment, P (A t |L t ) = 1 1+e −x ,where x = C 1 * ∆ map + C 2 * ∆ cvp + C 0 .
The DAG in Figure 8 depicts the causal structure of D o . Note that only observed covariates have arrows pointing into treatment (treatment is actually only a function of CVP, MAP, and ∆ P ), so that there is no unobserved confounding and g-computation may be applied.