Causal Contrastive Learning for Counterfactual Regression Over Time

Authors: mouad  el bouchattaoui
myriam  tami
benoit  lepetit
paul-henry  cournède
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 20:34:04

Abstract: Estimating treatment effects over time holds significance in various domains, including precision medicine, epidemiology, economy, and marketing. This paper introduces a unique approach to counterfactual regression over time, emphasizing long-term predictions. Distinguishing itself from existing models like Causal Transformer, our approach highlights the efficacy of employing RNNs for long-term forecasting, complemented by Contrastive Predictive Coding (CPC) and Information Maximization (InfoMax). Emphasizing efficiency, we avoid the need for computationally expensive transformers. Leveraging CPC, our method captures long-term dependencies within time-varying confounders. Notably, recent models have disregarded the importance of invertible representation, compromising identification assumptions. To remedy this, we employ the InfoMax principle, maximizing a lower bound of mutual information between sequence data and its representation. Our method achieves state-of-the-art counterfactual estimation results using both synthetic and real-world data, marking the pioneering incorporation of Contrastive Predictive Encoding in causal inference.

Main paper: 
 introduction: 
 It's vital in real-world applications to estimate potential responses, i.e., responses under hypothetical treatment strategies. Individuals show diverse responses to the same treatment, emphasizing the need to quantify individual response trajectories. This enables personalized interventions, enhancing decision-making efficacy. In medical contexts, precise response estimation enables tailored treatments for patients [2, 70, 49] . This paper focuses on counterfactual regression over time, estimating responses under hypothetical treatment plans based on individual records, including past covariates, responses, and treatment sequences up to the current prediction time [64, 62] . In addressing the challenges of this time-varying setting, we tackle: (1) Time-dependent confounding [55] : confounders influenced by past treatment, impacting subsequent treatments and responses. (2) Selection bias: imbalanced covariate distributions across treatment regimes in observational data, requiring time-aware handling beyond methods in static settings [64, 67, 41] . (3) Long-term dependencies: enduring interdependencies among covariates, treatments, and responses, enabling long-range interactions [15, 54] . Recent advancements in neural networks, such as Recurrent Marginal Structural Networks (RMSNs) [42] , Counterfactual Recurrent Networks (CRN) [7] , and G-Net [39] , have tackled these causal inference challenges. However, their reliance on RNNs limits their ability to capture long-term dependencies. Recent studies [48] propose integrating transformers to better represent temporal dynamics. Rather than viewing this as a limitation of RNNs, we see it as an opportunity to emphasize their strengths. We design specific architectures for counterfactual regression over large horizons, avoiding complex, hard-to-interpret models like transformers. Our approach leverages prediction quality. Additionally, previous works [62, 64, 41, 39, 48] did not consider the role of invertible representation in improving counterfactual regression. Here, we introduce an InfoMax regularization term to make our encoder easier to invert. Appendix C.1 provides a detailed overview of counterfactual regression models.InfoMax Principle The InfoMax principle aims to learn a representation that maximizes MI with its input [43, 5] . Estimating MI for high-dimensional data is challenging, often addressed by maximizing a simple and tractable lower bound on MI [32, 56, 40] . Another approach involves maximizing MI between two lower-dimensional representations of different views of the same input [3, 29, 75, 77] , offering a more practical solution. We adopt this strategy by dividing our process history into two views, past and future, and maximizing a tractable lower bound on MI between them. This encourages a "reconstructable" representation of the process history. To our knowledge, the only work applying an InfoMax approach to counterfactual regression, albeit in static settings, is [16] . They propose maximizing MI between an individual's representation and a global representation, aggregating information from all individuals into a single vector. However, the global representation lacks clarity and interpretability, raising uncertainties about its theoretical underpinnings in capturing confounders. Furthermore, there's a lack of theoretical analysis on how minimizing MI between individual and treatment-predictive representations could yield a treatment-invariant representation. As a novelty, we extend the InfoMax principle to longitudinal data, providing a theoretical guarantee of learning balanced representations. Appendix C.2 discusses self-supervision and MI, with all proofs in Appendix G. Setup In the framework of Potential Outcomes (PO) [65] , and following [63] , we track a cohort of individuals (units) i ∈ {1, 2, . . . , N } over t max time steps. At each time t ∈ {1, 2, . . . , t max }, we observe the following: (1) Discrete treatment W it ∈ W = {0, 1, . . . , K -1}, e.g., in medical contexts, W it may represent treatments like radiotherapy or chemotherapy. (2) Outcome of interest Y it ∈ Y ⊂ R, such as tumor volume. 
