Structured Learning of Compositional Sequential Interventions
Authors: Jialin Yu
Andreas Koukorinis
Nicolò Colombo
Yuchen Zhu
Ricardo Silva
Venue: NeurIPS 2024 poster (2024)
Keywords: Causality
sequential data
Primary area: causal_inference
TL;DR: We describe a way of identifying and learning the effect of combinations of interventions in a sequential setup, in the regime of sparse data where limited combinations are jointly observed.

Abstract: We consider sequential treatment regimes where each unit is exposed to combinations of interventions over time. When interventions are described by qualitative labels, such as "close schools for a month due to a pandemic" or "promote this podcast to this user during this week", it is unclear which appropriate structural assumptions allow us to generalize behavioral predictions to previously unseen combinations of interventions. Standard black-box approaches mapping sequences of categorical variables to outputs are applicable, but they rely on poorly understood assumptions on how reliable generalization can be obtained, and may underperform under sparse sequences, temporal variability, and large action spaces. To approach that,  we pose an explicit model for composition, that is, how the effect of sequential interventions can be isolated into modules, clarifying which data conditions allow for the identification of their combined effect at different units and time steps. We show the identification properties of our compositional model, inspired by advances in causal matrix factorization methods. Our focus is on predictive models for novel compositions of interventions instead of matrix completion tasks and causal effect estimation. We compare our approach to flexible but generic black-box models to illustrate how structure aids prediction in sparse data conditions.

