G-Transformer for Conditional Average Potential Outcome Estimation over Time

## Abstract

Abstract Estimating potential outcomes for treatments over time based on observational data is important for personalized decision-making in medicine. Yet, existing neural methods for this task either (1) do not perform proper adjustments for time-varying confounders, or (2) suffer from large estimation variance. In order to address both limitations, we introduce the G-transformer ( GT ). Our GT is a novel, neural end-to-end model which adjusts for time-varying confounders, and provides low-variance estimation of conditional average potential outcomes (CAPOs) over time. Specifically, our GT is the first neural model to perform regression-based iterative G-computation for CAPOs in the time-varying setting. We evaluate the effectiveness of our GT across various experiments. In sum, this work represents a significant step towards personalized decision-making from electronic health records.

## 1 Introduction

Causal machine learning has recently garnered significant attention with the aim to personalize treatment decisions in medicine (Feuerriegel et al., ). Here, an important task is to estimate conditional average potential outcomes (CAPOs) from observational data over time (see Fig. 1). Recently, such data has become prominent in medicine due to the growing prevalence of electronic health records (EHRs) (Allam et al., ; Bica et al., ) and wearable devices (Battalio et al., ; Murray et al., ).

Several neural methods have been developed for estimating CAPOs over time. However, existing methods suffer from one of two possible limitations (see Table 1): 1 Methods without proper adjustments for time-varying confounding (Bica et al., ; Melnychuk et al., ; Seedat et al., ) exhibit significant bias, as they do not target the correct estimand. Hence, these methods have irreducible estimation errors irrespective of the amount of available data, which renders them unsuitable for medical applications. 2 Methods that perform proper time-varying adjustments (Li et al., ; Lim et al., ) suffer from large variance. Here, the causal adjustments are based on the estimation of either the distributions of all time-varying covariates, or on inverse propensity weighting at several time steps in the future. While the former is impracticable when granular patient information is available, the latter suffers from strong overlap violations in the time-varying setting. To the best of our knowledge, there is no method that can address both 1 and 2.

To fill the above research gap, we propose the *G-transformer* (GT), a novel, neural end-to-end transformer that overcomes both limitations of existing methods. Our GT builds upon G-computation (Bang & Robins, ; Robins & Hernán, ). However, unlike existing neural models that perform G-computation (Li et al., ), our GT is based on an iterative *regression* scheme and does *not* require estimating any probability distribution. As a result, our GT has two clear strengths: it performs 1 proper adjustments for time-varying confounding, and it has 2 low variance.

Our contributions are three-fold:(^1^11Code is available at [https://github.com/konstantinhess/G_transformer](https://github.com/konstantinhess/G_transformer).)(1) We introduce the first neural end-to-end model for estimating CAPOs over time with 1 proper adjustments for time-varying confounding, while maintaining 2 a low estimation variance. (2) To the best of our knowledge, we are the first to leverage regression-based iterative G-computation for estimating CAPOs over time. (3) We demonstrate the effectiveness of our GT across various experiments. In the future, we expect our GT to help personalize decision-making from patient trajectories in medicine.

## 2 Related Work

Figure: Figure 1: Trajectories with outcomes under *observational* vs. *interventional* treatment sequences.
Refer to caption: x1.png

**Table 1: Overview of key neural methods for estimating CAPOs over time. Methods that perform proper adjustments for time-varying confounding target the correct causal estimand and, therefore, have no infinite-data bias. Methods with low estimation variance neither require estimation of full probability distribution nor inverse propensity weighting.**
|  | CRN | TE-CDE | CT | RMSNs | G-Net | GT (ours) |
| --- | --- | --- | --- | --- | --- | --- |
| 1 Proper adjustments for time-varying confounding | ✗ | ✗ | ✗ | ✓ | ✓ | ✓ |
| 2 Low estimation variance | ✓ | ✓ | ✓ | ✗ | ✗ | ✓ |

Estimating CAPOs in the static setting: Extensive work on estimating potential outcomes focuses on the *static* setting (e.g., Alaa & van der Schaar, ; Frauen et al., ; Johansson et al., ; Louizos et al., ; Melnychuk et al., ; Yoon et al., ; Zhang et al., )). However, observational data such as electronic health records (EHRs) in clinical settings are typically measured *over time* (Allam et al., ; Bica et al., ). Hence, static methods are not tailored to accurately estimate potential outcomes when (i) time series data is observed and (ii) multiple treatments in the future are of interest.

Estimating APOs over time: Estimating average potential outcomes (APOs) over time has a long-ranging history in classical statistics and epidemiology (Lok, ; Robins, ; Rytgaard et al., ; van der Laan & Gruber, ). Popular approaches are the G-methods (Robins & Hernán, ), which include marginal structural models (MSMs) (Robins & Hernán, ; Robins et al., ), structural nested models (Robins, ; Robins & Hernán, ) and the G-computation (Bang & Robins, ; Robins, ; Robins & Hernán, ). G-computation has also been incorporated into neural models (Frauen et al., ). However, these works do not focus estimating CAPOs. Therefore, they are not suitable for personalized decision-making.

Estimating CAPOs over time: In this work, we focus on the task of estimating the heterogeneous response to a sequence of treatments through *conditional average potential outcomes* (CAPOs).(^2^22This is frequently known as *counterfactual prediction*. However, our work follows the potential outcomes framework (Neyman, ; Rubin, ), and we thus use the terminology of CAPO estimation.)
There are some non-parametric methods for this task (Schulam & Saria, ; Soleimani et al., ; Xu et al., ), yet these suffer from poor scalability and have limited flexibility regarding the outcome distribution, the dimension of the outcomes, and static covariate data; because of that, we do not explore non-parametric methods further but focus on neural methods instead.(^3^33Other works are orthogonal to ours. For example, (Hess et al., ; Vanderschueren et al., ) are approaches for informative sampling and uncertainty quantification, respectively. However, they do not focus on the causal structure in the data, and are therefore *not* primarily designed for our task of interest.) Hence, we now summarize key neural methods that have been developed for estimating CAPOs over time (see Table 1). However, these methods fall into two groups with *important limitations*, as discussed in the following:

Limitation 1 proper adjustments: A number of neural methods for estimating CAPOs have been proposed that do not properly adjust for time-varying confounders (Bica et al., ; Melnychuk et al., ; Seedat et al., ). Therefore, they are *biased* as they do not target the correct estimand. Here, key examples are the counterfactual recurrent network (CRN) (Bica et al., ), the treatment effect neural controlled differential equation (TE-CDE) (Seedat et al., ) and the causal transformer (CT) (Melnychuk et al., ). These methods try to account for time-varying confounders through balanced representations. However, balancing was originally designed for reducing finite-sample estimation variance and *not* for mitigating confounding bias (Shalit et al., ). Hence, this is a heuristic and may even introduce another source of representation-induced confounding bias (Melnychuk et al., ). Unlike these methods, our GT performs proper adjustments for time-varying confounders.

Limitation 2 variance: Existing neural methods with proper causal adjustments require estimating full probability distributions at several time steps in the future, which leads to large variance. Prominent examples are the recurrent marginal structural networks (RMSNs) (Lim et al., ) and the G-Net (Li et al., ). Here, the RMSNs leverage MSMs (Robins & Hernán, ; Robins et al., ) and construct pseudo outcomes through inverse propensity weighting (IPW) in order to mimic data of a randomized control trial. However, IPW is particularly problematic in the time-varying setting due to *severe overlap violations*. Hence, IPW can lead to extreme weights and, therefore, large variance. Further, the G-Net (Li et al., ) uses G-computation (Robins, ; Robins & Hernán, ) to adjust for confounding. G-computation is an adjustment that marginalizes over the distribution of time-varying confounders under an interventional sequence of treatments (see Supplement A). However, G-Net proceeds by estimating the *entire distribution of all confounders at several time-steps in the future*. Hence, when granular patient information is available, G-Net is subject to the curse of dimensionality and, therefore, also suffers from large variance. Different to these methods, our GT makes use of regression-based G-computation, which leads to low estimation variance.

Research gap: None of the above neural methods leverages G-computation (Bang & Robins, ; Robins, ) for estimating CAPOs through iterative regressions. Therefore, to the best of our knowledge, we propose the first neural end-to-end model that 1 properly adjusts for time-varying confounders through regression-based iterative G-computation. Hence, our GT yields estimates of CAPOs over time that have 2 low estimation variance.

## 3 Problem Formulation

Setup:
We follow previous literature (Bica et al., ; Li et al., ; Lim et al., ; Melnychuk et al., ) and consider data that consist of realizations of the following random variables: (i) outcomes $Y_{t}\in\mathbb{R}^{d_{y}}$, (ii) covariates $X_{t}\in\mathbb{R}^{d_{x}}$, and (iii) treatments $A_{t}\in\{0,1\}^{d_{a}}$ at time steps ${t\in\{0,\ldots,T\}\subset\mathbb{N}_{0}}$, where $T$ is the time window that follows some unknown counting process. We write ${U_{t:t+\tau}=(U_{t},\ldots,U_{t+\tau})}$ to refer to a specific subsequence of a random variable $U_{t}\in\{Y_{t},X_{t},A_{t}\}$. We further write $\bar{U}_{t}=U_{0:t}$ to denote the full trajectory of $U$ including time $t$. Finally, we write ${\bar{H}_{t+\delta}^{t}=(\bar{Y}_{t+\delta},\bar{X}_{t+\delta},\bar{A}_{t-1})}$ for $\delta\geq 0$, and we let ${\bar{H}_{t}=\bar{H}_{t}^{t}}$ denote the collective history of (i)–(iii).

Estimation task:
We are interested in estimating the *conditional* average potential outcome (CAPO) for a future, interventional sequence of treatments, given the observed history. For this, we build upon the potential outcomes framework (Neyman, ; Rubin, ) for the time-varying setting (Robins & Hernán, ; Robins et al., ). Hence, we aim to estimate the potential outcome ${Y_{t+\tau}[a_{t:t+\tau-1}]}$ at future time $t+\tau$, $\tau\in\mathbb{N}$, for an interventional sequence of treatments $a=a_{t:t+\tau-1}$, *conditionally* on the observed history $\bar{H}_{t}=\bar{h}_{t}$. That is, our objective is to estimate

$$ $\displaystyle\mathbb{E}\left[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h} _{t}\right].$ (1) $$