problem formulation: 
 Y≤t(3) Time-varying context X it ∈ X ⊂ R dx , containing information about the individual that may influence treatment decisions and outcomes. X it is a d x -dimensional vector of confounders, such as health records or clinical measurements. (4) Static confounders V ∈ V ⊂ R dv , such as gender, which remain constant over time. (5) Partially observed potential outcomes Y it (ω i,≤t ), representing the outcomes that would have been observed for individual i at time t under treatment sequence ω i,≤t = (ω i,1 , . . . , ω i,t ) ∈ W t . We define the history process up to time t + 1 as H t+1 = [V, X ≤t+1 , W ≤t , Y ≤t ], capturing all information prior to the assignment of treatment W t+1 . This history is illustrated in the causal graph shown in Figure 1 .Goal Given a training dataset {H i,t+1 , i = 1, . . . , N } sampled from the empirical distribution P Ht+1 , we address the following causal inference problem: Given a history process H t+1 , how can we efficiently estimate counterfactual responses up to time t + τ (where τ ≥ 1 is the prediction horizon) for a potential treatment sequence ω t+1:t+τ = (ω t+1 , . . . , ω t+τ )? The goal is to estimate the causal quantity: E(Y t+τ (ω t+1:t+τ ) | H t+1 ), i.e., the expected outcome at time t + τ , given the history H t+1 and a sequence of treatments ω t+1:t+τ . We identify this causal quantity from observational data using the assumption of sequential ignorability [62, 42, 39, 48] , which is implicitly assumed in Figure 1 and explicitly discussed in Appendix B.1. This allows us to express the counterfactual as:E(Y t+τ (ω t+1:t+τ ) | H t+1 ) = E (Y t+τ | H t+1 , W t+1:t+τ = ω t+1:t+τ ) .5 Causal CPCWe now prove the statement in the core paper, which asserts that the conditional entropyH(H t | (C h t , C f t )) ≥ 0 is minimized if H t is almost surely a function of (C h t , C f t ).The proof is adapted from [17] .Proposition G.2. If H(A | B) = 0, then A = f (B) almost surely.Proof. For simplicity, suppose A and B are discrete random variables. Assume, by contradiction, that there exists b 0 and two distinct values a 1 and a 2 such that p(a 1 | b 0 ) > 0 and p(a 2 | b 0 ) > 0. Then, the conditional entropy is given by:H(A | B) = - b p(b) a p(a | b) log p(a | b).In particular, we have:H(A | B) ≥ p(b 0 ) (-p(a 1 | b 0 ) log p(a 1 | b 0 ) -p(a 2 | b 0 ) log p(a 2 | b 0 )) > 0.Since -t log t ≥ 0 for 0 ≤ t ≤ 1 and is strictly positive for t not equal to 0 or 1, the conditional entropy H(A | B) = 0 if and only if A is almost surely a function of B.To ensure identifiability in the latent space, we leverage recent advances in causal and disentangled representation learning. Suppose the true data-generating process is given by H t = g(z t ), where z t represents the true latent factors. In the sequential context, we assume that the same function g generates two historical subsequences:H f t = g(z f t ), H h t = g(z h t ).We assume a general dependency of the form:p(z f t | z h t ) = Q(z f t ) Z(z h t ) exp(-d(z f t , z h t )).Here, Φ is an encoder, and we use the InfoMax regularization term as follows:L (Inf oM ax) (Φ, d, B) := -E B log exp(-d(Φ(H f t ), Φ(H h t ))) |B| l=1 exp(-d(Φ(H f l,t ), Φ(H h t ))).According to [47] , under certain conditions, if the encoder f minimizes L (Inf oM ax) , then h = g • f is a scaled permutation matrix. This result suggests that when the encoder achieves a minimizer for L (Inf oM ax) , the encoder function f closely approximates an invertible transformation of g. From a causal inference perspective, if Y it (ω it ) ⊥ ⊥ W it | H it and H it = g(Z it ), then an invertible function g • f ensures that:Y it (ω it ) ⊥ ⊥ W it | g • f (H it ). Thus, Y it (ω it ) ⊥ ⊥ W it | g(C it) and since g is invertible, we have:Y it (ω it ) ⊥ ⊥ W it | C it .This demonstrates that the representation C it retains the essential independence structure, facilitating accurate counterfactual inference.To prove the Theorem 5.4, we first prove the following lemma and proposition.Lemma G.3. Let Φ be a fixed representation function. Given that q(W t+1 | Φ(H t )) is the conditional likelihood of observing the treatment W t+1 , denote the probability of observing each treatment value as q j = q(Φ(H t )) := q(W t+1 = j | Φ(H t )) for j ∈ {0, 1, . . . , K -1}. Then, the optimal treatment prediction function is such thatEQUATIONProof. For a fixed representation Φ, finding the optimal treatment probabilities amounts to solving the following optimization problem:EQUATIONFirst, we write the likelihood q(W t+1 | Φ(H t )) using the conditional probabilities q j (Φ(H t )).q(W t+1 | Φ(H t )) = K-1 j=0 q j (Φ(H t )) 1 {W t+1 =j}Then, the treatment likelihood can be written asE P(Φ(Ht),Wt+1) [log q(W t+1 | Φ(H t ))] = E P(Φ(Ht),Wt+1) K-1 l=0 log(q l (Φ(H t )))1 {Wt+1=j} = K-1 l=0 log(q l (Φ(H t ))1 {Wt+1=j} p(W t+1 | Φ(H t ))p(Φ(H t ))dW t+1 dΦ(H t ) = K-1 l=0 log(q l (Φ(H t ))p(W t+1 = l | Φ(H t ))p(Φ(H t ))dΦ(H t ) = K-1 l=0 log(q l (Φ(H t ))p(Φ(H t ) | W t+1 = l)p(W t+1 = l)dΦ(H t )Let's denote α l = p(W t+1 = l), the marginal probability of observing the l-th treatment regime, andp Φ l (H t ) = p(Φ(H t ) | W t+1 = l) with a corresponding probability distribution P Φ l . We intend to maximize point-wise the objective in Eq. ( 21). Plugging the latter formulation of the conditional likelihood in Eq. ( 21) and writing the Lagrangian function, we getEQUATIONComputing the gradient w.r.t q l (Φ(H t )) for l ∈ {0, 1, . . . , K -1} and setting to zero, we haveEQUATIONThen, by the equality constraint, we find that λ = -K-1 l=0 α l p Φ j (H t ). Proposition G.4. Let Φ be a fixed representation function. The I CLU B objective when the treatment prediction function is optimal (i.e. q = q * ) has the following form:EQUATIONProof. First, recall thatI CLUB (Φ(H t ), W t+1 ; q * ) = E P (Φ(H t+1 ),W t+1 ) [log q * (W t+1 | Φ(H t+1 ))]-E P Φ(H t+1 ) E P W t+1 (log q * (W t+1 | Φ(H t+1 ))) I CLUB (Φ(H t ), W t+1 ; q * ) = A -B Let's detail A and B separately, A = K-1 j=0 α j log(q l, * (Φ(H t ))p Φ j (H t )dΦ(H t ) = K-1 j=0 α j log( α j p Φ j (H t ) K-1 l=0 p Φ l (H t )α l )p Φ j (H t )dΦ(H t ) = K-1 j=0 α j log( p Φ j (H t ) K-1 l=0 p Φ l (H t )α l )p Φ j (H t )dΦ(H t ) + log(α j )α j = K-1 j=0 α j D KL (P Φ j || K-1 l=0 α l P Φ l ) + K-1 j=0 log(α j )α jFinally, we can writeEQUATIONFor the remaining term B, we haveB = E P Φ(H t ) E P W t+1 (log q * (W t+1 | Φ(H t+1 ))] = K-1 j=0 E P Φ(H t ) E P W t+1 log(q j (Φ(H t )))1 {Wt+1=j} = K-1 j=0 E P Φ(H t ) α j log(q j (Φ(H t ))) = K-1 j=0 α j log α j p Φ j (H t ) K-1 l=0 p Φ l (H t )α l p(Φ(H t ))dΦ(H t ) = K-1 j=0 α j log p(Φ(H t )) K-1 l=0 p Φ l (H t )α l p(W t+1 = j | Φ(H t )) p(W t+1 = j) p(Φ(H t ))dΦ(H t ) -H(W t+1 ) = K-1 j=0 α j log p(Φ(H t )) K-1 l=0 p Φ l (H t )α l =0 p(Φ(H t ))dΦ(H t ) + K-1 j=0 α j log p(W t+1 = j | Φ(H t )) p(W t+1 = j) p(Φ(H t ))dΦ(H t ) -H(W t+1 )The final form of B is thereforeEQUATIONThe proposition follows immediately from Equations ( 25) and (26) .Proof. (Theorem 5.4) Since by lemma G. In this section, we present a detailed overview of the training procedure for Causal CPC. Initially, we train the Encoder using only the contrastive terms, as outlined in Algorithm 1. Our primary objective is to ensure that, for each time step t, the process history H t is predictive of future local features Z t . However, calculating the InfoNCE loss for a batch across all possible time steps t = 0, . . . , t max can be computationally demanding.To address this, we adopt a more efficient approach by uniformly sampling a single time step t per batch. Subsequently, the corresponding process history H t is contrasted. The sampled H t is then employed as input for the InfoMax objective and randomly partitioned into future H f t and past H h t sub-processes. The decoder is trained while taking the encoder as input (Algorithm 2), utilizing a lower learning rate compared to the untrained part of the decoder. It is trained autoregressively and without teacher forcing. This implies that for each time step t, our GRU-based decoder should predict the future sequence of treatments Ŷt+1:t+τ with its hidden state initialized to the representation Φ t of the historical process up to time t.Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
representation learning: 
 Contrastive Predictive Coding We employ contrastive learning to efficiently represent the process history H t . Causal forecasting over multiple time horizons requires representations that capture variability in H t . For short-term predictions, local features and smooth signal variations are critical, while long-term predictions rely on capturing global structures and long-term dependencies, as shared information between history and future points diminishes. To achieve this, we learn a representation of H t that predicts future components over multiple time steps. For each horizon j = 1, . . . , τ , future components are defined as U t+j = [V, X t+j , W t+j-1 , Y t+j-1 ] ∈ U ⊂ R (dv+dx+K+1) . First, local features are extracted by encoding [V, X t , W t-1 , Y t-1 ] into Z t = Φ θ1 ([V, X t , W t-1 , Y t-1 ]). Then, the full process history H t is summarized into a context representation C t , given by an autoregressive model: Φ ar θ2 (Z ≤t ) = C t , where Φ ar θ2 is implemented with a GRU [14] . This results in the representation function: Φ θ1,θ2 (H t ) = C t . To train the model, we use a contrastive loss that encourages the context C t to predict future local features Z t+1 , . . . , Z t+τ while distinguishing them from the features of other individuals. This is done by minimizing the InfoNCE loss L (Inf oN CE) j for each horizon j:EQUATION)where B is a batch containing individual histories, and Γ j is a weight matrix. The discriminator T j (., .) classifies the correct future feature among negative samples from other individuals:EQUATIONIn practice, C t predicts Ẑt+j , and prediction quality is measured by the dot product Z ⊤ t+j Ẑt+j . Minimizing the InfoNCE loss L (Inf oN CE) j provides a lower bound on the MI between the context and future features I(U t+j , C t ) [52] :EQUATIONFor multiple forecasting horizons j = 1, 2, . . . , τ , we learn long-term dependencies by minimizing the InfoNCE loss across all horizons:EQUATIONThus, minimizing L CP C maximizes the shared information between the context and future components as shown in Eq. ( 5), pushing the model to capture the global structure of the process over large horizons-crucial for counterfactual regression across multiple time steps:EQUATIONInfoMax Principle We introduce a regularization term to make the context representation of the process history H t "reconstructable." We leverage the InfoMax principle to maximize the MI between H t and the context C t . However, we avoid computing the contrastive loss between C t and H t for two main reasons. First, H t is a high-dimensional sequence, making the loss computation very demanding. Secondly, we are still interested in incorporating inductive bias toward capturing global dependencies, this time by pushing any subsequence to be predictive of any future subsequence within H t . Hence, we divide the process history into two non-overlapping views, H h t := U 1:t0 , H f t := U t0+1:t representing a historical subsequence and a future subsequence within the process history H t , with t 0 randomly chosen per batch. We then maximize the MI between the representations of these views, C h t and C f t , resulting in a lower bound to the InfoMax objective as formulated below:Proposition 5.1. I(C h t , C f t ) ≤ I(H t , (C h t , C f t )). We provide an intuitive discussion of the inequality by providing an exact writing of the gap in 5.1: Theorem 5.2.EQUATIONBoth terms on the RHS of Eq. ( 6) are positive, providing an alternative proof to Proposition 5.1.When equality holds, it impliesI(H t ; C f t | C h t ) = 0, indicating that H t is independent of C f t given C h t .This suggests that C h t retains sufficient information from H t that is predictive of C f t . The symmetry of MI also leads to the occurrence of the second term on the RHS when conditioning on C f t . The equality in Proposition 5.1 impliesP C h t |ht = P C h t |c f t, suggesting that C f t efficiently encodes its subsequence while sharing maximum information with C h t .By considering the proposed variant of the InfoMax principle, we can compute a contrastive bound to I(C h t , C f t ) more efficiently, as the random vectors reside in a low-dimensional space thanks to the encoding. We define a contrastive loss using InfoNCE similar to Eq. (1):EQUATIONWe use a non-linear discriminator T η parametrized by η (detailed in Appendix I). The representation of the past subsequence C h t is mapped to a prediction of the future subsequence Ĉft := F η (C h t ) and T η = C f ⊤ tt . Theorem 5.2 and Proposition 5.1 justify using the loss in Eq. ( 7) by showing that our InfoMax simplification provides a valid lower bound. Thus, the contrastive loss in Eq. ( 7) is valid, as:EQUATIONThe "mental model" behind our regularization term comes from the MI,I(H t , (C h t , C f t )) = H(H t )- H(H t | (C h t , C f t )), which can be written using entropy. Since the entropy term is constant and parameter-free, minimizing the conditional entropyH(H t | (C h t , C f t )) ≥ 0 ensures that H t is almost surely a function of (C h t , C f t ) (Appendix G.4, Proposition G.2).When MI is maximized, the theoretical existence of such a function suggests that the learned context C t can decode and reconstruct H t . Beyond the idea of reconstruction, it was shown that the InfoNCE objective implicitly learns to invert the data's generative model under mild assumptions [86] . Recent works [18, 44] extend this insight to multi-modal settings, which can reframe our InfoMax problem: H h t and H f t can be seen as two coupled modalities, allowing us to identify latent generative factors up to some mild indeterminacies (e.g rotations, affine mappings). We plan to extend multi-modal causal representation learning to the longitudinal setting, where we anticipate minimizing our InfoMax objective, in the limit of infinite data, will effectively invert the data generation process up to a class of indeterminacies that we conjecture to be broader and under weaker assumptions than those in current causal representation learning literature, given our focus is on causal inference rather than the identification of causal latent variables. We initiate a formal basis for this claim in Appendix G.5.Motivation Our goal is counterfactual regression, specifically estimating E(Y t+τ (ω t+1:t+τ ) | H t+1 ). For simplicity, with τ = 1, we estimate the potential outcome for a given treatment W t+1 = ω t+1 , where W t+1 ∈ {0, 1, . . . , K -1}, expressed as E(Y t+1 (ω t+1 ) | H t+1 ), which under standard assumptions is identified as:E(Y t+1 (ω t+1 ) | H t+1 ) = E(Y t+1 | H t+1 , W t+1 = ω t+1 ).The RHS can be estimated from data as E(Y t+1 | H t+1 , W t+1 ) = f (H t+1 , W t+1 ). Since only one treatment is observed per individual at each time step, W i,t+1 = ω i,t+1 , our model f generates counterfactual responses by switching treatments f (h i,t+1 , ω ′ t+1 ), where ω ′ t+1 ̸ = ω t+1 (e.g., chemotherapy vs. radiotherapy). The challenge is that H t+1 and W t+1 are not independent, introducing potential bias in counterfactual estimation [61] , leading to covariate shift or selection bias. To address this, we learn a representation Φ(H t+1 ) that enforces distributional balance during decoding.Setup To mitigate selection bias, we leverage the context representation C t of H t and introduce two sub-networks: one for response prediction and one for treatment prediction, both using a mapping of the context representation:Φ t = SELU(Linear(C t )) = Φ θ R (H t ),where SELU represents the Scaled Exponential Linear Unit [37] , and θ R denotes all parameters of the representation learner, i.e., θ R = [θ 1 , θ 2 ]. Following [7, 48] , our objective is to learn a representation that accurately predicts outcomes while remaining distributionally balanced across all possible treatment choices W t = 0, 1, . . . , K -1. To achieve this, we frame the problem as an adversarial game: one network learns to predict the next treatment from the representation, while a regularization term ensures that the representation is non-predictive of the treatment. L Y (θ R , θ Y ) = -log p θ Y (y t+1:t+τ | Φ t , ω t+1:t+τ ) = - τ j=1 log p θ Y (y t+j | y t+1:t+j-1 , Φ t , ω t+1:t+j ).We denote I j t := [Y t+1:t+j-1 , Φ t , W t+1:t+j ] and assume a Gaussian distribution for the conditional responsesY t+j | I j t ∼ N (G Y (I j t ), σ 2 ), where G Y models the mean of the conditional response (see right side of Figure 2 ). We set σ = 0.05 throughout our experiments. The response sequence is estimated autoregressively using a GRU-based decoder without teacher forcing [81] to ensure model training's consistency with testing in real-world scenarios (Figure 2 and Algorithm 2 in Appendix H).Treatment prediction We learn a treatment prediction sub-network parameterized by θ W that takes as input the representation Φ t+1 and predicts a distribution q θ W (ω t+1 | Φ t+1 ) over the treatment W t+1 by minimizing the negative log-likelihood, L W = -log q θ W (ω t+1 | Φ t+1 ). To assess the quality of the representation in predicting the treatment, the gradient from L W only updates the treatment network parameters θ W and is not backpropagated through the response of the parameters for the representation Φ t+1 (Algorithm 2, Appendix H).Adversarial learning To create an adversarial game, we update the representation learning parameters, and in the next step, the treatment network q θ W (• | Φ t+1 ) with adverse losses such that the representation Φ t+1 becomes invariant with respect to the assignment of W t+1 . Different from SOTA models (as highlighted in related work) and in line with our information guidelines principles, learning a balanced representation Φ t+1 amounts to ensuring Φ t+1 ⊥ ⊥ W t+1 , which is equivalent to I(Φ t+1 , W t+1 ) = 0. Hence, we minimize the MI as a way to confuse the treatment classifier. Specifically, we minimize an upper bound on I(Φ t+1 , W t+1 ), namely the CLUB of MI [13] .EQUATIONWe use the objective in Eq. ( 9) to update the representation learner Φ(.) [9, 32] . This update aims to minimize the discrepancy between the conditional likelihood of treatments for units sampled from P (Ht,Wt+1) and the conditional likelihood of treatments under the assumption of independent sampling from the product of marginals P Ht+1 ⊗ P Wt+1 . In practice, we generate samples from the product of marginals by shuffling the treatment W t+1 across the batch dimension similar to [9, 32] . When minimizing L W , q θ W (ω t+1 | Φ t+1 ) gets closer to the true conditional distribution p(ω t+1 | Φ t+1 ), and, in this case, the objective in Eq. ( 9) provides an upper bound of the MI between representation and treatment. We formalize the intuition by adapting the result of [13] :Theorem 5.3. [13] Let q θ W (Φ t+1 , ω t+1 ) := q θ W (ω t+1 |Φ t+1 )p(Φ t+1 ) be the joint distribution induced by q θ W (ω t+1 |Φ t+1 ) over the representation space of Φ t+1 . If:D KL (p(Φ t+1 , ω t+1 )||q θW (Φ t+1 , ω t+1 )) ≤ D KL (p(Φ t+1 )p(ω t+1 )||q θW (Φ t+1 , ω t+1 )), then I(Φ t+1 , W t+1 ) ≤ I CLU B (Φ t+1 , W t+1 ; q).Based on Theorem 5. ) and achieving balance. We theoretically formulate such behavior by proving in the following theorem that, at the Nash equilibrium of this adversarial game, the representation is exactly balanced across the different treatment regimes provided by W t+1 .Theorem 5.4. Let t ∈ {1, 2, . . . , t max }, Φ = Φ θ R and q = q θ W are, respectively, any representation and treatment network. Let P Φ(Ht) be the probability distribution over the representation space and P Φ(Ht)|Wt+1 its conditional counterpart. Then, there exist Φ * and q * such that:EQUATIONSuch an equilibrium holds if and only if P Φ(Ht)|Wt+1=0 = P Φ(Ht)|Wt+1=1 = . . . = P Φ(Ht)|Wt+1=k-1 , almost surely.Since we target multi-timestep forecasting, covariate balancing in the representation space extends beyond t + 1. For simplicity, we presented it for t + 1, but in practice, the adversarial game applies the balancing across all forecasting horizons (Algorithm 2). The theorem also holds for other horizons by replacing Φ(H t ) with Φ(H t+j-1 ) and W t+1 with W t+j , for 2 ≤ j ≤ τ .The Causal CPC model is trained in two stages: (1) Encoder pretraining:We first learn an efficient representation of the process history by minimizing loss:L enc = L CP C (θ 1 , θ 2 , {Γ j } τ j=1 ) + L (Inf oM ax) (θ 1 , θ 2 , γ).(2) Decoder training: After pretraining, we fine-tune the encoder by optimizing the factual outcome and treatment networks in the adversarial game from Theorem 5.4. Formally:min θR,θY L dec (θ R , θ Y , θ W ) = L Y (θ R , θ Y ) + I CLUB (Φ θR (H t ), W t+1 ; q θW ), min θW L W (θ W , θ R ) = -E Φθ R (Ht) [log q θW (W t+1 | Φ θR (H t ))] .Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: Code is provided in the supplementary material at submission.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We compare Causal CPC with SOTA baselines: MSMs [64] , RMSN [41] , CRN [7] , G-Net [39] , and CT [48] . All models are fine-tuned via a grid search over hyperparameters, including architecture and optimizers. Model selection is based on mean squared error (MSE) on factual outcomes from a validation set, and the same criterion is used for early stopping. Further details on hyperparameters and training procedures are provided in Appendices J and D.Tumor Growth We use the PharmacoKinetic-PharmacoDynamic (PK-PD) model [24] to simulate responses of non-small cell lung cancer patients, following previous works [41, 7, 48] . We evaluate our approach on simulated counterfactual trajectories, varying the confounding level via the parameter γ (Appendix E.1). Unlike [48] , who used larger datasets (10,000 for training, 1,000 for testing), we use a smaller, more challenging dataset (1,000 for training, 500 for testing) to reflect real-world data limitations. For long-horizon forecasting, we set the prediction horizon to 10 and evaluate two training sequence lengths: 60 and 40, with covariates of dimension 4. MSM is excluded due to high prediction errors.We tested all models on the cancer simulation data across three confounding levels, γ = 1, 2, 3. Figures 3a, 3b and 3c show the evolution of Normalized Root Mean Squared Error (NRMSE) over counterfactual tumor volume as the prediction horizon increases. Causal CPC consistently outperforms all baselines at larger horizons, demonstrating its effectiveness in long-term predictions. This confirms the quality of C t in predicting future components across multiple time steps, capturing the global structure of the process as discussed in Eq. ( 5).Extended results are provided in Appendix E.2.1. In the more challenging case where the maximum sequence length is 40 (Figures 3d, 3e and 3f), the error evolution remains similar to Figure 3a , 3b and 3c. Our model maintains its advantage, outperforming most baselines in long-term forecasting. However, Causal CPC does not outperform other models in short-term forecasting, a consistent limitation across experiments. Still, the model's superior long-term performance highlights its potential in applications requiring long-term accuracy. At higher levels of confounding, the model does not always outperform SOTA models at certain time steps. This may be due to the low dimensionality of the time-varying components and static covariates,U t = [V, X t , W t-1 , Y t-1 ],which have only four dimensions. Our model leverages contrastive learning-based regularization to excel on datasets with higher confounding dimensions, as demonstrated on MIMIC-III where U t has 72 dimensions. In this setting, our model consistently outperforms baselines at longer prediction horizons. The occasional underperformance of Causal CPC at the final horizon is due to τ = 10 being the last contrasted horizon, not an issue specific to τ = 10. To support this, we reran all models with a sequence length of 60, τ = 15, and γ = 2.As shown in Figure 4 , Causal CPC still outperforms SOTA for horizons beyond τ = 10 due to the encoder's retraining, where the InfoNCE loss is computed over all 15 time steps. The last prediction error remains close to SOTA, suggesting that training over larger horizons than initially intended may be beneficial.Semi-synthetic MIMIC-III We used a semi-synthetic dataset constructed by [48] based on the MIMIC-III dataset [35] , incorporating both endogenous temporal dependencies and exogenous dependencies from observational patient trajectories, as detailed in Appendix F.1. The , where K is the number of possible treatments. This is particularly relevant when generating multiple treatment plans, such as minimizing tumor volume. Table 2 shows the models' complexity (number of parameters) and running time, split between model fitting and prediction. Causal CPC is highly efficient during prediction due to its simple 1-layer GRU, similar to CRN (1-layer LSTM), while providing better ITEs estimation. In contrast, CT is less efficient due to its transformer architecture and teacher forcing, which requires recursive data loading during inference. G-Net also has longer prediction times due to Monte Carlo sampling. Overall, Causal CPC strikes a strong balance between accuracy and efficiency, making it well-suited under constrained resources. Real MIMIC-III Data We evaluated our model on real MIMIC-III data, where counterfactual trajectories cannot be assessed due to the absence of observed counterfactual responses. However, performance can still be measured by forecasting factual (observed) responses over time. Our model estimates responses for each individual based on their observed treatment trajectory. As shown in Figure 6 , Causal CPC consistently outperforms all baselines, especially at larger horizons, demonstrating its robustness and effectiveness in real-world settings.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Models' hyperparameter search range is provided in Appendix J, the selection method is provided at the beginning of Section 6, and remaining details about training and testing are provided in the experimental protocol (Appendix D).Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
discussion: 
 Why does Causal CPC outperform SOTA at large horizons? Our context C t is designed to capture shared information across future representations, particularly covariates, by minimizing the InfoNCE loss over multiple time steps (Eq. 1). As shown in Eq. ( 5), minimizing L CP C maximizes shared information between the context and future components, helping capture the global structure of the process. This is especially beneficial for counterfactual regression over long horizons, explaining the model's superior performance. However, it may not always outperform SOTA in shorter-term predictions due to its focus on long-term dependencies.Short-term Counterfactual Regression While our model is designed for long-term predictions, it may not consistently outperform SOTA for short-horizon tasks. However, the use of contrastive loss, particularly InfoNCE (Eq. 4), suggests potential adaptability to balance both short-and long-term predictions without retraining. A trade-off could be achieved by adjusting the contrastive term weights across time steps in Eq. ( 4), which we leave for future work. Ablation study We examined the model's performance in various configurations-full model, without CPC, and without InfoMax. Table 3 and Figure 7 show that removing either term reduces counterfactual accuracy across all horizons, underscoring their significance. Additionally, replacing our ICLUB objective with CDC loss [48] or removing balancing increases errors. We also tested different MI lower bounds like NWJ [50] and MINE [4] Estimating MI bounds for high-dimensional variables is challenging and expensive [59, 57] , often limited to low-dimensional inputs or Gaussian assumptions.In MI-constrained models, batch size is crucial. As shown in Eq. ( 3), increasing the batch size B tightens the lower bound via log(|B|), and to balance memory and performance, we chose batch sizes of 256 for the encoder and 128 for the decoder. While these bounds may not be perfectly tight, mutual information and self-supervision biases significantly enhance performance (Appendix C.2), as confirmed by ablation studies. Other MI estimators like NWJ and MINE [50, 4] did not improve performance; our initial setup consistently performed better (Appendix F.2.1).Extending Causal CPC to Continuous Treatment Our approach could be extended to continuous treatments by replacing the treatment classifier with a regressor. Since the method maximizes likelihood, the equilibrium in Theorem 5.4 remains valid. However, in practice, continuous treatments will be represented by a single dimension, unlike discrete treatments with K-dimensional one-hot encoding. This risks losing important treatment information in counterfactual predictions. A simpler adaptation to our model could involve discretizing continuous treatments.We introduced a novel approach to long-term counterfactual regression, combining RNNs with CPC to achieve SOTA results without relying on complex transformer models. Prioritizing computational efficiency, we incorporated contrastive loss-based regularization guided by mutual information (MI). Our method consistently outperforms existing models on both synthetic and realworld datasets, marking the first application of CPC in causal inference. Future work could focus on improving interpretability by integrating Shapley values into the causal framework [30] . Additionally, developing uncertainty-aware models tailored for longitudinal data is crucial for enhancing the reliability of predictions in our causal framework [19, 33, 84] .Our paper seeks to advance the field of Trustworthy Machine Learning by focusing on the accurate estimation of counterfactual trajectories. This capability holds significant potential to enhance decision-making processes across various domains, particularly in healthcare, where clinicians can leverage models designed to mitigate bias and promote fairness. Additionally, by focusing on efficiency, our contributions extend beyond traditional machine learning considerations to address environmental concerns associated with energy consumption. By advocating for the prudent use of computational resources, especially in training complex models deployed in real-world scenarios, we aim to promote sustainability in developing and applying machine learning solutions.In this section, we detail the assumptions used for the identifiability of the counterfactual responses E(Y t+τ (ω t+1:t+τ ) | H t+1 ). As briefly stated in Section 3, we follows similar assumptions to [62, 64, 7, 48] , namely Assumption B.1 (Consistency). For every time step t and given any manner by which a unit i receives the sequence of treatment ω i,≤t , we always observe the potential outcome Y it (ω i,≤t ). Formally:W i,≤t = w i,≤t =⇒ Y it = Y it (w i,≤t ).Assumption B.2 (Sequential Ignorability). Given any time step t, we have the conditional independence:Y it (ω it ) ⊥ ⊥ W it |H it = h it ∀(ω it , h it )Assumption B.3 (Overlap/positivity). Given any time step t, and for any possible historical context h t , the probability of observing any of the possible treatment regimes is strictly positive but not deterministic:p(h t ) ̸ = 0 =⇒ 0 < p(W t = ω t |h t ) < 1The three assumptions are sufficient for the identification of the counterfactual responses from observational data, which we formulate in the following proposition. E(Y t+τ (ω t+1:t+τ ) | H t+1 ) = E (Y t+τ | H t+1 , W t+1:t+τ = ω t+1:t+τ ) Proof. See [62]We repeat the causal graph introduced in Figure 1 to explain the data generation process. here, all of the past observed data encompassed in H t+1 confounds future treatments and responses, W t+1 , W t+2 , . . . , W tmax and Y t+1 , Y t+2 , . . . , Y tmax , which create long-term dependencies. The fact that post-covariates are affected by past treatments creates time-dependent confounding. The static covariates are assumed to be affecting all of the time-varying variables. Since we suppose sequential ignorability, there are no possible exogenous noises affecting both treatments and responses. However, such noise may possibly affect responses, time-varying covariates, and response variables.In the figure, for simplicity, we represent past treatments as W ≤t such that each element in that sub-sequence confounds the next treatment and response W t+1 and Y t+1 . Idem for Y ≤t and X ≤t .The static covariates V are assumed to be affecting all the time-varying variables. We omit the representation of exogenous noise for simplicity. Interactions between W ≤t , X ≤t , and Y ≤t were also omitted for simplicity. In this section, we give a brief overview of models included in our experiments: MSMs [62] , RMSN [41] , CRN [7] , G-Net [39] , and CT [48] . To delineate the differences between these models and Causal CPC, we detail in Table 5 the main design differences between all these models. Our work relies on the assumption of sequential ignorability, yet several alternative models operate under different assumptions, often addressing the presence of unobserved confounders. Some of these models are rooted in deconfounding theory [45, 60, 80] , which has been extended to time-varying settings. Deconfounding involves imposing a factor model on treatment assignment, where each treatment becomes conditionally independent given latent variables that act as proxies for unobserved confounders. Examples of this approach include [8, 28, 10] . Other models assume the presence of proxy variables, inferring a representation of unobserved confounders through probabilistic models based on these proxies [79, 12, 38] .Y ≤t W ≤t X ≤t Y t+1 W t+1 X t+1 V C ExtendedIn contrast to our setting, which is governed by the three causal assumptions in Appendix B.1, many models assume a data-generating process similar to [73, 72, 58] . These methods, often non-or semiparametric, tend to either ignore static covariates or handle them linearly, leading to computational inefficiencies and scalability issues. Nevertheless, some non-or semi-parametric approaches-such as [68, 69, 19, 31] -align with our causal assumptions but extend them to continuous time, treating sequential ignorability in a continuous setting. Additionally, models like [34] incorporate continuous-time and assume interactions between units, where an individual's outcome depends on both their treatment and the treatments of others. [6] , focusing on binary treatment sequences, requires a stronger version of sequential ignorability-conditional on current covariates-whereas our model assumes a weaker version, conditioning on the entire history of covariates to account for long-lasting confounding effects. Furthermore, [11] focuses solely on binary treatments and targets the estimation of the average treatment effect on the treated (ATT). The authors assume a specific treatment regime, where individuals enter a post-treatment state after a defined point in time. This assumption is restrictive compared to our framework, which allows for complex, individualized treatment assignment mechanisms and non-binary treatments, where treatment values can fluctuate over time. As a result, [11] is incompatible with our causal assumptions. Other methods, like [82] , address high-dimensional counterfactual generation based on time-varying treatment plans under the same sequential ignorability assumption. However, they are not designed for causal forecasting over multiple time steps, as required in our setting. Similarly, [22] focuses on estimating the average causal effect and is not suited for predicting individual treatment effects or conditional counterfactual responses, as it targets marginal counterfactual expectations via gcomputation.C.2 Mutual Information and Self-Supervision Self-Supervised Learning and Mutual Information In self-supervised learning, Deep InfoMax [32] [85] encourages representations capable of discriminating between individual examples in the dataset, our preference for CPC arises from its adaptability in processing sequential features in an ordered and autoregressive manner, which aligns with the requirements of our specific context, especially when dealing with counterfactual regression over time.Mutual Information and Inductive Bias. Mutual information (MI) estimation success relies not only on MI's properties but also on the inductive biases from feature representation choices and MI estimator parameterization [77] . Experimental evidence shows that, although MI remains invariant under homeomorphisms, maximization with an invertible encoder during random initialization enhances downstream performance. While higher-capacity critics yield tighter MI bounds, findings consistent with [59] suggest that simpler critics provide better representations, even with looser MI bounds. Accordingly, we selected a simple bilinear critic function for contrastive losses. In vision tasks, augmentations and contrastive loss properties are crucial for representation efficiency [1, 76, 27] , and [66] highlights that inductive bias, via function class representation and optimizers, significantly affects downstream performance, offering theoretical, non-vacuous guarantees on representation quality.The estimation of MI faces inherent challenges, particularly within variational lower bounds. These bounds often degrade as MI increases, creating a delicate trade-off between high bias and high variance. To address this, methods that utilize upper bounds on MI have been developed, attempting to mitigate challenges associated with variational bounds. One strategy for MI maximization involves computing gradients of a lower MI bound concerning the parameters of a stochastic encoder. This computational approach potentially eliminates the need for direct MI estimation, providing a more tractable solution. However, estimating MI from samples remains challenging, and traditional approaches encounter scalability issues in modern machine-learning problems. It's crucial to note that higher estimated MI between observations and learned representations does not consistently translate to improved predictive performance in downstream supervised learning tasks. CPC is an example, exhibiting less variance but more bias, with estimates capped at log |B|.Strategies to reduce bias, such as increasing the batch size, introduce higher computational complexity, requiring additional evaluations for estimating each batch with the encoding function.In our empirical approach, we adopt a specific sampling strategy for sequences, considering a onetime step per batch. This facilitates computing the InfoNCE between local summary features at time t and the future prediction of local features, leading to a reduction in algorithmic complexity for contrastive loss computation. Empirical observations demonstrate non-decreased representation quality and improved prediction of factual and counterfactual outcomes.Other MI lower bounds. The Mutual Information Neural Estimator (MINE) [4] leverages the relationship between MI and the Kullback-Leibler (KL) divergence. MI can be expressed as the KL divergence between the joint distribution and the product of marginals: I(X; Z) := D KL (P (X,Z) ∥ P X ⊗ P Z )MINE employs the Donsker-Varadhan representation [20] of the KL divergence:EQUATIONHere, the supremum is taken over all functions T where the expectations exist. For a specific class of functions F, potentially represented by a class of neural networks, we obtain the lower bound:EQUATIONIn practice, we maximize ÎMINEγ (P ∥ Q) = E P [T γ ] -log E Q [e Tγ ] ,where T γ is a discriminator parameterized by γ, representing neural network parameters. The MINE estimator is a strongly consistent estimator of the true MI (Theorem 2, [4] ). Alternatively, the f-divergence representation of D KL [51] allows us to derive another MI lower bound, known as the Nguyen, Wainwright, and Jordan (NWJ) estimator [50] :EQUATIONThis results in the estimator:ÎNWJ γ (P, Q) = E P [T γ ] -log E Q [e Tγ -1 ] .Unlike the InfoNCE estimator, which exhibits high bias and low variance, the NWJ estimator has a low bias but high variance [57] .All models were implemented using PyTorch [53] and PyTorch Lightning [21] . In contrast to the approach in [48] , we employed early stopping for all models. The stopping criterion was defined as the Mean Squared Error over factual outcomes for a dedicated validation dataset. Specifically, for the Causal CPC encoder, the stopping criterion was determined by the validation loss of the encoder. While all models in the benchmark were trained using the Adam optimizer [36] , we opted for training Causal CPC (encoder plus decoder without the treatment subnetwork) with AdamW [46] due to its observed stability during training. Similar to the common practice in training GAN discriminators, the treatment subnetwork was optimized using SGD with momentum [74] .The CT employed the Exponential Moving Average (EMA) [83] of parameters to enhance training stability. However, this technique was not applied to Causal CPC, as experimental evidence suggested only marginal improvements. Weight decay was set to zero for all models.For each experiment, the models were trained over five different seeds, and the reported performance metrics include the mean and standard deviation of the results. The counterfactual trajectories are generated following two strategies:• Single sliding treatment [7, 48] : Trajectories are generated with a single treatment per trajectory while the treatment slides over the forecasting range to generate multiple trajectories.Similar to [7] , we apply such a generation scheme to cancer simulation data. • Random trajectories: Trajectories are generated such that at each time step, treatment is generated randomly. We apply random trajectories to semi-synthetic MIMIC data.For the falsifiability test on MIMIC III datset, we mask two confounders from the inputs of the benchmark models, namely sodium and glucose measurements.We present a tumor growth simulation model, focusing on the PharmacoKinetic-PharmacoDynamic (PK-PD) model as discussed in [24] , a recent approach to predicting treatment responses in non-small cell lung cancer patients. This simulation models the evolution of tumor volume, denoted by V (t), in discrete time, where t represents the number of days since diagnosis:V (t) =     1 + Λ log K V (t -1) Tumor Growth -κ c C(t) Chemotherapy -κ rd Rd(t) + υRd(t) 2 Radiation + e t Noise     V (t -1) (14)Here, the model parameters Λ, K, κ c , κ rd , υ are sampled for each patient based on prior distributions from [24] . Additionally, Rd(t) represents the radiation dose applied at time t, and C(t) denotes the drug concentration. We introduce confounding into the assignment of radiotherapy/chemotherapy treatment by making it dependent on past tumor volume evolution. Treatment is simulated using a Bernoulli distribution with probability σ(π t ), where:EQUATIONIn this context: D(t) represents the average tumor diameter over the last 15 days, D max = 13 cm is the maximum tumor diameter, δ is set to δ = D max /2. The parameter γ controls the level of time-dependent confounding, with a higher γ value assigning more weight to the history of tumor diameter in treatment assignment.We report in this section detailed counterfactual errors for Causal CPC and baselines over the cancer simulation dataset, which are responsible for Figure 3 .Table 6 : Results on the synthetic data set with sequence length 60: mean±standard deviation of NRMSEs. The best value for each metric is given in bold: smaller is better. γ = 1 Model τ = 1 τ = 2 τ = 3 τ = 4 τ = 5 τ = 6 τ = 7 τ = 8 τ = 9 τ =We detail here the results of the ablation study conducted on the cancer simulation dataset (Table 3 ). The (full) Causal CPC model, as presented in the core paper, gives, in most cases, better results than any ablation configuration.Table 7 : Results on the synthetic data set with sequence length 40: mean±standard deviation of NRMSEs. The best value for each metric is given in bold: smaller is better.γ = 1 Model τ = 1 τ = 2 τ = 3 τ = 4 τ = 5 τ = 6 τ = 7 τ = 8 τ = 9 τ = 10Causal CPC (ours) 1.21±0.07 1.13±0.12 1.22±0.12 1.27±0.17 1.35±0. 19 F Experiments on semi-synthetic data: DetailsIn this section, we provide a concise overview of the simulation model built upon the MIMIC III dataset, as introduced by [48] . Initially, a cohort of 1,000 patients is extracted from the MIMIC III data, and the simulation proposed by [48] extends the model of [68] . Let d y be the dimension of the outcome variable. In the case of multiple outcomes, untreated outcomes, denoted as Z j,(i) tfor j = 1, . . . , d y , are generated for each patient i within the cohort. The generation process is defined as follows:EQUATIONwhere: the B-spline B-spline(t) is an endogenous component, g j,(i) (•) is sampled independently for each patient from a Gaussian process with a Matérn kernel and f j Z (•) is sampled from a Random Fourier Features (RFF) approximation of a Gaussian process. To introduce confounding in the assignment mechanism, current time-varying covariates are incorporated via a random function f l Y (X t ) and the average of the subset of the previous T l treated outcomes, ĀT l Y t-1 . For d a binary treatments A l t , where l = 1, . . . , d a , the assignment mechanism is modeled as:p A l t = σ γ l A ĀT l Y t-1 + γ l X f l Y (X t ) + b l , A l t ∼ Bernoulli p A l t .Subsequently, treatments are applied to the untreated outcomes using the following expression:EQUATIONThe final outcome combines the treatment effect and the untreated simulated outcome: Proof. In the following, we draw inspiration from the proof of [52] . The InfoNCE loss corresponds to the categorical cross-entropy of classifying the positive sample U t+j correctly, given the context C t , with a probability:EQUATIONexp(T j (U t+j , C t )) |B| l=1 exp(T j (U l,t+j , C t )).The positive sample U t+j is one element in the batch B, where the remaining elements serve as negative samples. Let pos ∈ {1, . . . , |B|} be the indicator of the positive sample U t+j . The optimal probability is given by: .p(Index = pos | B, C t ) = p(u pos,t+j | C t )For the score exp(T j (U t+j , C t )) to be optimal, it should be proportional to EQUATIONThe approximation in the third equation, Eq. ( 19), becomes more precise as the batch size increases.We now prove Proposition 5.1, which states that:I(C h t , C f t ) ≤ I(H t , (C h t , C f t )).Proof. This follows from two applications of the data processing inequality [17] , which states that for random variables C h t Φ θ 1 ,θ 2 •trunc f ← ---------H t Φ θ 1 ,θ 2 •trunc h ---------→ C f t , which is Markov equivalent to: C h t Φ θ 1 ,θ 2 •trunc f ---------→ H t Φ θ 1 2 •trunc h ---------→ C f t .By the data processing inequality, this results in I(C To begin, we split the process history into two non-overlapping views (Figure 2 ): H h t := U 1:t0 and H f t := U t0+1:t , representing a historical subsequence and a future subsequence within the process history H t , respectively. We then computed representations of these two views denoted C h t and C f t , respectively. This naturally gives rise to the Markov chain, as in showed in the proof of proposition 5.1:C h t ← -H t - → C f twhich is Markov equivalent to:C h t - → H t - → C f tFollowing this Markov chain, we can show that [71] :I(C f t , C h t ) = I(H t , C h t ) -E ht∼P H t E c f t ∼P C f t |h t D KL [P C h t |ht ||P C h t |c f t ]On the other hand, by applying the chain rule of the mutual information to I(H t ; (C h t , C f t )) we get:I(H t ; (C f t , C h t )) = I(H t , C h t ) + I(H t ; C f t | C h t )Combining these equations, the tightness of our bounds can be written as:I(H t ; (C f t , C h t )) -I(C f t , C h t ) = I(H t ; C f t | C h t ) + E ht∼P H t E c f t ∼P C f t |h t D KL [P C h t |ht ||P C h t |c f t ] G.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [No] Justification: It is computationally demanding to compute errors bars for all neural network models in our benchmark. However, we reported the mean and standard deviation of metrics for each experiment, computed from multiple runs.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]• The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: The paper presents a novel approach to counterfactual regression over time, i.e., long-term predictions. It combines RNNs with CPC and InfoMax. This approach aims to capture long-term dependencies in data without relying on complex transformer models. There is a general   consensus among the reviewers that the paper contains interesting ideas and will make non-negligible contributions to the field. The reviewers raised several important concerns, most of which were addressed in the authors' rebuttal. It should be possible for the authors to address the remaining concerns in the camera-ready version of the manuscript.

Hence, I recommend that the paper be accepted as a poster.

Reviews: review: paper_summary: paper_summary: The paper presents a novel approach to counterfactual regression over time, particularly focusing on long-term predictions. It introduces a method that leverages RNNs combined with CPC and InfoMax. This approach aims to capture long-term dependencies in data and improve computational efficiency without relying on complex transformer models.
strength_weakness: strength_weakness: 1. Complexity: While the method avoids transformers, the combination of CPC and InfoMax introduces its own complexity, which may pose challenges for implementation and replication. Also reported in Table 2 and Table 11, where the proposed method requires more training time with fewer parameters compared to the Causal Transformer and G-Net. 
2. Scope of evaluation: Although the evaluation is comprehensive, it could be expanded to include more diverse datasets to further validate the robustness and generalizability of the method. For example, 1) what is the total trajectory length in the synthetic experiment? I wonder about the performance and an additional experiment showing the performance under different trajectory lengths would be appreciated. 2) what is the dimensionality of the variables? I wonder about the performance when dealing with the high-dimensional case.
questions: questions: 1. In the problem formulation, the treatment of W is assumed to be discrete, is it possible to extend it to be continuous and would such an extension violate the theoretical result?
2. What is the performance of the approach when dealing with time series under different trajectory lengths and high-dimensional cases? 
3. How should the counterfactual treatment be chosen in the experiment? Randomly chosen from the support or do the authors have some rules?
4. Could the authors explain why the proposed method 'excels at large horizon predictions, it does not outperform SOTA models on short-term predictions'?
limitations: limitations: The authors claimed in the Checklist.
paper_summary: paper_summary: This paper proposes a new method for counterfactual outcome prediction over time, with the goal of avoiding highly complicated models and expensive computation while achieving state-of-the-art prediction accuracy for time series with long-range dependencies. To this end, the presented method employs RRNs for long-range forecasting, rather than more complicated network architecture like Transformers, together with contrastive predictive coding (CPC) to capture long-range dependencies in the presence of time varying confounders and information maximization (InfoMax) to help with learning invertible representation for addressing the identification issue. Experiments with synthetic and semi-synthetic data are conducted to demonstrate the efficacy and superiority of the proposed method.
strength_weakness: strength_weakness: 1. The efficacy of the proposed method for long-range predictions is only demonstrated by experiments to certain extend, and from the results, it is not clear how the method would perform when $\tao$ is greater than 10 steps. In particular, from the middle and right diagrams in Fig. 3, the advantage of the  proposed method over the baselines start to drop and some of the baselines start to outperform the proposed method. The paper would be much stronger if some theoretical evidence of the performance of the proposed method could be provided.
2. There are some doubts with the experiment results or performance of the proposed method (please see the Questions section for more details).
questions: questions: 1. Why a mapping of $\textbf{C}_t$, instead of $\textbf{C}_t$ itself is used as input of the decoder part of the proposed architecture?
2. Could you provide any explanation on why the performance advantage of the proposed method starts to drop when $\tao$ becomes bigger (as shown in the middle and right diagrams in Figure 3, and Table 10)? Moreover, why the baselines in general over outperform the proposed method with small horizon predictions?
3. Table 2: the results shown are for the tumor growth data only and $\gamma=1$. How the number of model parameters, training and testing time when different $\gamma$ values are used?
4. Why MSM is not used/reported in all the experiments? Having poor performance may not be a convincing reason for excluding it.
limitations: limitations: The paper has indicated some limitations of the proposed method, e.g. in Fig. 3 and Section 6.4, but more explicit discussions on the limitations in terms of decreased performance/superiority when the prediction horizon is small or large and when there are unobserved confounders should be given.
paper_summary: paper_summary: This paper presents causal CPC, a framework for predicting counterfactual responses under time-varying treatments. The proposed method consists of two components: an encoder that leverages contrastive predictive coding and infomax principle to learn a representation of the history H, and a decoder that leverages a loss function based on adversarial game to predict outcome and treatment with balanced representations to account for selection bias. The concept of utilizing MI-based principles for casual representation learning in the longitudinal setting is interesting, whose efficiency has been demonstrated by the SOTA. However, I believe the paper could benefit from more streamlined writing with a clearer justification of how the use of MI-based principles enables efficient learning of long-term dependency.
strength_weakness: strength_weakness: 1. Even thought the paper demonstrates the benefits of using CPC and infomax for learning long-term dependency with RNN, it still remains unclear to me, both in intuition and in methodological details, how this benefit was attained. The paper can benefit from some more insights, particularly in section 5.1, that links the proposed information theoretic methods to their actual benefits. The theoretical analysis seems a bit detached from the context and not so helpful in providing insights for the effectiveness of the proposed method.