Structured Learning of Compositional Sequential Interventions
Jialin Yu
Andreas Koukorinis
Nicolò Colombo²
Yuchen Zhu'
Ricardo Silva
University College London
Royal Holloway, University of London
Abstract
We consider sequential treatment regimes where each unit exposed to combina- tions of interventions over time. When interventions are described by qualitative labels, such as "close schools for month due to a pandemic" or "promote this podcast to this user during this week", it is unclear which appropriate structural assumptions allow us to generalize behavioral predictions to previously unseen combinations of interventions. Standard black-box approaches mapping sequences of categorical variables to outputs are applicable, but they rely on poorly understood assumptions how reliable generalization can be obtained, and may underperform under sparse sequences, temporal variability, and large action spaces. To approach that, we pose explicit model for composition, that is, how the effect of sequential interventions can be isolated into modules, clarifying which data conditions allow for the identification of their combined effect different units and time steps. We show the identification properties of our compositional model, inspired by advances in causal matrix factorization methods. Our focus on predictive models for novel compositions of interventions instead of matrix completion tasks and causal effect estimation. We compare our approach to flexible but generic black-bor models illustrate how structure aids prediction in sparse data conditions.
1 Contribution
Many causal inference questions involve treat- ment sequence that varies over time. In the (discrete) time-varying scenario, an unit n is exposed sequence D1, D D3 D of tions, or treatments. potentially controllable by an external agent [22] Ch. 19] intervening on them. Each D' can be interpreted as a cause of and its modeling amenable to the tools of causal infer- ence [e.g..B35 36
D
D2
D
D
X
the subsequent behavior X X x
49). For instance, each D can take values in the space of par- ticular drug dosages which can be given to an in-patient, or combinations of items to be pro- moted to a user of recommender system. Here, the units, and X denotes
Figure 1: Within unit n, actions D interact with (latent) random effect parameters 8, to produce behavior X represented as a dense graphical model with square vertices denoting interventions do(d)::) [33] [6]. Further assumptions will be required for the identifiability of the impact of in- terventions and their combination, including how temporal impact takes shape and the number of independent units observation.
patients and users a measure of their state or behavior under time t. Dependencies can result in dense causal graph, as in Fig. T Further background on this classical problem provided in Appendix A
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Scope. We introduce an approach to predict the effects of actions within combinatorial and sparse space of categorical intervention sequences. That is, for an unit n and given past behavioral sequence x x := [x n.... ), manipulable variables D1:t-1 := [D
........ D and time-fixed covari- ates Zn. we predict future behavior X X : X under the (partially hypothetical) interven- R2
tional regime do(d) The latter denotes controlling the action variables by external intervention
33 For mean-squared error losses, prediction means learning E[ _./ii-1,, 2n. do d):++))), A. We consider the case where sequential ignorabilityy holds, by randomization or adjustment
33
22
in order to focus on predictive guarantees for previously unseen action sequences.
Challenge. Much of the sequential intervention literature in AI and statistics considers generic black-box or parametric model to map sequences D == [D D ..... D to behavior X and beyond, sometimes exploiting strong Markovian assumptions or short sequences and small action spaces [1] 43 An unstructured black-box for sequences, such as recurrent neural networks and their variants [e.g. 13 45), can learn to map a sequence to a prediction. This is less suitable when actions are categorical and sparse, in the sense that most entries correspond to some baseline category of "no changes". That is, there is no obvious choice of smoothing or interpolation criteria to generalize from seen to unseen strings DA:T, with the practical assumption being that we observe enough variability to reliably perform this interpolation. There is an increased awareness of problems posed by large discrete spaces among pre- pre-treatment covariates [47 48), which is of different nature as losses can be minimized with respect to an existing natural distribution, as opposed to the interventional problem where the action variables
allowed to be chosen and so sparsity conditions can easily result in very atypical test cases of interest. Several papers address large intervention spaces (mostly for non-sequential problems) without particular concerns about non-smooth, combinatorial identifiability A few more recent papers address explicit assumptions of identifiability directly in combinatorial categorical spaces by sparse regression 2 or energy functions (8) but without a longitudinal component. More classically, tools like the do-caleulus (3) and their extensions [5] can be used with directed acyclic graphs (DAGs) infer combinatorial effects, although they often restrict each manipulated variable to directly affect only a small set of variables (c.g. (5]).
25
32 38
39
With strong assumptions and computational cost, this may include latent variables [50). Alternative representation learning ideas for carefully extrapolating to unseen interventions are presented by (37 in the non- -equeeuuet,,,,, non-combinatorial setting.
In this paper, we consider the conservative problem of identification guarantees for the effect of sequential combinations of interventions which may never have co-occurred jointly in training data, formalizing assumptions that allow for the transfer and recombination of information learned across units. These are not of the same nature as identifiability of causal effects from observational data (33] or causal discovery [41], but of causal extrapolation: even in the randomized (or unconfounded) regime, is not a given that the distribution under a particular regime do (o(d);) can be reliably inferred from limited past combinations of action sequences without an explicit structure on the causal generative model. We will consider the case where actions imply non-stationary transitions: once unit n is exposed to D' = d",, then X will have transition dynamics that can be affected not only by the choice of d', but also the difference A, resulting in dense dependency structures such as the one in Figure I As multiple treatments are applied to the same unit n. there will be composition of effects that will depend not only on the labels of the actions, but also on the order by which they take place. If a substantive number of choices of treatments at time i is available, a large dataset may have no two sequences DI:T and D1:T taking the same values, even when considering only the unique treatment values applied to and n', regardless of ordering and time stamps.
Problem statement. Let each individual unit n be described by time-invariant features Zn and a pair of time-series (X D!:T). ). Here, /41t22 denotes a slice of -ime time-series of variables V from 4 to 12, inclusive. That observations for all units stop at time T is not requirement, but just a way of simplifying presentation. Each D`, is categorical action variable which, when under external control, represents a potential intervention, encoded by values in N. We assume data is available under sequentially ignorable selection of actions B3 22]. either by randomization or consequences of an adjustment. D describes a (possibly hypothetical) intervention which causes changes to the distribution of Xt::0, The special value D', = 0 denotes "idle" "default" treatment level that can be interpreted as "no intervention". We are given a training set of units
We include here models which explicitly parameterize causal contrasts in longitudinal interventions, such as the blip models of (36).
2

---

X The goal is to predict future sequence of behavior X for (training or test) unit n' under hypothetical intervention sequence do(d),++)) This process takes place under the stable unit treatment value assumption (SUTVA) [24], meaning that action variables Dj only affect unit n.
Each intervention on D* is allowed to have an impact on the whole future series t, f + 1, 1, ...... It can correspond to instantaneous shock ("give five dollars of credit time (²) or an action that takes place over an extended period of time ("promote this podcast from time / to time L 1+ +-):: we just assume the meaning to be directly encoded within arbitrary category labels 0, 2, 2...... In a real-world scenario where a same action can be applied more than once, different symbols in N should be used for each instance (e.g., 1 arbitrarily standing for "give five dollars", and 2 for "give five dollars second time" and so on).
In Section 2 we provide account of what we mean by compositionality of interventions, with explicit assumptions for identifiability. In Section 3] we describe an algorithm for likelihood-based learning, along with approaches for predictive uncertainty quantification. Further related work is covered in Section before we perform experiments in Section 5 highlighting the shortcomings of black-box alternatives.
2 A Structural Approach for Intervention Composition
For simplicity, we will assume scalar behavioral measurements X as the multivariate case readily follows. Consider the following conditional mean model for the regime do(d)."),
B[X² xn 1=--1, 2n. do(d,,"") = (o')) (S @ (i
Σ Ou
(I)
where @ is the elementwise product, with dh Bn and via defined below. We postpone explaining the motivation for this structure to Section 2. First, we describe this structure in more detail.
Using p(-) to denote generic probability mass/density functions for random variables given by context, assume p(af,, | do(d)ht )) p(a do(d);")) for any t' (future interventions do not affect the past), and p do(d).") )) = p(zn) for all t. We do not explicitly condition on D.dtt in most of our notation, adopting a convention where intervention indices do(d):²) are always in agreement with the (implicit) corresponding observed D dia We define & == (c(x11111 2n as evaluations of basis functions or(-). For now, we will assume such (b(-) functions to be known. Each unit has individual- level coefficients 3n. Causal impact, as attributed to do(d):)) is given by
U == II (t
(2)
where function (d, f', captures a trajectory in time for the effect of Dj = d, for l' < t. Compositionality will be described in the sequel in two ways: first, as way of combining the effects of actions in a type of recursive functional composition (compositional recursion); second, by interpreting this function composition at time t in terms of an operator warping baseline unit-level feature vector 3n with an interventional-level embedding vector (ddt, t, t) (compositional warping).
Effect families. We consider two variants for v--). The first is a time-boue bounde variant which, for d > 0, defines edddd f', .t) == Wat,e-e', where Wat,e('' is a free parameter for 0 < s-- -t' ka- otherwise Wait-e' -a-- The intuition is that the effect of intervention level d has shape constraints but cannot affect the past and it must settle to a constant after a chosen time window hyperparameter ka > O. The default level d
has effect and no free parameters, with (0, f'', t) = 1. The second variant is motivated by real-world phenomena where causal impacts We consider time-unbounded
diminish their influence over time and result in new equilibrium 9. shape with exponential decay, parameterized as
(t) == o(wia)) t- x W2dl + W3dl.
(3)
where (-) is the sigmoid function, so that 0 < a(w1)) < 1. As t grows, aWwa)) goes to 0. Wya the stationary contribution of an intervention at level d. In particular, for d = 0, we define wiot w200 = 0 and w3or = 1, while for d > O we have that wiat W2at and wzat are free parameters to be learned. Let us interpret the previous ha hyperparameter as the dimensionality of intervention level d. As the model of Eq. ( is given by vectors Wid, Wad, Wad, we here define ka- = 3.
3

---

2.1 Representation Power: Compositional Recursion and and Compositional Warping
Eq. T is reminiscent of tensor factorization and its uses in causal modeling 7 where the primary motivation comes from the imputation of missing potential outcomes [24] taking place on a pre-determined period in the past. From a predictive perspective, can be motivated from known results in functional analysis (such as Proposition 1 of (25]) that allows us to represent a function fa) by first partitioning its input into (Fa; Tb) and controlling the approximation error given by an inner product of vector- vector-valued functions,
(F)) R (FW)))(()))))
which we call the factorization trick. In practice, we choose the dimensionality of such vectors (-) fo(-)) by a data-driven approach. Although those results apply typically to smooth functions with continuous inputs (such as the Taylor series approximation, which relies on vectors of monomials with coefficients constructed from derivatives), discrete inputs such as D1:T" at a fixed dimensionality T, can be separated from the other inputs
['n, d'.T I.TT .2) = fJ (dd x I ,2).
(5)
Here, similarly to classical ANOVA models, fa(+) is one-hot encoding vector for the entire (exponentially sized) space of possible d!-7 trajectories. That is, (-) is binary with exactly one non-zero entry, corresponding to which of the possible trajectories d! took place. The codomain of (tns(t:;))) must be a set of vectors of corresponding (exponential in 7) length.
In what follows, we will start from the premise of further partitioning ((n,d22 ,,1.1" Zn, do (d:2) }] by applying the factorization trick of Eq. (4) repeatedly, justifying the generality of the right- side (RHS) of Eq. (T). We will accomplish this by first introducing more structure into the respective fa(d':T) based on recursive composition of functions. We will then show in Proposition|] hhow the RHS of Eq. ① follows. We also provide a complementary theoretical result in Proposition 2 which, by assuming bounded T, describes the choice of dimensionality in the application of Eq. 4 as controlling a formal approximation error for general class of functions.
x11 E[X
factorization trick with compositional recursion. As we do not want to restrict ourselves to fixed T (and N), not to mention an exponential cost of representation, we will resort to embeddings of sequences and identities. In particular, we design fa(d to have a finite and tractable output dimensionality ry for all t. Inspired by sequential hidden representation models such as recurrent neural networks, for each entry / E 1,2, .y in the output of fa(-), we define fa lt) via the following recursive definition:
Main results:
fail (dd tt,t hilgiddiit-
- 1, d' tttt,, 0 octst,
(6)
hu(v,d, d, = hy, (u) hig (d)
The base case for L' is gi (dd=0 0,1) 1 for all I, with d2:0 undefined. Moreover, we define he(u, 0, t) = for any (v, t) that is, intervention level D does nothing.
In order to define hu d, f',') f) for d > O, we will apply the same trick of splitting the input into two subsets, and defining vector representation on each. The representational choice hu(v, d, 1', () == hi, (T)
hi. (d, t', t) requires choice of dimensionality Th, for the inner vectors hi, hıg- Recursively applying (6) in this product form can still be done relatively efficiently, but choosing Th; = 1 for all / simplifies matters. Finally, we apply the factorization trick one last time to define fnee(n, 1:T within Eq. 5 to get to factorization form analogous to Eq. (I). More formally:
Proposition 1 Let (i) f(n, dI.T xx.Txx..)2)) :) == :====fl (d (n, ,2), where function sequences
fT( and frzz (n, x: I.T 2) are defined for all E N+ and have codomain R (ii) fa ( be given as in Eq. (6) with Tby 1: (iii) the 1-th entry of (n, ,2) be given by Inau(1, (..T
== un(n)T mn( 2). where both u((:) and u((.-)) have codomain RR';; (iv) vi(xiit, ;)) be defined analogously to fat(d::t) as in Eq. T carrying out the fixed-size a as an extra argument. Then there exists some integer and three functions a, b and e with codomain
R' so that f(n, x ,2) Ei ar(a!:" ,2) x bi(n) x af (T,"),)) Moreover ai ,T) =
II mild', t, I) for some function my : 143 + R.

---

Proofs of this and next results are given in Appendix B The above shows that if we start with an assumption that the ground truth function takes some desired form, then there is some finite r for which our combination of basis functions can parameterize the ground truth function for arbitrary T. analogous to Eq. ( where (a, b, c) plays the role of (Ø, 13, v²). On the other hand, we show a companion result below: that if we hold T fixed, there is some fixed value for r, which depends polynomially on T, such that any measurable ground truth function can be approximated by our combination with some appropriate choice of basis functions.
Proposition 2 Given fixed value of TI, suppose we have measurable function f -Ui Xlit-1 x -'_ D' x Z + Rd, S T, where X1it-1 and Z are compact, D is finite. Letx E U2= xli-11 and let d E U!==1 D'. For d = =1,... ,dmaz lel ma
-i=11 D' - {1.... .T} with ma(d) = t if s.t. d, = d else O. Then, for any > 0, there exist r = O(Tdmaot1) and measurable
functions di : x x B + Rd and thi : :D' - R. and real numbers 3r, such that
dman
sup "PPxEU"_ deU deU?_ =2 f(x, Σ (xxx 31
II d=1
the (ma (1)) Sc. (7)
The theoretical results above are reminiscent of the line of work matrix factorization (e.g.. Appendix A of [4). Eq. ( and Eq. 2 are motivated by the success of approaches for representation learning of sequences, which assume that we can compile all the necessary information from the past using a current finite representation at any time step t. However, Propositions |
Notes.
and 2 further imply that, for f(n,, di:T .T..T 2) == E[X," x -.ITT1 2n, do(d):)) ). the representation of functions in a multilinear combination of adaptive factors (Eq. (T)) is motivated by more fundamental results in functional analysis. Our assumption that value of (other than the baseline 0) can be used more than once has mostly an empirical motivation, as this would imply unrealistic assumptions about the same intervention having the same effect when applied multiple times. Ultimately, we assume that the factors hNg (d, L', 1) leading to Eq. (2) are time-translation invariant, being a function of d and I - l' only. While it is possible generalize beyond translation invariant models, done in the synthetic controls literature [ that inspired causal matrix factorization ideas 7 3 4]. this would complicate matters further by requiring fourth factor into the summation term of Eq. (I) encoding absolute time. Models for time-series forecasting with parameter drifts can be tapped into and integrated with our model family, but we leave these out as future work.
compositional warping. At the start of the process, where all units are assumed to be at their "natural" state (assuming D 0), can interpret 7h1 = Bn a finite basis representation, with respect to 2/. of B[X; n Zm, do(d", 0)] ($)) Tm/. Each T/u can be seen as a latent feature of unit n. Effect vector vit detines warping of n/ that describes how treatments affect the behaviour of an individual in feature space, with the resulting n/h := of- @ vh. The modifier vi,, can be interpreted as reverting, suppressing or promoting particular latent features that are assumed to encode all information necessary to reconstruct the conditional expectation of X from the chosen basis (notice that intervention level 0 implies 1//- ./2 since we define (0, L' 1). This presents a sequential warping view of compositionality of interventions. Our framing also conveniently reduces the problem of identifiability of conditional effects to the problem of identifying baseline vectors 3, and warp function embeddings (d, i', t), as we will see next.
Interpretation:
2.2 Identification and Data Assumptions
Assume for now that functions é are known. As commonly done in the matrix factorization literature, assume also that we have access to some moments of the distribution. In particular, for given N data points and T time points, we have access to E[X; x 1it-1 En. do do(d).)) )]] We will impose conditions on the realizations of X1:T and values chosen for di.1" well as values of N and T as a function of r, in order to identify each Sn, and the parameters of (/d L', ()) for all intervention levels d of interest. The theoretical results presented, which describe how parameters are identifiable from population expectations and known function spaces determined by prescribed basis 6 of known and finite dimensionality r, will then provide the foundation for a practical learning algorithm in the sequel.
Assumption For a given individual n, we assume there is an initial period To., 2 r with "no interventions" i.e.. D = for all 1 < t < To..
5

---

To. - is given by Assumption 2 Let An be (To.. - 1) x r matrix, each row t 1.2 6x (1rh" .2,) realizations under regime do(d.."). Let Ag be full rank with left pseudo-inversy A
The first assumption can be interpreted as allowing for "burn-in" period for unit under only the default action. The second assumption ensures enough diversity among realized features on so that a least-squares projection can be invoked to identify Bn, formalized in the following proposition.
Let by be a (To, 1) 1 vector with entries [[ a 2n do(d) Under Assumptions 7 and Sn is identifiable from An and by-
Proposition 3
The next assumption is the requirement that any action level d of interest is carried out across large enough number of units. Moreover, they must be sufficiently separated in time from follow-up non-default action levels that their parameters can be identified. Assumption 4 is a counterpart to Assumption 2
ensuring linear independence of feature trajectories.
Assumption 3 For a given intervention level d, assume vhr(d, 1', (1) + 0 for all 1, and that there at least one time index I, and one sel Na containing all unils n where d d, such that: (i) Na = |Nal 2 r: (ii) vn Na, =.... =d d intervention level d.
O, where ka is the dimensionality of the
Assumption 4 Let n1, (1,,,22,--- m12....,,N, index the elements of a set of units Na. For 1,, ki let Aae be Na X r matrix where each row i ,......u Na given by on, (2 (2littt'-1, )n.)) realized under regime do (o(d2t+t'-] We assume that each Ady is full rank with left pseudo-inverse Ada.
This leads to the main result of this section:
Theorem Assume we have a dataset of N units and T time points (d), ,. 1::T 2n) generated by model partially specified by Eq. (7). Assume also knowledge of the conditional expectations E[X; n an Zn; do(d" []]] and basis functions du (=:) for all 1,2,.... and all S S and 1 < < T. Then, under Assumptions 4 applied to all individuals and all intervention levels d that appear our dataset, we have that all Bn and all (i, ..)) are identifiable.
Notice that nothing above requires prior knowledge of all intervention levels which will exist, and could be applied rolling basis new interventions to start synchronously at = 1: the framing of the theory assumes so to simplify presentation, with the only requirement being that each unit is given a "burn-in" period of at least To steps and that each new intervention level dl is applied to at least units which have not been perturbed recently by relatively novel interventions. Ultimately, the gap between a novel intervention level and the next one should be sufficiently large according to the complexity of the model, as indexed by This makes explicit that do not get free identification: the more complex the domain requirements, the more information we need, both in terms of the time waited and the number of observations.
invented. There no need for all units
3 Algorithm and Statistical Inference
This section introduces learning algorithm, as well as ways of quantifying uncertainty in prediction.
We also allow for the learning of adaptive basis functions Φ. As Eq. (T) does not define a generative model, which will be necessary for multiple-steps ahead prediction, for the remainder of this section we will assume the likelihood implied by X x -n, do do(d22))
[x =n, df:1), )²))) Here, the conditional mean J(.:5.) :) is given by Eq. (1) As o² not affected by D, it can be easily shown to not require further identification results. In general, if parameters in our likelihood are implied by finite set of estimating equations, we can parameterize
the multilinear form analogous to (T) and repeat the analysis of the previous section. In practice, the functionals (such as the conditional expectations of X₂) used an estimating equation are unknown. Matrix factorization methods can be used directly by first applying smoothing method to the data to get plug-in estimates of these functionals 42), but we will adopt instead likelihood-bart -aseed approach.
3.1 Algorithm: CSI-VAE
We treat each 3n as random effects vector, giving each entry But an independent zero-mean Gaussian prior with variance oz. Along with all hyperparameters, optimize the (marginal) log-

---

likelihood by mudieent-based optimization. With Gaussian likelihood, the posterior and filtering distribution of each 5, can be computed in closed form, but in our implementation we used a black-box amortized variational inference framework 227 30 34 anyway that can be readily put together without specialized formulas and is easily adaptable to other likelihoods. We use a mean- field Gaussian approximation with posterior mean and variances produced by gated recurrent unit (GRU) model [14] composed with a multilayer perceptron (MLP). Therefore, the approximate posterior at É time steps follows from Mq.8,n == MLP(GRU, ni.1 (d. ,a/h"" )nn),, and log 0,8,n ;= MLP(GGUU, n0.2 (dh ==1 z,)). Prediction for X(+1:++ is done by sampling M trajectories, where for each trajectory we first sample a new Br, from the mean-field Gaussian approximation. We set each forward X to be the corresponding marginal Monte Carlo average. Each basis vector of parameterized via another GRU model such that d = MLP(GRU, & (x]]t-1 z,)). Although the theory in the previous section relies on known basis functions, we could train this GRU up to a threshold time point of To and freeze from that point, each 3n being defined as the unit-level coefficient vector under this basis. In practice, we found that doing end-to-end learning provides a modest improvement, and this will be the preferred approach in the experiments. We do find that it is sometimes more stable to condition only up to time To for the computation of the variational posterior of f. Hence, we adopt this pipeline for any results reported in this paper, unless otherwise specified. We call our method Compositional Sequential Intervention Variational Autoencoder (CSI-VAE).
3.2 Distribution-free Uncertainty Quantification
provides prediction intervals with coverage guarantees. The intervals are computed using a calibration set of labeled samples and include the future samples with non-asymptotic bourde probability. We consider Split Conformal Prediction in two setups.
Conformal prediction (CP) 46 [9]
1. Hold-out predictions. We have a set of historical users, n I...., N, whose behavior has been observed until time I + A. The task is to predict the behaviour at time | + A of a new user, n = NN + 1, who has been observed up to time t, i.c. to predict X X++A given N 11 41
1:t
and DV If we assume we have used the history of the new user X2= = 1,..., N + 1, train the model, calibration and test samples are exchangeable.
2. Next-intervention predictions. We have observed the behavior of all users, n ,,...,N
1. up to time and aim to predict the effects of the next intervention, which happens at time t 1 for all users, i.e. D + O, holding D = O. The task is to predict X under do(d):*+4) for n = L.....NN N + 1. but in what follows we will drop the explicit do(.) indexing to keep notation lighter, referring explicitly to past observed D only (as its sampling distribution matters), and assuming sequential ignorability by assumption or randomization. Calibration and test are not exchangeable, because i) the joint distribution after time f,, Pr (X++A, X1a, Di:t) may be different from the one before time i, Pc ~ (X, (
n x X!:t' D ), t' <t - A and ii) we only used X D t' = 1....f, f, for training.
CP algorithms are applied on top of given the cross-sectional sample space of any X,. In our case, the underlying model is f. which predicts the expected user behaviour time L A given the user history up to time L, i.e. X++A := X1 1
models, where we will use X' to denote
Zn ~ E| n. 1 Zn. do d .0* 0++22++ }]] for a chosen implicit d'+1 and with 0 2--2++2 denoting the default action 0 being taken t ++2 2,....t A.
f (D).
Setup
Calibration and test scores, (Sn = |x++ f(D X1:,, Zn)ll and X - f(DN+ X Zn)l are exchangeable, i.c. Prob(S , Sy + 1) = Prob(Sa(1)...... (I)N where a is any permutation of [L....,N N + 1}. The Quantile Lemma, e.g. Lemma 1 of [44], implies the prediction interval
1
SN+1
C= (xr E 2 X < Qat
N
Qa inf
5 x na na = [(1 + +N((1 - a]1
is valid in the sense it obeys
Prob (XM+A E EC) 21 21aa. a,
(9)
where the probability is over the calibration and test samples.
7

---

S The prediction intervals defined in may not be valid, i.. 9 may not hold, because the calibration and test samples, S = x n ,,..,, .N, are not exchangeable. Theorem provides a bound on the coverage gap, i.e. measure of violation of (9), under the assumption that the distribution shift is controlled by perturbation parameter, E 0.
f(D x 7n)l with t' I and t' > t,
2
Theorem Assume we have N calibration samples,
= - " Za)l, t' t 1=1....,,N
(10)
where in the time user n experienced the last intervention before t. Assume there exists > 0 such that, for all n,
pr pr = c)pcc)p (po(St+A),
(11)
where pr and pe are the (unknown) densities of the test and calibration distributions, ps a bounded arbitrary shift density, and Pmain = min (c > 0. Then,
Prob (X EC) 2 -
(12)
Pmnnn 1 1-6
In the proof, we use the likelihood setiihood-ratio-oe vatio-weighting approach of 44 to obtain the empirical test distribution from the calibration samples and bound its quantile from below. The statement follows from standard inequalities on the mean and variance of the c-perturbed distribution. We prefer this approach to conformal prediction adaptive models for time series, e.g. [20] or (6), for two reasons: i)
adaptive schemes have asymptotic coverage guarantees that can not be used to estimate the uncertainty
on single time step and ii) optimized density estimates byproduct of our prediction model.
Further Related Work
9 leverages Bayesian structural time-series models to estimate causal effects, and motivates our exponential decay model in Eq. (3). Unlike [9], who focus single interventions, our model explicitly addresses the complexity arising from sequential interventions, taking a more detailed perspective the dynamic interplay of treatments
time. Several approaches for conformal
prediction in causal inference have emerged in recent years [e.g..[44, 28), including matrix completion 21] and synthetic controls [12.. Our focus has been individual nor average treatment effects, but directly on expected potential outcomes, similarly to the work on causal matrix completion 42]. Unlike the matrix completion literature, we focused on prediction problems, out-of- sample for both units and time steps. Moreover, the matrix completion literature is usually framed terms of marginal expectations E[X; (d)), as opposed to conditional expectations (notice that considers covariates, but those are included as part of the generative model). Marginal models have some advantages when modeling multiple-step-ahen computational considerations that we leave for future work. Compositionallity based on additivity instead of the factorization trick is discussed by e.g. [29]. but additivity lacks the connection to explicit principle for universal function approximation. Finally, there is rich literature on confounding adjustment where exogeneity Dú, cannot assumed, see 22 for a textbook
effects [18), but they also involve complex
treatment. We can rely on standard approaches of sequential ignorability [36 to justify our method in the absence of randomization.
5 Experiments
We run number of synthetic and semi-synthetic experiments to evaluate the performance of the DSI approach. In this section, we summarize our experimental set-up and main results. The code for reproducing all results and figures is available o In Appendix we provide a detailed description of the datasets and models. Appendix [ we present further analysis and more results. Finally, in Appendix E we present illustration of uncertainty quantification results.
"The code is available https://gitubb.cm/jjaiinnnnuuuuu VAE

---

Table 1: Main experimental results, averaged mean squared root error over five different seeds.
Full Synthetic
Semi Semi-Synthetic Spotify
Model
T+1
T+2
T+3
T+4
T+5
T+1
T+2
T+3
T+4
T+5
CSI- ESI-VAE-1
36.33
41.46
41.73
41.12
41.32
68.23
82.94
83.53
81.97
79.63
CSI-VAE-2
97.80
118.25
117.79
127.25
135.03
253. 85
312.53
305.08
303.68
302.83
CSI-V
138.78
164.02
141.71
132.59
125.55
757.94
937.07
800.55
704.66
634.72
GRU-O
229.72
269.66
220.95
208.30
188.43
215. 42
260.65
193,41
137.20
117.06
GRU-1
230.76
270.83
220.93
208.33
184.92
223.60
269.69
205.91
141.53
126.36
GRU-2
93.73
101.03
118.01
88.53
132.28
154.18
187.42
177.96
133.36
127.58
L.TTM
114.71
126.65
137.12
105.22
137.19
130.35
156.02
133.28
94.35
85.92
Transformer
111.66
122.08
150.57
175.84
87.89
133. 42
157.66
154.66
164.70
158.03
Table 2: P-values from tao-sample t-tests against the null hypothesis that models perform no better on average than its counterpart. Significance at 0.05 indicated with one asterisk, 0.001 with two.
Syddeed
Serui Syrthetic Speciy
Mekd
C3E-
CS--AE--
CRU-4
GRU-1
GRIF2
15134
CSE-AD--
CSE-VE--
GRR50
GRUSS
LSTM
Tranffera:
GRI GREWAFT
0:00
<0100
<0000
<0
CAL-AB-2
0.897
0.24
0.37
(1.33
0.086
1..41
< -RaY"
GR3-4
3350
< <0a00m
4.733
(L223
< <R03"
(LUB)
G
0.555
0.255
<.055"
LSTM
0.755
Datasets and oracular simulators. The first step to assess intervention predictions is to build set of proper ground truth test beds, where interventions. We build two sets of oracular simulators. (1) The Fully-synthetic simulator is constructed primarily based random models following our parameterization in Section [ (2) The Semi-synthetic simulator is constructed based on a real-world dataset from Spotify which aims to predict skip behavior of users based on their past interaction history. See details in Appendix C For each type of simulator, we generate 5 different versions with different random seeds. From simulator (1), construct a simulated dataset of size 50, 000, containing 5 different interventions happening to the units at any time after an initial burn in period of To 10, although only maximum 3 of these can be observed for any given unit. We set (2), we construct simulated datasets of size 3, 000, again containing 5 different interventions, an initial period To- = 25, a maximum of different interventions per unit, and r 10. The task is to predict outcomes for interventions not applied yet within any given unit (i.e., at least from the 2 options left). In simulator (2), parameters & and 3 are learned from artificial, but inspired by the process of showing different proportions of track types to an user in Spotify-like environment. For both setups, we use a data ratio of 0.7, 0.1, 0.2 for training, validation and test, respectively. We report the final results in Table [ with another constructed holdout set (5, 000 points for the fully-synthetic case, and 1, 000 for the semi-synthetic one).
control different levels of combinations of
as the true dimensionalitty of the model. For simulator
data. Interventions are
Compared models. We implemented three variations: (1) CSI-VAE-1 follows exactly setup in Section 2 (2) CSI-VAE-2 can be considered as an ablation study, which relaxes the product form of Eq. and replace it with a black-box MLP applied directly to (dh, Gn, vit)) (and hence may not guarantee identifiabilityy): (3) CSI-VAE-3 E-3 is another ablation study, where the equation for 1á (Eq. (2)) is replaced by black-box function taking the sequence of actions D standard temeseries. We compare our model against: (1) GRU-0, a black-box gated recurrent unit (GRU, (B) composed with MLP,, using only the past history of X1# and 7n; (2) GRU-1, another GRU composed with MLP that takes into account also the latest intervention D and (3) GRU-2, which uses not only n but also the entire past history Di:t just like CSI-VAE. In general, GRU-2 can be considered
very strong and generic black-box baseline model. In addition, we conduct experiments comparing other popular black-box baseline models: (4) LSTM, 233 and (5) Transformer, (45). For those, we use the same input setup as in the case of the GRU-2 model.
Results. Each experiment was repeated 5 times, using Adam (26) a learning rate of 0.01, with 50 epochs in the fully-synthetic case and 100 for the semi-synthetic, which was enough for convergence. We selected the best iteration point based on a small holdout set. The main results are presented in Tables T and which show the superiority of our model against strong baselines. We also observed
httpp:://ppe..sp0iff..cm

---

200
Frver (467))
200
150
1413
C313
DXL4
CRILD
DHULS
DELL
INTR
0313
0M410
18.2
DOIN
Error (MRE))
202
Figure 2: Top: 5-run evaluation of test mean squared error on the fully-synthetic (left) and semi- synthetic cases (case). CSI-3 was removed the right due very high errors. Bottom: how change as training sizes are increased, CSI-1 vs. GRU-2 (left: fully-synthetic, right: semi-synthetic).
that the identifiability results and compositional interactions of intervention effects are both critical, as evidenced by the drop performance for CSI-VAE 2 and Table
In Appendix D] we provide the following further experiments: (1) different choices of r (summary: using r less than the true value gracefully underfits, while there is evidence of some overfitting for choices which overshoot the true value mitigated by regularization); (2) different sizes of model consistently outperforms
the training data (summary: even with more data provided, the black-box models by significant margin; we show that a generic black-box cannot solve this problem by simply feeding more data); and (3) demonstration conformal prediction that allow us to better calibrate the predictive coverage compared to vanilla model-based prediction.
Conclusion
We introduced an approach for predictions of sequences under hypothetically controlled actions, with a careful accounting of when extrapolation to unseen sequences of controls is warranted.
Findings. Embeddings important given sparse categorical sequential data 45 However, large combinatorial interventional problems benefit from models that carefully lay down conditions for the identification of such embeddings. Naïve sequential models, however flexible, cannot fully do the heavy lifting of generalizing in the absence of structure. Information has to come from somewhere.
Limitations. Unlike the traditional synthetic control literature [ we assume model for time effects based on autoregression and truncated or parametric time/intervention interactions. While it is possible to empirically evaluate the predictive abilities of the model using validation sample, high-stakes applications (such as major interventions to counteract the effect of a pandemic) should take into consideration that uncontrolled distribution shifts may take place, and careful modeling of such shifts should be added to any analytical pipeline to avoid damaging societal implications.
Future work. Besides allowing for an explicit parameterization of drifts, accounting for unmeasured confounding that may take place among past actions states necessary to increase the applica- bility of the method. Moreover, when each action level d is itself combination of cross-sectional actions, causal energy-based models such as can combined with the compositional factorization idea to also generalize to previously unseen cross-sectional combinations of actions.
10

---

Acknowledgements
We thank Zhenwen Dai and Georges Dupret from Spotify for numerous helpful comments and suggestions. JY and RS were partially funded by the EPSRC Open Fellowship The Causal Continuum:: Transforming Modelling and Computation support of the UKRI Al programme, and the Engineering and Physical Sciences Research Council, for EPSRC CHAR-EPSRC Al Hub for Causality in Healthcare Al with Real Data (grant number EP/Y028856/1). The authors would also like to thank the four anonymous reviewers for fruitful suggestions on how to improve the presentation of this paper.
Causal Inference, EP/W024330/1. RS acknowledges
References
[1] A. Abadie. Using synthetic controls: Feasibility, data requirements, and methodological aspects. Journal of Economic Literature, 59:391-425, 2021.
[2] A. Agarwal, A. Agarwal, and S. Vijaykumar. Synthetic combinations: a causal inference framework for combinatorial interventions. Advances in Neural Information Processing Systems 36 (NeurIPS 2023), pages 19195-19216, 2023.
[3] A. Agarwal, M. Dahleh, D. Shah, and D. Shen. Causal matrix completion. Proceedings of Thirty Sixth Conference on Learning Theory, pages 3821-3826, 2023.
[4] A. Agarwal, D. Shah, and D. Shen. Synthetic interventions. arXiv:2006.0769/ (econ), 2023.
[5] V. Aglietti, X. Lu, A. Paleyes, and J. González. Causal Bayesian optimization. Proceedings of the Twenty Third International Conference on Artificial Infelligence and Statistics, 34:3155-
3 3164, 2020.
[6] A. Anastasios, E. Candès, and R. Tibshirani. Conformal PID control for time series prediction. Advances Neural Information Processing Systems, 36, 2024.
[7] S. Athey, M. Bayati, N. Doudchenko, G. Imbens, and K. Khosravi. Matrix completion methods for causal panel data models. Journal of the American Statistical Association, 116:1716-7700, 2021.
[8] G. Bravo-Hermsdorff, D. Watson, J. Yu, J. Zeitler, and R. Silva. Intervention generalization: view from factor graph models. Advances in Neural Information Processing Systems 36 (NeurIPS 2023), pages 43662-43675, 2023.
[9] K. Brodersen, F. Gallusser, J. Koehler, N. Remy, and S. Scott. Inferring causal impact using Bayesian structural time-series models. Annals of Applied Statistics, appear, 2015.
[10] B. Brost, R. Mehrotra, and Jehan. The music streaming sessions dataset. In The World Wide Web Conference, pages 2594-2600, 2019.
[11] B. Chakraborty and E. Moodie. Statistical Methods for Dynamic Treatment Regimes: Reinforce- ment Learning, Causal Inference, and Personalized Medicine. Springer, 2013.
[12] V. Chernozhukov, K. Wuthrichr, and Y. Zhu. An exact and robust conformal inference method for counterfactual and synthetic controls. Journal Of The American Statistical Association, OCT 2 2021.
116(536):
[13] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio. Learning phrase representations using RNN encoder-decoder for statistical machine translation. Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP 2014), pages 1724-1734, 2014.
[14] Chung, C. Gulcehre, K. Cho, and Bengio. Empirical evaluation of gated recurrent neural networks on sequence modeling. NIPS 2014 Workshop on Deep Learning. 2014.
[15] J. Correa and E. Bareinboim. A calculus for stochastic interventions:causal effect identification and surrogate experiments. Proceedings of the AAAI Conference on Artificial Intelligence (AAAI 2020), 34: 101110000010000 2020.

---

[16] P. Dawid. Decision-theoree foundations for statistical causality. Journal of Causal Inference, 9(56):39-77, 2021.
[17] V. Didelez, A. P. Dawid, and S. Geneletti. Direct and indirect effects of sequential treat- ments. Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence (UAI2006), 2006.
[18] R. Evans and V. Didelez. Parameterizing and simulating from causal models. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 2024.
[19] M. Fontana, G. Zeni, and S. Vantini. Conformal prediction: A unified review of theory and new challenges. Bernoulli, 29(1):1-23, FEB 2023.
[20] Isaac Gibbs and Emmanuel Candes. Adaptive conformal inference under distribution shift. Advances in Neural Information Processing Systems, 34: 34:1660-1672, 2021.
[21] Y. Gui, R. Barber, and C. Ma. Conformalized matrix completion. Advances in Neural Informa- tion Processing Systems, 36:4820-4844, 2023.
[22] M. Hernán and Robins J. Causal Inference: What If. Chapman & Hall/CRC, 2020.
[23] Hochreiter and J. Schmidhuber. Long short-term memory. Neural Computation, 9:1735-1780
1997.
[24] G. W. Imbens and D. B. Rubin. Causal Inference in Statistics, Social, and Biomedical Sciences: An Introduction. Cambridge University Press, 2015.
[25] J. Kaddour, Y. Zhu, Q. Liu, M. J. Kusner, and R. Silva. Causal effect inference for structured treatments. Advances in Neural Information Processing Systems 34 (NeurIPS 2021), pages 24841-24854 2021.
[26] D. Kingma and J. Ba. Adam: A method for stochastic optimization. 3rd International Conference for Learning Representations, 2015.
[27] D. Kingma and M. Welling. Auto-encoding variational Bayes. International Conference on Learning Representationss 2014.
[28] L. Lei and Candès. Conformal inference of counterfactuals and individual treatment effects. Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(5):911-99 1-998, 2021.
[29] M. Lotfollahi, Susmelj, C. DeDonno, L. Hetzel, Y. Ji, I. Ibarral, Srivatsan, M. Naghipour- far, R. Daza, B. Martin, J. Shendure, J. McFaline-Figueroa ad P. Boyeau, F. A. Wolf, N. Yakubova, S. Günnemann, C. Trapnell, D. Lopez-Paz, and F. Theis. Predicting cellu- lar responses to complex perturbations in high-throughput screens. Molecular Systems Biology, 19:e11517, 2023.
[30] A. Mnih and K. Gregor. Neural variational inference and learning in belief networks. In International Conference on Machine Learning, pages 1791-1799. PMLR, 2014.
[31] S. Murphy. Optimal dynamic treatment regimes. Journal of the Royal Statistical Society Series B, 6(2):331-355, 2003.
[32] R. Nabi, T. McNutt, and I. Shpitser. Semiparametric causal sufficient dimension reduction of multidimensional treatments. The 38th Conference on Uncertainty in Artificial Intelligence, pages 1445-1455, 2022.
[33] J. Pearl. Causality: Models, Reasoning and Inference, 2nd edition. Cambridge University Press, 2009.
[34] D. Rezende, S. Mohamed, and D. Wierstra. Stochastic backpropagation and approximate inference deep generative models. International conference on machine learning, pages 1278-1286 PMLR, 2014.
12

---

[35] J. Robins. A new approach to causal inference in mortality studies with sustained exposure period a control of the healthy worker survivor effect. Mathematical modelling. 19:11212339-1512, 1986.
[36] J. Robins. Causal inference from complex longitudinal data. In Latent variable modeling and applications to causality, pages 69-117. Springer, 1997.
[37] S. Saengkyongam, Rosenfeld, P. Ravikumar, N. Pfister, and Peters. Identifying representations for intervention extrapolation. International Conference on Learning Representations (ICLR 2024), 2024.
[38] Y. Saito, R. Qingyang, and T. Joachims. Off-policy evaluation for large action spaces via conjunct effect modeling. Proceedings of the 39th International Conference on Machine Learning (ICML 2023), pages 29734 29759, 2023.
[39] Y. Saito, J. Yao, and T. Joachims. POTEC: Off-policy learning for large action spaces via two-stage policy decomposition. Proceeding of the 40th International Conference on Machine Learning (ICML 2024), 2024.
[40] A.M. Schäfer and H.-G. Zimmermann. Recurrent neural networks universal approximators. International Journal of Neural Systems, 17(04):253228 2007.
[41] P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction and Search. Cambridge University Press, 2000.
[42] C. Squires, D. Shen, A. Agarwal, D. Shah, and Uhler. Causal imputation via synthetic interventions. First Conference on Causal Learning and Reasoning (CLeaR 2021), pages 688-77 688-711, 2022.
[43] R. Sutton and A. Barto. Reinforcement Learning: An Introduction. MIT Prsss, 2018.
[44] R. Tibshirani, R. F. Barber, E. Candès, and A. Ramdas. Conformal prediction under covariate shift. Advances in neural information processing systems, 32, 2019.
[45] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. Gomez, t.. Kaiser, and I. Polo- sukhin. Attention is all you need. Advances in Neural Information Processing Systems 31 (NIPS 2017), pages 6000-6010, 2017.
[46] V. Vovk, A. Gammerman, and Shafer. Algorithmic Learning in a Random World. Springer, 2005.
[47] L. Yao, S. Li, Li H. Xue, J. Gao, and A. Zhang. On the estimation of treatment effect with text covariates. Proceedings of the Twenty-Eight International Joint Conference on Artificial Intelligence (IJCAI-19), pages 4 1106-4113, 2019.
[48] Z. Zeng, S. Balakrishnan, Y. Han, and E. Kennedy. Causal inference with high-dimensional discrete covariates. arXiv:2405.06 0000000001188 18 (math), 2024.
[49] J. Zhang and F. Bareinboim. Designing optimal dynamic treatment regimes: causal rein- forcement learning approach. Proceedings of the 37th International Conference on Machine Learning (ICML 2020). 2020.
[50] J. Zhang, K. Greenewald, C. Squires, A. Srivastava, K. Shanmugam, and C. Uhler. Identifiability guarantees for causal disentanglement from soft interventions. Advances Neural Information Processing Systems, 06:5054450299 2023.
13

---

Appendix
A
General Background
Much of traditional causal inference aims at estimating the outcome of a single-shot treatment that is not intertwined within sequence with uncontrolled intermediate outcomes [24] However, many causal questions involve sequential treatments. For example, we might be interested estimating the causal effects of medical treatments, lifestyle habits, employment status, marital status, occupation exposures, etc. Because these treatments may take different values for a single individual over time, they are often refereed as time- vareivaying treatments (22]
In many designs, is assumed that a time-fixed treatment applied to individuals in the population happens at the same time (i.e. there is shared global clock for everyone). A more relaxed setting (and more practical, as described in this paper) is that (categorical) interventions can happen to different individuals different time steps, potentially resulting in pair of identically distributed time-series. The purpose now is to estimate the causal effect or to predict expected outcomes under time-varying treatments as a time-series, rather than estimating/imput imputing missing potential outcome snapshots from a fixed time interval that has taken place in the past.
One typical way of designing experiments under a binary treatment (0 or 1) treatments is to have dg = (0, o,... 0) ("never treated") and di = (1, 1..... :) ("always treated") as the two options in the design space. This allows to answer what the average expected outcome contrasted between these two options. A fine-grained estimation at every time step will require further
35 36 of time-varying
levels of contrast for other regimes, such as d = (0, 1, .................... 1) and d = (1,0, ... O). However, the number of combinations grows exponentially on T, the number of time steps, even when only binary treatments are at the stake.
In this paper, we did not consider identification strategies that may be necessary when confounding is an issue. Conditional ignorability holds when we are able to block unmeasured confounders using observable confounders only. Similarly, in time-varying setting. we achieve sequential conditional ignorability by conditioning on the time-varying covariates, including past actions ((22) and Chapter 4 of 33 ). In our manuscript, we will have little to say about the uses of conditional sequential ignorability and other adjustment techniques for observational studies, as these ideas are well-understood and can be applied on top of our main results. Hence, we decided to focus on treatments that act as-if randomized so that we can explore in more detail our contributions.
B Proofs
(n, x. 2), where function sequences Proposition 1 Let (i) f'n, d..T x...T (d. (nnz(n,
f and fnaz (n, x. I.T 2) are defined for all T E N+ and have codomain R. (ii) ( be given as in Eq. (6) with Thy 1; (iii) the 1-th entry of (nz I:T. 2) be given by (n, ; = on (i 22). where both ui(-) and vx(.,.) have codomain R²"; (iv) vi(x' .2) be defined analogously to fa(d as in Eq. T.. carrying out the fixed-izze z as an extra argument. Then there exists some integer r and three functions a, b and e with codomain
I" so that f(n, d':T x.:T .2) = Ei= a1 .2) x by(n) (d ,T). Moreover a (dd TT)
Ih= mi(d" ,1, T) for some function my : N3 + IR.
Proof. The proof is tedious and the result intuitive, but they help formalize the main motivation for Eqs. ① and (3).
fJ (d=:1) ) x I:T ,2) = Ex fay (dd:T) ( (n, x 2)
Ei Sar (dd:T) Ete we (n) ,2)
= EM Ei' up ( (n) (x (.,T,,2)faa aav
= (x 1:T ,2) x (1(n)) x (dd"T) for 1 f
where r r= -=r,>y Xr and (x 22) = urp urpa(x (x
I": (n(m)) == upp (n) for I = P' I":
and ((ddT)) far (d for l' [1/ro/, where [-] is the ceiling function.
14

---

Finally, each 7 is by definition some (d) From Eq. (6) and the assumptions, fav (d gr (d T, I), and gr (dtit, tt,t) == hy (gr (dd !t' - 1, 1, df f,,). For t' 0 and d > O, we have that by is given by by, (g ,t' 1, 1ug (d' f',t). For = O, (e t) = 1, gr (dlit'-1, implies (ax:2) = II -e-= so hig (d' ,U,T) and the result follows. C
and for t' > 0 and d = 0, hy (v, (o,d,tt (,,)
1,t). By recursive application, this
Proposition 2 Given a fixed value of T, suppose have measurable function :: .Ui_ X1:--1x U D' x Z + RR, S T, where X X1a-- and Z are compact, D' is finite. Let E X1a-1 and let d U Dt. For d = 1, ,dmax. let ma Ut=1 D' - (1..... _T} with ma(d)
t if s.t. d, = d else O. Then, for any > 0. there exist r O(T and measurable functions di x1a-1 x E + Bdd and thas D' - IR and real numbers 3r, such that
dmao
s"P S"PxEUT x deU? EE (x), - Σ x1
II
vhaa (ma Se. (13)
d=1
Proof. First, let
dmaz
f(x, d, 2) :=
Σ (x 2) B1
II
die (m (d)).
Notice that we can construct equivalent definition of D for some fixed / as follows. We for now call the new object Dr:
Rdmux na E (0, ... ,1)vd; va # ea whenever d' * d and oa + of
D' = v E
(14)
Note that there is one-to-one correspondence between De and D', and denote the bijection connect- ing them by g D + D Since both sets are finite, we can approximate with f. if and only if we can approximate f o x g x i by foi x g x i. So, without loss of generality, just consider De as D'. The set of all sequences of treatments is the union of DD over ( = 1,..., T, for some known T:
DJ U D
known T.
(15)
Note that
|D'| < (1 1)dmex
(16)
D (t + 1) dman
(17)
therefore the size of DT is polynomial in T.
Fixing the ground truth function f at D DT gives us function , X1:t-1 x 2 + X.
Analogously, fixing at some E DT also gives a new function the same function space. Since DT is finite, we can enumerate its elements. Thus, for every D(i) DT - [ write down f fixed at (+) as:
2)
(x () (12())))) ... Viidmar (D
(18)
We can then express the vector field as:
f(,DD),)
31 I. dman (D)))
B, Amas Vir (D)2)
(1(()))
f(. D(.DD)))).))
B1 Hmmes dn (n((D (D((DT1)
Br Ilmm the (DD(IDT)))
(.())))
(19)
15

---

By choosing r = |DT| and appropriate values for 3 and (D()) we can make W invertible.
Now, consider the vector field given by the ground truth function fixed at all values of D DT
f(., f(,D00,)))
(20)
f(, ((,DDDD
and consider its transformation under W
f(. f(.D(),))
W-1
(21)
(',, S(,,D(DD)) D..)
Since recurrent neural networks universal approximators 40) for every 8 > we can choose (1," .*. &r such that
(1 (x,,)
f(x, D(),2) 2) : (xx, (,2)
2
SUPXEUP XLi- 2E2
- W-
S.. 6.
(22)
(x(xx)) 2)
Since W is finite matrix, its operator norm well-defined. Then
Φ (xxx))
f(x, D(1,
2
SuPxEU?_, ,i sE2
W
-
(23)
(x,
f(x, D(l
1)
22)
Φ1 (x(x))) 2)
f(x, D(1)
SUUPXEUR N1:-11, sE2 W/12
W-1
(24)
Φr 2)
f(x,D(1 2)
2
Φ (x, 2) :
f(x, D(W)
< 11p sup 211-1,,
W-1
(25)
(x,,))
f(x,D((PT)) 2)
< [ 1
(26)
However,
sup 2 aeDT sezl/(x,d) d, - fxxx d,
(27)
1/2
SUPXEUF_ XIa mag f(x, - f(xx d, z)13)
(28)
1/2
SuPxEU",, pii-1, sE2 Σ
1)) dd,
(29)
dEDT
1/2
[1 2)
5(x, ((1,)))))))
SUPxEU? gli--1
W
(30)
f(x, D(OTT) 2)/
(x, 2)
(31)
Therefore, choosing 8 = = gives us the desired result.
Proposition Let bn be a (To., 1) x vector with entries E[x 1 1: =n, 2n, do(d):t1)) Under Assumptions and 2 3n is identifiable from An and by.
16

---

Proof. Let b, be a (To.. - x) vector where each entry t is given by B[x++] 2 do(0 }] where 0lit+1 denotes all actions kept at the default level of from time time t + 1. We know from standard least-square squares theory, and the given assumptions, that as long as the number To 2 r of observations larger than or equal to the dimension and the matrix full rank, that the system (over-)determined and Bn b4 M
Theorem 1 Assume we have a dataset of N units and time points (dh Ey 2n)) generated by model partially specified by Eq. 47. Assume also knowledge of the conditional expectations 1
E[X; n a--1 Zn; (lo(),")) and basis functions (i(;;)) :)) for all 1,2, .,,,...;? ,
and all S S N and
1 SSSS T. Then, under Assumptions 7 applied to all individuals and all intervention levels d that appear in our dataset, we have that all Bn and all th (i(d,, are identifiabl.
Proof. By Assumptions and 2 and Proposition 3 all 8n are identifiable.
For the next step, let
Tha But
II
1=
and
A.
II
d (,,,,,,t)
0d
Consider now the first time i where intervention level d > 0 is assigned to someone in the entire dataset. By Assumption there exists a dataset Ny with the given properties where, for all n E Na. we have Tha = .unt x Ie= (d t) = Put. since all interventions prior to f have been at level O.
Let ba be N₂ x vector where each entry i is given by E[X aist+t'-1,, An,, do(d);t+ )}, for n;
= 1,2..... hd 1. Consider the Lime-bounded case with free parameters war.1 I 1,2,..... r.
Given Assumption 4 we can solve for 7//+ = A ba, the vector with entries [++
1+1 = nha x As 7/ identifies wa,1- Also by 3 at time point 2-+2 we have that tha n = "/y th 1,
1,t
Vnl ++1 x wat,1, with non-zero Mht by Assumption 3 this also
for all units in Na. which by analogy to the previous paragraph, identifies vi(dd t + 1, + 2) and consequently wat22 As this goes all the way to t + ka- 1, this identifies all of s/y(d,
)).
If intervention d level is the time- unbounddd model of Eq. B by Assumption 3 we also identify
ag a and a. as they are defined by functions 1 (d, t + 1. +') (d, 1, for t' = 1, 2, 3 that we established identified by the reasoning the previous paragraph. As a (Wiat) W2d! of (w)a)3 x W24t ald 3 and
(wia)) W2nt ) x w2dt
2
dddl a1
( W2dt a (wid))) ) x W2dl
a
a
solving for the above results in
( 222 a'll
0-
Wid!
a a
0 0
W2d!
o (Wi1)) o (Wiw )3
Widl = a ( x Wadl
So far we have established identification of the (dd. functions for the first intervention level d that appears in the dataset (the value d is not unique, as it is possible to assign different intervention levels d' in parallel at time t, but to different units). We now must show identification for the remaining interventions d > 0 levels that take place in dataset after time L.
Assume the induction hypothesis that lg the s-th unique time intervention level >> 0 is assigned to some unit in the dataset, and all intervention levels d' that appeared up prior to ta have their
17

---

functions 4d(d",, identified. We showed this is the case for 1 t, which is our base case. We assume this be true for t, and we will show that this also holds for 1-+1-
Let be an intervention level assigned at 1x+1 and let Nu the corresponding subset of data assumed to exist by Assumption For any E NNa let 8- be the last time non-zero intervention d_ has been assigned to n prior to + 1. Also by Assumption 3 this must have taken place at least ka_ steps in the past. By the induction hypothesis, (kd -)) is identifiable, and therefore so the case for all vir (ka.."?)
) for intervention levels dx taking place for the first time prior to t.. This also implies that Vnt is identified, and the rest of the argument proceeds as in the base case.
Theorem 2
Assume we have N calibration samples, = -(D Xxittn Zn)l, t' t n=
(32)
where in is the time user experienced the last intervention before I. Assume there exists >> 0 such for all n,
Pr Pr( = (1 c)1ccppcs
cps(St
(33)
where pr and pC are the (unknown) densities of the test and calibration distributions, på is a bounded arbitrary shift density, and Pmain min ( St > 0. Then,
Prob (X N+1 EC) 2 1 - a Pmnnn 1-e"
(34)
Proof. Let Pr and Po be the joint distributions of (X,,, D for > t and t' S t and w/l pr(x .. For any t' > t, we can estimate Pr by weighting the calibration samples with why
(x )
n I,...,N, N, i.e.
N Px(X Σ l' ((X, (X D)),, Z
(35)
n=
Conditional on the calibration samples, the (1 - a--th empirical quantile of (r D) is
Qr(1 a²) =infg
7
" (Sn |
1
(36)
inl,
Σ
W
1
+
1(Sn 5
2 1 a
(37)
N
N
N
N
1
infa
Σ
1(Sm <g)
> 21- a - Σ
n1
1(Sn Sg)
(38)
N
Z
N
n=1
n=1
N
1
2
2 inf,
Σ
(Sm 5a)
2 21 VN
Σ
(39)
N
N
n=1
= Qc((1
gap).
(40)
BlweT with Elwx"] and V[w"] being the empirical mean and variance of w² estimated the N calibration samples. Under the theorem
vN where gap- 7 v En (w E[w 2/2 [(w"'
N
Y/w"
A(85') assumptions, we have wh = =(1 (1- + P
Then, E[w"] = 1 "ectBB8ss"
Bia(s" '1l "c(s
e and
(i(w=] = Elw" EEw S 2722 8(s"") The statement follows from
2 (mins (3* v2(8(5*)) S (mins PC( y)z
pc((*))
Prob Prob(Z<a)> a) Prob(Z
if 2 b,
(41)
and the Quantile Lemma on N which implies Prob N411 < Qo(1 -ax)) 2 Os, where a4 obeys
42 "
E
1
a, = a + gap < a +
(42)
E [
mins (C 1-6
18

---

C
Simulators
We design two simulators: (1) fully-synthetic simulator; and (2) semi-synthetic simulator.
C.1
General Setting
The simulators serve as fully controllable oracles to allow us to test the performance of our predictive causal inference problems. In particular, we have the following parameters:
N: the total number of training users.
M: the total number of testing users.
the total number steps in a time series.
To: the number of steps in time series when the intervention is the "default" one (D 0).
K: the number of different interventions (so that D E (0, K 11).
r: the dimensionality the feature space.
Adim: the number of time-invariant covariates.
/: the maximum number of non-zero intervention levels in time series.
Whenever possible, set the same random seeds of 2, 5 to aid reproducibility of results. For the fully-synthetic simulator, a different seed indicates that is a different simulator (as given by random draws of simulator parameters based this sced). This also reflects the randomness coming from neural network parameter initialization, and data splitting when training models. For the semi-synthetic simulator, the simulator parameters and 3 are learned from real world data, but we will need to also draw synthetic parameters for tb. In both cases, using different random seeds can be considered drawing problems where each has unique parameters.
C.2 Fully-Synthetic Simulator
The first simulator is purely synthetic, where all parameters are randomly sampled from some pre- defined distribution. The overall data generation process exactly follows the model we specified in Section
We generate synthetic dataset for training with the following parameters: 5, T = 20, To 10. K 5, N 50, 000, and Adi
We generate 200, 000 samples first place, but only use the first 50, 000 samples for training. To test our performance, we generate additional non- con-overlapping M = 5, 000 samples with the same parameters, this time with T 25. We make predictions on the last 5 time steps based on the first 20.
The intervention effect parameters sampled from the following distributions, each with size of K x For d = 1,2:
Wid U(1,2)
U(1,2)
W3d ~ 2d(,
where ((aa)) b) is the uniform distribution in (a, b). For d 3, 4:
Wid ~ 1-(-,, -1)
Wed (d(22-11
W3d (((110))
We set wid = 0, W2d = O and Way = to reflect the "idle" or "default" intervention level d = 0. Values are designed symmetric around the default intervention value. We randomly generate time-series [DL based on the identification results in Section 2.2 with the following rules:
1. have enough default-intervention points (Ta 10 r, Assumption
2. for each unique unit, non-zero action levels are sampled without replacement. We addition- ally assume least / > K/2 (in this case, 1 3) of the interventions show up in each unit-level time series.
19

---

3. once a d > 0 action happens, there must be at least two consecutive d = 0 actions (see
Theorem D.
4. in the test dataset, at prediction time, we assume unseen intervention in the first 20 time series shows up at exactly step 21, and then the time-series has intervention level of 0 until
the end. For instance, the training sequence is D [0. 0, ,,,11,0,, 0, ... 0, 0] with T = 20, and test sequence is D = [0, 0, ... .4.0,
0, 3, 0j + 02,000,0,0 with T 25.
We generate time series { }, based on Eq. Q using the corresponding Wid; W2d, W3d values. We generate r iid entries BBn ~ N (Pa, on) for each user, as follows:
Nn. ~ NN 1)
On exp exp(pX( (O( 1)).
For each user n. we generate Zdim iid unit-specific covariates {Zn} as:
Zn 30.
The same is used to set up an initial X (the initial starting point of the time series) with random- parameter MLP with input dimension of Zdim and output dimension of 1. Finally, to create the synthetic time series { (Xx-T } (0 comes from the initial starting point), we randomly generate a neural network consisting of single long-short term memory (LSTM) layer followed by a MI.P to generate é at each time step. We further compose & with a sigmoid function so that it is bounded between [- We iteratively sample the (X,) points with the conditional mean from Eq (based on X Zn. and D): plus a standard Gaussian additive error term at each step.
The outcome of this simulator gives us the following data: time-series (X,:") intervention sequence {D [T:T and covariates Zn.
C.3 Semi-Synthetic Simulator
The second simulator is semi-synthetic, in the sense that the simulator parameters are learned based on real world data from the WSDM Cup, organised jointly by Spotify, WSDM, and CrowdA The dataset comes from Spotify 5 which is an online music streaming platform with over 190 million active users interacting with a library of over 40 million tracks. The purpose of this challenge is to understand users' behaviours based on sequential interactions with the streamed content they are presented with, and how this contributes to skip track behaviour. The latter is considered as an important implicit feedback signal.
For a detailed description of the dataset, please refer to 10 We created our simulator based on the "Test_ Set.tar.gz (14G)" file. After unzipping the test data, we randomly chose file named prrggppreesttry__ee cettittor__21800990 20180918_0000 .00000000000.csx"" as our main data source. In our case, we only use the columns named: "session_id" "d"", "session "session_position". "session_lengee "skip_3". Here, each unique "session_id" indicates a different sequence of interactions from an unknown user. The length of the interaction is denoted by "session_length where each step is denoted by "session position". The column "track id clean" indicates the particular track user listens at each position in the process. "skip 3" a boolean value indicating whether most of the track was played by the user. We assign the number of 1 if is "Truc" and O otherwise. We use this value later create X
"track_id "rrack_id_elean,,, and
To process the data, we start by looking at the existing sequence lengths, and select the most frequent sequence length in the data (20) (the actual corresponding session sequence is 10). This accounts for approximately 58% of the original file. After filtering, there are around 228, 460 unique sessions and 301, 783 unique tracks in the data. To further reduce the computational cost of constructing a very large session-track matrix, we randomly sample a sub-population of 6, 000 unique session ids, which then bring entry value is the number of counts a particular track has been listened to over a particular session based on this sub-popun uub-population. We apply singular-value decomposition (SVD) 6 to this session-trad track matrix to get session embeddings and track embeddings. For numerical stability reasons, we save
26, 707 unique tracks. We build the session-track matrix where the
hattpe://open.spotiffy.. hutp://wwwnnyyooocooouttrmeeeeem
20

---

only the top 50 singular values and further normalize them using its empirical mean and standard deviation. For track embeddings, we take the top 10 singular values after the nornalization and assign them into 10 clusters, using constrained K-means to to encourage equal-sized cluster?" Using these clusters, we create a new categorical variable F and its continuous counterparth Φ F). which respectively represent the corresponding cluster and its centroid vector based on the value of its associated "track_id lii__dd_elii _clean" entry. Next, we set X as the cumulative sum of the corresponding binary "skip_3" value over each unique "session idd with added zero-mean Gaussian noise with standard deviation 1/2. We keep the first singular values for the session embedding normalization and consider them to be the session specific covariate Z. We save X, F, ((F) and Z separate matrices as the final output, which have shapes of 6, 000 x 10, 6, 0.0 x 10 x 10, 6, 000 x 10 x 10 and 6, 000 x 5, respectively.
To build the simulator, we further define the following two processes:
P(P. Fla-1 X1:-11 Zn) = softmax (f (f(F x1#-11 Zn)),
and
= = B, +
where ((F)) is the 10-dimensional embedding of categorical variable F2 given by the cluster centroid, and c, ~ (0,02).
Function .(.).) is modeled with deep neural network consisting of single layer GRU and MLP, trained on the real F with cross entropy loss (as a likelihood based model). The second equation is trained on the real X with an ordinary least square (OLS) regression with ridge regularization term (a = 0.01). Error variance is estimated by setting d to be equal to the variance of the residuals. This defines semi-synthetic model calibrated by our pre-processing of the real data.
We now need to define interventions synthetically, as they are not present on the real data. We choose K = 5 as the size of the space of possible interventions. We will define interventions in a way to capture the notion of changing the exposure of tracks in particular cluster for particular user n at time f. Since the real data does not contain information regarding how interventions can influence the behaviour of Spotify users, similarly to the simulator setup for the fully-synthetic one, we use 5 random seed to indicate different possible changes to user behaviour in the Spotify platform.
The intervention effect parameters are sampled from the following distributions, each with size of K x For d = 1
Wid ~ 1.5)
W2d ~ 1.5)
W'd ~ 2.0.,,
For = 3,
Wid ~ (.115, -)
W24 (11.5,,1111111 -1)
W3d 0.5).
We set wha = 0, w244 = 0 and waz = to reflect the "idle" or "default" intervention level d The values are designed to symmetric around the default intervention value. We randomly generate time-series (D*,} based on the identification results in Section 2.2 and the rules we adopted in the full-synthetic simulator in Appendix C.2 Notice that the above does not correspond exactly to the model class we discuss in Section 2 as the product of intervention features v in Eq. (2) is pipeline
with the softmax operation of f(P x Zn)) to get Et which is then composed with & to get (FF)).
We generate a synthetic training dataset with the following parameters: = 10, T = 35, To 25. N = 3, 000, and Zuim = (notice that we have the first 5, 000 for training, but only use the first 000 samples for the actual baseline training, as we have additional experiments to check the influence of more training data points). To test our performance, we generate additional noe-overlapping
M = 1, 000 samples with the same parameters, but this time with 7 = 40 for making predictions on the last time-series steps based on the first 35.
21

---

To sample data from this simulator, we take the initial state of X and F from the pre processed data. We use the Zn vector which comes from the real-world user embedding. Then, we iteratively generate the next Fú via categorical draws from the softmax output based on user embedding 7n all previous Fn and X,, Once a categorical value of Fn drawn, we replace it with embedding form @ Fn), which is the centroid of its cluster. We use this & Fx) to generate X,, plus an additive Gaussian noise with zero mean and homoscedastic variance learned from data.
The outcome of this simulator gives us the following data: time-series (X,²') intervention sequence {D.') and covariates Zn
D
More Results
In this Section, we present more experimental results for our method. This includes: (1) the effect of choosing r; (2) further examples visualization and empirical results.
D.1 Effect of the Choice of
The effect of different choices of presented Fig. 3 with being the ground truth parameter from which the data is generated. We have also plotted the performance of our method under and r 10. We observed that when r smaller than the ground truth, model tend to under-fit the data with slightly higher MSE error. When is bigger than the actual value, the MSE error can be further reduced, but resulting in very large variance. We also show that this can be mitigated with further regularization, such applying Lı and La penalty terms during the likelihood learning process.
The insight we draw from this is that, when applying our approach to real-world problems, we could potentially choose a higher without much concern about what an exact value should be. Using large leads to higher variance but in general performs reasonably well compare to knowing the true r (here, 5). We also notice that applying regularization techniques, even when we know the right dimension of r, can be beneficial. This leads to more stable performance (that is, with lower variance), but does not change much of the mean estimate shown in Fig. 3 This maybe partially caused by the gradient optimization process when training the deep neural network. Hence, claim that the main price to paid for a large choice of is that this makes data requirements stricter (such as the number of steps prior to the first intervention) in order to guarantee identification.
120
100
Error (MSE)
60
Figure 3: Effect of changing for the fully-synthetic dataset.
22

---

D.2 Additional Figures
R.S. (Bost) Recons Peth for user
Pool V.S. (Bes0 Rocons Post for User
Real VSS. (Busti Recors Path user
Neon
Rocan
100 U.S. (Bost) Recons Peth Urrr
Seal (Bect Recons Port for User
Real V.S. (Best) Recore Path user
Roll
Tere
Figure Examples of reconstruction of training data for CSI-VAB-1 model.
Rea VSS (Bee)) Reco to Path for User
Real V.S. (Beet) Recons Part for User
Real VSS (Beet) Reconss Fath for User
ovar Bitchio
V.S. (Best) Roco Path User
Real'sS (Bert) Recons Pach for User
Roal (Best) Rocons Rath User
en and 20 hed Toxe reeall
-en Halls
Roun boo0
Rou book
Time
Figure 5: Model-based uncertainty quantification of reconstruction for CSI-VAE-I.
Time Scries of for User
Time Series of for User
Time Series of X User
Cenrol Perh
Contrll Rat?
-10
Control Safh
25
Falt
0C
Fath
Fath
1090
Figure 6: Demonstration of prediction, in the synthetic data case, for CSI-VAE-I.
Histogram at
Histogram Erior at
Hetooram of Emor T
300
132
Vake
VA
-40
Figure 7: Residual distribution examples for CSI-VAE-1. In general, we observe a very long tail effect across model predictions.
E Uncertainty Quantification
For the sake of illustration, we generate another new dataset with a random seed of 42 using the synthetic simulator, based on the following specifications: 5, T 20, To 10, K 5.
23

---

N = 50, 000, Adim 5, 1 = 3. We also generated two further datasets of M = 2, 2,000, one for calibration and one for testing.
E.1 Plug-in Model-based Uncertainty Quantification
We first present the model-based uncertainty quantification at level of a = 0.95, in Fig. 8 Given estimates of the conditional mean and homoscedastic error variance o² learned from data, we can calculate predictive intervals as (f - 1.96 x a. - 1.96 x g). For purely model based uncertainty quantification with plug-in in parameter estimates, we do not need an additional calibration dataset, but it is an obviously problematic one as it does not take into account estimation error and sampling variability. The coverage rate is 57. 75% on the test dataset without further calibration.
Modal Prad coen Interveis and Nasss a Suboet instoncus, Coverage 57.755
rendution rnterval
Figure 8: Plug-in model-based predictive interval and coverage.
E.2 Conformal Prediction
We then present the model-free uncertainty quantification at level of a = 0.95 (based on the Sctup in Section 3.2), shown in Fig. B2 We use the absolute value of the predictive residual
conformity score function, ly - f(x) We calibrate our model output based on the calibration dataset and then apply it on the test dataset. We obtained coverage rate of 95.30%. This is significant improvement
upon what we see in Fig. X and demonstrates the effectiveness of the conformal prediction approach.
Moder Precection Interveas and thue Valiass for Subortt of 20 Text Instancers Coverage 1a 2%%
300
Foe valae
-20
300
WT-
m
Figure 9: Conformal prediction predictive interval and coverage.
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We explain the particular scope in terms of examples believe to be intuitive, why we see this is a challenge and highlight what to expect from off-the-shelf solutions. We point out empirical gains particular regime of data ("sparse data conditions") to avoid overgeneralizing where we add value, and the gains are explicitly presented in the empirical section of the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
The experimental section outlines scenarios where we do not do well, partic- ularly on low sample size regimes and undershooting the right dimensionality of the true process.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model vell-specifiiattin, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
25

---

3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Section 2 in particular lists assumptions explicitly.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
Data sources and references are properly accounted for, no-frills view of the main algorithm is given in Section 3 that directly allows for implementing the required steps using off-the-shelf libraries. Appendices delineate in full detail how reproduce the necessary and important preprocessing steps.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
26

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Yes, Jupyter notebooks with data preprocessing pipelines and running of baselines is provided.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cct
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Yes, the main text provides the high level choices made on running the algorithms and the details are in the appendices.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Statistical tests for pairwise comparison of methods are provided, as well as full boxplots for descriptive summarization of results.
Guidelines:
27

---

The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report I-sigma error bars, but one should state The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
As described in the Experiments section and appendices.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: httttp:://eeuipss wwwpeewwww.pppueeeeeeeeeeeeeeeee cc/puullic/thiiccsuueeellll..eee
Answer: [Yes]
Justification:
We have read and complied with the Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdictio))
10. Broader Impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The introduction lines up why the problem is important, while the Conclusions mentions briefly possible misuses to be avoided.
28

---

Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considcrations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning conrectly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses such risks beyond any of the usual issues that generic predictive model may bring.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer:
[Yes]
Justification:
A link to the main empirical data source used is provided in the Experiments
section.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
29

---

The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the palaaets.ppperssttt paperssii thcode. com con/datases package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The assets (code) are provides as Jupyter notehooks.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
30

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
31


Reviews:
[{"review_id": "m6tQ4YGXoh", "reviewer": "Reviewer_Ni2f", "summary": "The authors propose CSI-VAE, a method that solves the problem of forecasting potential outcomes in multiple interventions. In particular, the authors seem to propose a solution to the problem of a large (combinatorial) space of future treatment plans using the controlled past treatment sequences for inference.", "questions": "I find the paper to be quite well done, my only remark is wrt the benchmark settings which I would like to see addressed in a rebuttal.", "limitations": "see above", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "STRENGTHS\n\n* The paper is well structured, written, and generally easy to follow. \n\n* Including an automated (distribution free) uncertainty quantification using CP is a nice to have and will definitely improve adoption into practice. in particular in medical applications\n\n* Focusing on identification in this tricky area is very welcome and I would encourage other authors to also include at least a discussion on identification, thank you for this", "weaknesses": "WEAKNESSES\n\n* The problem sounds fairly general which would allow other treatment effects over time models to also be relevant to discuss and in particular benchmark against. Currently, it seems the authors only benchmark against ablations of CSI-VAE and naive methods using GRU? Am I mistaken?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review, and for finding that our paper is **“quite well done”** and **“generally easy to follow”**.\n\nThe question about the benchmark is a good opportunity for further clarifications. In order to keep the paper focused, we stripped away complementary discussions about dealing with confounding to directly address the value of taking a structured approach to causal extrapolation. In the benchmark, this means asking directly whether a black-box method can learn to smooth the gaps in the data by itself, without the structure we build. Several methods could be used. We focused on the GRU because it is a well established method for learning sequential predictions out of categorical input sequences - more flexible than RNNs, less complex than LSTMs or transformers, which would be major overkills: we already overfit with GRUs, a LSTM or beyond would not help (and did not, in preliminary experiments). Autoregressive linear models were also attempted in simpler initial versions of the dataset, but they were already underfitting badly and we did not report them.\n\nWe added further methods (LSTMs, transformers) as requested to the shared rebuttal box. Code has been updated (not uploaded yet, we don't think we are allowed to), and it basically replaces calls to GRU with calls to other methods within the PyTorch library. As we anticipated and had seen in preliminary runs, these more complex methods didn't really add much of substance to our central thesis. \n\nThanks again for your feedback, and happy to take further questions!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "EZmN8YhyjB", "reviewer": "Reviewer_yF1i", "summary": "This paper studies treatment effect estimation under sequential, discrete interventions.\nIt is assumed that all treatments are fully independent and impact all future outcomes.\nIt studies the problem of \"causal extrapolation\" where some combinations of sequential interventions may not have been observed during training.\nParametric assumptions on the causal mechanisms are made and used for generalising to unseen interventions.\nA matrix factorisation approach is developed to fill in the gaps for the unseen interventions.\nThe approach is empirically tested on synthetic and semi-synthetic data.", "questions": "- Eq. 1: Where does this come from? Can any treatment effect that follows the causal structure in Fig. 1 be expressed like this? Or is this a simplifying modelling assumption?\n- L136: What are ANOVA models?\n- L123: What does the last sentence in this paragraph mean?\n- Equation in L135: What does this function correspond to? The equation seems not to hold for general functions: e.g. if the dimension of the two RHS functions is 1, then this is certainly not a statement that is true in general. Under what conditions can you decompose a function like this? What are the terms on the RHS? How do they fit into the initial problem setting?\n- Prop. 3: Does this assume that equation 1 holds?", "limitations": "- Fig. 1: The assumed causal structure severely limits the applicability of the presented approach. There is no allowed confounding and, more importantly, the interventions are assumed to be independent. The latter assumption would require full randomisation during data generation. It does simplify things a lot, since it avoids having open paths through conditioning on colliders. But, in essence, it is assumed that the data comes from an RCT. Therefore, the contribution is applicable to RCTs that do not cover the full range of combinatorial interventions.\n- The identifiability Assumption 1 means that we essentially have full access to probe each unit (we need to have the \"no intervention\" applied for several time steps). That is, we essentially can perform an RCT on all units. In that case, why not keep the units in the lab, rather than predicting what will happen under future interventions?", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "- This paper studies an important and, in my opinion, an understudied problem in causality: causal extrapolation, as it is called in the paper. Many papers look at assumptions that can be made on the causal structure alone, whereas here, the functional form of the structural equations is considered as well.\n- The paper studies a way to generalize over unseen combinations of interventions. That is an important concern, especially in the time series context.", "weaknesses": "- The mathematical formalisation is very difficult to follow. Wild statements are made without properly introducing the objects, motivating the statement or clarifying what assumptions are made in order for the statement to hold (see questions below). It is unclear which parts follow from the assumption of the causal relationships in Fig. 1 and what additional assumptions or approximations are made.\n- The assumptions on the data generating process are too restrictive for the method to have much hope of being applied to real data (see limitations).\n\n**Minor:**\n\n- L21 and L 86: Typo \"an unit\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your review, and for the much appreciated point that the problem is **“important and, … an understudied problem in causality”**, which is one of the main messages we wanted to convey! We address all clarification questions below.\n\n**Eq. 1: Where does this come from?** \n\nThis is an excellent point, which is why we dedicated the entire Section 2.1 of our paper solely about the representation power of Eq. 1. For a detailed answer to your question, please consult it. Here is a summary.\n\nEq. 1 takes a standard tensor decomposition format. It can be motivated by a variety of starting points, e.g. Taylor series approximations. Its hyperparameter is $r$, for which there are standard ways of choosing it, see e.g. [1]-[4]. \n\nWe address it from two perspectives. *First*, when time is unbounded: there we assume that the true model can be constructed from the exact form in line 149, as motivated by results such as Proposition 1 of [25], which can then be transformed into the structure of Eq. 1. The reasoning is formalized as Proposition 1 in our paper.  *Second*, when time is bounded: we can allow the true data generating functions to belong to a broad class of functions, showing that there exists a finite $r$ that controls the approximation error of Eq. 1 to any a priori degree. This is formalized as Proposition 2 in our paper. \n\n**L136 : ...ANOVA...?**\n\nAnalysis of variance. They are the workhorse of analysis of experiments and widely taught across applied sciences.\n\n**L123: … last sentence in this paragraph…**\n\n$k_d$ is the number of parameters associated with a particular level $d$. The symbol was already used in line 114 to describe the time-bounded model. It does not appear explicitly in the time-unbounded model. Since it has three parameters, we define $k_d = 3$, so that we can refer to it later in e.g. Assumption 3 without having to mention which of the intervention models is used.\n\n**Equation in L135**\n \nIt is a standard way of writing a regression function with categorical inputs. For instance, if we have two categorical inputs $d_1 \\in \\{0, 1, 2}; d_2 \\in {0, 1}$, we can have the following (overcomplete) ANOVA model with parameters $\\theta_{00}, \\theta_{01}, \\theta_{10}, \\theta_{11}, \\theta_{20}, \\theta_{21}$, with indicator functions $f_{(d_1’, d_2’)}(d_1, d_2) \\equiv I(d_1 = d_1’, d_2 = d_2’)$:\n\n$$f(d_1, d_2) = f_{(0, 0)}(d_1, d_2) \\times \\theta_{00} + \\dots + f_{(2, 1)}(d_1, d_2) \\times \\theta_{21}$$\n\nWhen there are other inputs beyond the categorical variables of interest, then each $\\theta$ can be a function of those.\n\nThis appears in many places, e.g. energy functions, exponential families etc., but also all sorts of machine learning methods - as we don’t want to span all exponentially many combinations of the discrete space. In a regression tree with categorical inputs, for instance, we don’t use all combinations, but only those that correspond to paths from the root to a leaf in some tree, with $\\theta$ being the piecewise constant expected value of the output at inputs mapped to that leaf. This is a standard interpretation that appears in e.g. Hastie et al.’s *Elements of Statistical Learning*. One interpretation of what we do in Eq.1 is to use a differentiable embedding of the $d$ vector with dimension $r$, instead of $r$ indicator functions given by a tree of $r$ leaves.\n\n**Prop. 3: …equation 1…?**\n\nIt does assume that Eq. 1 holds (what would $\\beta$ mean otherwise?) The entire Section 2.2. refers to symbols introduced in Eq. 1. \n\nWe therefore argue that the clarifications above fully address weakness #1, “The mathematical formalisation...”\n\n**Fig. 1: The assumed causal structure severely limits the applicability … RCT on all units.**\n\nThanks for raising these points, which quite possibly could come from a few other interested readers too, so it’s very useful to put them to rest. \n\nThe structure in Fig. 1 is far more flexible than the standard structure for sequential interventions: it allows interventions to have direct effects indefinitely into the future. \n\nThe structure in Fig. 1 does allow confounding: it is only depicting what happens when each $D_t$ is controlled. This is the standard graphical way of denoting an intervention in a random variable $D_t$ (e.g, [22, 32]): wipe the edges into it, adopt a different symbol (here, a square) to indicate that’s a fixed index, not a random variable anymore. This structure is compatible with a non-manipulated graph where the entire past causes $D_t$. It is ill-posed to say that “interventions are assumed to be independent”: intervention variables are not random variables, the concept of probabilistic independence does not apply to them.\n\nWe do allow for observational data: the method is introduced *as if* we are allowed to *interpret* each $D_t$ as controlled. To quote line 39, \n\n*“We consider the case where each $D_n^t$ behaves **as if** it was randomized…”*\n\nFor avoidance of doubt, in lines 81 and 321:\n\n*\"sequentially unconfounded with the system by randomization* **or assumption**\"\n\n*“We can rely on standard approaches of sequential ignorability [35] to **justify our method in the absence of randomization**.”*\n\nSee also lines 526-532 (Appendix A).\n\nThat is, *nothing in our results change at all* if confounding can be blocked by observables, e.g., if in the non-intervened graph we have background variables $Z$ pointing to all $D$, and past $D$ variables pointing to all future $D$ variables (we do not explicitly consider the case of variables $Z^t$ in between variables $D^t$, $D^{t + 1}$, but they can be absorbed into $\\phi$ and the do-calculus is still there for us to use.) That is, we can still identify functionals such as Eq. 1, and the conformal prediction results in Section 3.2 do not assume $D$ is randomized. We can make these points more explicitly in the paper, this is a good idea.\n\nWe therefore argue that the clarification above fully addresses weakness #2, “The assumptions….”\n\nThanks again for the questions!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "ObyFcL3OBF", "reviewer": "Reviewer_2QJb", "summary": "This paper considers a special case of a series of interventions where the interventions are categorical and sparse. And interventions can effect later timestamp. This is a form of causal extrapolation. Authors propose to study this using a conditional mean model that utilizes basis functions and subsequently develop corresponding algorithms.", "questions": "1. In equation (1), do you need the lower case x to be a potential outcome? \n2. Is the notation changed in 2.1? f does not appear in equation (1).\n3. For CSI-VAE, why use GRU instead of more modern arch like attention or even LSTM?", "limitations": "The only limitation is how practical these methods are.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. The problem is interesting and well-motivated. \n2. As far as I can tell, the theory is sound.", "weaknesses": "1. Section 2.1 is really hard to parse. \n2. The assumptions in section 2.2 is not explained intuitively and I don’t know how necessary they are. \n3. The experiments are synthetic and semi-synthetic. Although since the nature of the paper is mostly theoretically, I wouldn’t consider that to be a huge issue.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your review, and for agreeing that **“the problem is interesting and well-motivated”**! We would like to see more of that in the community.\n\n**Section 2.1 is really hard to parse**\n\nWe hope the following helps. Eq. 1 takes a standard tensor decomposition format. It can be motivated by a variety of starting points, for instance Taylor series approximations. Its main hyperparameter is $r$, for which there are standard ways of choosing it - see for instance [1]-[4]. \n\nIn Section 2.1, we address it from two perspectives. *First*, when time is unbounded: there we assume that the true model can be constructed from the exact form in line 149, as motivated by results such as Proposition 1 of [25], which can then be transformed into the structure of Eq. 1. The reasoning is formalized as Proposition 1 in our paper.  *Second*, when time is bounded: we can allow the true data generating functions to belong to a broad class of functions, showing that there exists a finite $r$ that controls the approximation error of Eq. 1 to any a priori degree. This is formalized as Proposition 2 in our paper. \n\nIn particular, separating treatment variables from the rest, as in the equation in line 135, can be motivated as follows. First, by recognizing that this is a standard way of writing a regression function with categorical inputs. For instance, if we have two categorical inputs $d_1 \\in \\{0, 1, 2\\}, d_2 \\in \\{0, 1\\}$, we can have the following (overcomplete) ANOVA model with parameters $\\theta_{00}, \\theta_{01}, \\theta_{10}, \\theta_{11}, \\theta_{20}, \\theta_{21}$, with indicator functions $f_{(d_1’, d_2’)}(d_1, d_2) \\equiv I(d_1 = d_1’, d_2 = d_2’)$:\n\n$$f(d_1, d_2) \\equiv f_{(0, 0)}(d_1, d_2) \\times \\theta_{00} + \\dots + f_{(2, 1)}(d_1, d_2) \\times \\theta_{21}$$\n\nWhen there are other inputs beyond the categorical variables of interest, we can make each $\\theta$ a function of those.\n\nSimilar constructions appear in many places, e.g. energy functions, exponential families, but also all sorts of machine learning methods - as we don’t want to span all exponentially many combinations of the discrete space. In a regression tree with categorical inputs, for instance, we don’t use all combinations, but only those that correspond to paths from the root to a leaf in some tree, with $\\theta$ being the piecewise constant expected value of the output at inputs mapped to that leaf. This is a standard interpretation that appears in e.g. Hastie et al.’s *Elements of Statistical Learning*. One interpretation of what we do in Eq.1 is to use a differentiable embedding of the $d$ vector with dimension $r$, instead of $r$ combinatorial paths down a tree of $r$ leaves.\n\nWe are very happy to take follow up clarification questions on the above.\n\n**The assumptions in section 2.2 is not explained intuitively and I don’t know how necessary they are.**\n\nAssumption 1 allows a large enough window of time so that we can identify $\\beta$ first by least squares. For least-squares to be well-posed we need that the corresponding rows in the observation matrix are linearly independent. This is what Assumption 2 states (pseudo-inverses basically mean the matrix that provides the least-squares projection).\n\nAssumption 3 states that, when we need to learn about intervention level $d$, then we need a sufficient number of units getting that assigned, and that we leave these units unperturbed by enough time so that we don’t have further interventions getting conflated with $d$. Here, “enough time” is formalized in terms of the number of intervention parameters $k_d$.\nAssumption 4 is analogous to Assumption 2, now applied to the identification of $\\psi$: having identified $\\beta$ and fixing $\\phi$, this will become another least-squares problem where again we need an explicit assumption about the rank of the matrix of observations. \n\nLet us know whether the above is helpful.\n\n**The experiments are synthetic and semi-synthetic. Although since the nature of the paper is mostly theoretically, I wouldn’t consider that to be a huge issue.**\n\nThat was our reasoning too. We believe that the main point is to call the attention that off-the-shelf black-boxes may not be the wisest option here (a point we want more people to appreciate), so keeping the benchmark to be as uncluttered and controllable as possible was one of our main objectives. We nevertheless provide updated experiments in the shared rebuttal box.\n\n**In equation (1), do you need the lower case x to be a potential outcome?**\n\nAt least implicitly it is necessary, since past $X$s must be exposed to the levels of past treatments as the future $X$s. We are not against other pieces of notation, such as having a single $do(d_n^{1:t})$ operator to indicate a single-world set of interventions. Sometimes it is slightly more convenient to have the potential outcome notation, as in the text of Assumption 4 where we refer to $\\phi$ as explicitly generated at the specific intervention level $d_n^{1:t + t’ - 1}$, but it has some disadvantages too (such as a heavy syntax).\n\n**Is the notation changed in 2.1? f does not appear in equation (1).**\n\nIt is changed to emphasize that in 2.1 we are talking about purely generic function representations - but yes, it is motivated to eventually have the shape of Eq. 1. In an eventual camera-ready version, we will modify the start Section 2.1 to provide a short summary of what will follow to facilitate reading.\n\n**For CSI-VAE, why use GRU instead of more modern arch like attention or even LSTM?**\n\nGRU was actually introduced several years after the LSTM as a way of simplifying it for applications where a full-blown LSTM was an overkill. We considered GRU to be better suited to our benchmarking, as we were not dealing with very long sequences and it was already overfitting - so little to be gained (and, actually, opening up the possibility of doing worse) by falling back into a LSTM, not to mention transformers.\n\nThank you again for your time and feedback!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "ptSYYGsNWQ", "reviewer": "Reviewer_MJr7", "summary": "Estimating the causal effect of a sequence of interventions on another sequence is a central task in causal inference / treatment effect estimation. However, for large discrete spaces, canonical assumptions such as Markovian assumptions, short sequences, and so forth will not apply, while very general black-box models may present poor generalizability due to the sparsity of the observed treatment sequences. This work presents a model where the effect of a sequence of interventions on the target sequence can be decomposed as a product of individual time-point-wise effects. This parameterization is fairly general as it parameterizes functions of a specific form exactly, and any measurable function approximately. Identification results are given for the parameter of the proposed parameterization. This motivates a VAE-based fitting approach, together with predictive intervals coming from conformal prediction. Experiments are given, showing competitive performance of the method.", "questions": "1) I struggle to understand what these mean, can you develop a bit?\n\na) l.38-39 : \"including the bursts of non-stationarity that the application of an intervention brings in\"\n\nb) l.106-108 : \"We do not explicitly condition on $D_{1:t}^n$ in most of our notation, adopting a convention where potential outcome indices $D_{1:t}^n$ are always in agreement with the (implicit) corresponding observed $D_{1:t}^n$\".\n\n2) l.173 : how exactly does $\\mathbb{E}[X_n^t(d_n^{1:t}) | x_n^{1:t-1}(d_n^{1:t-1}), z_n]$ depend on $x^{1:T}$ exactly? Which components of the latter are used?\n\n3) How stringent are Assumptions 3 and 4?\n\n4) Section 3.1 :\na) Is identifiability preserved with variational inference in CSI-VAE?\nb) What is $\\mu$ in the equation of l.250, and how does this incorporate the functional form of Equation 1?", "limitations": "Limitations have a dedicated paragraph, and social impact is no different than any other work on treatment effect estimation more broadly.\n\n**EDIT (2024/08/11)** : increased my score following authors' rebuttal and discussion.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The work is original as it shows that coming with a specialized parameterization can help, both theoretically (identifiabiity) and empirically (better performance), than a general black-box model. It is also significant as the parameterization is actually not particularly restrictive, as Propositions 1 and 2 show that it can recover exactly or approximately general types of functions. Thus the method is generally very principled and should be very applicable in the field.", "weaknesses": "1) IMO the main drawback of the work is in clarity.\n\n1a) It took me significant time to parse what the field is missing and what is the exact contribution of the method. The contribution is specified in the \"Scope\" paragraph l.37-42 and later in l.64-78, and the previous work with its drawbacks is scattered in l.43-63, Section 4 l.304-322 and Appendix A. I think that these sections should be completely re-organized, notably all the previous work should be put together and the contributions outlined clearly in contrast to such previous work ; especially as the contribution consists actually in specifying a more restricted functional form that's more tailored to a specific type of intervention sequence.\n\n1b) It is not clear what the exact estimand is, or what the exact estimands are, as they are not stated in the problem statement l.79-90. The functional form $\\mathbb{E}[X_n^t(d_n^{1:t}) | x_n^{1:t-1}(d_n^{1:t-1}), z_n]$ is also confusing, as it suggests that it does not depend on covariate values but only on $n$, $t$ and interventions, while from l.173 it also depends on covariates. It is only in Section 4 that I understood that I understood that only predictions were to be made conditional on previous covariates. Thus, I would be happy if authors could state in their rebuttal what estimand(s) are we looking for, and if they could write it here.\n\n1c) It is also not clear how the functional form of Equation 1 incorporates \"compositionality\" of interventions. While it can be guessed from Equation 2, it is only in l.184-192 that it becomes clear this functional form for $\\psi$ incorporates compositionality (if I understood it well?)\n\n2) While experiments generally suggest the relevance of the method, I'd have a few concerns :\n\n2a) There seems to be rather few baselines : GRU-0 and GRU-1 seem rather redundant as they look like mutilated versions of the GRU-2 baseline ; the latter seems like the only actually fair comparison against CSI-VAE as only that baseline takes the whole intervention history into account. The bottom part of Figure 2 indirectly confirms this as only GRU-2 is evaluated. Thus, all of this might leave GRU-2 as the only real baseline. Could another one be implemented?\n\n2b) Conformal prediction is incorporated in the paper, but it is only evaluated in the Appendix and only on the method. Could an analysis include baselines and be moved to the main part?\n\nNote : I am aware that it might not easy to obtain new empirical results in the rebuttal period, I am happy to discuss all of this with other reviewers and ACs", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and your remarks that **“The work is original…”** and that is its **“generally very principled and should be very applicable in the field”**! The many comments about presentation are very useful.\n\n**IMO the main drawback of the work is in clarity. 1a.: …time to parse…1b. It is not clear what the exact estimand is…**\n\nThe goal is \"behavioral forecasting under hypothetical future interventions\" (lines 41-42), where \"behaviour\" is defined in line 25.  We want (line 87) to predict $X_{n^\\star}^{T + 1:T + \\Delta}(d_{n^\\star}^{1:T + \\Delta})$ using whatever is observed up to time $T$, extending the past $D_{n^\\star}^{1:T}$ with a hypothetical plan of future actions $D_{n^\\star}^{T + 1:T + \\Delta}$.\n\nWe recognize that, even though this description fits standard machine learning formulations, it is typical for a causal modeling paper to have equations of estimands presented explicitly, even if this is not natural when the problem is predictive. To address this, we propose bringing predictive potential outcomes in terms of functionals earlier on in the write-up.  Under the usual mean-squared error loss, this commonly boils down to the expected value of each $X_{n^\\star}^{T + 1}$, $X_{n^\\star}^{T + 2}$, etc. given the past up to $T$. We originally didn’t want to prescribe this because one may want to use a different loss function other than mean squared error, where plain expectations wouldn’t be the answer. So, it’s not that straightforward to introduce a single equation such as e.g. $\\mathbb E[X_n^{T’}(d_n^{1:T’}) |  x_n^{1:T}(d_n^{1:T}), z_n]$ for $T’ > T$ as the target estimand, because *we are not particularly committed to any given functional of the predictive distribution*: we only want to deal with (potential) observables, not with functionals (that’s why we focus on conformal prediction instead of confidence intervals, and predictive potential outcomes as opposed to cross-world causal effects).\n\nNevertheless, we should anticipate that some readers (like the reviewer) may find it jarring not to have an explicit estimand equation laid out right at the beginning. Given that Eq. 1 is presented just for the predictive mean, we can use this as the family of estimands. So in a revised manuscript we will further emphasize the predictive aspect within the context of means (motivated by mean squared losses) earlier on. \n\n(we have one question for the reviewer: we didn’t quite get what the reviewer meant by “... it suggests that it does not depend on covariate values but only on $n$, $t$,  and interventions, …” - the functional form listed seems to be the same as in l.173?)\n\n**1c) It is also not clear how the functional form of Equation 1 incorporates \"compositionality\" of interventions...  (if I understood it well?)**\n\nThe reviewer got it correctly, but as mentioned it won’t hurt to anticipate early on in the text where we are going with compositionality, including the function composition of Eq. 4.\n\n**Experiments**\n\nWe focused on the GRU because it is a well established method for learning sequential predictions out of categorical input sequences - more flexible than RNNs, less complex than LSTMs or transformers, which would be major overkills: we already overfit with GRUs, a LSTM or beyond would not help (and did not, in preliminary experiments). We nevertheless provide updated experiments in the shared rebuttal box.\n\nGRU-0 and GRU-1 should be seen as ablation studies. GRU-0 quantifies the total contribution of the signal coming from $D$, showing that results are not just artifacts of modeling the evaluation of the $X$ series. GRU-1 excludes long-term histories, focusing only on the latest iteration of intervention, and illustrates that a very strong Markovian structure cannot emulate long-term direct effects. Doing this in the context of GRU shows that even a flexible model that potentially overfits nevertheless benefits from long-term, non-Markovian contributions of past interventions.\n\nRegarding conformal prediction, as the results are done in (semi)synthetic experiments we expect theory to agree with, and so they are in less need to be thoroughly empirically assessed. But we will do our best to include them in the main text in the eventual acceptance of the paper, which would allow us to have an extra page.\n\n**1a) l.38-39**\n\nWhen an intervention is applied, we expect a system to behave in a non-stationary way in the short-term, eventually settling down to a possibly new equilibrium. See for instance Fig. 1 of [9].\n\n**1b) l.106-108**\n\nWe do not allow conditioning of $X^{1:t}$ generated under different $D^{1:t}$ other than the one used in forecasting future trajectories i.e. no cross-world counterfactuals.\n\n**2. l. 173**\n\nThey are used in the definition of $\\phi_{nl}^t$ (see line 109).\n\n**3. Stringency Assumptions 3 and 4.**\n\nAssumption 4 is relatively benign, basically that we don’t have linear dependency between (features of) past trajectories.\n\nAssumption 3 essentially requires that there are enough units of time for a large enough group of individuals to be exposed to a particular treatment level $d$ prior to being exposed to something else. It will depend on how smooth a treatment effect is. Citing again Fig. 1 of [9], realistic effect shapes in many domains can be described with relatively few parameters (effect smoothly going up, then down, then settling), which should be amenable to this assumption holding in practice.\n\n**4. Section 3.**\n\nThe VAE still uses a likelihood (decoder) where the conditional mean is given by Eq. 1, and the conditional variance is homoscedastic and easily shown to be identifiable. Identifiability of the parameters of the encoder is not fundamental (it just represents a posterior distribution, it does not specify a causal structure) and its variance will go to zero as the number of time points increases.\n\nMany thanks again for the useful questions!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:
We take this opportunity to once again thank all reviewers for their time and suggestions! Individualized answers have been provided to each of you.

We use this space to report an update on experiments. To summarize the context, we chose the GRU family merely as an illustration of a modern black-box method that does sequential prediction with representation learning. We didn't think the message would be materially any different if we used a classical RNN, or LSTM, or transformers: the GRU itself was designed to be a good compromise between RNNs and more complex architectures.

In any case, it may be simpler to just run the experiments than argue the above (although we do stress that the main point is structured vs black-box, we never cared whether it was GRU or something else). Since it is straightforward to replace GRUs with other related methods in PyTorch, we did just this. Below are updates including LSTMs and transformers. Autoregressive linear models were also attempted in simpler initial versions of the dataset, but they were already underfitting badly and we do not report them.

*  Fully Synthetic

| Model        | T+1      | T+2      | T+3      | T+4      | T+5      |
| :----------- | :------: | :------: | :------: | :------: | -------: |
| CSI-VAE-1    |  $36.53$ |  $41.46$ |  $41.73$ |  $41.12$ |  $41.32$ |
| CSI-VAE-2    |  $97.80$ | $118.25$ | $117.79$ | $127.25$ | $135.03$ |
| CSI-VAE-3    | $138.78$ | $164.02$ | $141.71$ | $132.59$ | $125.55$ |
| GRU-0        | $229.72$ | $269.66$ | $220.95$ | $208.30$ | $188.43$ |
| GRU-1        | $230.76$ | $270.83$ | $220.93$ | $208.33$ | $184.92$ |
| GRU-2        |  $93.73$ | $101.03$ | $118.01$ |  $88.53$ | $132.28$ |
| LSTM         | $114.71$ | $126.65$ | $137.12$ | $105.22$ | $137.19$ |
| Transformer  | $111.66$ | $122.08$ | $150.57$ | $175.84$ | $87.89$  |

* Semi-Synthetic Spotify

| Model        | T+1      | T+2      | T+3      | T+4      | T+5      |
| :----------- | :------: | :------: | :------: | :------: | -------: |
| CSI-VAE-1    |  $68.23$ |  $82.94$ |  $83.53$ |  $81.97$ | $79.63$  |
| CSI-VAE-2    | $253.85$ | $312.53$ | $305.08$ | $303.68$ | $302.83$ |
| CSI-VAE-3    | $757.94$ | $937.07$ | $800.55$ | $704.66$ | $634.72$ |
|  GRU-0       | $215.42$ | $260.65$ | $193.41$ | $137.20$ | $117.06$ |
| GRU-1        | $223.61$ | $269.69$ | $205.91$ | $141.53$ | $126.36$ |
| GRU-2        | $154.18$ | $187.42$ | $177.96$ | $133.36$ | $127.58$ |
| LSTM         | $130.35$ | $156.02$ | $133.28$ |  $94.35$ |  $85.92$ |
| Transformer  | $133.42$ | $157.66$ | $154.61$ | $164.70$ | $158.03$ |

The attached pdf with box plots provides a visualisation of the above.

---

## Code modifications

Our code will also be updated accordingly (the changes are very localized, basically replacing calls for GRU with other methods already are available in PyTorch).  Changing from GRU to LSTM boils down to commenting a line and adding another:

```
#self.rnn_z_x = nn.GRU(self.z_dim+1, hidden_dim, self.num_layers, batch_first=True)
self.rnn_z_x_d = nn.LSTM(self.z_dim+1+hidden_dim, hidden_dim, self.num_layers, batch_first=True)
```

For transformers, we need to first define a transformer encoder with masking.

```
import torch
import math

class PositionalEncoding(nn.Module):
    def __init__(self, hidden_dim, max_len=5000):
        super().__init__()
        self.hidden_dim = hidden_dim

        pe = torch.zeros(max_len, hidden_dim)
        position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1)
        div_term = torch.exp(torch.arange(0, hidden_dim, 2).float() * (-math.log(10000.0) / hidden_dim))
        pe[:, 0::2] = torch.sin(position * div_term)
        pe[:, 1::2] = torch.cos(position * div_term)
        pe = pe.unsqueeze(0).transpose(0, 1)
        self.register_buffer('pe', pe)

    def forward(self, x):
        # x of size B, T, H
        B, T, _ = x.size()
        # print(x.size())
        # print(self.pe.size())
        pe = torch.permute(self.pe, (1, 0, 2))
        # print(self.pe.size())
        x = x + pe[:, :T, :].expand(B, -1, -1)
        # print(x.size())
        # x += self.pe
        return x