Identifiability:
In order to estimate the causal quantity in Eq. (1) from observational data, we make the following identifiability assumptions (Robins & Hernán, ; Robins et al., ) that are standard in the literature (Bica et al., ; Li et al., ; Lim et al., ; Melnychuk et al., ; Seedat et al., ): (1) *Consistency:* For an observed sequence of treatments $\bar{A}_{t}=\bar{a}_{t}$, the observed outcome $Y_{t+1}$ equals the corresponding potential outcome $Y_{t+1}[\bar{a}_{t}]$. (2) *Positivity:* For any history $\bar{H}_{t}=\bar{h}_{t}$ that has non-zero probability $\mathbb{P}(\bar{H}_{t}=\bar{h}_{t})>0$, there is a positive probability $\mathbb{P}(A_{t}=a_{t}\mid\bar{H}_{t}=\bar{h}_{t})>0$ of receiving any treatment $A_{t}=a_{t}$, where $a_{t}\in\{0,1\}^{d_{a}}$. (3) *Sequential ignorability:* Given a history $\bar{H}_{t}=\bar{h}_{t}$, the treatment $A_{t}$ is independent of the potential outcome $Y_{t+\delta}[a_{t:t+\delta-1}]$, that is, ${A_{t}\perp Y_{t+\delta}[a_{t:t+\delta-1}]\mid\bar{H}_{t}=\bar{h}_{t}}$ for all $a_{t:t+\delta-1}\in\{0,1\}^{\delta\times d_{a}}$.

Figure: Figure 2: During inference, future time-varying confounders are *unobserved* (here: $(X_{t+1},Y_{t+1})$). In order to estimate CAPOs for an interventional treatment sequence without time-varying confounding bias, proper causal adjustments such as G-computation are required.
Refer to caption: x2.png

G-computation:
Estimating CAPOs without confounding bias poses a non-trivial challenge in the time-varying setting. The issue lies in the complexity of handling future time-varying confounders. In particular, for $\tau\geq 2$ and ${1\leq\delta\leq\delta^{\prime}\leq\tau-1}$, future covariates $X_{t+\delta}$ and outcomes $Y_{t+\delta}$ may affect the probability of receiving certain treatments $A_{t+\delta^{\prime}}$, as we illustrate in Fig. 2. Importantly, the time-varying confounders are *unobserved* during inference time, which is generally known as *runtime confounding* (Coston et al., ). Therefore, in order to estimate the direct effect of an interventional treatment sequence, one needs to adjust for the time-varying confounders. That is, it is in general insufficient to only adjust for the history (Frauen et al., ) via

$$ $\displaystyle\mathbb{E}\left[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h} _{t}\right]\neq\mathbb{E}\left[Y_{t+\tau}\mid\bar{H}_{t}=\bar{h}_{t},A_{t:t+ \tau-1}=a_{t:t+\tau-1}\right].$ (2) $$

As a side note, the problem of *time-varying confounding does not arise for one-step ahead predictions* (i.e., $\tau=1$). Here, under assumptions (i)–(iii), conditioning on the observed history is equivalent to backdoor-adjustments in the static setting.

One way to adjust for time-varying confounders is IPW (Robins & Hernán, ; Robins et al., ), which is leveraged by RMSNs (Lim et al., ). However, as we show in Supplement C, IPW is subject to large variance.
Instead, we leverage G-computation (Bang & Robins, ; Robins, ; Robins & Hernán, ), which provides a rigorous way to account for the time-varying confounders through proper adjustments. Formally, G-computation identifies the causal quantity in Eq. (1) via

$$ $\displaystyle\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h}_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots\mathbb{E}\big{\{} \mathbb{E}[Y_{t+\tau}\mid\bar{H}_{t+\tau-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}] \mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\big{\}}$ (3) $\displaystyle\quad\quad\quad\ldots\bigl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_{ t:t+1}\bigg{]}\bigl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}.$ $$

A derivation of the G-computation formula for CAPOs is given in Supplement A. However, due to the nested structure of G-computation, estimating Eq. (3) from data is challenging.

So far, only G-Net (Li et al., ) has used G-computation for estimating CAPOs in a neural model. For this, G-Net makes a Monte Carlo approximation of Eq. (3) through

$$ $\displaystyle\int_{\mathbb{R}^{d_{x}\times\tau-1}\times\mathbb{R}^{d_{y}\times \tau-1}}\mathbb{E}[Y_{t+\tau}\mid\bar{H}_{t+\tau-1}^{t}=\bar{h}_{t+\tau-1}^{t} ,A_{t:t+\tau-1}=a_{t:t+\tau-1}]$ $\displaystyle\times\prod_{\delta=1}^{\tau-1}p(x_{t+\delta},y_{t+\delta}\mid \bar{h}_{t},{x}_{t+1:t+\delta-1},{y}_{t+1:t+\delta-1},{a}_{t:t+\delta-1}) \mathop{}\!\mathrm{d}({x}_{t+1:t+\tau-1},{y}_{t+1:t+\tau-1}).$ (4) $$

However, Eq. (D) requires estimating the *entire distribution* of all time-varying confounders at several time steps in the future, which leads to large variance. In particular, *all moments* of a $(\tau-1)\times(d_{x}+d_{y})$-dimensional random variable need to be estimated, which leads to estimation of nuisance. We provide more details in Supplement D.

In contrast, our GT does not rely on high-dimensional integral approximation through Monte Carlo sampling. Further, our GT does not require estimating any probability distribution. Instead, it performs *regression-based iterative G-computation* in an end-to-end transformer architecture. Thereby, we perform proper adjustments for time-varying confounding, while providing low-variance estimates of Eq. (3).

## 4 G-transformer

In the following, we present our G-transformer. Inspired by (Bang & Robins, ; Robins, ; Robins & Hernán, ) for APOs, we reframe G-computation for CAPOs over time through recursive conditional expectations. Thereby, we precisely formulate the training objective of our GT through iterative regressions. Importantly, existing approaches for estimating APOs do not estimate potential outcomes on an individual level for a given history $\bar{H}_{t}=\bar{h}_{t}$, because of which they are not sufficient for estimating CAPOs. Therefore, we proceed below by first extending regression-based iterative G-computation to account for the heterogeneous response to a treatment intervention. We then detail the architecture of our GT and provide details on the end-to-end training and inference.

### 4.1 Regression-based iterative G-computation for CAPOs

Our GT leverages G-computation as in Eq. (3) and, therefore, properly adjusts for time-varying confounders in Eq. (1). However, we do not attempt to integrate over the estimated distribution of all time-varying confounders. Instead, one of our main novelties is that our GT performs *iterative regressions* in a neural end-to-end architecture. This allows us to estimate Eq. (1) with low estimation variance.

We reframe Eq. (3) equivalently as a recursion of conditional expectations. Thereby, we can formulate the iterative regression objective of our GT. In particular, our approach resembles an *iterative pseudo-outcome regression*. For this, let

$$ $\displaystyle g_{t+\delta}^{a}(\bar{h}_{t+\delta}^{t})=\mathbb{E}[G_{t+\delta+ 1}^{a}\mid\bar{H}_{t+\delta}^{t}=\bar{h}_{t+\delta}^{t},A_{t:t+\delta}=a_{t:t+ \delta}],$ (5) $$

where the *pseudo-outcomes* are defined as

$$ $\displaystyle G_{t+\tau}^{a}=Y_{t+\tau}$ (6) $$

and

$$ $\displaystyle G_{t+\delta}^{a}=g_{t+\delta}^{a}(\bar{H}_{t+\delta}^{t})$ (7) $$

for $\delta=0,\ldots,\tau-1$. By reformulating the G-computation formula through recursions, the nested expectations in Eq. (3) are now given by

$$ $\displaystyle G_{t+\tau-1}^{a}=\mathbb{E}[Y_{\tau}\mid\bar{H}_{t-1}^{t},A_{t:t +\tau-1}=a_{t:t+\tau-1}],$ (8) $\displaystyle G_{t+\tau-2}^{a}=\mathbb{E}\Big{[}\mathbb{E}[Y_{\tau}\mid\bar{H} _{t-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}]\mid\bar{H}_{t-2}^{t},A_{t:t+\tau-2}= a_{t:t+\tau-2}\Big{]},$ (9) $\displaystyle\ldots$ (10) $$

Hence, the G-computation formula in Eq. (3) can be rewritten as

$$ $\displaystyle g_{t}^{a}(\bar{h}_{t})=\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid \bar{H}_{t}=\bar{h}_{t}].$ (11) $$

We show in the following proposition that our iterative pseudo-outcome regression recovers the CAPOs and thus performs proper adjustments for time-varying confounding.

###### Proposition 1 .

Our regression-based iterative G-computation yields the CAPO in Eq. (1).

###### Proof.

See Supplement B.1.
∎

To further illustrate our regression-based iterative G-computation, we provide two examples in Supplement B.3, where we show step-by-step how our approach adjusts for time-varying confounding.

In order to correctly estimate Eq. (2) for a given history $\bar{H}_{t}=\bar{h}_{t}$ and an interventional treatment sequence $a=a_{t:t+\tau-1}$, all subsequent pseudo-outcomes in Eq. (7) are required.
However, the ground-truth realizations of the pseudo-outcomes $G_{t+\delta}^{a}$ are *not available in the data*. Instead, only realizations of $G_{t+\tau}^{a}=Y_{t+\tau}$ in Eq. (6) are observed during the training. Hence, when training our GT, it alternately generates predictions $\tilde{G}_{t+\delta}^{a}$ of the pseudo-outcomes for $\delta=1,\ldots\tau-1$, which it then uses for learning the estimator of Eq. (5).

Therefore, the training of our GT completes two steps in an iterative scheme: First, it runs a A *generation step*, where it generates predictions of the pseudo-outcomes Eq. (7). Then, it runs a B *learning step*, where it regresses the predictions $\tilde{G}_{t+\delta}^{a}$ for Eq. (7) and the observed $G_{t+\tau}^{a}=Y_{t+\tau}$ in Eq. (6) on the history to update the estimator for Eq. (5). Finally, the updated estimators are used again in the next A *generation step*. This procedure resembles an iterative pseudo-outcome regression. Thereby, our GT is designed to simultaneously A *generate* predictions and B *learn* during the training. Both steps are performed in an end-to-end architecture, ensuring that information is shared across time and data is used efficiently. We detail the architecture as well as training and inference of our GT in the following sections.

###### Proposition 1 .

###### Proof.

### 4.2 Model architecture

Figure: Figure 3: Neural end-to-end architecture of our G-transformer.
Refer to caption: x3.png

We first introduce the architecture of our GT. Then, we explain the iterative prediction and learning scheme inside our GT, which presents one of the main novelties. Finally, we introduce the inference procedure.

Our GT consists of two key components (see Fig. 3): (i) a *multi-input transformer* $z_{\phi}(\cdot)$, and (ii) several *G-computation heads* $\{g_{\phi}^{\delta}(\cdot)\}_{\delta=0}^{\tau-1}$, where $\phi$ denotes the trainable weights. First, the multi-input transformer encodes the entire observed history. Then, the G-computation heads take the encoded history and perform the iterative regressions according to Eq. (5). We provide further details on the transformer architecture in Supplement F. For all $t=1,\ldots,T-\tau$ and $\delta=0,\ldots,\tau-1$, the components are designed as follows:

(i) *Multi-input transformer:* The backbone of our GT is a multi-input transformer ${z_{\phi}(\cdot)}$, which consists of three connected encoder-only sub-transformers $z_{\phi}^{k}(\cdot)$, $k\in\{1,2,3\}$. At time $t$, the transformer $z_{\phi}(\cdot)$ receives data $\bar{H}_{t}=(\bar{Y}_{t},\bar{X}_{t},\bar{A}_{t-1})$ as input and passes them to one corresponding sub-transformer. In particular, each sub-transformer $z_{\phi}^{k}(\cdot)$ is responsible to focus on one particular $\bar{U}_{t}^{k}\in\{\bar{Y}_{t},\bar{X}_{t},\bar{A}_{t-1}\}$ in order to effectively process the different types of inputs. Further, we ensure that information is shared between the sub-transformers, as we detail below. The output of the multi-input transformer are hidden states $Z_{t}^{A}$, which are then passed to the (ii) G-computation heads.

We build upon a state-of-the-art multi-input transformer (Melnychuk et al., ). For each $\bar{U}_{t}^{k}\in\{\bar{Y}_{t},\bar{X}_{t},\bar{A}_{t-1}\}$, the sub-transformer $z_{\phi}^{k}(\cdot)$ first performs a masked *multi-headed self-attention* (Vaswani et al., ) on $\bar{U}_{t}^{k}$, $k\in\{1,2,3\}$, and outputs $Q_{t}^{k}$. Further, our GT builds upon masked *in-between cross-attentions* in order to allow the different sub-transformers to share information. For this, our GT uses $Q_{t}^{k}$ as a query to keys and values $K_{t}^{l}=V_{t}^{l}\in\{Q_{t}^{l}\mid l\in\{1,2,3\}\setminus\{k\}\}$ from both of the other two sub-transformers, respectively. Thereby, we ensure that information is shared between each of the sub-transformers. Then, the outputs of these cross-attention mechanisms are passed through feed-forward networks and averaged. Finally, this yields a hidden state $Z_{t}^{A}$, which encodes all the information that is required for estimating the CAPO. We provide more details in Supplement F.

(ii) *G-computation heads:* The *G-computation heads* $\{g_{\phi}^{\delta}(\cdot)\}_{\delta=0}^{\tau-1}$ are the read-out component of our GT. As input at time $t+\delta$, the G-computation heads receive the hidden state $Z_{t+\delta}^{A}$ from the above multi-input-transformer. Recall that we seek to perform the iterative regressions in Eq. (5) and Eq. (2), respectively. For this, we require estimators of $\mathbb{E}[G_{t+\delta+1}^{a}\mid\bar{H}_{t+\delta},\bar{A}_{t+\delta}]$. Hence, the G-computation heads compute

$$ $\displaystyle\hat{\mathbb{E}}[G_{t+\delta+1}^{a}\mid\bar{H}_{t+\delta},{A}_{t+ \delta}]=g_{\phi}^{\delta}(Z_{t+\delta}^{A},A_{t+\delta}),$ (12) $$

where

$$ $\displaystyle Z_{t+\delta}^{A}=z_{\phi}(\bar{H}_{t+\delta})$ (13) $$

for $\delta=0,\ldots,\tau-1$. As a result, the G-computation heads and the multi-input transformer together give the estimators that are required for the regression-based iterative G-computation. In particular, we thereby ensure that, for $\delta=0$, the last G-computation head $g_{\phi}^{0}(\cdot)$ is trained as the estimator for the CAPO as given in Eq. (2). That is, for a fully trained multi-input transformer and G-computation heads, our GT estimates the CAPO via

$$ $\displaystyle\hat{\mathbb{E}}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h }_{t}]=g_{\phi}^{0}(z_{\phi}(\bar{h}_{t}),a_{t}).$ (14) $$

For this, we present a tailored training procedure in the following.

### 4.3 Iterative training and inference time

We now introduce the iterative training of our GT, which consists of a A *generation step* and a B *learning step*. Then, we show how inference for a given history $\bar{H}_{t}=\bar{h}_{t}$ can be achieved. We provide pseudocode in Supplement G.

Iterative training: Our GT is designed to estimate the CAPO $g_{t}^{a}(\bar{h}_{t})$ in Eq. (2) for a given history $\bar{H}_{t}=\bar{h}_{t}$ and an interventional treatment sequence $a=a_{t:t+\tau-1}$ via Eq. (14). Therefore, the G-computation heads in Eq. (12) require the pseudo-outcomes $\{{G}_{t+\delta}^{a}\}_{\delta=1}^{\tau}$ from Eq. (7) during training. However, they are only available in the training data for $\delta=\tau$. That is, we only observe the factual outcomes $G_{t+\tau}^{a}=Y_{\tau}$.

As a remedy, our GT first predicts the remaining pseudo-outcomes $\{{G}_{t+\delta}^{a}\}_{\delta=1}^{\tau-1}$ in the A *generation step*. Then, it can use these generated pseudo-outcomes and the observed $G_{t+\tau}^{a}$ for learning the network weights $\phi$ in the B *learning step*. In the following, we write $\{\tilde{G}_{t+\delta}^{a}\}_{\delta=1}^{\tau-1}$ for the generated pseudo-outcomes and, for notational convenience, we also write $\tilde{G}_{t+\tau}^{a}={G}_{t+\tau}^{a}$.

A *Generation step:* In this step, our GT generates $\tilde{G}_{t+\delta}^{a}\approx{g}_{t+\delta}^{a}(\bar{H}_{t+\delta}^{t})$ as substitutes for Eq. (7), which are the pseudo-outcomes in the iterative regression-based G-computation. Formally, our GT predicts these via

$$ $\displaystyle\tilde{G}_{t+\delta}^{a}=g_{\phi}^{\delta}(Z_{t+\delta}^{a},a_{t+ \delta}),$ (15) $$

where

$$ $\displaystyle Z_{t+\delta}^{a}=z_{\phi}(\bar{H}_{t+\delta}^{t},a_{t:t+\delta-1 }),$ (16) $$

for $\delta=0,\ldots,\tau-1$. For this, all operations are *detached* from the computational graph. Hence, our GT now has pseudo-outcomes $\{\tilde{G}_{t+\delta}^{a}\}_{\delta=0}^{\tau}$, which it can use in the following B *learning step*. Of note, these generated pseudo-outcomes will be noisy for early training epochs. However, as training progresses, the G-computation heads perform increasingly more accurate predictions, as we explain below.

B *Learning step:* This step is responsible for updating the weights $\phi$ of the multi-input transformer $z_{\phi}(\cdot)$ and the G-computation heads $\{g_{\phi}^{\delta}(\cdot)\}_{\delta=0}^{\tau-1}$. For this, our GT learns the estimator for Eq. (5) via

$$ $\displaystyle\hat{\mathbb{E}}[G_{t+\delta+1}^{a}\mid\bar{H}_{t+\delta}^{t},A_{ t:t+\delta}]=g_{\phi}^{\delta}(Z_{t+\delta}^{A},A_{t+\delta}),$ (17) $$

where

$$ $\displaystyle Z_{t+\delta}^{A}=z_{\phi}(\bar{H}_{t+\delta})$ (18) $$

for $\delta=0,\ldots,\tau-1$. In particular, the estimator is optimized by backpropagating the squared error loss $\mathcal{L}$ for all $\delta=0,\ldots,\tau-1$ and $t=1,\ldots,T-\tau$ via

$$ $\displaystyle\mathcal{L}=\frac{1}{T-\tau}\sum_{t=1}^{T-\tau}\left(\frac{1}{ \tau}\sum_{\delta=0}^{\tau-1}\left(g_{\phi}^{\delta}(Z_{t+\delta}^{A},A_{t+ \delta})-\tilde{G}_{t+\delta+1}^{a}\right)^{2}\right).$ (19) $$

Then, after $\phi$ is updated, we can use the updated estimator in the next A *generation step*.

Here, it is important that for $\delta=\tau$, the pseudo-outcome $\tilde{G}_{t+\tau}^{a}=Y_{t+\tau}$ is *available in the data*. By estimating $Y_{t+\tau}$ with $g_{\phi}^{\tau-1}(Z_{t+\tau-1}^{A},A_{t+\tau-1})$, it is ensured the last G-computation head $g_{\phi}^{\tau-1}(\cdot)$ is learned on a ground-truth quantity. Thereby, the weights of $g_{\phi}^{\tau-1}(\cdot)$ are gradually optimized during training. Hence, the predicted pseudo-outcome

$$ $\displaystyle\tilde{G}_{t+\tau-1}^{a}=g_{\phi}^{\tau-1}(Z_{t+\tau-1}^{a},a_{t+ \tau-1})$ (20) $$

in the next A *generation step* become mores accurate. Therefore, the G-computation head $g_{\phi}^{\tau-2}(\cdot)$ is learned on a more accurate prediction in the following B *learning step*, which thus leads to a better generated pseudo-outcome $\tilde{G}_{t+\tau-2}^{a}$, and so on. As a result, the optimization of the G-computation heads gradually improves from $g_{\phi}^{\tau-1}(\cdot)$ up to $g_{\phi}^{0}(\cdot)$.

Inference at runtime: Finally, we introduce how inference is achieved with our GT. Given a history $\bar{H}_{t}=\bar{h}_{t}$ and an interventional treatment sequence $a=a_{t:t+\tau-1}$, our GT is trained to estimate of Eq. (1) through Eq. (2). For this, our GT computes the CAPO via

$$ $\displaystyle\hat{g}_{t}^{a}(\bar{h}_{t})=\hat{\mathbb{E}}[G_{t+1}^{a}\mid\bar {H}_{t}=\bar{h}_{t},A_{t}=a_{t}]=g_{\phi}^{0}(z_{\phi}(\bar{h}_{t}),a_{t}).$ (21) $$

We summarize this in the following proposition.

###### Proposition 2 .

Our GT estimates the G-computation formula as in Eq. (2) and, therefore, performs proper adjustments for time-varying confounders.

###### Proof.

See Supplement B.2.
∎

###### Proposition 2 .

###### Proof.

## 5 Experiments

We show the performance of our GT against key neural methods for estimating CAPOs over time (see Table 1). Further details (e.g., implementation details, hyperparameter tuning, runtime) are given in Supplement H. We report ablation studies of our GT in Supplement E.1.

### 5.1 Synthetic data