2. Following on the second point, in the ablation study, it seems that the removal of the InfoNCE loss is non-significant. This seems somewhat curious to me given the amount of emphasis the authors put on the importance of CPC.

3. Some recent work in counterfactual outcome prediction with time-varying treatments should be discussed in the related work:
 - Berrevoets et al, Disentangled counterfactual recurrent networks for treatment effect inference over time 

 - Chen et al, A Multi-Task Gaussian Process Model for Inferring Time-Varying Treatment Effects in Panel Data 

 - Wu et al, Counterfactual Generative Models for Time-Varying Treatment 

 - Frauen et al, Estimating average causal effects from patient trajectories
questions: questions: 1. How to decide the relative weighting for the encoder and decoder during fine-tuning?
2. Related to Weakness #2, how to better illustrate the role of L_infoNCE given the relatively mild drop in performance after removing it?
limitations: limitations: Limitations were mentioned in the appendix, but not much in the main text. Adding more discussion of model limitations in the Conclusion section might help.
paper_summary: paper_summary: The paper introduces a novel algorithm for long-term counterfactual forecasting over time through representation learning of historical information and balanced representation learning that predicts the outcome given balanced treatments. The representation learning process combines Contrastive Predictive Coding and Information Maximization, while the training of balanced representation learning involves an adversarial game in the actual outcome and treatment networks.
strength_weakness: strength_weakness: 1. Is it possible to conduct experiments on real datasets?
2. It appears there is only one synthetic dataset and one semi-synthetic dataset. In the synthetic experiment, is the standard deviation of the errors calculated from experiments with 5 different seeds instead of random settings of the data generating process? Would it be feasible to conduct experiments on multiple synthetic datasets with random settings to obtain error bars, rather than relying on just one synthetic dataset?
3. There don't seem to be any comparisons of time consumption with baselines.
4. The loss functions are constructed based on the lower or upper bound of the mutual information, but there is no discussion or measurement regarding the tightness of the bound.
questions: questions: 1. In equation 1, $U_{t+j}$ in the numerator has one subscript, while $U_{l,t+j}$ in the denominator has two subscripts. Could you please make them consistent? Is $U_{t+j}$ information for one specific individual or an aggregate for all individuals? Is InfoNCE calculated for each individual or as the expectation over all individuals? Given that in equation 1, the subscript is $\mathcal{B}$, is it the individual index
2. Could you elaborate on how Theorem 5.2 influences the construction of the loss function in Equation 7? Specifically, does Theorem 5.2 impact the formulation or design of the loss function described in Equation 7? The authors discuss what we can infer if equality holds from lines 189 to 192, but will equality actually hold? The description in this paragraph suggests that equality can be achieved, yet Proposition 5.1 and Theorem 5.2 indicate positive gaps between the two mutual information measures without specifying the magnitude of these gaps, and the construction of the contrastive loss does not appear to be influenced by Theorem 5.2, as long as Proposition 5.1 holds.
3. Could you explain what is $\gamma$ in equation 7? Is $\gamma$ a similar definition of horizons?
4. In line 202, is there any proof provided regarding "The theoretical existence of such a function"
5. Since the goal is to estimate the counterfactual responses given in the equation at the bottom of page 3, could you explain why we still need to predict the treatments as the conditional expectation is given $ W_{t+1:t+\tau} $? Is this training needed because we want balanced learning? If so, could you explain why selection bias should be avoided given the goal of learning the expectation conditioning on potential treatments?
6. In Theorem 5.3, could you explain why or when the “if” condition holds? 
7. Does the balanced learning only work for $W_{t+1}$? Why is it not necessary to extend it to $W_{t+\tau}$ for $\forall \tau\geq 1$?
8. In section 5.3, is $\theta_{R}=[\theta_1, \theta_2]$?
limitations: limitations: The authors have discussed the broader societal impacts.
paper_summary: paper_summary: This paper leverages Contrastive Predictive Coding (CPC) with RNN for counterfactual regression over time to provide a compelling alternative to transformer-based approaches (which are challenging to interpret).  By leveraging CPC to capture long-term dependencies and InfoMax for "reconstructable" representation, the method achieves state-of-the-art results in counterfactual estimation.
strength_weakness: strength_weakness: - Experiments: performance in random trajectories for the pharmacokinetic-pharmacodynamic model of tumor growth is missing. 

