Estimating counterfactual treatment outcomes over time in complex multiagent scenarios

The multiagent observational data is denoted as the following $X_{t},A_{t},Y_{t}$ at time stamp $t$. Examples of these variables are shown in Fig. 1. Let $X_{t}$ be the time-dependent covariates of the observational data such that $X_{t}=\{x^{(1)}_{t},...,x^{(n)}_{t}\}$, where the $x^{(i)}_{t}$ denotes the covariates for $i$-th multiagent sample with $K$ agents, and $n$ denotes the number of samples. The relationships between agents are represented by the theory-based computation and GVRNN described in Section III. Although some related papers [5, 7] consider the static covariates $C$, which do not change over time, here we do not explicitly consider $C$ because we can easily formulate and implement our methods to add $C$ by conditioning. The treatment (or intervention) assignments are denoted as $A_{t}=\{a^{(1)}_{t},a^{(2)}_{t},...,a^{(n)}_{t}\}$, where $a^{(i)}_{t}$ denotes the treatments assigned in the $i$-th sample. We consider $a^{i}_{t}\in\{0,1\}$, where $1$ is considered treated whereas $0$ is the control (i.e., a binary treatment setting), and we estimate the effect of the treatment assigned at time stamp $t$ on the outcomes $Y_{t+1}=\{y^{(1)}_{t+1},y^{(2)}_{t+1},...,y^{(n)}_{t+1}\}$ at time stamp $t+1$. Note that in observational data, a multiagent sample can only belong to one group (i.e., either a treated or control group), thus the outcome from the other group is always missing and referred to as counterfactual. To represent the historical sequential data before time stamp $t$, we use the notation $\overline{X}_{t}=\{X_{1},X_{2},...,X_{t-1}\}$ to denote the history of covariates observed before time stamp $t$, and $\overline{A}_{t}$ refers to the history of treatment assignments. Combining all covariates and treatments, we define $\mathcal{H}^{(i)}_{t}=\{\overline{x}^{(i)}_{t},\overline{a}^{(i)}_{t}\}$ as all the historical data collected before time stamp $t$.

We adopt the potential outcomes framework (e.g., [17]) and extended by [18] to account for time-varying treatments. The potential outcome $y^{(i)}_{a_{t}=a,t+1}$ of the $i$-th sample given the historical treatment can be formulated as $y^{(i)}_{a_{t}=a,t+1}=\mathbb{E}[y|x^{(i)}_{t},\mathcal{H}^{(i)}_{t},a_{t}=a]$, where $a=\{0,1\}$. Then the ITE on the temporal observational data is defined as:

Here, the observed outcome $y^{(i)}_{a_{t}=a,t+1}$ under treatment $a$ is called factual outcome, while the unobserved one $y^{(i)}_{a_{t}=1-a,t+1}$ is the counterfactual outcome.

Our estimation of ITE is based on the following standard assumptions [19, 4, 20], and we further extend the assumptions in our scenario to include time-varying observational data.

Assumption II.2 means that, for each time $t$, each treatment has a non-zero probability of being assigned. Besides these two assumptions, much of the existing work is based on the strong ignorability assumption as follows:

Definition 1 means that there are no hidden confounders, i.e., all covariates affecting both the treatment assignment and the outcomes are present in the observational dataset. However, this condition is difficult to be guaranteed in practice especially in real-world observational data (in other words, it is not testable in practice [19, 21]). In this paper, we relax such a strict assumption by acknowledging potential hidden confounders. Our proposed methods can learn the representations of the hidden confounders and eliminate the bias between the treatment assignments and outcomes at each time stamp.

In our approach, the learned representations (denoted by $Z_{t}=\{z^{(1)}_{t},z^{(2)}_{t},...,z^{(n)}_{t}\}$) can be leveraged to infer the unobserved confounders and act as substitutes of hidden confounders. That is, we extend the strong ignorability assumption by considering the existence of hidden confounders $Z_{t}$, which influence the treatment assignment $A_{t}$ and potential outcomes $Y_{t+1}$. Given the hidden confounders $Z_{t}$, the potential outcome variables are independent of the treatment assignment at each time stamp. We aim to learn the representations of hidden confounders $Z_{t}$ for bias elimination based on the following assumptions [8]:

Based on the premise, we study the identification of ITE:

Here we explain VRNN [15] and GNN [16] used in the following section III.

Here, as a main approach in Fig. 3A, we extend a GVRNN [13] for local multiagent locations $x^{l}_{t}$ (i.e., specific for each agent) with theory-based computation. As its variant (e.g., used for the ablation study), a pure data-driven model combining GVRNN and VRNN [15] for global variables $x^{g}_{t}$ (i.e., common for all agents) is also considered. Since the global variables do not usually have the graph structure, VRNN without GNN is suitable. Here we describe the representation learning of hidden confounders.

Our methods predict time-varying covariates, potential outcomes, and treatment for balancing by combining data-driven and theory-based approaches. Here our contributions are to propose the theory-based computation of global variables and the prediction of time-varying covariates. In this subsection, we describe the theory-based computation and the prediction of time-varying covariates, potential outcomes, and treatment for balancing.

Since we consider the binary treatment in this paper, we use a cross-entropy loss for the treatment prediction as follows:

One critical consideration is balancing the representations of treated and control groups, which helps reduce the confounding bias [5, 7, 8] and minimize the upper bound of the outcome inference error [29]. In this paper, we incorporate the adversarial learning-based balancing method using a gradient reversal layer [30] according to the related studies [5, 8]. Figure 4 illustrates our adversarial training process. Let $\rm{GNN_{enc}(\cdot;\theta_{g})}$ be the GNN encoder with parameters $\theta_{g}$. Also, let $\rm{MLP_{a}(\cdot;\theta_{a})}$ and $f^{y}_{theory}(\cdot;\theta_{y})$ be the MLP for treatment and outcome prediction with parameters $\theta_{a}$ and $\theta_{y}$ (if $f^{y}_{theory}$ is a neural network). $\gamma$ and $\lambda$ are the hyperparameters (or weights) of loss functions as described below. We add the gradient reversal layer before the treatment classifier $\rm{MLP_{a}(\cdot;\theta_{a})}$ to ensure the confounder representation distribution of the treated and that of the controlled are similar at the group level [5, 8]. The gradient reversal layer does not change the input during the forward-propagation, but in the back-propagation, reversing the gradient by multiplying by a negative scalar (i.e., $-\lambda\frac{\partial\mathcal{L}_{a}}{\partial\theta_{g}}$ in Fig. 4). Although the model will minimize the loss of the treatment prediction, the adversarial learning process will contribute to the distribution balancing during the prediction of the potential outcome and the covariates.

The loss function for our method is defined as

where $\mathcal{L}_{y}$ is the factual outcome prediction loss, $\mathcal{L}_{gvrnn}$ is the ELBO in GVRNN, $\mathcal{L}_{x}$ is the covariate prediction loss, and $\mathcal{L}_{a}$ is the treatment prediction loss. $\alpha$, $\gamma$, and $\lambda$ are hyperparameters to balance the loss function. In a pure data-driven model, we add $\mathcal{L}_{vrnn}$ to $\mathcal{L}$. To prevent over-fitting, we select the best-performing model using the loss function on the validation set. The sensitivity analysis in the hyperparameters is presented in Appendix I.

To verify our method, we compared the performances to infer the causal effect with those in various baselines using two synthetic datasets with ground truth. We used two types of simulations: autonomous vehicle and biological agent (Boid) simulations. All simulations are rule-based rather than learning-based as described below. As for performance metrics, we first adopted commonly used potential outcome and covariate prediction errors. Potential outcome prediction errors were computed as an absolute error of all simulated outcomes: $L_{outcome}^{(a_{t},t,i,t^{\prime})}=|\hat{y}_{a_{t},t}^{(i,t^{\prime})}-y_{a_{t},t}^{(i,t^{\prime})}|$, where a treatment $a_{t}=\{0,1\}$, time $t=[T_{b},T+1]$, simulated temporal variation of the treatment timing $t^{\prime}\in T_{i}$, and the burn-in period in RNN $T_{b}$. Covariate prediction errors were computed as $L_{2}$ prediction error of all simulated covariates: $L_{covariates}^{(a_{t},t,i,t^{\prime})}=\|\hat{x}_{a_{t},t}^{(i,t^{\prime})}-x_{a_{t},t}^{(i,t^{\prime})}\|_{2}$. We took the average over all variables for evaluation.

