Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Authors: Xin Wang
Shengfei Lyu
Chi Luo
Xiren Zhou
Huanhuan Chen
Venue: ICML 2025 poster (2025)
Keywords: causal inference
counterfactual reasoning
counterfactual target achievement
Primary area: general_machine_learning->causality
TL;DR: We introduce a novel counterfactual target achievement problem in temporal systems and propose VCIP to find optimal interventions by modeling achievement probability.

Abstract: A key challenge in personalized healthcare is identifying optimal intervention sequences to guide temporal systems toward target outcomes, a novel problem we formalize as counterfactual target achievement. In addressing this problem, directly adopting counterfactual estimation methods face compounding errors due to the unobservability of counterfactuals. To overcome this, we propose Variational Counterfactual Intervention Planning (VCIP), which reformulates the problem by modeling the conditional likelihood of achieving target outcomes, implemented through variational inference. By leveraging the g-formula to bridge the gap between interventional and observational log-likelihoods, VCIP enables reliable training from observational data. Experiments on both synthetic and real-world datasets show that VCIP significantly outperforms existing methods in target achievement accuracy.

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Xin Wang' Shengfei Lyu? Chi Luo
Xiren Zhou'
Huanhuan Chen'
Abstract
paradigm has advanced understanding of intervention ef- fects, decision-makert often face critical inverse problem of identifying optimal intervention sequences that guide systems toward desired target outcomes.
A key challenge in personalized healthcare is iden- tifying optimal intervention sequences to guide temporal systems toward target outcomes, novel problem we formalize as counterfactual target achievement. In addressing this problem, directly adopting counterfactual estimation methods face compounding errors due to the unobservability of counterfactuals. To overcome this, we propose Variational Counterfactual Intervention Planning (VCIP), which reformulates the problem by mod- eling the conditional likelihood of achieving tar- get outcomes, implemented through variational inference. By leveraging the g-formula to bridge the gap between interventional and observational likeeiiiisoos VCIP enables reliable training from observational data. Experiments on both synthetic and real-world datasets show that VCIP significantly outperforms existing methods in tar- get achievement accuracy.
This paper introduces the counterfactual target achievement problem, a novel sequential decision makinn problem that aims to find intervention sequences guiding system outputs toward specific target outcomes based on historical observa- tions. Although similar to the sequential decision-making framework in dynamic treatment regimes (DTRs) (Laber et al., 2014), this problem differs fundamentally in its ob- jective. While DTRs aim to maximize population-level outcomes, decision-makers in many contexts may focus on achieving specific target outcomes for individual patients. For example, in diabetes management, physicians adjust in- sulin interventions maintain each patient's 's blood glucose within a personalized target range, rather than optimizing glycemic control for the population (Bergenstal, 2018; Pe- terson et al., 2007).
An intuitive solution is to leverage counterfactual estima- tion methods (c.g.. CRN (Bica al., 2020)) to estimate potential outcomes under different intervention sequences and select interventions by comparing these outcomes to the target outcome. However, as shown in the left panel of Figure 1, while sequence 2 is predicted to be closer to the target (blue dashed line), the true outcome (blue solid line) deviates significantly. This illustrates a key challenge: since counterfactuals are unobservable in practice, predic- tion errors cannot be accurately assessed, leading biases in estimated target distances and misleading intervention selections. Indeed, sequence I achieves better true perfor- mance (red solid line) despite its predicted trajectory (red dashed line) appearing less favorable.
1. Introduction
Causal inference in time series data plays a crucial role in healtheare and clinical decision support (Wu et al., 2022; Zheng al., 2020; prosperi et al., 2020). With the widespread adoption of electronic health records (EHRs), practitioners face the challenge of evaluating long-term treat- ment effects on patient outcomes, making accurate causal analysis essential for informed clinical decisions (Lim et al., 2018; Jensen et al., 2012; Chakraborty et al., 2022).
Existing research in temporal causal inference primarily focuses the "forward prediction" problem of predicting future system trajectories given historical observations and intervention sequences (Lim et al., 2018; Bica et al., 2020; Melnychuk et al., 2022; Wang et al., 2024). While this
To address these limitations, we reformulate the counterfac- tual target achievement problem by directly modeling the probability of reaching the target outcome, conditioned on the observed history, as shown in the right panel of Figure I. By introducing latent variables to capture temporal system state evolution, we propose the Variational Counterfactual Intervention Planning (VCIP) framework, which avoids the pitfalls of counterfactual estimation methods. VCIP first constructs an evidence lower bound (ELBO), a tool used for optimizing likelihoods under interventional distributions,
Universsity of Science and Technology of China NNayynng Technological University. Correspondence to: Huanhuan Chen <hchen @ www.e.uuuuusss.......... Shengfei Lyu <shengfei.lyue tbengtei.lyuentuucduug......
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Based on Prodisted Pokataal Outcomes
Based on Likethood Prodiction
6.0
6M
Historical data
Target
5.0
Prodisted
Truc
p(YLE J== Your! (M)
Iistorical data
Predicted 2
Discance Target:
3.
3.5
Sequence 1:
GH
Acmal Best
Trearment
Histertcal
No Treatment
Sequence
Proficted
25
25
20
Time
Time
Figure 1. Illustration of two prediction methods for counterfactual target achievement. Left: The counterfactual estimation method directly estimates potential outcomes under different intervention sequences, showing both predicted (dashed lines) and true (solid lines) trajectories, with distance indicators demonstrating each sequence's 's closeness to target outcome. Right: VCIP evaluates intervention sequences by modeling the probability of achieving the target outcome, conditioned on the observed history, and optimizing the likelihood of reaching the target.
and leverages the g-formula (Robins, 1986) to establish connections between interventional and observational likeli- hoods. VCIP circumvents error accumulation that typically arises from relying on explicit predictions and offers a robust training framework, ensuring reliability with observational data. Our main contributions are summarized as follows:
ever, these methods imposed strong prior assumptions on model structures, limiting their practical applicability.
Recent advances primarily leverage deep neural networks, including RMSN (Lim al., 2018), CRN (Bica al., 2020), G-Net (Li et al., 2021). CT (Melnychuk et al., 2022), and ACTIN (Wang et al., 2024). Architecturally, while RMSN, CRN, and G-Net build upon basic LSTM networks, CT adopts the more powerful Transformer architecture (Vaswani et al., 2017) to address limitations in capturing long-range dependencies among time-varying confounders in longitudinal data (Hochreiter al., 2001). ACTIN further proposes dual-module architecture to balance effectiveness and efficiency.
We formulate the novel temporal counterfactual target achievement problem of identifying intervention se- quences that guide dynamic systems toward specified target outcomes.
We propose Variational Counterfactual Intervention Planning (VCIP), framework that directly models target achievement likelihood through variational in- ference, with theoretical guarantees for training on observational data.
These methods employ time-varying adjustment techniques to mitigate confounding bias. RMSN builds dual propensity networks with IPTW scoring. G-Net integrates g- g-formula methods. CRN, CT, and ACTIN learn balanced represen- tations through gradient reversal, reverse KI. divergence minimization, and distribution-bace adversarial training re- spectively. While these approaches advance temporal coun- terfactual estimation, they have limitations for the target achievement task.
Through extensive experiments on both synthetic and real-world datasets, we demonstrate that VCIP con- sistently outperforms existing approaches in terms of target achievement accuracy.
2. Related Work
2.2. Dynamic Treatment Regimes
2.1. Temporal Counterfactual Estimation
Dynamic treatment regimes (DTRs) provide a formal frame- work for personalizing treatments in sequential medical decision-making while optimizing the population-level out- comes (Murphy, 2003). Most existing work in DTR opti- mization focuses offline learning from observational data (Zhang & Bareinboim, 2019). Recent works have explored reinforcement learning (RL) approaches for DTR optimiza- tion (Luckett et al., 2020; Raghu et al., 2017; Zhang & Bareinboim, 2019), addressing challenges such as policy evaluation biases from confounded data (Luo et al., 2024)
Estimating temporal counterfactual outcomes is closely re- lated to but fundamentally different from the proposed target achievement task. Early work in this field focused on statis- tical methods including g-formula, structural nested mod- els, and marginal structural models (Robins, 1986; Robins et al., 2000; Robins & Hernán, 2009). Due to the limi- tations of linear regression in handling complex temporal dependencies (Mortimer et al., 2005), research gradually shifted towards Bayesian nonparametric frameworks (Xu al., 2016; Soleimani al., 2017; Roy al., 2017). How-
2

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
and continuous state spaces (Tsirtsis & Gomez-Rodriguez, 2023). However, these methods primarily focus find- ing optimal treatment policies at the population level while maintaining individual-leve adaptivity. In contrast, our work aims to predict and optimize treatment sequences for specific individuals based on their historical trajectories, differing from DTRs in two key aspects: (1) we focus on individual-level optimization rather than population-level policy learning, and (2) directly maximize the probabil- ity of achieving patient-specific target outcomes instead of optimizing expected outcomes across the population.
7
Figure 2. Causal model for temporal dynamics. The latent factor Z. at time s represents the system state, which encapsulates the historical information H. Specifically, Z. is encoded from his- torical information H., and the outcome Yo+r is decoded from the final latent state Ze+r- The dashed arrows indicate the omitted intermediate steps in the temporal sequence.
3. Problem Formulation
system state. We assume that Z. encodes all relevant histori- cal information contained in Hg. thus serving as a sufficient representation of the system's 's temporal evolution.
Let D represent a longitudinal observational dataset com- prising records from N subjects, structured as D
[xf]) at Yi Ji=1 [((())] ;-1" Each subject's trajec- tory, denoted by index i, consists of sequential measure- ments over (1) time points, where x!") E X' represents time-varying covariates, A E [0, 1] or {0, 13 denotes continuous-valuee or binary treatments, respectively, and
(0) (0
4.1. Evidence Lower Bound (ELBO) Derivation
The optimal intervention sequence for Equation 1 is ob- tained by maximizing:
== lo aa.r)= Y He);
(2)
y E y indicates the measured outcomes. Time-invariant characteristics are captured in V()) E V. For simplicity in notation, we will drop the subject indicator (i) when the context is clear.
where denotes the parameters of the generative model. Following the causal model in Figure 2, we introduce vari- ational distribution go(2 (Zer++ H,, Yr-r,ārr)) to derive the ELBO for optimizing O:
The temporal progression of subject can be captured through their historical trajectory H, = (X,, Ar-1, Y, V), where the bar notation represents the sequence up to time t: X, = (X1,.. X.), Y, (Y,,..., YY), and A,- (A1,-.. 11
ELBO; EEgg Dog Po Y t+T, Zt,7+1 HA, d)(ā,,)))) (3) go (Z. Hr, Y YetT, Ar.r)
Under the potential outcomes framework (Rubin, 1978), let āc, (a.. a++r-1) denote sequence of treatments from time t to t + r -- 1. We denote Y[ā,,z] as the poten- tial outcome at time t + that would be observed under treatment sequence āz,= (Robins et al., 2000). For coun- terfactual target achievement, given a desired future state Y Trgget& at time + v, our research question focuses on deter- mining the optimal sequence of treatments that maximizes the likelihood of achieving this target:
where the do(.) operator (Pearl, 2009) represents atomic interventions that set variables to specific values, breaking the causal link from Z, to a, in the original causal model. As a practical approximation, due to the inaccessibility of the intervention distribution, we substitute with the obser- vational distribution. Therefore, we introduce a variational distribution 90 (Ze.r+1 H,, at,r, approximate the posterior distribution under intervention, where the ef- fects of interventions are partially captured in the observed outcomes Y E+T.
a² argmax, a PYY(a Y Larget H.).
(1)
However, in our causal model, the intervention assignment
For causal identification with observational data, we rely on standard assumptions including consistency, sequential ig- norability, and positivity (see Appendix A). For comparison of counterfactual larget achievement with counterfactual estimation and DTRs, refer readers to Appendix C.
depends on latent factors, making infeasible to directly obtain the interventional distribution and optimize ELBO,, Instead, we can optimize the following log- lig-likelihood objec- tive based on observational data:
O' = -log po Y x++
He, ar,7)
(4)
4. Methodology
Similarly, the ELBO for optimizing O' can be derived as:
In this work, we adopt a causal model as illustrated in Fig- ure 2. The model incorporates latent factors Z, at each timestep (s = t,..., ,t T) to represent the underlying
Pe Y ltr: Ze, He ae.r) (Z., He, Y Y ac.r)
ELBO, Eyg. (log
(5)
3

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
To bridge the gap between interventional and observational log-likelihoods (O and O'), we establish the following the- orem:
Algorithm Optimize Intervention Sequence
Require: Historical information He, target Y target trained
VCIP model, learning rate maximum epochs M Ensure: Optimal intervention sequence a
Theorem 4.1. Given the causal model illustrated in Fig-
ure 2, assume there exists constants 61, 62 > 0 such that:
I: Initialize at.r randomly
for epoch to M do
O ELBO1 Sc,, O' ELBO2 S c2.
Calculate CELBO using Equation 19
Then optimizing O can be approximated by optimizing:
āt,7 at,r a a., LELBO
5: end for
6: return ā7,- = āt,7
ELBO2-
Σ
Beg (a. Z.)] + le He) (6)
with the error bounded by 61 62-
the interventions are continuous with values bounded in [0, 1], we use Beta distributions:
The key proving Theorem 4.1 lies in leveraging the G- formula (Hedeker & Gibbons, 2006), which states; Property 1 (g- (g-formula). The expected potential outcome under intervention sequence ät, can be expressed as:
(a |)) Beta (Z.), (Z.))
(12)
Po (ar, |) IT Beta(af (H, ar.s). 8 ar,))) (13)
Po (Y = YotT H.)
When the interventions are binary, we use Bernoulli distri- butions:
polY YYtr tov HH, Z1, ) x (Z) H.)x
T++
a.
(a, 1-) Bern (Z.)).
(14)
Po | IT Bern (x) (H, a..)).
(15)
For the complete proof and detailed decomposition of the ELBO terms, refer readers to the Appendix.
where au, So, and To are implemented using MLPs, while on B1, and कi are implemented using LSTM networks capture the temporal dependencies in interventions.
4.2. Model Architecture
Generative model. Our causal model in Figure 2 assumes the following factorization of the generative model:
Inference model. We use an inference model qo (Zu:+r He Y Y++r) to approximate the posterior distribution of the latent variables, which can be factorized as:
Po Y t+r Ze,++ 1 H., ar,)) ) = pa (6 (Y,+ IZ+++)
11 Po (Z,+1 1 Zo )..r.) ) x po (Ze (H:),
go (Zucr He, Yiso, ac.)) = (Z Ht)
where 7, = r +1 -8 and each conditional distribution is parameterized as a Gaussian with neural network-derived parameters:
IT 20 (Z.41 Ye+r.a,.7. Z.);
(16)
where T, = I Each of the corresponding factors described as:
(4+ 1-) = NN 2, (Z+++)),
(9)
po 1) = =N(H, (Z., a., 0, (Z., ae, )),
(10)
2 (z( (He,, of (H))), qo (2.+++l))= (Y+, Z.), o (Y & Z,)), (11)
(17)
(Z. 1-) = -NV(L (H.), 0 (H.)).
(18)
Specifically, we employ a multilayer perceptron (MLP) for My and oy to decode the final system state Ze++ into the target output Y.+r. Given the sequential nature of histori- cal information H,, we utilize Long Short-Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) architectures for Mz, and oz to encode the initial system state Z₁ Simi- larly, the state transition functions (z, and o2 which map the current state and intervention the next state, are also implemented using LSTMs capture the iterative dynamics of the system.
where μ4, o are implemented by LSTM networks. Note that Z, depends on the descendant latent factors after time s. To implement this in practice, representation of intervention sequences a,, Specifically, at each time step 8, encode the current action a, and sub- sequent interventions through an LSTM network, producing the latent representation of the intervention sequence.
an LSTM to learn the
4.3. Traning and Inference
For the intervention distributions in Equation 6, we model them differently based on the type of interventions. When
To solve the optimal sequence in Equation 1, we first train VCIP maximize ELBO, under observational data. Ac-

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
cording the previous section, we define:
relevance, extend the binary intervention space to con- tinuous domain, capturing the dose-dependent nature of chemotherapy and radiotherapy treatments. The detailed data generation protocol is presented in Appendix.
++T
CELDO-
Σ
KL. (Z.l-)) (o((..)))
5=1
Bgg pog s+T Z[++]]
For real-world evaluation, we leverage the Medical Informa- tion Mart for Intensive Care III (MIMIC-III) database (John- son et al., 2016). which contains de-identified electronic health records from intensive care unit patients. The tar- get outcome variable is diastolic blood pressure under two concurrent interventions: vasopressor administration and mechanical ventilation support. The feature set comprises 25 time-varying physiological measurements and 3 static patient characteristics, following the preprocessing proto- col established in recent works (Hatt & Feuerriegel, 2021; Kuzmanovic et al., 2021; Melnychuk al., 2022).
(+T-1
+ A
Σ
Bag (og
Z.)] (o..r H₂) (19)
The first term on the RHS minimizes the KL divergence to align the inference model ge (Z.l)))) with the generative model pa (Zal-): the second term maximizes the probability of achieving the target output Y r- According to Theorem optimizing these two terms is equivalent maximizing the observational ELBO2- When À 1, the third term bridges the gap between interventional and observational log- oogliiklihoood.
5.1. Ranking-based Evaluation
During VCIP training, we split observational data into his- torical and future 7-step information to optimize param- eters ē and 0 by minimizing the evidence lower bound CELBO. After VCIP training, we fix all model parameters and optimize the intervention sequence as learnable param- eters. The specific process is shown in Algorithm 1 1. For cnatterfactual-sstmaatii methods, we employ a similar op- timization algorithm, but the objective is to minimize the distance between the predicted output and the target.
To evaluate the model's ability in ranking intervention se- quences, we design ranking- based evaluation framework. Specifically, for each test sample (He at,r: Yurga), where Y target Y [{,r,]] represents the output corresponding to the ground truth sequence, we randomly generate k candidate sequences {a{)) ji=1, with some generated by adding ran- dom perturbations to the ground truth sequence. For details, please refer to Appendix D. This results in evaluation set containing k + 1 intervention sequences including the ground truth sequence.
5. Experiments
In this section, conduct comprehensive experiments to evaluate the performance of VCIP against several state-of- the-art baselines. Specifically, we validate our model on both simulated tumor and real-world healthcare datasets to demonstrate its effectiveness in intervention ranking and counterfactual inference.
Based on this evaluation set, we obtain sequence rankings through different models. For baseline models capable of predicting conditional expectations Y [@r.-] TE| (( H.). we rank the intervention sequences in ascending order based on the deviation between predicted and target values Y [|["]] arggeel,, denoted as Fped A higher ranking position indicates that the intervention sequence is more likely to achieve the target output Y Yurget- For the VCIP model, we rank sequences based on the ELBO loss LELBO. where a lower value indicates a higher conditional likelihood H.).
Baselines. We compare our method against state-of-the- art models for counterfactual estimation over time: RMSN (Lim et al., 2018), CRN (Bica et al., 2020), CT (Melnychuk et al., 2022), and ACTIN (Wang et al., 2024). To ensure fair comparison, we perform hyperparameter tuning for all baselines (see details in Appendix.)
log Po (Y [(a,, (i) Y urrget
Based on the predicted ranking Trred we propose the Ground Truth Ranking Position (GRP) metric to evaluate the nor- malized ranking of the ground truth sequence:
Datasets. The tumor dataset simulates lung can- treatment dynamics through pharmacokinetic- pharmacodynamic framework (Geng et al., 2017). This dataset serves as a benchmark for evaluating causal infer- ence in sequential decision (L.im et al., 2018; Bica et al., 2020). A key advantage of this bio-mathematical model its confounding parameter y that modulates treat- ment history's 's impact on subsequent intervention decisions. Higher values of r correspond to stronger historical con- founding effects in treatment allocation. To enhance clinical
1: +
GRP
(20)
k
where & denotes the position of the ground truth sequence in the predicted ranking (from I to k + 1). A higher GRP value indicates better ranking position of the ground truth sequence, reflecting better model performance.
To comprehensively evaluate the model's 's ranking ability, we further propose the Ranking Correlation Score (RCS) metrie. Specifically, we first obtain the true ranking Tiree
Cooee available at: Bangxin0126/VVII..o
https://github. com/
5

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
r=2
VCP
ACTIV
RMSy
ORN
0.00
VCP
ACTBV
CT
RMON
CRM
1000
ACTIN
CT
RMSN
CRN
VCB
ACTB
RMSN
CRR
Figure Comparison of GRP across different models MIMIC-III dataset with varying prediction horizons 2, Higher GRP indicates better ranking performance.
==2
==4
=
10-
GR0
0.4
0.2
0.2
0.2
0.2
0.0
0.0
0.0
0.0
0.5
0.5
RCS
010
0.-
0.0
0.0
-0.5
0055
-0.-
-0.-
-10
-1.0
-1.0
-10
VCIP
ACTIN
CT
CRN
RMSN
VOIP
ACTIN
CT
OEN
RAIN
VCP
ACTI
CT
RMSN
VOP
ACTIV
C7
CEN
RMMN
Figure 4. Evaluation results on the tumor simulation dataset with varying prediction horizons (7 = 2, 2, 6, 8) and fixed confounding bias level (7 4). The upper shows GRP and the lower presents the RCS that the correlation between predicted and ground truth rankings (ranging from -1 to 1), where higher values indicate better performance in both metrics.
by sorting all sequences in ascending order based on their actual target distances Y | murge ||. RCS is defined as the Spearman correlation coefficient between the predicted ranking r Pped and the true ranking Fuse (Spearman, 1961):
of prediction errors. In contrast, VCIP, which directly com- putes ELBO using target outputs and intervention sequences, not only maintains but improves its ranking accuracy with larger т. This improvement may be attributed to the model's ability to leverage richer intervention information contained in longer sequences.
RCS prodd )
(21)
where -) denotes the Spearman correlation coefficient. The RCS value ranges from [-,, with values closer indicating better alignment between the predicted and true rankings, reflecting more accurate assessment of interven- tion sequence quality by the model.
On the simulated tumor dataset with fixed confounding bias (y 4), we compare the GRP and RCS metrics across different models, with results shown in Figure 4. The eval- uation demonstrates that VCIP significantly outperforms baseline models on both GRP and RCS metrics across all prediction horizons (T = 2, 4, 6, 8). Notably, the mean RCS maintains high level above 0.7, indicating that VCIP's predicted intervention sequence rankings strongly corre- late with ground truth rankings and accurately reflect how different intervention sequences approach the target state. Furthermore, as the prediction horizon increases, VCIP ex- hibits an increasing performance advantage over baseline models, demonstrating robust predictive stability.
Results. We first evaluate different models' GRP perfor- mance on the real-world MIMIC-III dataset, with results shown in Figure 3. Overall, our proposed VCIP model significantly outperforms baseline methods across all pre- diction horizons т, demonstrating its superiority in ranking intervention sequences. From the perspective of prediction horizon, baseline methods based on counterfactual estima- tion (ACTIN, CT, RMSN, and CRN) tend to suffer from per- formance degradation with larger т, due the accumulation
6

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Table 1. The long-range prediction results on the tumor dataset (7 = 4) under identical intervention strategies for both training and sets are reported as target distance (mean + std over five runs).
6
10
11
ACTIN
0.42
0014
0.71
0.13
1,0
= 0.21
1.30
+ 0.13
1.47
1.55
0.13
1.59 0.16
1.63 1.18
1.68
+ 0.21
CT
0.55 = 0.20
0.88
+ 0.25
143
0.36
1.69
4 0.27
1,87
1 0.39
2,01
+ 0.43
2.04
- 0.46
2
2.14
= 0.59
CRN
0.38 = 0.09
0.60
1 0.08
092
0.08
1.19
1 0.18
1.33
0.21
1,40
1 0.24
1.49
1 0.32
1.59
0.1
1.62
1 047
RMSN
0.30
0.45
1 0.10
0.75
0.16
0.98
4: 0.22
1.15 0.23
1.22
4 0.25
1.28
1 0.34
1.43
0.32
1.47
+ 0.35
VCIP
0.29
0.08
0.42 1:
0.60 0.15
1
0.92 0.24
0.95
0.24
0.99
1
1.44
0.27
1: 0.29
Table 2. The prediction results on the tumor dataset (7 = 4) with distinct intervention strategies applied to training and test sets are reported target distance (mean + std five runs).
2
7 6
T
9
10
T
11
r 12
ACTIN
0.45 -0133
0.78
1 0.09
1.15 1: 0.08
1.49
1 0.07
1.77 1: 0.23
1,90
1 0.25
1.94
1 0.29
2.06
1: 0.27
2.10 1 0.30
CT
= 0.19
0.3
0.21
1.56
- 0.25
1,98
4
0.37
2.34
043
2.33
0:45
2
0.40
4:
CRN
0.47 0.11
0.900.000
1.00
0,16
1.30
: 0.20
1.64
0.34
1.73
0.36
1.80
+ 0.24
1,91
0.29
2.05
+ 0.31
RMSN
0.37 0.11
0.53
+ 0.07
081
0.14
1.04
+ 0.24
1.34
134::00 0.30
1.48
= 0.31
1.59
+ 0.22
1.72
+ 0.28
1.83
= 0.22
VCIP
0.29 = 0.04
0.44 + 0.06
0.71 = 0.09
0.87 + 0.10
1.08 + 0.22
1.24 = 0.26
1.28 + 0.25
1.30
+ 0.18
1.36 + 0.26
5.2. Intervention Sequence Optimization
1.0
Target Distano -Y
BLBU Craso
0.8
shedoed Distance (V Yogel
When interventions are continuous-valued. the approach of selecting from a finite set of ranked intervention candi- dates may prove inadequate for capturing the full spectrum of possible interventions. To address this limitation, we propose gradient-based intervention sequence optimization algorithms for both the baseline methods and VCIP, In this section, we evaluate the models' capability to optimize in- tervention sequences for achieving desired outcomes under continuous intervention scenarios. Specifically, given a tar- get trajectory Y Yurget for an individual, we first obtain the optimal intervention sequence a using Algorithm 1, then compute the target distance || |||ar_] - uurger| between the potential outcome and the target trajectory. A smaller distance indicates that the obtained optimal sequence leads to potential outcomes closer to the target output.
Aherioe titrr (T J Vll
0.6
Nendlard
0.4
0.2
0.0
20
60
80
130
Surtal Sumpli Index
Figure 5. Individual- predictions under 100 candidate inter- ventions on the tumor dataset (r = 6, y 4) for VCIP and ACTIN, Each point represents normalized prediction score for specific intervention option.
Results. Due to the unobservable nature of counterfactu- als, we conduct evaluations exclusively on the simulated tumor dataset. We compare the performance of different models across varying confounding levels (7 1, 3) and prediction horizons (r = 1 to 6), with results illus- trated in Figure 6. The experimental results demonstrate that VCIP achieves comparable performance when confounding is minimal (? = 1) and shows superior performance as con- founding intensities increase. While the prediction errors of all models increase with both the confounding strength and prediction horizon, VCIP maintains clear advantage over baseline methods under stronger confounding (especially for 7 = 2, 3), exhibiting better stability and lower prediction errors in long-range predictions.
Case study. To investigate why baselines perform poorly in RCS ranking, we conduct a case study by selecting repre- sentative samples. As shown in Figure 5, the experimental results reveal a fundamental limitation of counterfactual methods based on predicted potential outcomes for rank- ing. Specifically, there exists a significant discrepancy be- tween the predicted target distance | (ai.)) Y mrges || and the true target distance [Y line). This discrepancy stems from the prediction error [| - Y a I. which cannot be evaluated in practice due to the unobservability of the counterfactual outcome Y| In contrast, VCIP's ELBO loss better preserves the monotonic relationship with the true target distance, resulting in more reliable ranking.
Y carget/ (grey dashed
To further investigate the impact of prediction horizon on model performance, we present the results for different pre- diction horizons
1 to 7 12) with ? 4 in Table 1.
7

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
1=1
0.12
VETP AUTIN
U.5
10-
VCBP
ACTINN
-0-
ACTIN
Mean distarce
0.80-
CT
0.08-
RMSN
RMSN
0-8
0.2
0.2
2
Figure 6. Comparison of target distances across different models on the tumor dataset under varying confounding levels. The x-axis represents different time steps, and the y-axis shows the mean target distance. Each subplot corresponds to different confounding levels.
Table 3. Ablation study results comparing model performance with and without confounding adjustment.
GPR
RCS
Target distance
r
7 >= 2
T
Y
RMSN
0.863 1: 0.181
0.796
1: 0.267
0.400
1 0.311
0.251
1: 0.299
0.078
0.017
0.301
0.131
0.599
1 0.204
0.985
1. 0.215
RMSN w/o adjustment
0.797 +: 0.318
0.747
0.348
0.461
+ 0.460
0,213
0.407
0.087
0.020
0.280
0.144
0.694
0.201
1.132 4: 0.090
VCIP
0.944
0.091
0.972
0.063
0.772
1: 0.206
0.156
0.192
0.382
0.085
0.746 1. 0.198
VCIP w/o adjustment
0.791 + 0.334
0.796
0.356
0.566 + 0.508
0.595
0.518
0.092 + 0.080
0.284 + 0.4
0.756 + 0.299
0.912 + 0.194
As shown in the table, the target distance increases with longer prediction horizons across all models. VCIP consis- tently achieves the lowest target distance across all predic- tion horizons, outperforming all baselines by substantial margin. Specifically, for short-term predictions (r 1), VCIP achieves a target distance of 0.29, showing modest improvement of 3.3% compared to RMSN (0.30). The per- formance gap widens significantly for long-term predictions. When - = 12, VCIP maintains a target distance of 1.09, outperforming the best baseline RMSN (1.47) by 25.9%, demonstrating VCIP's superior capability in handling long- range predictions under strong confounding (7 4).
5.3. Abalation Study
Previous research has shown that addressing confounding bias through proper adjustment is crucial for accurate coun- terfactual estimation in longitudinal settings. In this work, we employ g-formula for adjustment, specifically using the last term in Equation (19) to bridge interventional and ob- servational log- logliikliioooss. To evaluate the effectiveness of adjustment, we conduct ablation studies by setting À = O in Equation (19) and design corresponding ablations for RMSN. Table presents the results for both ranking tasks (at 7 = 4) and sequence optimization tasks (with r 6 evaluated across different 7 values).
When applying identical intervention strategies for both training and test sets, VCIP demonstrates superior perfor- mance as shown above. To further evaluate the generaliza- tion capability of different models under more challenging scenarios, we conduct experiments where the intervention strategies differ between training and test sets (see Appendix for details of strategy generation). The results are presented in Table 2. Even under this more challenging setting, VCIP maintains its superior performance across all prediction hori- zons. For short-term predictions (r 1), VCIP achieves a target distance of 0.29, outperforming the best baseline RMSN (0.37) by 21.6%. The advantage becomes more pro- nounced for long-term predictions, where at 5 = 12, VCIP achieves target distance of 1.36, showing significant improvement of 25.7% compared RMSN (1.83). These results show VCIP's strong generalization capability when facing unseen intervention strategies.
The ranking results demonstrate that adjustment signifi- cantly enhances the performance of both models. For VCIP, incorporating adjustment improves the GPR metric to 0.944 and 0.972 at 7 = 2 and す = 4 respectively (compared to 0.791 and 0.796 without adjustment), while RCS shows stronger ranking capability (0.772 and 0.869 vs. 0.566 and 0.595). Similarly, RMSN's GPR improves from 0.797 and 0.747 0.863 and 0.796 respective horizons, accompa- nied by enhanced RCS scores. These results indicate that adjustment enables more accurate evaluation and ranking of intervention sequences.
From the target distance perspective, adjustment helps both models achieve more stable performance under strong con- founding. Notably, at high confounding levels (? =3, 3, 4), VCIP with adjustment achieves target distances of 0.382 and 0.746. substantially outperforming its counterpart with- out adjustment (0.756 and 0.912). RMSN exhibits simi- lar trend, with adjustment yielding better results at = 4
8

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
(0.985 vs. 1.132). This confirms that adjustment effectively handles confounding effects, enabling the models identify superior intervention sequences. Interestingly, the perfor- mance gains from adjustment are less pronounced under low confounding conditions (? = 1,
391(10128):1334-1333 2018.
Bica, I., Alaa, A. M., Jordon, J., and van der Schaar, M. Estimating counterfactual treatment outcomes over time through adversarially balanced representations. Inter- national Conference on Learning Representations. 2020.
6. Conclusion
Chakraborty, S., Bhattacharya, S., Garg, S., and Subrama- nian, D. A survey on causal inference in healthcare using electronic health records. ACM Computing Surveys, 55 (1):1-35, 2022.
This paper introduces the novel counterfactual target achievement problem, which aims to identify intervention sequences that guide systems toward desired target out- comes, addressing a critical need in many real-world applica- tions like personalized healthcare. To address this problem, we propose the VCIP framework that models the likelihood of achieving target states through intervention sequences. Unlike traditional counterfactual estimation methods that accumulate prediction errors, VCIP provides more robust solution by learning target achievement probability. Exten- sive experiments on both synthetic and real-world datasets validate VCIP's effectiveness.
Geng, C., Paganetti, H., and Grassberger, C. Prediction of treatment response for combined chemo-and radiation therapy for non-small cell lung cancer patients using a bio-mathematical model. Scientific reports, 7(1):13542, 2017.
Hatt, T. and Feuerriegel, S. Sequential deconfounding for causal inference with unobserved confounders. arXiv preprint arXiv:2104.09323, 2021.
Hedeker, and Gibbons, R. D. Longitudinal data analysis. John Wiley & Sons, 2006.
Acknowledgements
This research is supported by National Key R&D Program of China (No. 202IZD0111700 1170), National Nature Science Foundation of China (No. 62137002, 621762,5), Special Foundation for Science and Technology Innovation and En- trepreneurship of CCTEG NO. 2020-2-TDD-CXY006, Key Science and Technology Special Project of Anhui Province (No. 202103a07020002). We thank the anonymous review- for their constructive comments that help improve the manuscript.
Hochreiter, S. and Schmidhuber, J. Long short-term memory.
Neural computation, 1735-178,, 1997.
Hochreiter, S., Bengio, Y., Frasconi, P., Schmidhuber, J., et al. Gradient flow in recurrent nets: the difficulty of learning long-term dependencies, 2001.
Jensen, P. B., Jensen, L. J., and Brunak, S. Mining electronic health records: towards better research applications and clinical care. Nature Reviews Genetics, 13(6):395-405. 2012.
Impact Statement
Johnson, A. E., Pollard, T. J.,, Shen, L.,, Lehman, L.-w. H., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., An- thony Celi, L., and Mark, R. G. Mimic-iii, a freely ac- cessible critical care database. Scientific data, 3(1):1-9, 2016.
This research advances temporal causal inference by intro- ducing a framework for identifying intervention sequences that guide systems toward specific target outcomes, with ap- plications in healthcare and economic mikiin.. While method can help practitioners make more informed de- cisions about sequential interventions to achieve desired targets (e.g., normalizing physiological indicators or stabi- lizing economic metrics), we emphasize that should serve as a decision support tool rather than replace expert judg- ment. The framework's 's effectiveness depends on training data quality and representativenesss and careful consider- ation must be given to potential biases and the need for ongoing validation before deployment in critical real-world applications.
Kuzmanovic, M., Hatt, T., and Feuerriegel, S. Deconfound- ing temporal autoencoder: estimating treatment effects over time using noisy proxies. In Machine Learning for Health, pp. 143-155. PMLR, 2021.
Laber, E. B., Lizotte, D. J., Qian, M., Pelham, W. E., and Murphy, S. A. Dynamic treatment regimes: Technical challenges and applications. Electronic journal of statis- tics, 8(1): 1225, 2014.
Li, R., Hu, S., Lu, M., Utsumi, Y., Chakraborty, P., Sow, D. M., Madan, P., Li, J., Ghalwash, M., Shahn, Z. and Lehman, L.-w. G-net: recurrent network approach to g-computation for counterfactual prediction under a dynamic treatment regime, In Proceedings of Machine
References
Bergenstal, R. M. Continuous glucose monitoring: trans- forming diabetes management step by step. The Lancet,

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Learning for Health, volume 158 of Proceedings of Ma- chine Learning Research, pp. 282-299. PMLR, 04 Dec 2021.
Robins, J. M., Hernan, M. A., and Brumback, B. Marginal structural models and causal inference in epidemiology. Epidemiology, pp. 550-560, 2000.
Lim, B., Alaa, A., and van der Schaar, M. Forecasting treatment responses time using recurrent marginal structural networks. Advances in neural information pro- cessing systems, 31, 2018.
Roy, J., Lum, K. J., and Daniels, M. J. A bayesian non- parametric approach to marginal structural models for point treatments and continuous or survival outcome. Biostatistics, 18(1):32-47, 2017.
Luckett, D. J., Laber, E. B., Kahkoska, A. R., Maahs, D. M., Mayer-Davis, E.,, and Kosorok, M. R. Estimat- ing dynamic treatment regimes in mobile health using v-learning. Journal of the American Statistical Associa- tion, 115(530):692-766, 2020.
Rubin, D. B. Bayesian inference for causal effects: The role of randomization. The Annals of statistics, pp. 34-58, 1978.
Soleimani, H., Subbaswamy, A., and Saria, S. Treatment- response models for counterfactual reasoning with continuous-time continuous-valued interventions. In 33rd Conference on Uncertainty in Artificial Intelligence, UAI 2017. AUAI Press Corvallis, 2017.
Luo, Z., Pan, Y., Watkinson, P., and Zhu, T. Position: Rein- forcement learning in dynamic treatment regimes needs critical reexamination. Proceedings of the 41st Inter- national Conference on Machine Learning, 2024.
Spearman, C. The proof and measurement of association between two things. 1961.
Melnychuk, V., Frauen, D., and Feuerriegel, S. Causal trans- former for estimating counterfactual outcomes. In Inter- national Conference on Machine Learning, pp. 15293- 15329. PMLR, 2022.
Tsittsis, S. and Gomez Rodriggz, M. Finding counterfactu- ally optimal action sequences in continuous state spaces. Advances in Neural Information Processing Systems, 2023.
Mortimer, K. M., Neugebauer, R., Van Der Laan, M., and Tager, B. An application of model-fitting procedures for marginal structural models. American Journal of Epidemiology, 162(():382-388, 2005.
Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. At- tention is all you need. Advances in neural information processing systems, 30, 2017.
Murphy, S. A. Optimal dynamic treatment regimes. Jour- nal of the Royal Statistical Society Series B: Statistical Methodology, 65(2):331-355 2003.
Wang, S., McDermott, M. B., Chauhan, G., Ghassemi, M., Hughes, M. C., and Naumann, T. Mimic-extract: A data extraction, preprocessing, and representation pipeline for mimic-iii. In Proceedings of the ACM conference on health, inference, and learning, pp. 222-235, 2020.
J. Causality. Cambridge university press, 2009.
Peterson, K., Silverstein, J., Kaufman, F., and Warren- Boulton, E. Management of type 2 diabetes in youth: an update. American family physician, 76(5):658 664, 2007.
Wang, X., Lyu, S., Yang. L., Zhan, Y., and Chen, H. A dual- module framework for counterfactual estimation over time. In Forty-first International Conference on Machine Learning, 2024.
Prosperi, M., Guo, Y., Sperrin, M., Koopman, S., Min, J. S., He, X., Rich, S., Wang, M., Buchan, I. E., and Bian, J. Causal inference and counterfactual prediction in machine learning for actionable healthcare. Nature Machine Intelligence, 2(7):369-375, 2020.
Wu, J. Q. Horeweg, N., de Bruyn, M., Nout, R. A., Jürgenliemk-Sehhn I. M., Lutgens, L. C., Jobsen, J. J., Van der Steen- Baans-Ba,,s,,, E. M., Nijman, H. W., Smit, T., al. Automated causal inference in application ran- domized controlled clinical trials. Nature Machine Intel- ligence, 4(5):436 444,444444, 2022.
Raghu, A., Komorowski, M., Celi, L. A., Szolovits, P., and Ghassemi, M. Deep reinforcement learning for sepsis treatment. arXiv preprint arXiv: /7/ 1.09602, 2017.
Robins, J. A new approach to causal inference in mortality studies with a sustained exposure crrriddaaplicattiion to control of the healthy worker survivor effect. Mathemati- cal modelling,
Xu, Y., Xu, Y., and Saria, S. A bayesian nonparametric ap- proach for estimating individualized treatment- responss curves. In Machine learning for healthcare conference, pp. 282-300. PMLR, 2016.
1986.
Robins, J. M. and Hernán, M. A. Estimation of the causal ef- fects of time-varying exposures. CRC Press, Boca Raton, FI., 2009.
Zhang. J. and Bareinboim, E. Near-optimal reinforcement learning in dynamic treatment regimes, arXiv preprint arXiv: / 1905 059505 2019.
10

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Zheng, C.,, Dai, R., Gale, R. P., and Zhang, M.-J. Causal inference in randomized clinical trials. Bone marrow transplantation, 55(1):4-8, 2020.

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
A. Assumptions
For accurate estimation treatment effects using observational data, we consider three fundamental assumptions that form the theoretical foundation of our work. These assumptions align with established frameworks in causal inference literature (Lim al., 2018; Bica al., 2020; Li al., 2021; Melnychuk al., 2022; Wang et al., 2024):
Assumption A.1 (Consisteney). The observed outcome 141 matches the potential outcome Y [ar] when treatment a, applied at time expressed as Y 4+1 Y++1 (a.).
Assumption A.2 (Sequential Overlap). For any treatment a, at time there exists non-zero probability of assignment conditioned the historical trajectory h; Mathematically, 0 < P(A, at H. h.) Va, A when P(H, h₂) > 0.
Assumption A.3 (Sequential Ignorability). Treatment assignment at time exhibits conditional independence from potential outcomes at + 1, given the observed history: At 1 Y [ar] H,, Va, E A. This ensures no unmeasured confounding affects the weatment-outconn relationship.
B. Proof
We first derive the Evidence Lower BOund (ELBO) for both the log interventional likelihood and log observational likelihood.
For the log interventional likelihood:
log Po Y [ar.r] Y target H,)
= log
pe Y Yarget
H,, ac,r, Ze,r-1) po(Ze H. s0t1 Pe Z a as-1) 1)d2e,411
He, Y target, az,7 pol Y target, Ze, He, (a))))))) dZz,,41
log
2ol (Z (,,41 H,, Y targers a1,7)
Po Y target; Ze, HH, do (at.r)) = ELBO
Egg [log
(22)
qo (Z, He, Y Larget, ac,7)
For the log observational likelihood, we derive two different formulations. The first formulation is:
log Y -Yuagee HH, at.-)
=log| (Yuggr Ze, Hy, 2e,7) dZt.s-1
log / qos (o(Zerri Y target Rt.7
Pe Ytarget: Ze, H ag,7) 2o (Z [1, H,, Y targett ar,)) dZz,T4
( target: Ze. H )
Egg
log qol (Z,,+1 H,, Y Yimrgats ar,r)
ELBO2
(23)
12

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
The alternative formulation of the log observational likelihood is:
log po Y t+r Y targer He, at.7)
log po YY Turget, H,, āt,-) log po (H,, ag,r)
log
po Y larget H1, a Zt, ) 11 log po (He, ās.7)
log
po uurget
Hr, r )po Hi, r, Zt, 11, ddZ log Po as,7)
log
pol Y urget
Re, Re, z1, Po (H, H₂) 11 8-141 Po(Zs Z,-1 a u1 (a Zo)dZe,7+4
log (H,,
log
po Y uugat Zt He, do(a),)))) po (Ht) I (a Z,)dZi, - log po (He, aH,)))
/ [9o (Z, He. Y Larget: ā,,7)
po Y target. Zt, He, do(ā)r)))
log
pe (H,) 11 (a Zz)dZe, - log pe (H,, ar,))
9ol [t,
H, Y Yuuger, a
Po Po(Yunge Iugget: Z1,7+1 H., do (ar,)))
Bgg [log
(H,))(H,,)
log po 2.)] log po H,,
qo (Zi r+1 He, Y usrget: )
s=t
Po Y target:
Z Z Hy, do(ār,))
Egg [log
po
2.)] - -o H,)
Tol Zz, r+1
He, Y Lugget: a
ELBO +
Σ
Eg log (a, [.)] - log pa (āter
Hc)
(24)
s=t
The fifth equation follows from Property 1.
Theorem B.1. Given the causal model illustrated in Figure 2, assume there exists constants E1, 62 > 0 such that:
O ELBO, S1, O' - EELOO S 62
Then optimizing O can be approximated by optimizing:
++r-1
ELBO₂
Σ
Egg [log Zz]] + log pap (āt,r H.)
(25)
with the error bounded by 61 + E2-
Proof. From the alternative formulation of the observational distribution, we have:
log (Y target
H, ār,+) 2 ELBO
Eg [log Z.)] - log (ar.r H,)
log Po (Y target
H, ār,r)
2 ELBO2
By the theorem assumptions:
O ELBO & 61,
E1,
O' ELBO2 S 62
Substituting the first inequality into the first line of the observational distribution expression:
O < ELBO, +
t+ Egg [log Pa Z.)] log pe (ar,r H₂) + E1
< log pa Y target ar,r)
H+
(ELBO2 + (2)
D
gg [log (a
Z.)] log po (ar,r H.) +
13

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Therefore:
H Σ
O (ELBO2
Eg Po
Z.)] + pa H,)) S 61 + c2
This establishes that optimizing the expression on the right-hand side approximates optimizing O, with an error bounded by
6 + C2
C. Problem Comparisons
To highlight the distinct characteristics of our proposed counterfactual target achievement problem, compare it with two widely-adopted problems in sequential decision decision-making: counterfactual estimation and dynamic treatment regimes (DTRs). While all three approaches address temporal treatment planning, they differ fundamentally in their objectives, methodologies, and levels of analysis, as summarized in Table 4.
C.1. Counterfactual Estimation
The counterfactual estimation framework focuses predicting potential outcomes under specific treatment sequences. Given patient's history H; = (X,, A.- Y,, V), the goal is to estimate the expected outcome:
B [a.- H.)
(26)
This framework enables individual-lever outcome prediction but does not directly optimize treatment sequences.
C.2. Dynamic Treatment Regimes
DTRs provide formal framework for personalizing treatments in sequential decision-making while optimizing the population-level outcomes. Unlike traditional fixed treatment protocols, DTRs adapt treatment decisions based on each patient's evolving characteristics and treatment history, yet aim to maximize the expected outcome across the entire patient population.
A DTR is defined as structural causal model M = (U, V, P(u)), where U represents unobserved exogenous variables, and V = {Kk, Sk. Y} comprises treatment variables Xx = (Xi...., Xx} over K stages, time-varying covariates S = .S....., Sx1, ), and a primary outcome Y E [0, 1]. The structural functions F determine how variables evolve: state transitions follow S₁ - Th (āk-1, Sk-1, u). treatment assignments are given by X₃ + fu (5k. øx-1, u). and the outcome is generated as +- r (Tx- , SK: u). The distribution P(u) specifies the underlying randomness in the system.
1K treatment policy consists of sequence of decision rules (Th === where each T& Sx Xx-1 P(X) maps the history of states and treatments distribution over the next treatment. When executed, the policy induces interventional distribution:
(x Bx, = Pax (v|5k)
II k=0
Ps. (54411 (k41,5a)
(27)
The core objective in DTR optimization is to find a policy x² that maximizes the expected outcome across the entire population:
argmax F V2 (MM) an
(28)
While DTRs have been extensively studied in the medical decision-making literature, our work addresses fundamentally different problem. There are two key distinctions between DTRs and our proposed problem:
Optimization Objective: DTRs aim to learn optimal treatment policies that maximize population-level outcomes (Y)) In contrast, our problem focuses on identifying treatment sequences for specific individuals that maximize the probability of achieving target state, formulated as a" argmax,22 (Y Y [ā,rr] Y Yurge: H.).
14

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Table Comparison of key characteristies different temporal treatment planning prohlems. Our counterfactual target achievement problem differs from existing approaches in its objective function, output form, analysis level, and decision timeline.
Counterfactual Target
Counterfactual Estimation
Dynamic Treatment
Aspect
Achievement (ours)
Regimes
Objective
max V+(M)
maxa p Y (
=Ytagge Larget
H.)
E Y ++ [+x.]]
H.
Function
max B, [Y]
Potential outcomes under treatment sequence
Decision rules
Output
Complete treatment plan a
{T : Sk x 24-1 +
Form
a
P(R4)]
Analysis
Individual trajectory
Individual-level outcome prediction
Population-level
Level
optimization
optimization
Decision
Estimate future outcomes given treatment sequence
Plan full sequence in advance
Make decisions sequentially
Timeline
While DTR policies allow for individual-level adaptation, their optimization is fundamentally conducted the population level. Our approach operates purely at the individual level, making predictions and optimizing treatment sequences based each patient's unique historical trajectory H,.
Analysis Level:
These fundamental differences in objectives and methodology explain why traditional DTR evaluation metrics and com- parison frameworks are not directly applicable our work. Our problem more closely aligns with personalized trajectory planning, where success is measured by the accuracy of individual-level predictions and the achievement of patient-specific target states.
D. Datasets description
D.1. Synthetic tumor dataset
Research presented (Geng et al., 2017) employs a Tumor Growth (TG) simulator to forecast tumor volume progression over + 1 days following initial cancer detection, generating single-dimensional outputs. The framework encompasses two therapeutic interventions: radiation treatment (A,) and pharmaceutical therapy (Af). Our adaptation transforms the intervention variables from discrete continuous values between 0 and 1. The therapeutic impacts manifest differently: radiation demonstrates immediate effects d(t) on subsequent measurements, while pharmaceutical intervention exhibits extended influence C(t) across multiple time points, as expressed by:
K Y = 1 plog Y. - BoC(1) (a,d(t) (,d())))++++++++++++++)))))))+++ + Y..
(29)
The simulation parameters p. K, Bc, Dry, Br are predefined, with representing random variation drawn from N(0., 0.01 Patient-specific responses are characterized through parameters fc, Orr, Br. derived from three-component truncated normal mixture distribution. The mixture components serve as unchanging patient characteristics. Detailed parameter specifications are available in the anonymous repository. Time-dependent confounding affects both interventions through biased assignment protocols. The allocation probabilities for both treatments follow beta distributions:
A Aç ~ Beta(2n;; 20x))
(30)
where
=
D Dmmav
(Dis(Y,-1) - Dmaz/22)
(31)
15

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Here, a(-) denotes sigmoid activation, D Dmaz represents maximum tumor size, D1s (-1) indicates average tumor dimensions a 15-day window, and Y controls confounding intensity. When y equals zero, treatment assignments become purely random, while increasing values intensify time-varying confounding effects. In our implementation, d, and Ct are modeled through cubic spline transformations of A; and AG, respectively:
d(t) (A](,
(32)
C(I) (p),()),
(33)
where v, and the represent cubic spline transformations for the respective treatments. This nonlinear functional approach creates more realistic treatment-response relationships.
To ensure different strategies between testing and training phases, during the testing stage, interventions are randomly selected (with probability 17) to follow an independent strategy, where Af Aç Beta(a, B), with and | being constants independent of historical data.
The dataset comprises 1,000 training trajectories, 100 validation sequences, and 100 test cases across different confounding levels 7.. Individual trajectories extend up to 60 time steps, though early termination may occur due to patient outcomes.
Following established methodologies (Bica et al., 2020; Melnychuk et al., 2022), performance evaluation utilizes normalized target distance, calculated relative to the maximum tumor volume Vmaa = 1150 cubic centimeters.
D.2. MIMIC-III Clinical Dataset
The Medical Information Mart for Intensive Care (MIMIC-III) database (Johnson et al., 2016) constitutes a rich repository of de-identified clinical records from intensive care units. This comprehensive database encompasses diverse healthcare data elements including physiological measurements, medication administrations, laboratory results, clinical documentation, diagnostic classifications, and patient outcomes. To ensure systematic data processing and reproducibility, employ the MIMIC-Extract framework (Wang et al., 2020), which implements standardized preprocessing procedures for the MIMIC-III database.
In alignment with recent methodological developments (Hatt & Feuerriegel, 2021; Kuzmanovic et al., 2021; Melnychuk et al., 2022), our analysis incorporates 25 longitudinal physiological indicators and 3 time-invariant patient characteristics. These features, particularly the categorical variables, undergo one-hot encoding transformation to facilitate numerical analysis. The intervention space comprises two binary treatment decisions commonly encountered in critical care settings: the administration of vasopressors and the implementation of mechanical ventilation support. We designate diastolic blood pressure our primary outcome measure, given its clinical significance and sensitivity to both therapeutic interventions, thereby providing meaningful insights into patient trajectory management.
The experimental cohort consists of 5,000 patients with intensive care episodes spanning between 30 and 60 hours. This dataset undergoes a strategic partition into training (70%), validation (15%), and testing (15%) subsets. For temporal predictions extending steps ahead (r 2 2), we establish Tmax as an upper bound on the prediction horizon. We extract sub- bob-trajectories of minimum length Tmax + 1 using rolling windows, with temporal alignment achieved by removing initial observations up to (() Tiax 1, thus maintaining temporal causality. For immediate-horiin predictions (r = 1), utilize complete trajectory sequences from the test cohort.
D.3. Candidate Generation
We employ hybrid approach generate candidate sequences. Our framework creates random sequences (50%-80% of candidates) and perturbed ground truth sequences (20%-50% of candidates). The perturbation strategy is treatment-mode specific:
For discrete interventions: We randomly flip bits in the ground truth sequence with probability 0.2.
For continuous interventions: We apply context-award shifts where values are modified based on their magnitude (low values shifted up, high values shifted down, middle values shifted randomly).
16

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Table Model comparison results across different interference levels (γ) and prediction horizons (r), with same intervention strategies applied during training and testing phases, Results are reported as target distance values (mean + std over five runs).
7
== 2
==
==
T=
==
ACTIN
0.05++0.03
0.04+0.02
0.08+0000
0.10=0.07
0.12=0.11
0.08=0.04
CT
0,04+0003
0.05+0001
0,09+0.05
0.11=0.09 0.09=0.06
0.11=0.06
CRN
0.04.10.03 0.0410002
0.06110.3
0.08
0.09
0.08
RMSN
0.04.10.02 0.05:10.02
0.08
0.08
0.11=0.07
0.08=0.02
VCIP
0.02.10.01
0.04
0.07
0.07
0.10=0.08
0.10-0007
7 2
ACTIN
0.10::0.02 0.15+:0.07
0.22
0.27+0012
0.32=0.16
0.36=0.16
CT
0.13++0.05 0.21+0.19
0.27+0.18
0.38
0.42=0.29
0.50
CRN
0.10++0.04 0.16+00.88
0.23++0.10
0.28+0017
0.32-0.17
0.37-0.23
RMSN
0.07 0.11
0.17
0.20
0.31=0.23
0.30=0.13
VCIP
0.06110.02
0.0910006
0.1311000
0.1611.09
0.18-0.13
0.19-0.10
3
ACTIN
0.2410007 0.4010.11
0.51
0.6010.11
0.66=0.10
0.74
CT
0.36.10.12 0.5410.11
0.7010.15
0.8410.17
0.97
1.04_0.17
CRN
0.23 0.36
0.4710.10
0.5710.10
0.15
0.71
RMSN
0.16++0.05 0.25+
0.32
0.41
0.51=0.15
0.60=0.20
VCIP
0.15+0.04
0.18+0.06
0.29
0.34
0.38
0.38
7 7= 4
ACTIN CT
0.42 0.7110.13
0.91
1.05
1.19-0.10
1.10
0.55 0.85
0.18
0.3
1.58=0.22
1.69=0.27
CRN
0.38 0.60
0.79
0.92
1.11=0.23
1.19=0.18
RMSN
0.30:10.10 0.45
0.61
0.75
0.89=0.28
0.98=0.22
VCIP
0.29++0.08
0.42+00.13
0.53+
0.60
0.68-0.19
0.75-0.20
E. Appended Results
This section presents additional experimental results of different models on the tumor simulation dataset. The evaluation was conducted from main aspects focusing on model performance under identical and different intervention strategies.
Under identical intervention strategies (Table 5), the results demonstrate that the interference level y increases (from 1 to 4), prediction errors of all models show an upward trend, indicating that stronger interference reduces prediction accuracy. The VCIP model achieved optimal performance in most cases, particularly showing significant advantages at higher interference levels (y=3,4), As the prediction horizon extends, prediction errors increase across all models, which aligns with the intuition that longer-term predictions are more challenging.
Under different intervention strategies (Table 6), prediction errors increased compared identical intervention scenarios, reflecting the generalization challenges models face when encountering unseen intervention strategies. VCIP and RMSN models demonstrated better adaptability, especially in short-term predictions (7=1,2). At high interference levels (y=3,4), performance differences between models became more pronounced, with VCIP maintaining relatively stable predictive capabilities.
Figures 7 through 10 show that the subtle consistency between GPR and RCS metrics under both identical and distinct intervention strategies suggests that the ranking process itself may introduce variations in intervention patterns, potentially affecting the comparative assessment model performance across different evaluation metrics.
F. Hyperparameter Tuning
The hyperparameter settings used in our experiments are detailed in Table 7. For all baseline models including RMSN, CRN, CT, and ACTIN, we followed the hyperparameter optimization strategy and ranges consistent with those reported in (Wang et al., 2024). Specifically, hyperparameter optimization was conducted through random grid search. Readers can refer to (Wang et al., 2024) for detailed hyperparameter configurations of these baseline models.
17

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
Table Model comparison results across different interference levels (7) and prediction horizons (T), with distinct intervention strategies applied during training and testing phases, Results are reported as target distance values (mean + std over five runs).
== 2
==
==
==
==
ACTIN
0.15+0.09
0.18+0.110
0.27+0.20
0.33+0.22
0.40=0.28
0.45=0.33
CT
0.17+0.07
0.23+0.10
0.38+0.12
0.45
0.50=0.17
0.59=0.23
CRN
0.12.10.04 0.1510.04
0.1910.07
0.2710008
0.30-0.07
0.36=0.09
RMSN
0.10.10.03 0.14100.2
0.20.10.05
0.28
0.33
0.36=0.09
VCIP
0.10.10.05
0.14
0.1911.110
0.27
0.24
0.22-0.08
7 2
ACTIN
0.17+:0.05 0.28
0.41
0.47
0.60=0.18
0.60=0.16
CT
0.21+00.88 0.32+0.111
0.49
0.58+0.11
0.70=0.26
0.74
CRN
0.19+00.55 0.30++.111
0.4140.16
0.49-0011
0.62-0.28
0.60-0.21
RMSN
0.1310006 0.20
0.2910.07
0.39
0.53-0.23
0.50=0.17
VCIP
0.1110006
0.25
0.25
0.38
0.35-0.36
0.35-0.10
3
ACTIN
0.4410035 0,6710.45
0.9010.64
0.96
0.58=
1.10=0.58
CT
0.52.10.13 0.8110.21
1.1110.88
1.27
1.41 0.32
CRN
0.31 0.5
0.7110.22
0.82.0022
0.21
1.04
RMSN
0.22++0066 0.35
0.50
0.55
0.59=0.06
0.69=0.09
VCIP
0.17
0.39
0.42
0.49
0.49=0.16
0.62=0.25
7 7= 4
ACTIN CT
0,45.10.13 0.78
1.001..11
1.15
1.38
1.49-0.07
0,62#0.19 0.93+0221
1.30
0.56
1.81=0330
1.98=0228
CRN RMSN
0.47 0.69
0.90
1.03
1.25=0.13
1.30=0.20
0.37 0.53
0.68
0.81
0.96=0,22
1.04=0.24
VCIP
0.29++0.04
0.44
0.61+-0.11
0.71
0.77-0.11
0.87
Table 7. Specified ranges for hyperparameter tuning of VCIP across various datasets.
Hyperparameter
Range (tumor)
Range (MIMIC-III)
Learning rate Minibatch size Representation size hidden size (Generative model)
0.01, 0.001
0.01. 0.001
128, 256, 512, 1024
128, 256, 512, 1024 8, 16, 32
8, 12, 16
8, 12, 16
8, 16, 32
FC hidden units (Generative model)
8, 12, 16
8, 16, 32
hidden size (Inference model)
8, 12, 16
8, 16, 32
FC hidden units (Inference model)
8, 12, 16
8, 16, 32 0, 0.1, 0.2,
Dropout rate
O, 0.1, 0.2, 0.3
Random search iterations
30
30
Number of epochs
100
100
G. Limitations
Our approach hinges on the standard causal identification assumptions consistency, positivity and sequential ignorability. In particular, positivity requires that each candidate intervention has sufficient support in the observational data, and violations of this assumption that all confounders affecting both treatment assignment and outcomes are observed and correctly modeled; unobserved confounding remains an open challenge and may introduce bias into both estimation and downstream planning.
lead to severe degradation in optimization performance. Sequential ignorability assumes
Beyond these identification requirements, practical deployment poses additional challenges. The variational inference and sequential planning stages are computationally intensive, making real-time decision support and scaling to large cohorts or long horizons difficult without further algorithmic improvements. Clinical time-series data are often sparse, irregularly sampled and plagued by missing values; while imputation and generative augmentation can help, they introduce extra hyperparameters and potential sources of bias. Finally, personalized intervention policies must be accompanied by rigorous fairness, transparency and privacy aafeguards-tthice and governance issues that our current framework does not explicitly address but which are critical for roal-world applications.
18

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
CRR
0.0
11
1110
0110
0.7
0.5
0.5-
0.0
08.
05-
-10
ACTN
RMON
-10
VEPP
ACTM
RMO
-10-
VElP
ACTM
CT
RMOy
01
RMON
(a)
T=2
T24
7=6
T=8
1.0
131
141
3.0
GRR
0.4
0.4
1.0
3.0
0.00
61
105
-0.3
VCF
ACTIN
RMNN
VCB
ACTIN
RMo
VGP
ACTIN
RMO
YGB
ACTIN
RMN
(b) y
Figure 7. Evaluation results (Part 1) the tumor simulation dataset with varying prediction horizons 2, 4, 6,8)) under different confounding bias levels (7 = 2). The real intervention sequence in test set follows the identical intervention strategy as training set. The figures show GRP metrics that model performance, where higher values indicate better predictive capability.
19

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
T=2
10
165
11
015
CRP
0.2
02-
0.5
RCS
0.5
06
0.6
ACTNV
D25
RMNN
1DP
ACTN
Run
ACT
RMo
VOP
ACTM
RKDN
(a)
3
2=2
T=6
CRRP
0.-
38
D6-
0.4
0.4
0.0
000
0.5
0.5
RCS
0.09
0.0
300
008
1.0
10
1.
1.
ACTDN
RMa
veB
ACTEN
C
RMbo
ACTIN
Raoo
Yog
ACTDC
CT
Ratos
(b) 7
Figure 8. Evaluation results (Part 2) on the tumor simulation dataset with varying prediction horizons (r = 2,4,6, 8) under different confounding bias levels (7 3, 4).
20

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
GRR
0.0
11
0110
0110
0.7
0.5
0.5-
0.0
08.
05-
-10
VE3P
ACTN
C7
RMDN
-10
VEPP
ACTM
RMD
-10
VelP
ACTM
CT
RMOy
01
RMON
(a)
T=2
T24
T=6
T=8
1.0
131
140
3.0
GRRP
0.4
0.4
D01-
1.0
3.0
-0.
VCF
ACTIV
RMBV
VCB
ACTIN
RMo
VGP
ACTIN
RMo
YGB
ACTIN
RMNN
(b) y
Figure Evaluation results (Part 1) the tumor simulation dataset with varying prediction horizons 2, 4, 4,4,66 6, under different confounding bias levels (7 = 1, 2). The real intervention sequence in test set follows a different intervention strategy from training set. The figures show GRP metrics that model performance, where higher values indicate better predictive capability.
21