- While good results are achieved with fewer patient data, it is important to check if the results are consistent as the patient number increases. 

- Line 96: "*...the role of invertible representation in improving counterfactual regression. Here, we introduce an InfoMax regularization term to make our encoder easier to invert.*" The authors provide some intuitions (lines 204-213) in favor of invertibility (beyond reconstructable) but they are not sufficient. A formal proof is needed. 

- For a “reconstructable” representation of the process history $H_{t}$,  While a lower bound of InfoMax objective ( regularization) is optimized during pre-training of the encoder, there is no guarantee that finetuning of pre-training will retain the desired “reconstructable” property of the representation.
Hence, the role of invertible (or reconstructable) representation in improving counterfactual regression is not well understood even though ablation studies show improvement when $\mathcal{L}^{Infomax}$ is introduced.
questions: questions: -  Line 178-179: *$H_{t}$ is a sequence of high-dimensional covariates and the computation of such loss is computationally demanding.* 
- Line 165-167: *...where at each horizon the discriminator parameter matrix $\Gamma_{j}$ is of a small number of dimensions since being a map between two lower-dimensional representations*.

Lower bounds are motivated by the problem of high-dimensionality of $H_{t}$.
What is the dimension of $H_{t}, Z_{t}$ and $C_{t}$ in the performed experiments?
limitations: limitations: See Weaknesses.
paper_summary: paper_summary: the paper proposes a method based on CPC to estimate the causal effect of treatments over a period of time
strength_weakness: strength_weakness: - Not clear from an intuition perspective if the contrastive approach incentivises the model to learn what the authors want 
- The complexity of the method is significant , raising questions about applicability and reproducibility 
- The point about invertible processes is not a great one as invertible processes severely limit the function space available without offering much in return. They also increase computational complexity significantly
 - not the first use of a contrastive objective in causal learning , id suggest that this contribution is toned down a bit