We also adopted two widely-used evaluation metrics in causal inference: rooted precision in the estimation of heterogeneous effect (PEHE) [37] $\sqrt{\epsilon^{t}_{PEHE}}=\sqrt{\frac{1}{n}\sum_{i\in n}(\hat{\tau}_{t}^{(i)}-\tau_{t}^{(i)})^{2}}$ and mean absolute error of the average treatment effect (ATE) [73] to measure the quality of the estimated individual treatment effects at different time stamps: $\epsilon^{t}_{ATE}=|\frac{1}{n}\sum_{i\in n}\hat{\tau}_{t}^{(i)}-\frac{1}{n}\sum_{i\in n}\tau_{t}^{(i)}|$. We took the average over all time stamps for evaluation.

Regarding interpretability, we emphasize that there has been no previous work to visualize or interpret future covariates (i.e., multiagent trajectory), which can confirm under what circumstances the intervention is effective. We then illustrated trajectory prediction results in each domain. We compared the customized baselines to predict future covariates to demonstrate the superiority of our approach, which are also quantitatively shown as the covariate losses $L_{covariates}$.

We predict the safe driving distance of the ego car from the starting point as an outcome while driving without any collision with obstacles. That is, the treatment effect of the full observation model is defined as the difference between safe driving distances with and without interventions. As the local covariates $x^{l}_{t}$, we used the position, velocity, pose, and size information of the ego car and obstacles in the 2D map. The maximum number of dynamic obstacles was $120$, but since the obstacles related to the ego car were limited and the number changes over time, we used the nearest $10$ obstacles’ information as the covariates. As the global covariates $x^{g}_{t}$, we used the current driving distance $x_{t}^{dist}$ and (binary) collision $x_{t}^{coll}$ information of the ego car. The theory-based function $f^{x}_{theory}$ mathematically computed $\hat{x}_{t+1}$ using the predicted velocity and pose information, and added the learned positive velocity in intervention. $f^{y}_{theory}$ computed $\hat{y}_{t+1}$ using the predicted $\hat{y}^{\prime}_{t+1}$ as the output of the MLP such that $\hat{y}_{t+1}=(1-x_{t}^{coll})\hat{y}^{\prime}_{t+1}$ based on the definition of the safe driving distance.

Quantitative verification results are shown in Table I. Our full model and its ablated variants show lower prediction errors of the outcome and covariates, and $\sqrt{\epsilon^{t}_{PEHE}}$, and $\epsilon^{t}_{ATE}$ than other baselines. Totally, our full model and the variant without amortized variational inference (TG-CRN) show the best performances, suggesting that the necessity of complex modeling for long-term prediction using this dataset would be smaller than that using the following datasets. In our four models, we found that on this dataset the combination of theory-based computation and GNN worked well in the covariate prediction, but did not in the outcome prediction. Example results of our method are shown in Fig. 5, which can interpret the results. Our method (TGV-CRN) shows better counterfactual prediction with intervention than the baseline (GCRN$+$X). On average, the outcome values with intervention compared to the absence of the intervention (i.e., ITE) were $0.005\pm 0.000$ (km) in our full model, $-0.021\pm 0.000$ in the baseline (GCRN+X), and $0.013\pm 0.001$ in the ground truth. As Fig. 5 shows, the baseline did not model the intervention effects. We expect that this method can be applied to the effect of human intervention in Level 3 autonomous vehicle simulations, which need human intervention. This approach can examine the effect of autonomous control in Level 4 or 5 autonomous vehicle simulations (without human intervention) in a case when human intervention is necessary in Level 3.

In this model, $20$ agents are described by a 2-D vector with a $1$ m/s constant velocity in a 15 $\times$ 15 m boundary square. At each time stamp, a member will change direction according to the positions of all other members based on three zones. The first is the repulsion zone with radius $r_{r}=0.5$ m, in which individuals within each other’s repulsion zone try to avoid each other by swimming in opposite directions. The second is the orientation zone, in which individuals try to move in the same direction; here we set radius $r_{o}=1$ to generate swarming behaviors before the intervention. To generate torus behaviors, we change $r_{o}=4$, which is the intervention in this study. The third is the attractive zone (radius $r_{a}=7.5$ m), in which agents move towards each other and tend to cluster.