---

Variational Counterfactual Intervention Planning to Achieve Target Outcomes
T=2
104
155
035
CRZ
0.2
02-
0.2
0.5-
RCS
0.5
0.5
06
0.6
ACTNV
D25
RMNN
1DP
ACTN
RMn
RMo
VOP
ACTM
CR6
RKON
(a) 3
CRRP
0.-
D65
0.4
0.4
0.1
000
0.5
0.5
RCS
0.09
6.0
300
1.0
10
1.
1.
ACTDN
RMa
veB
ACTEN
C
RMbo
ACTIN
RMoo
Yog
ACTDC
CT
Rabos
(b) 7
Figure 10. Evaluation results (Part 2) on the tumor simulation dataset with varying prediction horizons (r = 2, 8) under different confounding bias levels (7 3, 4).
22


Reviews:
[{"review_id": "lc7wyu4ey7", "reviewer": "Reviewer_whTh", "summary": "This paper presents a new method for finding desirable intervention sequences for individual instances. \nFirst, the authors formulate the task of finding effective intervention sequences as an optimization problem that maximizes the likelihood of the target outcome after the intervention. \nThen, the authors propose a framework called VCIP, which uses variational inference to construct a surrogate function that approximates the likelihood of the target outcome. \nThe numerical experiments demonstrate that the proposed method could effectively find effective intervention sequences to achieve the desired outcomes.\n\n## update after rebuttal \n\nThank you for your response. I appreciate the authors' efforts to clarify the points I raised. Since my main concerns have been addressed, I maintain my evaluation.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thorough review and appreciation of our work. Below, we address your main concerns:\n\n**Explanation of GRP and RCS Metrics:**\n\nGRP focuses on how the model ranks a sequence $\\bar{a}\\_{t,\\tau}$ that can definitely achieve $Y\\_{target} = Y[\\bar{a}\\_{t,\\tau}]$. Ideally, GRP should be 1. This metric only requires knowledge of the potential outcome of sequence $\\bar{a}\\_{t,\\tau}$, which is available in observational data, making it applicable to real-world datasets (MIMIC-III dataset).\n\nRCS evaluates the model's performance across the entire candidate set, testing whether it can compare the likelihood of different sequences achieving $Y\\_{target}$. We use actual target distances for comparison (line 261). However, this requires calculating the true potential outcome for each candidate sequence for each individual, thus it can only be used with simulated datasets (tumor datasets).\n\nTherefore, an increase in GRP indicates the model can more accurately identify sequences that can achieve the target, while an increase in RCS indicates the model can better rank the entire candidate set by likelihood of achieving the target, with higher values showing better alignment with true outcomes.\n\n**Rationale for Setting $Y\\_{target} = Y[\\bar{a}\\_{t,\\tau}]$:**\n\nWe set $Y\\_{target} = Y[\\bar{a}\\_{t,\\tau}]$ to ensure that we can establish that sequence $\\bar{a}\\_{t,\\tau}$ can achieve outcome $Y\\_{target}$. In this scenario, an ideal model should rank $\\bar{a}\\_{t,\\tau}$ first among candidate sequences, giving a GRP of 1. Regarding whether $\\bar{a}\\_t$ is necessarily an optimal intervention sequence to obtain $Y\\_{target}$, indeed there may be multiple intervention sequences that can achieve $Y\\_{target}$, but this does not affect the principle that an ideal model should rank $\\bar{a}\\_{t,\\tau}$ first among candidate sequences (possibly tied with others).\n\n**Regarding Concerns about Theoretical Claims:**\n\nThank you for this valuable feedback on Theorem 4.1. We agree that the current claim could be stated more precisely. We appreciate your suggestion to modify the claim to explicitly state that the intervention sequence that minimizes Equation (6) is an $\\epsilon$-optimal solution to the problem that maximizes Y. \n\nWe will revise Theorem 4.1 in the updated manuscript to clarify this relationship and remove any ambiguity. The formal statement will be adjusted to more accurately reflect the theoretical guarantee provided by our approach.\n\n**Regarding Concerns about the Random Perturbations:**\n\nWe employ a hybrid approach to generate candidate sequences. Our framework creates **random sequences** (50%-80% of candidates) and **perturbed ground truth sequences** (20%-50% of candidates).\n\nThe perturbation strategy is treatment-mode specific:\n- **For discrete interventions**: We randomly flip bits in the ground truth sequence with probability 0.2.\n- **For continuous interventions**: We apply context-aware shifts where values are modified based on their magnitude (low values shifted up, high values shifted down, middle values shifted randomly).\n\n**Generalizability Considerations**. Our mixed-generation strategy ensures robust evaluation across different perturbation levels by testing against both arbitrary interventions and \"near-miss\" candidates that challenge the model's discrimination ability. This design ensures our evaluation is generalizable beyond specific perturbation patterns, as it comprehensively tests the model's ability to rank interventions across the similarity spectrum.\n\nFull implementation details are available in our codebase under `src/utils/helper/generate_perturbed_sequences`. We will add relevant details in the updated manuscript.", "claims_and_evidence": "Overall, the claims made in this paper are well-supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "For the proposed framework, the authors provide a clear motivation and rationale for the design of the VCIP framework. \n\nFor the numerical experiments, perhaps it may be just my misunderstanding, but I have a slight concern about the ranking-based evaluation in Section 5.1. \nIn Section 5.1, the authors set $Y_{\\\\text{target}} = Y[\\\\bar{a}\\_{t,\\\\tau}]$, and compared whether the model-based ranking of the intervention sequences is consistent with the ground-truth ranking. \nI think that \\\\(\\bar{a}\\_{t,\\tau}\\\\) is not necessarily an optimal intervention sequence to obtain \\\\(Y_{\\text{target}}\\\\), so I am not sure what the authors aim to validate through this comparison.", "theoretical_claims": "The claim of Theorem 4.1 is somewhat ambiguous. \nIt might be clearer to modify the claim of Theorem 4.1 to state that the intervention sequence that minimizes Equation (6) is an $\\epsilon_1+\\epsilon_2$-optimal solution to the problem that maximizes $\\mathcal O$.\nI  have checked the proof of Theorem 4.1.", "experimental_designs_or_analyses": "The experiments were generally conducted in a reasonable manner.\n\nHowever, I have a concern about the Ranking-based Evaluation. In this evaluation, random perturbations are added to the ground truth intervention sequences to generate $k$ new intervention sequences, but the details of the generation process are not clearly described.\nAs a result, it is unclear to what extent the reported results hold for different levels of perturbation, making the generalizability of this evaluation unclear.\n\nAdditionally, as mentioned in ''Methods and Evaluation Criteria'', there is a concern about whether Ranking-based Evaluation is appropriate for discussing the effectiveness of the proposed method.", "relation_to_broader_scientific_literature": "While existing studies have focused on deriving intervention sequences at the population level, this study proposes a method for optimizing intervention sequences at the individual level. The proposed approach enables learning models that effectively guide decision-making to achieve desirable outcomes. Providing such guidance is increasingly important, particularly in the context of Explainable AI, and is expected to have significant value in fields such as causal inference and algorithmic recourse.", "essential_references_not_discussed": "None in particular.", "other_strengths_and_weaknesses": "None in particular.", "other_comments_or_suggestions": "None in particular.", "questions_for_authors": "1. As mentioned in \"Methods and Evaluation Criteria\", could you clarify what you aim to demonstrate in Section 5.1? Additionally, why is it reasonable to consider a model desirable if GRP and RCS increase for intervention sequences that are not necessarily optimal or close to optimal for the target outcome?\n\n2. Why are the results of RCS not reported for the MIMIC-III dataset? Was the same evaluation performed for this dataset? If so, could you explain the results?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "03I4yxsIFX", "reviewer": "Reviewer_XzCZ", "summary": "The paper introduces an approach named \"variational counterfactual intervention planning (VCIP)\" to address the problem of optimal sequences of interventions selection towards a target outcome. The method is useful particularly in healthcare scenarios. Traditional counterfactual estimation methods suffer from compounding errors due to their inherent reliance on unobservable counterfactual outcomes. VCIP addresses this issue by reformulating the problem through a variational inference framework, directly modeling the conditional likelihood of achieving target outcomes, hence avoiding explicit prediction of counterfactuals. Experiments are conducted on both synthetic and real-world datasets, demonstrated superior performance compared to existing methods.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thoughtful review and for acknowledging our work. We will address your main concerns below.\n\n**Regarding Concerns about Claims and Evidence:**\n\nThe **Consistency** assumption is typically satisfied in clinical settings where treatments are well-defined and outcomes can be stably measured after treatment administration.\n\nTo test model performance under Positivity assumption violations, we conducted experiments on the Tumor dataset at $\\gamma$=4, setting the probability of receiving treatment to 0 when treatment values >0.5 (while simple, this effectively violates the Positivity assumption). Below are results under the Optimization experiment :\n\n| Method           | τ=2           | τ=4           | τ=6           |\n| ---------------- | ------------- | ------------- | ------------- |\n| RMSN (violated)  | $2.24\\pm0.90$ | $3.41\\pm1.26$ | $4.49\\pm1.17$ |\n| RMSN (satisfied) | $0.45\\pm0.10$ | $0.75\\pm0.16$ | $0.98\\pm0.22$ |\n| VCIP (violated)  | $1.32\\pm0.30$ | $2.10\\pm0.40$ | $2.81\\pm0.46$ |\n| VCIP (satisfied) | $0.42\\pm0.13$ | $0.60\\pm0.15$ | $0.75\\pm0.20$ |\n\nAs shown, violating the Positivity assumption degrades performance for both RMSN and VCIP, as models encounter interventions with no supporting observational data during optimization. However, VCIP's performance degrades notably less than RMSN's, likely because RMSN suffers more severely from compounding errors in such scenarios.\n\nTesting robustness to **Sequential Ignorability** violations requires specialized datasets with unobserved confounding, which we couldn't explore due to time constraints. Such violations introduce bias, and while approaches like Time Series Deconfounder exist, resolving this in our context remains an open question outside our current scope. We will gradually explore these directions in our subsequent research work.\n\nAdditionally, we acknowledge personalized healthcare applications face multifaceted challenges. Computationally, VCIP requires 2400s for running, and RMSN 4200s, demonstrating clinical feasibility while scaling remains challenging. Data sparsity issues can be addressed through multiple imputation, GAN augmentation, and transfer learning. We recognize the importance of ethical considerations including fairness, transparency, and privacy protection. The revised version will elaborate on these challenges, analyzing computational efficiency, methods for handling sparse data, and expanding ethical discussions to enhance practical applicability.\n\n**Regarding Concerns about Theoretical Claims:**\n\nWhen we **omit** $do(\\bar{\\mathbf{a}}\\_{t,\\tau})$, the variational distribution does not explicitly model the intervention process, but if $\\mathbf{Y}\\_{t+\\tau}$ is of high quality (in experimental settings, sufficient and accurate $\\mathbf{Y}\\_{t+\\tau}$ information can often be observed), the impact brought by the intervention is largely \"reflected\" in $\\mathbf{Y}\\_{t+\\tau}$. To illustrate this point, we compare by explicitly adding observed intervention information, namely using $q\\_{\\phi}(\\bar{\\mathbf{Z}}\\_{t,\\tau+1} \\mid \\bar{\\mathbf{H}}\\_t, \\mathbf{Y}\\_{t+\\tau},\\bar{\\mathbf{a}}\\_{t,\\tau})$\n\n1. **Optimization under $\\gamma=3$**  \n\n   |                                      | $\\tau=2$      | $\\tau=6$      | $\\tau=8$      | $\\tau=12$     |\n   | ------------------------------------ | ------------- | ------------- | ------------- | ------------- |\n   | omit $do(\\bar{\\mathbf{a}}\\_{t,\\tau})$ | $0.39\\pm0.26$ | $0.63\\pm0.33$ | $0.67\\pm0.33$ | $0.79\\pm0.33$ |\n   | with $\\bar{\\mathbf{a}}\\_{t,\\tau}$     | $0.37\\pm0.23$ | $0.60\\pm0.25$ | $0.64\\pm0.23$ | $0.78\\pm0.28$ |\n\n   As can be seen, compared to omitting $do(\\bar{\\mathbf{a}}\\_{t,\\tau})$, explicitly introducing $\\bar{\\mathbf{a}}\\_{t,\\tau}$ shows a minor improvement in model performance, but the difference is not significant.\n\n2. **Ranking under $\\gamma=4$**  \n\n   |                                      | GRP $\\tau=2$  | RCS $\\tau=2$  |\n   | ------------------------------------ | ------------- | ------------- |\n   | omit $do(\\bar{\\mathbf{a}}\\_{t,\\tau})$ | $0.94\\pm0.09$ | $0.77\\pm0.21$ |\n   | with $\\bar{\\mathbf{a}}\\_{t,\\tau}$     | $0.95\\pm0.09$ | $0.79\\pm0.21$ |\n\n   Here again, explicitly including $\\bar{\\mathbf{a}}\\_{t,\\tau}$ in the variational distribution shows a slight performance improvement, though omitting the intervention sequence still maintains performance very close to the optimal value.\n\nTherefore, the approximation bias from omitting $do(\\bar{\\mathbf{a}}\\_{t,\\tau})$ has minimal negative impact, primarily because $\\mathbf{Y}\\_{t+\\tau}$ carries effective intervention result information, allowing variational inference to partially capture intervention effects on latent variable $\\bar{\\mathbf{Z}}\\_{t,\\tau+1}$ while learning $\\mathbf{Y}\\_{t+\\tau}$.\n\nAdditionally, the positivity of $\\epsilon\\_1$ and $\\epsilon\\_2$ can be derived from the third inequality in Eq. 22 and Eq. 23, or more rigorously, these values are non-negative, which we will correct in the revised version.", "claims_and_evidence": "The following claims need further evidence or clarity to fully support the assertions in a convincing way:\n\n1. Robustness under violations of standard causal assumptions\nThe paper implicitly assumes consistency, positivity, and sequential ignorability without extensive empirical or theoretical exploration of sensitivity to these assumptions. Real-world scenarios often violate these assumptions due to unobserved confounding, missing data, measurement error, or non-random treatment assignment patterns. The paper should explicitly evaluate or at least discuss VCIP’s performance under potential assumption violations. (Perhaps, sensitivity analysis, or ablation studies?)\n2. Practical applicability in personalized healthcare scenarios\nAuthors claim significant potential impact in personalized healthcare but provide limited discussion of realistic challenges such as computational complexity, data sparsity, missingness, and ethical issues.", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "I carefully checked the statements of theoretical claims provided in the paper but did not go into the detailed verification of the proofs.\n\nIn line 145, the authors explicitly mention that they \"omit the intractable intervention sequence $do(\\bar{a}_{t,\\tau})$ as its effects are partially captured in $Y_{t+\\tau}$. This step, while practically understandable, might slightly weaken the theoretical rigor. A suspicious aspect here is whether this approximation significantly affects the theoretical guarantees or introduces unacknowledged biases. This step is not fully justified, which leaves room for questioning the precision or generality of theoretical guarantees.\n\nIn Theorem 4.1, why can we safely assume $\\epsilon_1$ and $\\epsilon_2$ are positive? Any justification?", "experimental_designs_or_analyses": "The experimental design is methodologically robust, clear, and aligned with standard practices. However, I identified some issues:\n\nThe experimental design implicitly relies on approximations for intractable intervention sequences. While reasonable, this decision is not empirically validated, raising minor suspicion about possible biases or implications of this simplification.\n\nAlso, the paper lacks of robustness analysis against violations of standard assumptions.", "relation_to_broader_scientific_literature": "The paper’s key contributions relate closely to several well-established strands within the broader scientific literature on causal inference and variational inference.  In contrast, this paper introduces a novel inverse formulation—the \"counterfactual target achievement\" problem—shifting the goal from prediction to actively selecting intervention sequences that drive outcomes toward specified targets.", "essential_references_not_discussed": "No.", "other_strengths_and_weaknesses": "None", "other_comments_or_suggestions": "None", "questions_for_authors": "None", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "DVkjyAZbPf", "reviewer": "Reviewer_g91d", "summary": "This paper addresses the problem of time varying treatment effect, aiming at finding the sequence of treatments that optimize a target outcome, instead of the typical problem of predicting potential outcomes. It uses the g-formula and a variational approach to estimate the conditional likelihood of achieving target outcomes. This is then used to find an optimal sequence of treatments.\n\n## Update after rebuttal\nAfter the substantive discussion with the authors, who tried to clarify my main points of contention, and in light of the interesting and important topic, with promising experimental results, I am willing to change my score from 2: Weak Reject, to 3: Weak Accept.\n\nI would like to thank the authors for their time and effort in this discussion.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review. Below we address your main concerns:\n\n**Regarding Concerns about Evaluation Criteria:**\n\nOur primary contribution is a novel \"counterfactual target achievement\" problem formulation, which differs fundamentally from counterfactual estimation addressed by models like RMSN. This difference in objectives explains why we didn't initially compare models on counterfactual estimation tasks.\n\nFollowing your valuable suggestion, we conducted additional multi-step estimation experiments at γ=2,4 to better illustrate how VCIP compares with other models on counterfactual prediction:\n\n|           | τ=4       | τ=6     | τ=8     | τ=10    |\n| --------- | --------- | --------- | --------- | --------- |\n| RMSN  γ=2 | 0.90±0.23 | 1.07±0.32 | 1.18±0.36 | 1.37±0.49 |\n| VCIP  γ=2    | 0.90±0.70 | 0.84±0.67 | 0.84±0.62 | 0.79±0.65 |\n| RMSN γ=4  | 1.34±0.21 | 1.61±0.31 | 1.85±0.32 | 2.09±0.43 |\n| VCIP  γ=2   | 1.83±0.61 | 1.75±0.66 | 1.74±0.59 | 1.79±0.64 |\n\nAs Figure 1 shows, using counterfactual predictions is suboptimal for counterfactual target achievement, as treating $Y\\_{target}$ as an intermediate variable causes compounding errors. Our method avoids this by directly incorporating $Y\\_{target}$ into the likelihood. The experiments reveal interesting insights:\n\n- Even when VCIP's counterfactual estimations sometimes underperform RMSN, it still shows advantages in the achievement problem, validating our framework's effectiveness.\n\n- Accurate counterfactual estimation helps the achievement problem: at τ=2, where RMSN outperforms VCIP in estimation, the performance gap in our problem is small (approximately 0.02), while at τ=12, where VCIP excels in estimation, the achievement performance gap widens (approximately 0.4).\n\n**Regarding concerns about Theoretical Claims**:\n\nC1)\n\nWhile $\\mathrm{ELBO}\\_1$ captures the true causal mechanism through consideration of interventions (do-operator), in practical training we typically only have observational data for maximum likelihood estimation. As long as the model structure has sufficient expressivity and the observational data reasonably approximates the causal process, $\\mathrm{ELBO}\\_1$ can be indirectly approximated during learning, making $\\epsilon\\_1$ relatively small. This is supported by our ablation studies (Table 3): even without adjustment, VCIP performs comparably to or better than RMSN, indicating that maximizing eq. 5 (where $\\mathrm{ELBO}\\_2$ can be directly optimized based on observational data) effectively drives improvements in eq. 2 ($\\mathrm{ELBO}\\_1$), thus approximating $\\mathcal{O}$ and preventing $\\epsilon\\_1$ from becoming \"arbitrarily large.\"\n\nFor intuitive interpretation of eq. 6, term A maximizes observational likelihood, while terms B and C serve as adjustment terms:\n\n$$\n\\text{(A)}\\ \\mathrm{ELBO}\\_2\\ \\ +\\ \\text{(B)}-\\ \\sum\\_{s=t}^{t+\\tau-1}\\mathbb E\\_{q\\_\\phi}[\\log p\\_\\theta(\\mathbf{a}\\_s\\mid \\mathbf{Z}\\_s)]\\ +\\ \\text{(C)}\\ \\log p\\_\\theta(\\bar{\\mathbf{a}}\\_{t,\\tau}\\mid \\bar{\\mathbf{H}}\\_t),\n$$\n\nIncreasing term B encourages the model to learn states that cannot accurately predict interventions, intuitively allowing interventions to \"break free\" from observational distribution relationships, thus mitigating confounding bias. Term C compensates for inherently reasonable action sequences by providing a bonus, avoiding excessive penalties in term B.\n\nThis demonstrates our core intuition of jointly optimizing these three terms to elegantly optimize $\\mathrm{ELBO}\\_1$ through $\\mathrm{ELBO}\\_2$ combined with appropriate penalties and bonuses for action sequences.\n\nTo verify the effects of B and C separately, we conducted more detailed ablation experiments.\n\n|        | GRP τ=2 | RCS τ=4 | Optimization τ=6 | τ=8  | τ=10 | τ=12 |\n| ------ | ------- | ------- | ---------------- | ---- | ---- | ---- |\n| Ours   | 0.94    | 0.87    | 0.75             | 0.92 | 0.97 | 1.08 |\n| w/o B  | 0.92    | 0.83    | 0.77             | 0.95 | 1.00 | 1.18 |\n| w/o C  | 0.91    | 0.61    | 0.74             | 0.91 | 0.92 | 1.01 |\n| w/o BC | 0.76    | 0.60    | 0.91             | 1.14 | 1.27 | 1.48 |\n\nFrom the results, we can see that both B and C improve performance on Ranking and Optimization tasks. However, considering both tasks comprehensively, incorporating both B and C simultaneously is the most appropriate approach.\n\nC2)\n\nPlease refer to our response to reviewer XzCZ's \"Regarding Concerns about Claims And Evidence\" section.\n\nC3)\n\nWe will revise our description of distribution $q_\\phi$ in the updated manuscript.\n\n**Regarding Concerns about Categorical Interventions**\n\nWe didn't design additional optimization algorithms for categorical treatments. Since categorical treatments are typically enumerable, one can use the ranking approach and enumerate to find the optimal intervention sequence. For more complex categories, techniques like Gumbel-Softmax could be used to design optimization processes similar to Algorithm 1.", "claims_and_evidence": "Proofs are offered for the main claims made in the paper. Some of them are, in my opinion, problematic. I deepen in that in the section of Theoretical Claims.", "methods_and_evaluation_criteria": "The idea of the rankings is consistent with the main objective of the paper. However, more direct comparisons with other benchmarks in traditional metrics like distance between counterfactual outcomes and ground truth counterfactuals would probably strengthen the paper.\n\nApart from this, it is a bit difficult to understand the content of the tables from the captions.", "theoretical_claims": "I Checked the theoretical claims and proofs, and I have several concerns:\n\nc1) I am not convinced by theorem 4.1, which is fundamental in this work. While the relations expressed in the equations are correct, I am not sure that optimizing the expression of eq. 6 amounts to optimizing eq. 2, which is the main claim. The problem that I see is that, if ELBO1 is not maximized, then \\epsilon_{1} can be arbitrarily large. Then, the error between the interventional loglikelihood (eq. 2) and the expression in eq. 6 will also be arbitrarily large. It would be very good if the authors could address this problem, as this is a major concern.\n\nc2) In 4.1, it is mentioned that a variational distribution is introduced to approximate the true posterior, but the term do(a) is omitted for practical considerations, as its effects are partially captured by Y. However, in light of the kind of problem that the paper addresses this seems that it is an important approximation that is not very discussed. Maybe the authors could try to better justify this approximation. \n\nc3) In the inference model, the claim that Z_{s} is obtained from its descendants does not seem very convincing. As mentioned in the last lines of this section, there is variable Z’_{s} (one from each step) obtained from the Z’_{s-1} and a_{s-1}. Then, those are the variables, for s={t+1,….t+\\thau}, used to obtain Z_{s}. Then, I think it would be better to say that Z_{s} is obtained from the descendant treatments after time s than from the descendants of latent factors.", "experimental_designs_or_analyses": "Overall, experiments look sound, but I didn’t check the code.\n\nIt would have been interesting to see how the proposed model compares to other models when estimating counterfactual outcomes, with distance metrics to ground truth counterfactuals. The estimated means of the output could be used to do that, and it would probably give a better idea of how trustable the model is.", "relation_to_broader_scientific_literature": "In general, previous important works like Causal Transformer, Counterfactual Recurrent Networks or ACTIN are properly discussed.", "essential_references_not_discussed": "All important papers are cited.", "other_strengths_and_weaknesses": "The paper addresses an important problem, and the idea seems interesting. However, there are some important issues that need to be clarified. Also, the paper is difficult to understand, and the authors could have done more efforts to improve clarity and give more explanations, for example on theorem 4.1 and the appendix proofs, which are fundamental for the paper despite being in the appendix. Also, it would be interesting if some intuitive interpretation of the terms in eq.6 could be given.", "other_comments_or_suggestions": "-", "questions_for_authors": "The treatment sequence is optimized with gradient descent. Does the method offer a solution for categorical treatments? Is it explained in the paper?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TEtBeaKKMx", "reviewer": "Reviewer_ixW5", "summary": "The  paper introduces Variational Counterfactual Intervention Planning (VCIP), a framework for determining optimal intervention sequences in personalized healthcare and other temporal decision-making systems.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thorough review and recognition of this work. Below we address your main concerns:\n\n**Regarding Concerns about Data:**\n\n\"Dependence on observational data quality\" is a common challenge in causal inference, especially in medical data analysis where EHR data usage can introduce quality issues. The VCIP framework mitigates error propagation by using variational inference to introduce latent variables that capture system state evolution, and directly modeling target achievement probability rather than relying on explicit predictions. We utilize the g-formula to establish connections between intervention and observational distributions, enabling reliable training on observational data. Nevertheless, we acknowledge the importance of data preprocessing and will explore data cleaning and specific medical data processing strategies in future work to enhance model robustness when facing imperfect data.\n\n**Regarding Concerns about Prediction Intervals:**\n\nThank you for this valuable comment. In practice, many clinical decisions indeed aim to maintain patient metrics within \"safe ranges.\" We acknowledge that our current model's focus on single target point optimization is a limitation, but a straightforward approach would be taking a value from the \"range of values\" as the target, such as the midpoint of the interval. For more precise characterization, extending VCIP to model the probability of falling within a target range would be necessary. While this is beyond the current paper's scope, it represents an interesting research direction that could effectively enhance our method's clinical applicability.\n\n**Regarding Concerns about High-Dimensionality:**\n\nThis paper validates using the widely-used real-world MIMIC-III dataset, which includes 25 time-varying features and a 2-dimensional intervention space, representing a relatively complex dataset in current medical causal inference research. At this scale, VCIP's complete runtime (including training and intervention ranking) is approximately 2700 seconds, which falls within an acceptable range. Validating VCIP's scalability in higher-dimensional intervention spaces requires newer, more complex datasets. We acknowledge this has important application value, but building such datasets in the medical field faces significant challenges. Future work will explore optimizing algorithm efficiency and testing performance in more complex medical decision scenarios to further validate the method's practicality.\n\n**Regarding Concerns about Unseen Interventions:**\n\nIn fact, to ensure identifiability of causal effects, this paper makes standard assumptions, including Assumption 2 (Sequential Overlap), which theoretically guarantees that any intervention has a possibility of being observed (note that \"unseen\" in this paper refers to using different intervention strategies in testing than in training, which doesn't conflict with Assumption 2). This assumption enables good generalization when observational data is sufficient. When this assumption is violated, meaning interventions not present in training data are encountered, we conducted experiments on the Tumor dataset. Specifically, with $\\gamma=4$, we set the probability of interventions greater than 0.5 to 0 (this may be simple but effectively violates Assumption 2). Below are the Optimization experiment results:\n\n| Method           | τ=2           | τ=4           | τ=6           |\n| ---------------- | ------------- | ------------- | ------------- |\n| RMSN (violated)  | $2.24\\pm0.90$ | $3.41\\pm1.26$ | $4.49\\pm1.17$ |\n| RMSN (satisfied) | $0.45\\pm0.10$ | $0.75\\pm0.16$ | $0.98\\pm0.22$ |\n| VCIP (violated)  | $1.32\\pm0.30$ | $2.10\\pm0.40$ | $2.81\\pm0.46$ |\n| VCIP (satisfied) | $0.42\\pm0.13$ | $0.60\\pm0.15$ | $0.75\\pm0.20$ |\n\nWhen Assumption 2 is violated, both models face unseen intervention sequences. As shown, both VCIP and RMSN performance deteriorates, but VCIP exhibits a smaller performance drop, indicating superior generalization capabilities.", "claims_and_evidence": "yes", "methods_and_evaluation_criteria": "they do", "theoretical_claims": "checked but not proofs in appendix", "experimental_designs_or_analyses": "seem reasonable", "relation_to_broader_scientific_literature": "builds upon literature for causal reasoning in personalised healthcare", "essential_references_not_discussed": "no", "other_strengths_and_weaknesses": "Pros : \n\n- Addresses a critical problem in personalized healthcare—predicting the best intervention sequences rather than merely estimating outcomes.\n- Reduces compounding errors common in standard counterfactual estimations by directly modeling target achievement probability.\n- Outperforms baseline models in both simulated and real-world datasets, particularly in ranking interventions.\n- Uses principled causal inference via the g-formula and variational inference, ensuring theoretically sound predictions.\n- Avoids over-reliance on counterfactual predictions, which are inherently unobservable and prone to error accumulation.\n\nCons:\n- Dependence on quality of observational data—errors in training data can propagate through the model.\n- In medicine usually there is no single target we want to optimise for, but rather a range of values, as such the existence or validity of the target outcome is in question \n- Handling high-dimensional intervention spaces or multiple simultaneous treatments may be computationally expensive.", "other_comments_or_suggestions": "see below", "questions_for_authors": "How does VCIP handle unseen interventions? Since it relies on observational data, does it generalize well when encountering interventions not present in training data?\n\nHow scalable is VCIP to high-dimensional, multi-treatment scenarios? The paper focuses on limited intervention sequences—can this approach scale efficiently to complex, real-world medical decisions?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: Short summary: This paper formalizes and studies the problem of “counterfactual target achievement” – in which goal is to find the optimal intervention sequences for a given target in a single instance, while mitigating errors from unobserved counterfactuals. 

All reviews recommended a “weak accept” for the paper. Some of the key merits of the paper as identified by reviews include, importance of the problem, novelty of the formulation and contribution, a generally robust experimental analysis (although some suggestions for improvements are made), The reviews also identify some limitations. These include concerns raised around theoretical claims (theorem 4.1), paper being somewhat hard to read and room for improved exposition (improving captions to Tables, etc.) 

Furthermore, the rebuttal also provides extensive arguments and explanations to allay many of the reviewer concerns. Hence my reasoning is to accept the paper (medium priority)
Decision: Accept (poster)