questions: questions: not many, mostly how susceptible is the method to violations of the assumptions the authors have set out ?
limitations: limitations: adequately discussed

Comments: title: Thanks for your detailed responses
Rebuttal acknowledgement
Correction of a typo in rebuttal
comment: Thank you very much for reviewing our rebuttal. We appreciate your consideration and positive assessment of our paper.
Thank you very much for your response; we’re happy our responses were helpful and your concerns addressed.
Thanks the authors for your detailed and helpful responses. I will maintain my positive rating.
Thank you very much for your thoughtful feedback and for raising the score. We appreciate your suggestion about the figures and will work specifically on improving them for clarity in the final version.
Thank you for reviewing our rebuttal and for maintaining your positive score. We appreciate your support and consideration.
I acknowledge that I have read the author's rebuttal. I maintain my score of acceptance.
I acknowledge that I have read the authors rebuttal, and I maintain my score of acceptance
Thank the authors for their clarification and now I have a better understanding of the details. Although I still have some concerns about the higher training time, it is acceptable to me now and I will raise my score. Also it could be better if the authors could polish the figures, as the current font size is somewhat small and difficult to read clearly.
Thank you for reviewing our responses; we're glad your concerns have been addressed. We appreciate your positive assessment and constructive feedback.
Thank you for reviewing our responses; we’re happy your concerns have been addressed. We appreciate your thoughtful consideration.
Thanks for the comments. My concerns are addressed and I will raise my score to 5.
Thank you for the detailed answers. My concerns have been addressed. I will maintain the score.
In response to Question 1, line 2:  "However, the results should be reported as $\Phi(\mathbf{H\_t}) \perp W\_{t+1}$ ...".