To simulate the treatment assignments, we generate factual $20800$ samples ($20000$ training, $400$ validation, and $400$ test datasets). We set $T=14$ and $T_{b}=9$ and randomly pick the intervention point during the intervention period $T_{i}=9,\ldots,13$. The outcome is defined as the mean angular momentum among individuals at time $t+1$. We also created a counterfactual dataset only in the test dataset. Here we created all combinations of treatment points during the intervention period $T_{i}$. As the local covariates $x^{l}_{t}$, we used position, velocity, and directional change of all agents. As the global covariates $x^{g}_{t}$, we used the current mean angular momentum. The theory-based function $f^{x}_{theory}$ mathematically computed $\hat{x}_{t+1}$ using the direction change $d$ at the next step, maximum turn angle $\beta$ as body constraints (for $d$ and $\beta$, see also Appendix J), and attraction rule when agents are far from the center of the group. In addition, we added the orientation rule when agents are in the orientation zone and not in the repulsion zone in $f^{x}_{theory}$. $f^{y}_{theory}$ was replaced with a MLP such that $\hat{y}_{t+1}={\rm{MLP_{y}}}([z^{l}_{t},x^{g}_{t},a_{t}])$.

The results are shown in Table II left. In addition to the four indices in the CARLA experiment, we investigated the estimation error of the best intervention timing $|\operatorname*{arg\,max}_{t^{\prime}}(\hat{y}^{(i,t^{\prime})}_{T+1})-\operatorname*{arg\,max}_{t^{\prime}}(y^{(i,t^{\prime})}_{T+1})|$. The results indicate that our model and its variants with theory-based computation show better prediction performances in covariates and the best intervention timing than all of the baselines. However, the outcome prediction errors in our models were worse than the causal inference baselines (DSW, GCRN, and GCFN+$X$), which may lead to degraded performances in $\sqrt{\epsilon^{t}_{PEHE}}$ and $\epsilon^{t}_{ATE}$ of our models than the baselines. In our four models, all combination of core three components (T, G, and V) did not work well and on this dataset, TG-CRN without VRNN is the best performing model. On average, the ITE values were $0.160\pm 0.010$ in our best model (TG-CRN), $0.035\pm 0.001$ in the baseline (GCRN+X), and $0.091\pm 0.026$ in the ground truth. From this viewpoint (both average and variance), our best model was closer to the ground truth than the baseline. Example interpretable results of our method are shown in Fig. 6. Our best model (TG-CRN) shows better counterfactual covariate prediction with intervention than the baseline (GCRN$+$X). Moreover, the counterfactual prediction of the outcome in our model had variation among various intervention times, whereas that in the baseline did not. We consider that these were important properties in our problem, and improved prediction of the potential outcome was left for future work. We expect that our approach can estimate the effect of an experimenter’s interventions on multi-animal behaviors even if they did not perform the experiment in missing conditions. Our approach is expected to improve the efficiency of experimental procedures for observing desired movements.