class TransformerTimeSeries(nn.Module):
    def __init__(self, hidden_dim, nhead, num_layers=1):
        super().__init__()
        self.hidden_dim = hidden_dim
        self.pos_encoder = PositionalEncoding(hidden_dim)

        encoder_layers = nn.TransformerEncoderLayer(hidden_dim, nhead, dim_feedforward=hidden_dim*4, batch_first=True)
        self.transformer_encoder = nn.TransformerEncoder(encoder_layers, num_layers)

    def generate_square_subsequent_mask(self, sz):
        mask = (torch.triu(torch.ones(sz, sz)) == 1).transpose(0, 1)
        mask = mask.float().masked_fill(mask == 0, float('-inf')).masked_fill(mask == 1, float(0.0))
        return mask

    def forward(self, x):
        '''
        input:
        x of size (B, T, hidden)
        output:
        x_T of size (B, T, hidden)
        '''
        x = self.pos_encoder(x)
        _, T, _ = x.size()
        tgt_mask = self.generate_square_subsequent_mask(T).to(x.device)
        # print(tgt_mask)

        h = self.transformer_encoder(x, mask=tgt_mask, is_causal=True)

        return h
```
 and then modify the call as in

```
# self.rnn_z_x = nn.GRU(self.z_dim+1, hidden_dim, self.num_layers, batch_first=True) 
self.rnn_z_x_d = TransformerTimeSeries(self.z_dim+1+hidden_dim, nhead=8)
```

Many thanks again!

Meta-review: The submission addresses learning a structural model for sequential treatment regimes where interventions are composed over time. Theoretical results, a VAE based algorithm and experimental validation are provided.

While most reviewers consider the contribution to be sound and valuable for the field, there has been concerns that the manuscript is hard to parse. Overall, the contribution is sound and significant, but highly technical, such that it is hard to follow beyond a community with strong expertise in statistical models for potential outcomes. Given the significance and strengths of the contribution, the AC suggests acceptance of the paper as a poster. However, the AC recommends to make an effort towards the general readership by giving more insights into their model in the camera ready version, e.g. by exposing early in the manuscript an example in a simplified setting.
Decision: Accept (poster)