First, we follow common practice in benchmarking for causal inference (Bica et al., ; Li et al., ; Lim et al., ; Melnychuk et al., ) and evaluate the performance of our GT against other baselines on fully synthetic data. The use of synthetic data is beneficial as it allows us to simulate the outcomes under a sequence of interventions, which are unknown in real-world datasets. Thereby, we are able to evaluate the performance of all methods for estimating CAPOs over time. Here, our main aim is to show that our GT is *robust against increasing levels of confounding*.

Setting: For this, we use data based on the pharmacokinetic-pharmacodynamic tumor growth model (Geng et al., ), which is a standard dataset for benchmarking causal inference methods in the time-varying setting (Bica et al., ; Li et al., ; Lim et al., ; Melnychuk et al., ). Here, the outcome $Y_{t}$ is the volume of a tumor that evolves according to the stochastic process

$$ $\displaystyle Y_{t+1}=\Big{(}1+\underbrace{\rho\log\left(\tfrac{K}{Y_{t}} \right)}_{\text{Tumor growth}}-\underbrace{\alpha_{c}c_{t}}_{\text{ Chemotherapy}}-\underbrace{(\alpha_{r}d_{t}+\beta_{r}d_{t}^{2})}_{\text{ Radiotherapy}}+\underbrace{\epsilon_{t}}_{\text{Noise}}\Big{)}\,Y_{t},$ (22) $$

where $\alpha_{c}$, $\alpha_{r}$, and $\beta_{r}$ control the strength of chemo- and radiotherapy, respectively, and where $K$ corresponds to the carrying capacity, and where $\rho$ is the growth parameter. The radiation dosage $d_{t}$ and chemotherapy drug concentration $c_{t}$ are applied with probabilities

$$ $\displaystyle A_{t}^{c},A_{t}^{r}\sim\text{Ber}\left(\sigma\left(\tfrac{\gamma }{D_{\text{{max}}}}(\bar{D}_{15}(\bar{Y}_{t-1}-\bar{D}_{\text{{max}}}/2\right) \right),$ (23) $$

where $D_{\text{{max}}}$ is the maximum tumor volume, $\bar{D}_{15}$ the average tumor diameter of the last $15$ time steps, and $\gamma$ controls the confounding strength. We use the same parameterization as in (Melnychuk et al., ). For training, validation, and testing, we sample $N=1000$ trajectories of lengths $T\leq 30$ each.

Table: Table 2: RMSE on synthetic data based on the tumor growth model with $\tau=2$. Our GT consistently outperforms all baselines. We highlight the relative improvement over the best-performing baseline. Reported: average RMSE $\pm$ standard deviation over five seeds.

We are interested in the performance of our GT for increasing levels of confounding. We thus increase the confounding from $\gamma=10$ to $\gamma=20$. For each level of confounding, we fix an arbitrary intervention sequence and simulate the outcomes under this intervention for testing.

Results: Table 2 shows the average root mean squared error (RMSE) over five different runs for a prediction horizon of $\tau=2$. Of note, we emphasize that our comparison is fair (see hyperparameter tuning in Supplement H.1). We make the following observations:

First, our GT outperforms all baselines by a significant margin. Importantly, as our GT performs proper adjustments for time-varying confounding, it is robust against increasing $\gamma$. In particular, our GT achieves a performance improvement over the best-performing baseline of up to $17.4\%$. Further, our GT is highly stable, as can be seen by low standard deviation in the estimates, especially compared to the baselines. In sum, our GT performs best in estimating the CAPOs, especially under increasing confounding strength.

Second, the 1 baselines that do not perform proper adjustments (i.e., CRN (Bica et al., ), TE-CDE (Seedat et al., ), and CT (Melnychuk et al., )) exhibit large variations in performance and are thus highly unstable. This is expected, as they do not target the correct causal estimand and, accordingly, suffer from the increasing confounding.

Third, the baselines with 2 large-variance (i.e., RMSNs (Lim et al., ) and G-Net (Li et al., )) are slightly more stable than the no-adjustment baselines. This can be attributed to that the tumor growth model has no time-varying covariates $X_{t}$ and to that we are only focusing on $\tau=2$-step ahead predictions, both of which reduce the variance. However, the RMSNs and G-Net are still significantly worse than the estimates provided by our GT.

### 5.2 Semi-synthetic data

Next, we study how our GT performs when (i) the covariate space is *high-dimensional* and when (ii) the *prediction windows* $\tau$ *become larger*. For this, we use semi-synthetic data, which, similar to the fully-synthetic dataset allows us to access the ground-truth outcomes under an interventional sequence of treatments for benchmarking.

Setting: We build upon the MIMIC-extract (Wang et al., ), which is based on the MIMIC-III dataset (Johnson et al., ). Here, we use $d_{x}=25$ different vital signs as time-varying covariates and as well as gender, ethnicity, and age as static covariates. Then, we simulate observational outcomes for training and validation, and interventional outcomes for testing, respectively. Our data-generating process is taken from (Melnychuk et al., ), which we refer to for more details. In summary, the data generation consists of three steps:
(1) $d_{y}=2$ untreated outcomes $\tilde{Y}_{t}^{j}$, $j=1,2$, are simulated according to

$$ $\displaystyle\tilde{Y}_{t}^{j}=\alpha_{s}^{j}\text{B-spline}(t)+\alpha_{g}^{j} g^{j}(t)+\alpha_{f}^{j}f_{Y}^{j}(X_{t})+\epsilon_{t},$ (24) $$

where $\alpha_{s}^{j}$, $\alpha_{g}^{j}$ and $\alpha_{f}^{j}$ are weight parameters, B-spline$(t)$ is sampled from a mixture of three different cubic splines, and $f_{Y}^{j}(\cdot)$ is a random Fourier features approximation of a Gaussian process. (2) A total of $d_{a}=3$ synthetic treatments $A_{t}^{l}$, $l=1,2,3$, are simulated via

$$ $\displaystyle A_{t}^{l}\sim\text{Ber}(p_{t}^{l}),\quad p_{t}^{l}=\sigma\left( \gamma_{Y}^{l}Y_{t-1}^{A,l}+\gamma_{X}^{l}f_{Y}^{l}(X_{t})+b^{l}\right)$ (25) $$

where $\gamma_{Y}^{l}$ and $\gamma_{X}^{l}$ are fixed parameters that control the confounding strength for treatment $A^{l}$, $Y_{t}^{A,l}$ is an averaged subset of the previous $l$ treated outcomes, $b^{l}$ is a bias term, and $f_{Y}^{l}(\cdot)$ is a random function that is sampled from an RFF (random Fourier features) approximation of a Gaussian process. (3) Then, the treatments are applied to the untreated outcomes via

$$ $\displaystyle Y_{t}^{j}=\tilde{Y}_{t}^{j}+\sum_{i=t-\omega^{l}}^{t}\frac{\text {min}_{l=1,\ldots,d_{a}}\mathbbm{1}_{\{A_{i}^{l}=1\}}p_{i}^{l}\beta^{l,j}}{( \omega^{l}-i)^{2}},$ (26) $$

where $\omega^{l}$ is the effect window for treatment $A^{l}$ and $\beta^{l,j}$ controls the maximum effect of treatment $A^{l}$.

We run different experiments for training, testing, and validation sizes of $N=1000$, $N=2000$, and $N=3000$, respectively, and set the time window to $30\leq T\leq 50$. As the covariate space is high-dimensional, we thereby study how robust our GT is with respect to estimation variance. We further increase the prediction windows from $\tau=2$ up to $\tau=6$.

Table: Table 3: RMSE on semi-synthetic data based on the MIMIC-III extract. Our GT consistently outperforms all baselines. We highlight the relative improvement over the best-performing baseline. Reported: average RMSE $\pm$ standard deviation over five seeds.

Results: Table 3 shows the average RMSE over five different runs. Again, we emphasize that our comparison is fair (see hyperparameter tuning in Supplement H). We make three observations:

First, our GT consistently outperforms all baselines by a large margin. The performance of GT is robust across all sample sizes $N$. This is because our GT is based on iterative regressions and, therefore, has a low estimation variance. Further, it is stable across different prediction windows $\tau$. We observe that our GT has a better performance compared to the strongest baseline of up to $26.7\%$. Further, the results show the clear benefits of our GT in high-dimensional covariate settings and for longer prediction windows $\tau$. Here, the performance gain of our method over the baselines is even more pronounced. In addition, our GT is highly stable, as its estimates exhibit the lowest standard deviation among all baselines. In sum, our GT consistently outperforms all the baselines.

Second, 1 baselines that do not perform proper adjustments (i.e., CRN (Bica et al., ), CT (Melnychuk et al., )) tend to perform better than baselines with large variance (i.e., RMSNs (Lim et al., ), G-Net (Li et al., )). The reason is that the former baselines are regression-based and, hence, can better handle the high-dimensional covariate space. They are, however, biased as they do not adjust for time-varying confounders and thus still perform significantly worse than our GT.

Third, baselines with 2 large variance (i.e., RMSNs (Lim et al., ), G-Net (Li et al., )) struggle with the high-dimensional covariate space and larger prediction windows $\tau$. This can be expected, as RMSNs suffer from overlap violations and thus produce unstable inverse propensity weights. Similarly, G-Net suffers from the curse of dimensionality, as it requires estimating a $(d_{x}+d_{y})\times(\tau-1)$-dimensional distribution.

Conclusion: In this paper, we propose the GT, a novel end-to-end transformer that adjusts for time-varying confounders, while providing low-variance estimates of CAPOs. For this, we propose a regression-based learning algorithm that sets our GT apart from existing baselines. Therefore, we expect our GT to be an important step toward personalized medicine with machine learning.

## Appendix A Derivation of G-computation for CAPOs

In the following, we provide a derivation of the G-computation formula (Bang & Robins, ; Robins, ; Robins & Hernán, ) for CAPOs over time. Recall that G-computation for CAPOs is given by

$$ $\displaystyle\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h}_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots\mathbb{E}\big{\{} \mathbb{E}[Y_{t+\tau}\mid\bar{H}_{t+\tau-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}] \mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\big{\}}$ (27) $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}.$ $$

The following derivation follows the steps in (Frauen et al., ) and extends them to CAPOs:

$$ $\displaystyle\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h}_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h}_{t}, {A}_{t}={a}_{t}]$ (28) $\displaystyle=$ $\displaystyle\mathbb{E}[\mathbb{E}\{Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t+1 }^{t},{A}_{t}={a}_{t}\}$ (29) $\displaystyle\qquad\qquad\qquad\qquad\quad\>\>\mid\bar{H}_{t}=\bar{h}_{t},{A}_ {t}={a}_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}[\mathbb{E}\{Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t+1 }^{t},A_{t:t+1}=a_{t:t+1}\}$ (30) $\displaystyle\qquad\qquad\qquad\qquad\quad\>\>\mid\bar{H}_{t}=\bar{h}_{t},{A}_ {t}={a}_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}[\mathbb{E}\{\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid \bar{H}_{t+2}^{t},A_{t:t+1}=a_{t:t+1}]$ (31) $\displaystyle\qquad\qquad\qquad\qquad\qquad\>\>\mid\bar{H}_{t+1}^{t},A_{t:t+1} =a_{t:t+1}\}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\>\>\mid\bar{H}_{t}=\bar{h}_{t},{A} _{t}={a}_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}[\mathbb{E}\{\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid \bar{H}_{t+2}^{t},A_{t:t+2}=a_{t:t+2}]$ (32) $\displaystyle\qquad\qquad\qquad\qquad\qquad\>\>\mid\bar{H}_{t+1}^{t},A_{t:t+1} =a_{t:t+1}\}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\>\>\mid\bar{H}_{t}=\bar{h}_{t},{A} _{t}={a}_{t}]$ $\displaystyle=$ $\displaystyle\ldots$ $\displaystyle=$ $\displaystyle\mathbb{E}[\ldots\mathbb{E}\{\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1} ]\mid\bar{H}_{t+\tau-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}]$ (33) $\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;\,\mid\bar{H}_{t+\tau-2}^{ t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;\,\mid\ldots$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;\,\mid\bar{H}_{t}=\bar{h}_ {t},A_{t}=a_{t}]$ $\displaystyle=$ $\displaystyle\mathbb{E}[\ldots\mathbb{E}\{\mathbb{E}[Y_{t+\tau}\mid\bar{H}_{t+ \tau-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}]$ (34) $\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;\,\mid\bar{H}_{t+\tau-2}^{ t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;\,\mid\ldots$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;\,\mid\bar{H}_{t}=\bar{h}_ {t},A_{t}=a_{t}],$ $$

where Eq. (28) follows from the positivity and sequential ignorability assumptions, Eq. (29) holds due to the law of total probability, Eq. (30) again follows from the positivity and sequential ignorability assumptions, Eq. (31) is the tower rule, Eq. (32) is again due to the positivity and sequential ignorability assumptions, Eq. (33) follows by iteratively repeating the previous steps, and Eq. (34) follows from the consistency assumption.

## Appendix B Regression-based iterative G-computation

### B.1 Unbiased estimand

###### Proposition 1 .

Our regression-based iterative G-computation yields the CAPO in Eq. (1).

###### Proof.

For the proof, we only need to apply the definition of the pseudo-outcomes $G_{t+\delta}^{a}$:

$$ $\displaystyle\mathbb{E}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h}_{t}]$ (35) $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots\mathbb{E}\big{\{} \mathbb{E}[Y_{t+\tau}\mid\bar{H}_{t+\tau-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}] \mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (36) $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots\mathbb{E}\big{\{} \mathbb{E}[G_{t+\tau}^{a}\mid\bar{H}_{t+\tau-1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau -1}]\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (37) $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots\mathbb{E}\big{\{}g_{t +\tau-1}^{a}(\bar{H}_{t+\tau-1}^{t})\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}= a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (38) $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots\mathbb{E}\big{\{}G_{t +\tau-1}^{a}\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (39) $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}\mathbb{E}\bigg{[}\ldots g_{t+\tau-2}^{a}(\bar {H}_{t+\tau-2}^{t})\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_{t:t+1} \bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (40) $\displaystyle=$ $\displaystyle\ldots$ (41) $\displaystyle=$ $\displaystyle\mathbb{E}\bigg{\{}G_{t+1}^{a}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{ t},A_{t}=a_{t}\bigg{\}}$ (42) $\displaystyle=$ $\displaystyle g_{t}^{a}(\bar{h}_{t}),$ (43) $$

where Eq. (B.1) holds due the G-computation formula (see Supplement A).
∎

###### Proposition 1 .

###### Proof.

### B.2 Target of our GT

###### Proposition 2 .

Our GT estimates G-computation formula as in and, therefore, performs proper adjustments for time-varying confounders.

###### Proof.

For the proof, we perform the steps as in Supplement B.1:

$$ $\displaystyle\hat{\mathbb{E}}[Y_{t+\tau}[a_{t:t+\tau-1}]\mid\bar{H}_{t}=\bar{h }_{t}]$ (44) $\displaystyle=$ $\displaystyle\hat{\mathbb{E}}\bigg{\{}\hat{\mathbb{E}}\bigg{[}\ldots\hat{ \mathbb{E}}\big{\{}\hat{\mathbb{E}}[Y_{t+\tau}\mid\bar{H}_{t+\tau-1}^{t},A_{t: t+\tau-1}=a_{t:t+\tau-1}]\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau -2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (45) $\displaystyle=$ $\displaystyle\hat{\mathbb{E}}\bigg{\{}\hat{\mathbb{E}}\bigg{[}\ldots\hat{ \mathbb{E}}\big{\{}\hat{\mathbb{E}}[\tilde{G}_{t+\tau}^{a}\mid\bar{H}_{t+\tau- 1}^{t},A_{t:t+\tau-1}=a_{t:t+\tau-1}]\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2} =a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (46) $\displaystyle=$ $\displaystyle\hat{\mathbb{E}}\bigg{\{}\hat{\mathbb{E}}\bigg{[}\ldots\hat{ \mathbb{E}}\big{\{}g_{\phi}^{\tau-1}(a_{t+\tau-1},z_{\phi}(\bar{H}_{t+\tau-1}, a_{t:t+\tau-2}))\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+\tau-2}=a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (47) $\displaystyle=$ $\displaystyle\hat{\mathbb{E}}\bigg{\{}\hat{\mathbb{E}}\bigg{[}\ldots\hat{ \mathbb{E}}\big{\{}\tilde{G}_{t+\tau-1}^{a}\mid\bar{H}_{t+\tau-2}^{t},A_{t:t+ \tau-2}=a_{t:t+\tau-2}\big{\}}$ $\displaystyle\quad\quad\quad\ldots\biggl{\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_ {t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (48) $\displaystyle=$ $\displaystyle\hat{\mathbb{E}}\bigg{\{}\hat{\mathbb{E}}\bigg{[}\ldots g_{\phi}^ {\tau-2}(a_{t+\tau-2},z_{\phi}(\bar{H}_{t+\tau-2},a_{t:t+\tau-3}))\ldots\biggl {\lvert}\bar{H}_{t+1}^{t},A_{t:t+1}=a_{t:t+1}\bigg{]}\biggl{\lvert}\bar{H}_{t} =\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (49) $\displaystyle=$ $\displaystyle\ldots$ (50) $\displaystyle=$ $\displaystyle\hat{\mathbb{E}}\bigg{\{}\tilde{G}_{t+1}^{a}\biggl{\lvert}\bar{H} _{t}=\bar{h}_{t},A_{t}=a_{t}\bigg{\}}$ (51) $\displaystyle=$ $\displaystyle g_{\phi}^{0}(a_{t},z_{\phi}(\bar{h}_{t})).$ (52) $$

∎

###### Proposition 2 .

###### Proof.

### B.3 Examples

To illustrate how regression-based iterative G-computation works, we apply the procedure to two examples. First, we show the trivial case for one-step ahead predictions, and then, for two-step ahead predictions. Recall that the following only holds under our standard assumptions (i) *consistency*, (ii) *positivity*, and (iii) *sequential ignorability*.

(1) One-step ahead prediction:

This is the trivial case, as there is *no time-varying confounding*. Instead, all confounders are observed in the history. Therefore, we can simply condition on the observed history and resemble the *backdoor-adjustment* from the static setting. Importantly, this is not the focus of our work, but we show it for illustrative purposes:

$$ $\displaystyle\mathbb{E}\big{[}Y_{t+1}[a_{t}]\mid\bar{H}_{t}=\bar{h}_{t}\big{]}$ (53) $\displaystyle\underbrace{=}_{\text{Ass. (ii)+(iii)}}\mathbb{E}\big{[}Y_{t+1}[a _{t}]\mid\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\big{]}$ (54) $\displaystyle\underbrace{=}_{\text{Ass. (i)}}\mathbb{E}\big{[}Y_{t+1}\mid\bar{ H}_{t}=\bar{h}_{t},A_{t}=a_{t}\big{]}$ (55) $\displaystyle\underbrace{=}_{\text{Def. }G_{t+1}^{a}}\mathbb{E}\big{[}G_{t+1}^ {a}\mid\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\big{]}$ (56) $\displaystyle\underbrace{=}_{\text{Def. }g_{t}^{a}}g_{t}^{a}(\bar{h}_{t}).$ (57) $$

(2) Two-step ahead prediction:

Two-step ahead predictions already incorporate all the difficulties that are present for multi-step ahead predictions. Here, we need to account for future time-varying confounders such as $(X_{t+1},Y_{t+1})$ as in Figure 2:

$$ $\displaystyle\mathbb{E}\big{[}Y_{t+2}[a_{t:t+1}]\mid\bar{H}_{t}=\bar{h}_{t} \big{]}$ (58) $\displaystyle\underbrace{=}_{\text{Ass. (ii)+(iii)}}\mathbb{E}\big{[}Y_{t+2}[a _{t:t+1}]\mid\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\big{]}$ (59) $\displaystyle\underbrace{=}_{\text{Law of total prob.}}\mathbb{E}\Big{[} \mathbb{E}\big{[}Y_{t+2}[a_{t:t+1}]\mid\bar{H}_{t+1}^{t},A_{t}=a_{t}\big{]} \mid\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\Big{]}$ (60) $\displaystyle\underbrace{=}_{\text{Ass. (ii)+(iii)}}\mathbb{E}\Big{[}\mathbb{E }\big{[}Y_{t+2}[a_{t:t+1}]\mid\bar{H}_{t+1}^{t},A_{t:t+1}=a_{t:t+1}\big{]}\mid \bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\Big{]}$ (61) $\displaystyle\underbrace{=}_{\text{Ass. (i)}}\mathbb{E}\Big{[}\mathbb{E}\big{[ }Y_{t+2}\mid\bar{H}_{t+1}^{t},A_{t:t+1}=a_{t:t+1}\big{]}\mid\bar{H}_{t}=\bar{h }_{t},A_{t}=a_{t}\Big{]}$ (62) $\displaystyle\underbrace{=}_{\text{Def. }G_{t+2}^{a}}\mathbb{E}\Big{[}\mathbb{ E}\big{[}G_{t+2}^{a}\mid\bar{H}_{t+1}^{t},A_{t:t+1}=a_{t:t+1}\big{]}\mid\bar{H }_{t}=\bar{h}_{t},A_{t}=a_{t}\Big{]}$ (63) $\displaystyle\underbrace{=}_{\text{Def. }g_{t+1}^{a}}=\mathbb{E}\big{[}g_{t+1} ^{a}(\bar{H}_{t+1}^{t})\mid\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\big{]}$ (64) $\displaystyle\underbrace{=}_{\text{Def. }G_{t+1}^{a}}=\mathbb{E}\big{[}G_{t+1} ^{a}\mid\bar{H}_{t}=\bar{h}_{t},A_{t}=a_{t}\big{]}$ (65) $\displaystyle\underbrace{=}_{\text{Def. }g_{t}^{a}}g_{t}^{a}(\bar{h}_{t}).$ (66) $$

## Appendix C Variance of inverse propensity weighting

In this section, we compare two possible approaches to adjust for time-varying confounders: G-computation and inverse propensity weighting (IPW) (Robins & Hernán, ; Robins et al., ), which is leveraged by RMSNs (Lim et al., ).

For a fair comparison of G-computation and IPW, we compare the *variance of the ground-truth pseudo-outcomes* that each method relies on – that is, the $G_{t+\delta}^{a}$ of our GT and the inverse propensity weighted outcomes of RMSNs. Importantly, larger variance of the pseudo-outcomes will directly translate into larger variance of the respective estimator. We find that IPW leads to larger variance, which is why we prefer G-computation in our GT.

###### Proposition 3 .

Pseudo-outcomes constructed via inverse propensity weighting have larger variance than pseudo-outcomes in G-computation.

###### Proof.

To simplify notation, we consider the variance of the pseudo-outcomes in the *static setting*. The analogue directly translates into the time-varying setting.

Let $Y$ be the outcome, $X$ the covariates and $A$ the treatment. Without loss of generality, we consider the potential outcome for $A=1$.

For G-computation, the variance of the pseudo-outcome $g^{1}(X)$ is given by

$$ $\displaystyle\text{Var}[g^{1}(X)]=\text{Var}[\mathbb{E}[Y\mid X,A=1]]$ (67) $\displaystyle=\mathbb{E}\Big{[}\mathbb{E}[Y\mid X,A=1]^{2}\Big{]}-\mathbb{E} \Big{[}\mathbb{E}[Y\mid X,A=1]\Big{]}^{2}$ (68) $\displaystyle=\mathbb{E}\Big{[}\mathbb{E}[Y\mid X,A=1]^{2}\Big{]}-\mathbb{E} \Big{[}Y[1]\Big{]}^{2}.$ (69) $$

For IPW, the variance of the pseudo-outcome is

$$ $\displaystyle\text{Var}\Big{[}\frac{YA}{\pi(X)}\Big{]}=\mathbb{E}\Big{[}\left( \frac{YA}{\pi(X)}\right)^{2}\Big{]}-\mathbb{E}\Big{[}\frac{YA}{\pi(X)}\Big{]}^ {2}$ (70) $\displaystyle=\mathbb{E}\Big{[}\mathbb{E}\Big{[}\frac{Y^{2}A}{\pi^{2}(X)}\mid X \Big{]}\Big{]}-\mathbb{E}\Big{[}Y[1]\Big{]}^{2}$ (71) $\displaystyle=\mathbb{E}\Big{[}\mathbb{E}\Big{[}\frac{Y^{2}\pi(X)}{\pi^{2}(X)} \mid X,A=1\Big{]}\Big{]}-\mathbb{E}\Big{[}Y[1]\Big{]}^{2}$ (72) $\displaystyle=\mathbb{E}\Big{[}\underbrace{\frac{1}{\pi(X)}}_{\geq 1}\mathbb{E }[Y^{2}\mid X,A=1]\Big{]}-\mathbb{E}\Big{[}Y[1]\Big{]}^{2},$ (73) $$

and with

$$ $\displaystyle\mathbb{E}[Y\mid X,A=1]^{2}+\underbrace{\text{Var}[Y\mid X,A=1]}_ {\geq 0}=\mathbb{E}[Y^{2}\mid X,A=1]$ (74) $$

we have that

$$ $\displaystyle\text{Var}\Big{[}\frac{YA}{\pi(X)}\Big{]}\geq\text{Var}[g^{1}(X)].$ (75) $$

Therefore, we conclude that G-computation leads to lower variance than IPW and, hence, our GT has lower variance than RMSNs.

∎

Remarks:

- •
The inverse propensity weight is what really drives the difference in variance between the approaches. Note that, in the time-varying setting, IPW relies on *products of inverse propensities*, which can lead to even more extreme weights for multi-step ahead predictions.
- •
IPW is particularly problematic when there are overlap violations in the data. However, as the input history $\bar{H}_{t}$ in the time-varying setting is very high-dimensional (i.e., $t\times(d_{x}+d_{y})$-dimensional), overlap violations are even more problematic.

###### Proposition 3 .

###### Proof.

## Appendix D Comparison to G-Net

In this section, we compare our iterative regression-based approach to G-computation to the version that is employed by G-Net (Li et al., ).

G-Net makes a Monte Carlo approximation of Eq. (3) through

$$ $\displaystyle\int_{\mathbb{R}^{d_{x}\times\tau-1}\times\mathbb{R}^{d_{y}\times \tau-1}}\mathbb{E}[Y_{t+\tau}\mid\bar{H}_{t+\tau-1}^{t}=\bar{h}_{t+\tau-1}^{t} ,A_{t:t+\tau-1}=a_{t:t+\tau-1}]$ $\displaystyle\times\prod_{\delta=1}^{\tau-1}p(x_{t+\delta},y_{t+\delta}\mid \bar{h}_{t},{x}_{t+1:t+\delta-1},{y}_{t+1:t+\delta-1},{a}_{t:t+\delta-1}) \mathop{}\!\mathrm{d}({x}_{t+1:t+\tau-1},{y}_{t+1:t+\tau-1}).$ (76) $$

For this, G-Net requires estimating the full distribution

$$ $\displaystyle\prod_{\delta=1}^{\tau-1}\mathop{}\!\mathrm{d}p(x_{t+\delta},y_{t +\delta}\mid\bar{h}_{t},{x}_{t+1:t+\delta-1},{y}_{t+1:t+\delta-1},{a}_{t:t+ \delta-1}).$ (77) $$

That is, for $\tau$-step ahead predictions, G-Net estimates a $(\tau-1)\times(d_{x}+d_{y})$-dimensional probability distribution.

We compare the approach of G-Net to to our regression-based G-computation in Table 4.

Table: Table 4: We compare the approach to G-computation of G-Net (Li et al., ) to our regression-based version. For this, we compare the *dimensions of the estimated moments* for each method, respectively. G-Net requires estimating the full distribution of all time-varying confounders in the future. This means that all moments of all time-varying confounders at all time steps in the future need to be estimated. In contrast, our GT *only requires estimation of the first moment of the lower-dimensional target variable*.

## Appendix E Additional results

### E.1 Additional results and ablations

In the following, we report the performance of two ablations: the (A) G-LSTM and the (B) Biased transformer (BT). For this, we show (C) additional results of our GT, the baselines and the two ablations.

(A) G-LSTM: Our first ablation is the G-LSTM. For this, we replaced the transformer backbone $z_{\phi}(\cdot)$ of our GT by an LSTM network. We find that our G-LSTM is highly effective: it outperforms all baselines from the literature while our proposed G-transformer is still superior. This demonstrates that our novel method for iterative regression-based G-computation is both effective and general.

(B) BT: Additionally, we implement a biased transformer (BT). Here, we leverage the same transformer backbone $z_{\phi}(\cdot)$ as in our GT, but we directly train the output heads on the factual data. Thereby, the BT refrains from performing G-computation. We can thus isolate the contribution of the iterative G-computation to the overall performance. Our results show that the BT suffers from significant estimation bias and, therefore, demonstrates that our proper adjustments for time-varying confounders are required for accurate estimates of CAPOs.

(C) Additional results: We report additional results on both (i) fully synthetic data as in Section 5.1 and on (ii) semi-synthetic data as in Section 5.2.

For (i) fully synthetic data, we report the performance of all methods for lower levels of confounding in Figure 4 and additional prediction windows up to $\tau=6$ for fixed level of confounding $\gamma=10.0$ in Figure 5.

For (ii) semi-synthetic data, we report additional prediction windows up to $\tau=12$ for $N=1000$ in Figure 6.

Figure: Figure 4: Synthetic data: We decrease the confounding strength ($\gamma=6,7,8,9$) for $\tau=2$. Additionally, we report previous results of the baselines with the new ablations: G-LSTM and BT. Notably, our G-LSTM has competitive performance, while BT suffers from significant bias. Our *GT remains the strongest method*. We see a similar picture as for Figure 5 and Figure 6: our methods perform the best due to our novel, iterative G-computation.
Refer to caption: x4.png

Figure: Figure 5: Synthetic data: We increase the prediction horizon up to $\tau=6$ for confounding $\gamma=10$. Our G-LSTM and our GT have the overall *best performance on all prediction windows*. The results coincide with our results in Figure 4 and Figure 6; our approach to G-computation leads to the lowest prediction errors. (Please note that decreasing prediction errors for increasing $\tau$ is due to the strong heteroscedasticity of the outcome variable; smaller $\tau$ means that we predict more samples in the test data for very small $t$, where variance is the highest.)
Refer to caption: x5.png

Figure: Figure 6: Semi-synthetic data: We increase the prediction horizon up to $\tau=12$ for $N=1000$ training samples. We further implement two ablations: our G-LSTM and the biased transformer (BT). As in Figure 4 and Figure 5, our G-LSTM almost consistently outperforms the baselines, while the BT has large errors. Our *GT remains the best for all prediction windows*. This shows that our novel approach for G-computation leads to accurate predictions, irrespective of the neural backbone. Further, it shows that proper adjustments are important for CAPO estimation.
Refer to caption: x6.png

### E.2 Quality of pseudo-outcomes

Finally, we provide more insights into the quality of the generated pseudo-outcomes $\tilde{G}_{t+\delta}^{a}$ in Figure 7. Here, we added increasing levels of constant bias to the pseudo-outcomes during training. Our results show that these artificial corruptions indeed lead to a significant decrease in the overall performance of our GT. We therefore conclude that, without artificial corruption, our generated pseudo-outcomes are good estimates of the true nested expectations. Further, this shows that correct estimates of the pseudo-outcomes are indeed necessary for high-quality unbiased estimates. Of note, the quality of the predicted pseudo-outcomes is also directly validated by the strong empirical performance in Section 5.

Figure: Figure 7: During training, we add artificial levels of noise to the pseudo-outcomes of our GT (prediction window $\tau=2$, confounding strength $\gamma=10$ on synthetic data). We see that performance quickly deteriorates. This is expected, as it implies that the pseudo-outcomes generated by our GT are meaningful and important for accurate, unbiased predictions.
Refer to caption: x7.png

## Appendix F Architecture of G-transformer

In the following, we provide details on the architecture of our GT.

Multi-input transformer:
The multi-input transformer as the backbone of our GT is motivated by (Melnychuk et al., ), which develops an architecture that is tailored for the types of data that are typically available in medical scenarios: (i) outcomes $\bar{Y_{t}}\in\mathbb{R}^{d_{y}\times t}$, covariates $\bar{X_{t}}\in\mathbb{R}^{d_{x}\times t}$, and treatments $\bar{A_{t}}\in\{0,1\}^{d_{a}\times t}$. In particular, their proposed transformer model consists of three separate sub-transformers, where each sub-transformer performs *multi-headed self-attention mechanisms* on one particular data input. Further, these sub-transformers are connected with each other through *in-between cross-attention mechanisms*, ensuring that information is exchanged. Therefore, we build on this idea as the backbone of our GT, as we detail below.

Our multi-input transformer $z_{\phi}(\cdot)$ consists of three sub-transformer models $z_{\phi}^{k}(\cdot)$, $k=1,2,3$, where $z_{\phi}^{k}(\cdot)$ focuses on one data input $\bar{U}_{t}^{k}\in\{\bar{Y}_{t},\bar{X}_{t},\bar{A}_{t-1}\}$, $k\in\{1,2,3\}$, respectively.

(1) Input transformations:
First, the data $\bar{U}_{t}^{k}\in\mathbb{R}^{d_{k}\times t}$ is linearly transformed through

$$ $\displaystyle Z_{t}^{k,0}=(\bar{U}_{t}^{k})^{\top}W^{k,0}+b^{k,0}\in\mathbb{R} ^{t\times d_{h}}$ (78) $$

where $W^{k,0}\in\mathbb{R}^{d_{k}\times d_{h}}$ and $b^{k,0}\in\mathbb{R}^{d_{h}}$ are the weight matrix and the bias, respectively, and $d_{h}$ is the number of transformer units.

(2) Transformer blocks:
Next, we stack $j=1,\ldots,J$ transformer blocks, where each transformer block $j$ receives the outputs $Z_{t}^{k,j-1}$ of the previous transformer block $j-1$. For this, we combine (i) *multi-headed self- and cross-attentions*, and (ii) *feed-forward networks*.

(i) *Multi-headed self- and cross-attentions:*
The output of block $j$ for sub-transformer $k$ is given by the *multi-headed cross-attention*

$$ $\displaystyle Z_{t}^{k,j}=\tilde{Q}_{t}^{k,j}+\sum_{l\neq k}{\textrm{MHA}}( \tilde{Q}_{t}^{k,j},\tilde{K}_{t}^{l,j},\tilde{V}_{t}^{l,j}),$ (79) $$

where $\tilde{Q}_{t}^{k,j}=\tilde{K}_{t}^{k,j}=\tilde{V}_{t}^{k,j}$ are the outputs of the *multi-headed self-attentions*

$$ $\displaystyle\tilde{Q}_{t}^{k,j}=Z_{t}^{k,j-1}+\textrm{MHA}({Q}_{t}^{k,j},{K}_ {t}^{k,j},{V}_{t}^{k,j}).$ (80) $$

Here, $\textrm{MHA}(\cdot)$ denotes the multi-headed attention mechanism as in (Vaswani et al., ) given by

$$ $\displaystyle\textrm{MHA}({q},{k},v)=(\textrm{Attention}(q^{1},k^{1},v^{1}), \dots,\textrm{Attention}(q^{M},k^{M},v^{M})),$ (81) $$

where

$$ $\displaystyle\textrm{Attention}(q^{m},k^{m},v^{m})=\textrm{softmax}\left(\frac {q^{m}(k^{m})^{\top}}{\sqrt{d_{qkv}}}\right)v^{m}$ (82) $$

is the attention mechanism for $m=1,\ldots,M$ attention heads. The queries, keys, and values $q^{m},k^{m},v^{m}\in\mathbb{R}^{t\times d_{qkv}}$ have dimension $d_{qkv}$, which is equal to the hidden size $d_{h}$ divided by the number of attention heads $M$, that is, $d_{qkv}=d_{h}/M$. For this, we compute the queries, keys, and values for the *cross-attentions* as

$$ $\displaystyle\tilde{Q}_{t}^{k,m,j}=\tilde{Q}_{t}^{k,j}\tilde{W}^{k,m,j}+\tilde {b}^{k,m,j}\in\mathbb{R}^{t\times d_{qkv}},$ (83) $\displaystyle\tilde{K}_{t}^{k,m,j}=\tilde{K}_{t}^{k,j}\tilde{W}^{k,m,j}+\tilde {b}^{k,m,j}\in\mathbb{R}^{t\times d_{qkv}},$ (84) $\displaystyle\tilde{V}_{t}^{k,m,j}=\tilde{V}_{t}^{k,j}\tilde{W}^{k,m,j}+\tilde {b}^{k,m,j}\in\mathbb{R}^{t\times d_{qkv}},$ (85) $$

and for the *self-attentions* as

$$ $\displaystyle{Q}_{t}^{k,m,j}={Z}_{t}^{k,j-1}{W}^{k,m,j}+{b}^{k,m,j}\in\mathbb{ R}^{t\times d_{qkv}},$ (86) $\displaystyle{K}_{t}^{k,m,j}=Z_{t}^{k,j-1}{W}^{k,m,j}+{b}^{k,m,j}\in\mathbb{R} ^{t\times d_{qkv}},$ (87) $\displaystyle{V}_{t}^{k,m,j}=Z_{t}^{k,j-1}{W}^{k,m,j}+{b}^{k,m,j}\in\mathbb{R} ^{t\times d_{qkv}}.$ (88) $$

where $\tilde{W}^{k,m,j},W^{k,m,j}\in\mathbb{R}^{d_{h}\times d_{qkv}}$ and $\tilde{b}^{k,m,j},\tilde{b}^{k,m,j}\in\mathbb{R}^{d_{q}kv}$ are the trainable weights and biases for sub-transformers $k=1,2,3$, transformer blocks $j=1,\ldots,J$, and attention heads $m=1,\ldots,M$. Of note, each *self- and cross attention* uses relative positional encodings (Shaw et al., ) to preserve the order of the input sequence as in (Melnychuk et al., ).

(ii) *Feed-forward networks:* After the *multi-headed cross-attention* mechanism, our GT applies a feed-forward neural network on each $Z_{t}^{k,j}$, respectively. Further, we apply dropout and layer normalizations (Ba et al., ) as in (Melnychuk et al., ; Vaswani et al., ). That is, our GT transforms the output $Z_{t}^{k,j}$ for transformer block $j$ of sub-transformer $k$ through a sequence of transformations

$$ $\displaystyle\textrm{FF}^{k,j}(Z_{t}^{k,j})=\textrm{LayerNorm}\circ\textrm{ Dropout}\circ\textrm{Linear}\circ\textrm{Dropout}\circ\textrm{ReLU}\circ \textrm{Linear}(Z_{t}^{k,j}).$ (89) $$

(3) Output transformation: Finally, after transformer block $J$, we apply a final transformation with dropout and average the outputs as

$$ $\displaystyle Z_{t}^{A}=\textrm{ELU}\circ\textrm{Linear}\circ\textrm{Dropout}( \frac{1}{3}\sum_{k=1}^{3}Z_{t}^{k,J}),$ (90) $$

such that $Z_{t}^{A}\in\mathbb{R}^{d_{z}}$

G-computation heads: The *G-computation heads* $\{g_{\phi}^{\delta}(\cdot)\}_{\delta=0}^{\tau-1}$ receive the corresponding hidden state $Z_{t+\delta}^{A}$ and the current treatment $A_{t+\delta}$ and transform it with another feed-forward network through

$$ $\displaystyle g_{\phi}^{\delta}(Z_{t+\delta}^{A},A_{t+\delta})=\textrm{Linear} \circ\textrm{ELU}\circ\textrm{Linear}(Z_{t+\delta}^{A},A_{t+\delta}).$ (91) $$

## Appendix G Algorithms for iterative training and inference time

In Algorithm 1, we summarize the iterative training procedure of our GT and how inference is achieved.

Figure: Algorithm 1 Training and inference with GT.

## Appendix H Implementation details

In Supplements H.1 and H.2, we report details on the hyperparameter tuning. Here, we ensure that the total number of weights is comparable for each method and choose the grids accordingly. All methods are tuned on the validation datasets. As the validation sets only consist of *observational data* instead of interventional data, we tune all methods for $\tau=1$-step ahead predictions as in (Melnychuk et al., ). All methods were optimized with Adam (Kingma & Ba, ). Further, we perform a random grid search as in (Melnychuk et al., ).

On average, training our GT on fully synthetic data took $13.7$ minutes. Further, training on semi-synthetic data with $N=1000/2000/3000$ samples took $1.1/2.1/3.0$ hours. This is comparable to the baselines. All methods were trained on $1\times$ NVIDIA A100-PCIE-40GB. Overall, running our experiments took approximately $7$ days (including hyperparameter tuning).

### H.1 Hyperparameter tuning: Synthetic data

Table: Table 5: Hyperparameter tuning for all methods on fully synthetic tumor growth data. Here, ${d_{yxa}=d_{y}+d_{x}+d_{a}}$ is the overall input size. Further, $d_{z}$ denotes the hidden representation size of our GT, the balanced representation size of CRN (Bica et al., ), TE-CDE (Seedat et al., ) and CT (Melnychuk et al., ), and the LSTM (Hochreiter & Schmidhuber, ) output size of G-Net (Li et al., ). The hyperparameter grid follows (Melnychuk et al., ). Importantly, the tuning ranges for the different methods are comparable. Hence, the comparison of the methods in Section 5 is fair.

### H.2 Hyperparameter tuning: Semi-synthetic data

Table: Table 6: Hyperparameter tuning for all methods on semi-synthetic data. Here, ${d_{yxa}=d_{y}+d_{x}+d_{a}}$ is the overall input size. Further, $d_{z}$ denotes the hidden representation size of our GT, the balanced representation size of CRN (Bica et al., ), TE-CDE (Seedat et al., ) and CT (Melnychuk et al., ), and the LSTM (Hochreiter & Schmidhuber, ) output size of G-Net (Li et al., ). The hyperparameter grid follows (Melnychuk et al., ). Importantly, the tuning ranges for the different methods are comparable. Hence, the comparison of the methods in Section 5 is fair.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: See weaknesses above.
soundness: 2
strengths: - The paper provides a clear, detailed description of model architecture.

- The authors conduct extensive experiments with synthetic and semi-synthetic datasets, demonstrating GT's effectiveness under different levels of confounding settings.

- The code is provided.
summary: The paper studies conditional average potential outcome estimation over time. The authors propose a model called G-Transformer for addressing two primary challenges in the field: adjustment for time-varying confounders and low variance in estimates. The authors validate their approach through experiments on synthetic and semi-synthetic data, demonstrating improvements over current baselines.
weaknesses: - The contribution of this paper is somewhat limited. 
  * From the methodology perspective, the proposed method's improvement over existing approaches appears marginal. The primary contribution is exploring the transformer backbone in treatment effect estimation, which has been widely investigated in existing work (static setting [1] and dynamic setting [2]).
  * From the application perspective, the paper claims that their model can help personalize decision-making from patient trajectories in medicine. However, it lacks discussion on how the model's estimates would concretely translate to real-world decision support. Furthermore, the assumptions on which the model relies (e.g., positivity, ignorability) are restrictive and may not hold robustly in practical, real-world scenarios.

- While the paper claims to reduce estimation variance, it lacks a theoretical analysis of variance bounds. Given GT’s objective to minimize estimation variance, a more rigorous theoretical justification or a derivation of variance bounds would significantly strengthen this claim.

-  The paper does not consider the existence of unobserved (and potentially time-varying) confounders, which can introduce biases and affect CAPO estimation in real-world EHR settings.

- The study uses a semi-synthetic dataset derived from MIMIC-III. It would be valuable to see the model's performance on fully real-world datasets to assess its utility in practical clinical decision-making problems.

[1] Zhang, Yi-Fan, et al. "Exploring transformer backbones for heterogeneous treatment effect estimation." arXiv preprint arXiv:2202.01336 (2022).

[2] Melnychuk, Valentyn, Dennis Frauen, and Stefan Feuerriegel. "Causal transformer for estimating counterfactual outcomes." International Conference on Machine Learning. PMLR, 2022.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: Table 2 and 3 have low readability due to small font sizes.

The reviewer thinks Supplementary E.2 regarding the quality of pseudo-outcome is worth being mentioned in the main text
soundness: 4
strengths: The study is motivated and the problem is well-defined by the research gap of the previous studies. Specifically, the authors address the limitations of the previous studies regarding G-computation to estimate CAPOs with neural networks by providing a viable solution of iterative pseudo-outcome regression, although the neural network architecture itself is similar to the existing study (Melnychuk et al., 2022).

The author provides sufficient and proper supplementary materials to explain the details of the study.
summary: This study proposes G-Transformer to estimate conditional average potential outcomes (CAPOs) over time while addressing two major limitations of the previous studies: (1) adjustments for time-varying confounders; (2) and large estimation variance. Particularly, the authors propose a viable solution to G-computation, iterative pseudo-outcome regression, using Transformer.
weaknesses: The reviewer did not find any significant weaknesses.

The reviewer thinks a more principled and diverse way to evaluate the quality of the generated pseudo-outcomes will be desirable. While the experiment conducted in E.2 could empirically show the learning step is sensitive to the noise of the generated pseudo-outcomes, it might not be sufficient to show the quality of the generated pseudo-outcomes as it is titled. An ablation study that separates each step (i.e., another step is frozen when one step is trained) might be useful for a controlled experiment to measure the sensitivity of the performance in terms of how much the model is trained to generate the pseudo-outcomes.

---

confidence: 5
contribution: 1
initial_rating: 3
presentation: 4
questions: My suggestions correspond to the two weaknesses above. I think comparison to prior work that already proposed the combination of transformer and temporal-difference learning/pseudo-outcome regression is absolutely crucial; literature review related to pseudo-outcome regression ([1][3][4] and other related work) should be added, and contribution of this work should be properly and correctly stated.
soundness: 3
strengths: 1. Paper is well written, problem formulation, architecture descriptions and the presentation of learning procedures are all clearly stated and very easy to follow.

2. The results look compeling. Table 2 and 3 showed noticable improvement of the proposed model over the compared benchmarks.
summary: This work proposed a framework to estimate conditional average potential outcome over time using a combination of transformer and iterative pseudo-outcome regression. The framework is tested in experiments on a synthetic dataset based on tumor growth model and a semi-synthetic dataset based on MIMIC-III, and is shown to outperform several benchmark models on these two datasets.
weaknesses: 1. The novelty of this work is the most important concern to me. The combination of transformer and temporal-difference learning (i.e. pseudo-outcome regression as called by the authors) are already proposed in [1] for the estimation of CAPO (on top of that, [1] also went one step further and proposed the adjustment for marginal causal parameters). Besides, the heterogeneous token transformer architecture in [1] appears to be a more clean and intuitive solution as opposed to the proposed architecture in my judgement. The proposed transformer architecture is rather trivial and does not significantly differ from [2]. At the very least, the author should compare their architeture to [1], show experimental results and properly state their contributions compared to prior work.

2. The paper has a very weak description of prior work related to pseudo-outcome regression in my opinion. Not only is the exact same procedure temporal-difference learning proposed before as stated above, it is widely used in causal inference and reinforcement learning known as sequential regression [3] and Q-learning [4]. The omission of these works by the authors made it look like pseudo-outcome regression is a novel contribution of their own.

[1] Shirakawa, Toru, et al. "Longitudinal Targeted Minimum Loss-based Estimation with Temporal-Difference Heterogeneous Transformer." International Conference on Machine Learning. PMLR, 2024.

[2] Melnychuk, Valentyn, Dennis Frauen, and Stefan Feuerriegel. "Causal transformer for estimating counterfactual outcomes." International Conference on Machine Learning. PMLR, 2022.

[3] van der Laan, M. J. and Gruber, S. Targeted Minimum Loss Based Estimation of Causal Effects of Multiple Time Point Interventions. The International Journal of Biostatistics, 8(1), 2012.

[4] Chebotar, Yevgen, et al. "Q-transformer: Scalable offline reinforcement learning via autoregressive q-functions." Conference on Robot Learning. PMLR, 2023.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: - It seems that the GT approach has a clear advantage over existing approaches when $\tau \ge 2$. Could you provide scenarios where estimating CAPO over time ($\tau \ge 2$) has clear clinical benefits?
- Could you quantify the variance of the baseline approaches, for example, using the coefficient of variation?
- Could you provide additional results on the precision in the estimation of heterogeneous effects, average treatment effects, and qualitative analysis measures of over- or underestimation of outcomes?
soundness: 2
strengths: - The paper addresses the important and impactful real-world problem of estimating CAPO over time under time-varying confounding and treatment from observational data.
- The proposed GT, iteratively generates pseudo-outcomes by reformulating the G-computation formula as recursive conditional expectations, which seems interesting.
- Experimental results on synthetic and semi-synthetic datasets demonstrate that the proposed approach outperforms the baselines in terms of RMSE, especially under conditions of high time-varying confounding and longer CAPO estimation horizons.
- The code and data are publicly available, making it easy to reproduce the results.
summary: The paper proposes the G-transformer (GT) approach, which combines the G-computation formula and a multi-input transformer to estimate conditional average potential outcomes (CAPO) over time, given time-varying confounding and treatment observational data. GT iteratively generates pseudo-outcomes by reformulating the G-computation formula as recursive conditional expectations. Experimental results on synthetic and semi-synthetic datasets demonstrate that the proposed approach outperforms the baselines in terms of average root mean squared error (RMSE).
weaknesses: - The novelty is quite limited; the proposed approach could be considered a simple extension of the causal multi-input transformer (Melnychuk et al., 2022) with pseudo-outcome predictions.

- While the paper claims to be the first to use pseudo-outcomes for CAPO estimation over time, several works have already explored the use of pseudo-observations in causal inference and survival analysis, including [1,2,3]. I encourage the author(s) to discuss these works and highlight key differences with their proposed approach.

- The paper asserts that related works, such as G-Net, which properly adjust for time-varying confounding, result in high estimation variance without providing substantial experimental or theoretical results to support this claim. Only RMSE values with error bars are provided (Tables 2 and 3), which seems insufficient given that this is a central claim of the GT approach.

- The experimental results seem underwhelming. Only one semi-synthetic dataset is considered, and RMSE is provided to quantify performance. I encourage the author(s) to explore evaluations tied to meaningful clinical measures, including precision in the estimation of heterogeneous effects, average treatment effects, and qualitative analysis measures of over- or underestimation of outcomes.

- Clarity: The notation and writing need improvement for clarity and readability. The use of both $U$ and $H$ to denote the random variables $\{Y_t, X_t, A_t\}$ does not seem necessary. The paper seems to refer to $H$ throughout, and it's unclear why $U$ was introduced. The definitions for $\bar{Y}$, $\bar{X}$, and $\bar{A}$ are not provided. The recursion could be made more explicit by combining Equations (8) and (9) and providing a simple example, e.g., $\tau = 3$. The paper's use of $\delta$, $\tau$, and $t$ for indexing time needs further simplification. Additionally, the definition for $\tau$ should be provided in the *Setup* section.


**References**
- [1] Andersen et al. (2017), "Causal inference in survival analysis using pseudo-observations", Statistics in medicine.
- [2] Andersen and Perme (2010), "Pseudo-observations in survival analysis", Statistical methods in medical research.
- [3] Chien-Lin et al. (2022), "Causal inference for recurrent event data using pseudo-observations", Biostatistics.

**Minor**
- Table 1: I encourage the author(s) to include citations to related works for improved readability.
- Notation overload:  $\phi$  is used to parameterize both functions $g(\cdot  )$ and $z( \cdot)$.
- Typos: Line 140 should be  $U_{t+\tau}$ remove the 'minus'.
- I encourage the author(s) to increase the font size for tables and figures to be consistent with the main text.
- The definition for $\delta$ in Line 244 is not consistent with that in Line 267.

Meta-review: In this paper, the authors propose a framework to estimate the CAPO in temporal settings based on a G-transformer and pseudo-outcome regression. The paper is well written and sound, providing synthetic and semi-synthetic experiments. The authors provided a revised version of their paper with extensive changes.

However, Reviewer i9bM mentioned a lack of novelty compared to [1]. The authors have argued that the goal of that paper is limited to APO, but the reviewer suggested a misunderstanding of [1], which also performs CAPO. The authors did not engage further. Based on my cursory overview of [1], reviewer i9bM is correct, and I believe that comparing with [1] is a key step in demonstrating the efficiency of the proposed method. In addition, I would suggest extending the experiments to include real-world datasets.

[1] Shirakawa, Toru, et al. "Longitudinal Targeted Minimum Loss-based Estimation with Temporal-Difference Heterogeneous Transformer." International Conference on Machine Learning. PMLR, 2024.

Decision: Reject