Finally, we examined the applicability of our methods to a real-world basketball dataset from the NBA. Data acquisition was based on the contract between the league (NBA) and the company (STATS LLC.), not between the players and us. They are top-level players and then the data was not anonymized. The company was licensed to acquire this data, and it was guaranteed that the use of the data would not infringe on any rights of players or teams. In this study, we used attack sequences from 630 games from the 2015/2016 NBA season (https://www.stats.com/data-science/), which contained the trajectories of 10 players and the ball. We extracted 47,467 attacks (i.e., offensive plays) as samples, which were subsampled at 5 Hz. We separated the dataset into 34,696 pre-training (458 games for training the following classifier of effective attack), 11460 training (154 games, 1/10 of that is used as validation), and 1,305 test samples (18 games) in chronological order. Since scoring predictions are difficult in general (e.g., [76, 11]), we define the attack effectiveness as the outcome by predicting whether the attack is effective or not in the future. This is because evaluating team movements based on scores alone may not provide a holistic view, due to factors such as the shooting skills of individual players. In addition, we predict the attack effectiveness at the next time stamp using the pre-training samples and a logistic regression. The details are provided in Appendix L.

We verified our methods using factual data and provided insights using counterfactual prediction. In model verification using factual data, we set $T=95$ and $T_{b}=85$. For counterfactual predictions, we set $T=105$ and $T_{b}=95$ (at the end of attacks, i.e., a shot or turnover) and predicted all combinations of the counterfactual timing during the intervention period $T_{i}=95,\ldots,98$. As the local covariates $x^{l}_{t}$, we used the position and velocity of all agents including the ball. As the global covariates $x^{g}_{t}$, we used the ball player’s information and areas (see also Appendix L), distances from the nearest defender (about the ball player and other attackers), successful shot probabilities of all attackers, and game and shot clock. The theory-based function $f^{x}_{theory}$ and $f^{y}_{theory}$ computed $\hat{x}^{g}_{t+1}$ and $\hat{y}_{t+1}$ in rule-based manners. In addition, to perform long-term prediction, $f^{x}_{theory}$ also computed the ball and the defender nearest to the next ball player during a counterfactual pass in rule-based manners.

Our verification results using factual data are shown in Table I right. We examined the counterfactual pass effect $\hat{\tau}^{CFP}=\max_{t^{\prime}\in T_{i}}(\max_{t=[T_{b},T]}(\hat{y}^{(i,t^{\prime})}_{T+1}))-y^{(i)}_{T_{b}+1}$ in addition to the outcome and covariate prediction errors. Our methods show better prediction performances in the factual covariates than the fully data-driven baselines. In the factual outcomes, our methods outperformed the most appropriate baseline (GCFN$+$X) and RNN, but only the method without GNN shows competitive performance with DSW and GCRN. One of the possible reasons may be the difficulty in modeling the relationship between the future outcome and current covariates, whereas GNN worked well in only the covariate prediction as the previous work [13]. As additional analysis, we indicate endpoint errors and distributions of long-term covariate prediction for basketball data in Appendix M. Again, the strength of our method is to model the covariates at the next timestep for interpretability of the model. Results in Appendix M that the tendency of the endpoint prediction error for all models is similar to the mean prediction error, in which our approach (TV-CRN and TGV-CRN) show the best performance. In addition, the distribution of the player velocity in the figure of Appendix M shows that our approach without theory-based computation had a wider distribution like ground truth than other models, and those with theory-based computation had less zero-velocity bins like ground truth than other models. Although our approach did not completely model the player’s velocity, we show that they were partially effective for covariate prediction using each component. In $\hat{\tau}^{CFP}$, since all models show positive values, it suggests that there may be a more promising shot opportunity by an extra pass in the shot situation. Figure 7 shows realistic counterfactual prediction from our model to demonstrate interpretable results. Our model completed the counterfactual pass and the nearest defender chase, but the baseline model failed and predicted unrealistic behaviors (e.g., the attackers moved toward the outside of the court). In practice, our approach can estimate the effect of the selection of passes in basketball shot scenarios. We expect that our approach can evaluate the decision-making skills of players in a competitive game. Again, we emphasize these important properties in our problem.

The codes and data we used are provided at https://github.com/keisuke198619/TGV-CRN. This experiment was performed on an Intel(R) Xeon(R) CPU E5-2699 v4 ($2.20$ GHz $\times$ 16) with GeForce TITAN X pascal GPU. For the training of the proposed and baseline models, we used the Adam optimizer [77] with an initial learning rate of $0.0001$ and $20$ training epochs. We set the batchsize to 256. For the hyper-parameters in the loss function, we set $\alpha=0.1,\lambda=0.1$ for all datasets and $\gamma=0.1$ in the Boid and CARLA dataset and $\gamma=1$ in the NBA experiment because the latter dataset was more difficult to predict than the Boid and CARLA experiments.

We compared the performances of our methods to infer ITE with those in the following baselines: a simple RNN using GRU [72], deep sequential weighting (DSW) [7], dynamic networked observational data deconfounder (but for clarity, we change the name into GCRN: graph counterfactual recurrent network) [8].