Structured Learning of Compositional Sequential Interventions

Authors: jialin  yu
andreas  koukorinis
nicolò  colombo
yuchen  zhu
ricardo  silva
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 16:51:12

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.

Main paper: 
 claims: 
 Question: Do the main claims made in the abstract introduction accurately reflect the paper's contributions and scope? Answer: [Yes] Justification: We explain the particular scope in terms of examples we believe to be intuitive, why we see this is a challenge and highlight what to expect from off-the-shelf solutions.We point out to empirical gains on a 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 be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
a structural approach for intervention composition: 
 For simplicity, we will assume scalar behavioral measurements X t n , as the multivariate case readily follows. Consider the following conditional mean model for the regime do(dEQUATIONwhere ⊙ is the elementwise product, with ϕ t n , β n and ψ t n defined below. We postpone explaining the motivation for this structure to Section 2.1. 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(x t n | do(d1:t ′ n )) = p(x t n | do(d 1:t n ))for any t ′ ≥ t (future interventions do not affect the past), and p(z n | do(d 1:t n )) = p(z n ) for all t. We do not explicitly condition on D 1:t n in most of our notation, adopting a convention where intervention indices do(d 1:t ′ n ) are always in agreement with the (implicit) corresponding observed D 1:t n = d 1:t n . We define ϕ t nl := ϕ l (x 1:t-1 n , z n ) as evaluations of basis functions ϕ l (•). For now, we will assume such ϕ l (•) functions to be known. Each unit n has individual-level coefficients β n . Causal impact, as attributed to do(d 1:t n ), is given byEQUATIONwhere function ψ l (d, t ′ , t) captures a trajectory in time for the effect of D t ′ n = d, for t ′ < t. Compositionality will be described in the sequel in two ways: first, as a 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 a baseline unit-level feature vector β n with an interventional-level embedding vector ψ(d t n , t, t) (compositional warping).Effect families. We consider two variants for ψ(•). The first is a time-bounded variant which, for d > 0, defines ψ l (d, t ′ , t) := w dl,t-t ′ , where w dl,t-t ′ is a free parameter for0 ≤ t -t ′ < k d , otherwise w dl,t-t ′ = w dl,k d -1 .The intuition is that the effect of intervention level d has no shape constraints but it cannot affect the past and it must settle to a constant after a chosen time window hyperparameter k d > 0. The default level d = 0 has no effect and no free parameters, with ψ l (0, t ′ , t) = 1. The second variant is motivated by real-world phenomena where causal impacts diminish their influence over time and result in a new equilibrium [9] . We consider a time-unbounded shape with exponential decay, parameterized asEQUATIONwhere σ(•) is the sigmoid function, so that 0 < σ(w 1d ) < 1. As t grows, σ(w 1d ) t-t ′ goes to 0. w 3d is the stationary contribution of an intervention at level d. In particular, for d = 0, we define w 10l = w 20l = 0 and w 30l = 1, while for d > 0 we have that w 1dl , w 2dl and w 3dl are free parameters to be learned. Let us interpret the previous k d hyperparameter as the dimensionality of intervention level d. As the model of Eq. ( 3) is given by vectors w 1d , w 2d , w 3d , we here define k d = 3.Eq. ( 1) is reminiscent of tensor factorization and its uses in causal modeling [7, 3] , 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, it can be motivated from known results in functional analysis (such as Proposition 1 of [25] ) that allows us to represent a function f (x) by first partitioning its input into (x a , x b ) and controlling the approximation error given by an inner product of vector-valued functions,EQUATIONwhich we call the factorization trick. In practice, we choose the dimensionality of such vectors {f a (•), f b (•)} 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 D 1:T n , at a fixed dimensionality T , can be separated from the other inputs asEQUATIONHere, similarly to classical ANOVA models, f d (•) is a one-hot encoding vector for the entire (exponentially sized) space of possible d 1:T n trajectories. That is, f d (•) is binary with exactly one non-zero entry, corresponding to which of the possible trajectories d 1:T n took place. The codomain off nxz (•, •, •) must be a set of vectors of corresponding (exponential in T ) length.In what follows, we will start from the premise of further partitioning f (n,d 1:T , x 1:T , z) := E[X T +1 n | x 1:T n , z n , do(d 1:T n )] by applying the factorization trick of Eq. ( 4) repeatedly, justifying the generality of the right-hand side (RHS) of Eq. ( 1). We will accomplish this by first introducing more structure into the respective f d (d 1:T ), based on a recursive composition of functions. We will then show in Proposition 1 how the RHS of Eq. ( 1) follows. We also provide a complementary theoretical result in Proposition 2 which, by assuming a bounded T , describes the choice of dimensionality in the application of Eq. ( 4) as controlling a formal approximation error for a general class of functions.Main results: 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 fixed-dimensionality embeddings of sequences and identities. In particular, we design f d (d 1:t ) to have a finite and tractable output dimensionality r g for all t. Inspired by sequential hidden representation models such as recurrent neural networks, for each entry l ∈ 1, 2, . . . , r g in the output of f d (•), we define f dl (d 1:t ) via the following recursive definition:EQUATIONThe base case for t ′ = 0 is g l (d 1:0 , 0, t) ≡ 1 for all t, with d 1:0 undefined. Moreover, we define h l (v, 0, t ′ , t) ≡ v for any (v, t ′ , t) -that is, intervention level D t ′ n = 0 does nothing. In order to define h l (v, d, t ′ , t) for d > 0, we will apply the same trick of splitting the input into two subsets, and defining a vector representation on each. The representational choice h l (v, d, t ′ , t) := h l1 (v) T h l2 (d, t ′ , t) requires a choice of dimensionality r h l for the inner vectors h l1 , h l2 . Recursively applying (6) in this product form can still be done relatively efficiently, but choosing r h l = 1 for all l simplifies matters. Finally, we apply the factorization trick one last time to define f nxz (n, x 1:T , z) within Eq. ( 5) to get to a factorization form analogous to Eq. (1). More formally:Proposition 1 Let (i) f (n, d 1:T , x 1:T , z) := f T d (d 1:T )f nxz (n, x 1:T , z), where function sequences f d (d 1:T ) and f nxz (n, x 1:T , z) are defined for all T ∈ N + and have codomain R rg ; (ii) f dl (d 1:T ) be given as in Eq. ( 6) with r h l = 1; (iii) the l-th entry of f nxz (n, x 1:T , z) be given by f nxzl (n, x 1:T , z) := u l (n) T v l (x 1:T , z), where both u l (•) and v l (•, •) have codomain R r s ′ ; (iv) v l (x 1:t , z) be defined analogously to f dl (d 1:t ) as in Eq. ( 6), carrying out the fixed-size z as an extra argument. Then there exists some integer r and three functions a, b and c with codomain R r so that f (n,d 1:T , x 1:T , z) = r l=1 a l (x 1:T , z) × b l (n) × c l (d 1:T , T ). Moreover c l (d 1:T , T ) = T t=1 m l (d t , t, T ) for some function m l : N 3 → 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. (1) , where (a, b, c) plays the role of (ϕ, β, ψ). 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 a fixed value of T , suppose we have a measurable function f :T t=2 X 1:t-1 × T t=1 D t × Z → R d , t ≤ T , where X 1:t-1 and Z are compact, D t is finite. Let x ∈ T t=2 X 1:t-1 and let d ∈ T t=1 D t . For d = 1, • • • , d max , let m d : T t=1 D t → {1, • • • , T } with m d (d) = t if ∃t s.t. d t = d else 0.Then, for any ϵ > 0, there exist r = O(T dmax+1 ) and measurable functions ϕ l : X 1:t-1 × Z → R d and ψ ld ′ : D t → R, and real numbers β l , such thatsup x∈ T t=2 X 1:t-1 ,d∈ T t=1 D t ,z∈Z f (x, d, z) - r l=1 ϕ l (x, z) × β l × dmax d=1 ψ ld ′ (m d ′ (d)) ≤ ϵ. (7)Notes. The theoretical results above are reminiscent of the line of work on matrix factorization (e.g., Appendix A of [4] ). Eq. ( 1) 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 1 and 2 further imply that, for f (n,d 1:T , x 1:T , z) := E[X T n | x 1:T -1 n , z n , do(d 1:T n )], the representation of functions in a multilinear combination of adaptive factors (Eq. ( 1)) is motivated by more fundamental results in functional analysis. Our assumption that no value of d (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 h l2 (d, t ′ , t) leading to Eq. ( 2) are time-translation invariant, being a function of d and t -t ′ only. While it is possible to generalize beyond translation invariant models, as done in the synthetic controls literature [1] that inspired causal matrix factorization ideas [7, 3, 4] , this would complicate matters further by requiring a fourth factor into the summation term of Eq. (1) 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.Interpretation: compositional warping. At the start of the process, where all units are assumed to be at their "natural" state (assuming D 1 n = 0), we can interpret η 1 n := β n as a finite basis representation, with respect toϕ 1 n , of E[X 1 n | Z n , do(d 1 n = 0)] = (ϕ 1 n ) T η 1 n . Each η 1nl can be seen as a latent feature of unit n. Effect vector ψ t n defines a warping of η t-1 n that describes how treatments affect the behaviour of an individual in feature space, with the resulting η t n := η t-1 n ⊙ ψ t n . The modifier ψ t n 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 t n from the chosen basis (notice that intervention level 0 implies η t n = η t-1 n since we define ψ l (0, t ′ , t) = 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 β n and warp function embeddings ψ l (d, t ′ , t), as we will see next.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 toE[X t n | x 1:t-1 n , z n , do(d 1:t n )].We will impose conditions on the realizations of X 1:T n and values chosen for d 1:T n , as well as values of N and T as a function of r, in order to identify each β n and the parameters of ψ l (d, t ′ , t) 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 a prescribed basis ϕ of known and finite dimensionality r, will then provide the foundation for a practical learning algorithm in the sequel.Assumption 1 For a given individual n, we assume there is an initial period T 0n ≥ r with "no interventions" i.e., D t n = 0 for all 1 ≤ t ≤ T 0n .Assumption 2 Let A n be a (T 0n -1) × r matrix, each row t = 1, 2, . . . , T 0n -1 is given by ϕ n (x 1:t n , z n ) T realizations under regime do(d 1:t n ). Let A n be full rank with left pseudo-inverse A + n .The first assumption can be interpreted as allowing for a "burn-in" period for unit n under only the default action. The second assumption ensures enough diversity among realized features ϕ n so that a least-squares projection can be invoked to identify β n , as formalized in the following proposition.Proposition 3 Let b n be a (T 0n -1) × 1 vector with entries E[X t+1 n | x 1:t n , z n , do(d 1:t+1 n )].Under Assumptions 1 and 2, β n is identifiable from A n and b n .The next assumption is the requirement that any action level d of interest is carried out across a large enough number of units. Moreover, they must be sufficiently separated in time from follow-up non-default action levels so 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 ψ l (d, t ′ , t) ̸ = 0 for all l, t, t ′ and that there is at least one time index t, and one set N d containing all units n where d t+1 n = d, such that: (i)N d := |N d | ≥ r; (ii) ∀n ∈ N d , d t+2 n = • • • = d t+k d -1 n = 0, where k d is the dimensionality of the intervention level d. Assumption 4 Let n 1 , n 2 , . . . , n N d index the elements of a set of units N d . For t ′ = 1, 2, . . . , k d -1, let A dt ′ be a N d × r matrix where each row i = 1, . . . , N d is given by ϕ ni (x 1:t+t ′ -1 ni , z ni ) T realized under regime do(d 1:t+t ′ -1 ni). We assume that each A dt ′ is full rank with left pseudo-inverse A + dt ′ . This leads to the main result of this section:Theorem 1 Assume we have a dataset of N units and T time points (d 1:T n , x 1:T n , z n ) generated by a model partially specified by Eq. ( 1). Assume also knowledge of the conditional expectationsE[X t n | x t-1 n , z n , do(d 1:t n )]and basis functions ϕ l (•, •) for all l = 1, 2, . . . , r and all 1 ≤ n ≤ N and 1 ≤ t ≤ T . Then, under Assumptions 1-4 applied to all individuals and all intervention levels d that appear in our dataset, we have that all β n and all ψ l (d, •, •) are identifiable.Notice that nothing above requires prior knowledge of all intervention levels which will exist, and could be applied on a rolling basis as new interventions are invented. There is no need for all units to start synchronously at t = 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 T 0 steps and that each new intervention level d is applied to at least r 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 r. This makes explicit that we 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.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, particularly 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 a 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 well-specification, 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 on a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to 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 fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
algorithm and statistical inference: 
 This section introduces a learning algorithm, as well as ways of quantifying uncertainty in prediction. We also allow for the learning of adaptive basis functions ϕ. As Eq. ( 1) 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 byX t n | x 1:t-1 n , z n , do(d 1:t n ) ∼ N (f (x 1:t-1 n , z n , d 1:t n ), σ 2 ).Here, the conditional mean f (•, •, •) is given by Eq. ( 1). As σ 2 is 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 a finite set of estimating equations, we can parameterize it in the multilinear form analogous to (1) and repeat the analysis of the previous section. In practice, the functionals (such as the conditional expectations of X t ) used in an estimating equation are unknown. Matrix factorization methods can be used directly by first applying a smoothing method to the data to get plug-in estimates of these functionals [42] , but we will adopt instead a likelihood-based approach.We treat each β n as a random effects vector, giving each entry β nl an independent zero-mean Gaussian prior with variance σ 2 β . Along with all hyperparameters, we optimize the (marginal) log-likelihood by gradient-based optimization. With a Gaussian likelihood, the posterior and filtering distribution of each β n can be computed in closed form, but in our implementation we used a black-box amortized variational inference framework [27, 30, 34] anyway that can be readily put together without specialized formulas and is easily adaptable to other likelihoods. We use a meanfield Gaussian approximation with posterior mean and variances produced by a gated recurrent unit (GRU) model [14] composed with a multilayer perceptron (MLP). Therefore, the approximate posterior at t time steps follows from µ q,β,n := MLP(GRU η β,1 (d 1:t n , x 1:t n , z n )), and log σ q,β,n := MLP(GRU η β,2 (d 1:t n , x 1:t n , z n )). Prediction for X t+1:t+∆ n is done by sampling M trajectories, where for each trajectory we first sample a new β n from the mean-field Gaussian approximation. We set each forward Xt+i n to be the corresponding marginal Monte Carlo average. Each basis vector ϕ t n is parameterized via another GRU model such that ϕ t n := MLP(GRU ϕ (x 1:t-1 n , z n )). Although the theory in the previous section relies on known basis functions, we could train this GRU up to a threshold time point of T 0 and freeze it from that point, each β n 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 T 0 for the computation of the variational posterior of β. 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).Conformal prediction (CP) [46, 19] 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 lower-bounded probability. We consider Split Conformal Prediction in two setups.1. Hold-out predictions. We have a set of historical users, n = 1, . . . , N , whose behavior has been observed until time t + ∆. The task is to predict the behaviour at time t + ∆ of a new user, n = N + 1, who has been observed up to time t, i.e. to predict X t+∆ N +1 given X 1:tN +1and D 1:t N +1 . If we assume we have used the history of the new user X 1:t n , n = 1, . . . , N + 1, to train the model, calibration and test samples are exchangeable. ) for n = 1, . . . , N + 1, but in what follows we will drop the explicit do(•) indexing to keep notation lighter, referring explicitly to past observed D 1:t 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 t, P T ∼ (X t+∆ n , X 1:t n , D 1:t n ) may be different from the one before timet, P C ∼ (X t ′ +∆ n , X 1:t ′ n , D 1:t ′ n ), t ′ < t -∆ and ii) we only used X t ′ n , D t ′ n , t ′ = 1, . . . t,for training. CP algorithms are applied on top of given point-prediction models, where we will use X to denote the cross-sectional sample space of any X t . In our case, the underlying model is f , which predicts the expected user behaviour at time t + ∆ given the user history up to time t, i.e.Xt+∆ n := f (D 1:t n , X 1:t n , Z n ) ≈ E[X t+∆ n | D 1:t n , X 1:t n , Z n , do(D 1:t , d t+1 , 0 t+2:t+∆ )]for a chosen implicit d t+1 and with 0 t+2:t+∆ denoting the default action 0 being taken at t + 2, . . . , t + ∆.n = |X t+∆ n -f (D 1:t n , X 1:t n , Z n )|} N n=1 and S N +1 = |X t+∆ N +1 -f (D 1:t N +1 , X 1:t N +1 , Z n )| are exchangeable, i.e.Prob(S 1 , . . . , S N + 1) = Prob(S σ(1) , . . . , S σ(N +1) ) where σ is any permutation of {1, . . . , N + 1}. The Quantile Lemma, e.g. Lemma 1 of [44] , implies the prediction intervalEQUATIONwhere the probability is over the calibration and test samples.Setup 2. The prediction intervals defined in (8) may not be valid, i.e. ( 9) may not hold, because the calibration and test samples,S t ′ n = |X t ′ +∆ n -f (D 1:t ′ n , X 1:t ′ n , Z n )| with t ′ ≤t and t ′ > t, n = 1, . . . , N , are not exchangeable. Theorem 2 provides a bound on the coverage gap, i.e. a measure of violation of (9), under the assumption that the distribution shift is controlled by a perturbation parameter, ϵ > 0.Theorem 2 Assume we have N calibration samples,EQUATIONwhere t n is the time user n experienced the last intervention before t. Assume there exists ϵ > 0 such that, for all n,EQUATIONwhere p T and p C are the (unknown) densities of the test and calibration distributions, p δ is a bounded arbitrary shift density, andEQUATION)In the proof, we use the likelihood-ratio-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 ϵ-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 a single time step and ii) optimized density estimates are a byproduct of our prediction model.Question: For each 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 crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a 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. 
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 on single interventions, our model explicitly addresses the complexity arising from sequential interventions, taking a more detailed perspective on the dynamic interplay of treatments over 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 not been on individual nor average treatment effects, but directly on expected potential outcomes, similarly to the work on causal matrix completion [7, 3, 4, 42] . Unlike the matrix completion literature, we focused on prediction problems, out-ofsample for both units and time steps. Moreover, the matrix completion literature is usually framed in terms of marginal expectations E[X t n (d)], as opposed to conditional expectations (notice that [4] considers covariates, but those are included as part of the generative model). Marginal models have some advantages when modeling multiple-step-ahead effects [18] , but they also involve complex computational considerations that we leave for future work. Compositionality based on additivity instead of the factorization trick is discussed by e.g. [29] , but additivity lacks the connection to an explicit principle for universal function approximation. Finally, there is a rich literature on confounding adjustment where exogeneity of D t n cannot be assumed, see [22] for a textbook treatment. We can rely on standard approaches of sequential ignorability [36] to justify our method in the absence of randomization.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: Data sources and references are properly accounted for, a 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 to reproduce the necessary and important preprocessing steps. 
experiments: 
 We run a number of synthetic and semi-synthetic experiments to evaluate the performance of the CSI-VAE approach. In this section, we summarize our experimental set-up and main results. The code for reproducing all results and figures is available online 2 . In Appendix C, we provide a detailed description of the datasets and models. In Appendix D, we present further analysis and more results. Finally, in Appendix E, we present an illustration of uncertainty quantification results. 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 1 with another constructed holdout set (5, 000 points for the fully-synthetic case, and 1, 000 for the semi-synthetic one).Compared models. We implemented three variations: (1) CSI-VAE-1 follows exactly our setup in Section 2; (2) CSI-VAE-2 can be considered as an ablation study, which relaxes the product form of Eq. ( 1) and replace it with a black-box MLP applied directly to (ϕ t n , β n , ψ t n ) (and hence may not guarantee identifiability); (3) CSI-VAE-3 is another ablation study, where the equation for ψ (Eq. ( 2)) is replaced by a black-box function taking the sequence of actions D as a standard time-series. We compare our model against: (1) GRU-0, a black-box gated recurrent unit (GRU, [13] ) composed with a MLP, using only the past history of X 1:t n and Z n ; (2) GRU-1, another GRU composed with a MLP that takes into account also the latest intervention D t n ; and (3) GRU-2, which uses not only D t n but also the entire past history D 1:t n just like CSI-VAE. In general, GRU-2 can be considered as a very strong and generic black-box baseline model. In addition, we conduct experiments comparing other popular black-box baseline models: (4) LSTM, [23] ; and (5) Transformer, [45] . For those, we use the same input setup as in the case of the GRU-2 model.. Each experiment was repeated 5 times, using Adam [26] at 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 1 and 2 , which show the superiority of our model against strong baselines. We also observed that the identifiability results and compositional interactions of intervention effects are both critical, as evidenced by the drop in performance for CSI-VAE 2 and 3 in Table 1 .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 of r which overshoot the true value -mitigated by regularization); (2) different sizes of the training data (summary: even with more data provided, our model consistently outperforms the black-box models by a significant margin; we show that a generic black-box cannot solve this problem by simply feeding in more data); and (3) a demonstration of 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 are 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 [1] , we assume a 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 a 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 and states is necessary to increase the applicability of the method. Moreover, when each action level d is itself a combination of cross-sectional actions, causal energy-based models such as [8] can be combined with the compositional factorization idea to also generalize to previously unseen cross-sectional combinations of actions. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper poses no such risks beyond any of the usual issues that a generic predictive model may bring. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 a best faith effort. 12. Licenses for existing assets Question: 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? 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 use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL.• The name of the license (e.g., CC-BY should be included for each asset.• For scraped data from a particular source (e.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 package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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 not available online, the authors are encouraged to reach out to the asset's creators. 
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 notebooks.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of 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 is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
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 not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS 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. 
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 crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) 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. 
: 
 1 Contribution D 1 n D 2 n D 3 n D 4 n β n X 1 n X 2 n X 3 n X 4 nn interact with (latent) random effect parameters β n to produce behavior X t n represented as a dense graphical model with square vertices denoting interventions do(d 1:t n ) [33, 16] . Further assumptions will be required for the identifiability of the impact of interventions and their combination, including how temporal impact takes shape and the number of independent units of observation.Many causal inference questions involve a treatment sequence that varies over time. In the (discrete) time-varying scenario, an unit n is exposed to a sequence D 1 n , D 2 n , D 3 n , . . . of actions, or treatments, potentially controllable by an external agent [22, Ch. 19 ] intervening on them. Each D t n can be interpreted as a cause of the subsequent behavior X t n , X t+1 n , . . . , and its modeling amenable to the tools of causal inference [e.g., 35, 36, 31, 17, 11, 49] . For instance, each D t n can take values in the space of particular drug dosages which can be given to an in-patient, or combinations of items to be promoted to a user of a recommender system. Here, patients and users are the units, and X t n denotes a measure of their state or behavior under time t. Dependencies can result in a dense causal graph, as in Fig. 1 . Further background on this classical problem is provided in Appendix A.Scope. We introduce an approach to predict the effects of actions within a combinatorial and sparse space of categorical intervention sequences. That is, for an unit n and given past behavioral sequence X1:t-1 n := [X 1 n , . . . , X t-1 n ], manipulable variables D 1:t-1 n := [D 1 n , . . . , D t-1 n ] and time-fixed covariates Z n , we predict future behavior X t n , X t+1 n , . . . , X t+∆ n under the (partially hypothetical) interventional regime do(d 1:t+∆ n ). The latter denotes controlling the action variables by external intervention [33] . For mean-squared error losses, prediction means learning E[X t+δ n | x 1:t-1 n , z n , do(d 1:t+δ n )], 0 ≤ δ ≤ ∆. We consider the case where sequential ignorability 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 models 1 to map sequences D 1:t n := [D 1 n , D 2 n , . . . , D t n ] to behavior X t n and beyond, sometimes exploiting strong Markovian assumptions or short sequences and small action spaces [11, 43] . An unstructured black-box for sequences, such as recurrent neural networks and their variants [e.g. 23, 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 D 1:T n , 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-treatment covariates [47, 48] , which is of a different nature as losses can be minimized with respect to an existing natural distribution, as opposed to the interventional problem where the action variables are 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 [25, 32, 38, 39] . 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-calculus [33] and their extensions [15] can be used with directed acyclic graphs (DAGs) to infer combinatorial effects, although they often restrict each manipulated variable to directly affect only a small set of variables (e.g. [5] ). 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-sequential, 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, it is not a given that the distribution under a particular regime do(d 1:t n ) 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 t n = d t n , then X t+∆ n will have transition dynamics that can be affected not only by the choice of d t n but also the difference ∆, resulting in dense dependency structures such as the one in Figure 1 . As multiple treatments are applied to the same unit n, there will be a 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 t is available, a large dataset may have no two sequences D 1:T n and D 1:T n ′ taking the same values, even when considering only the unique treatment values applied to n and n ′ , regardless of ordering and time stamps.Problem statement. Let each individual unit n be described by time-invariant features Z n and a pair of time-series (X 1:T n , D 1:T n ). Here, V t1:t2 denotes a slice of a discrete-time time-series of variables V t from t 1 to t 2 , inclusive. That observations for all units stop at time T is not a requirement, but just a way of simplifying presentation. Each D t n is a categorical action variable which, when under external control, represents a potential intervention, encoded by values in N. We assume data is available under a sequentially ignorable selection of actions [33, 22] , either by randomization or consequences of an adjustment. D t n describes a (possibly hypothetical) intervention which causes changes to the distribution of X t:∞ n . The special value D t n = 0 denotes an "idle" or "default" treatment level that can be interpreted as "no intervention". We are given a training set of units(Z n , X 1:T n , D 1:T n ).The goal is to predict future sequence of behavior X n ⋆ for a (training or test) unit n ⋆ under a hypothetical intervention sequence do(d 1:T +∆ n ⋆). This process takes place under the stable unit treatment value assumption (SUTVA) [24] , meaning that action variables D t n only affect unit n. Each intervention on D t n is allowed to have an impact on the whole future series t, t + 1, t + 2, . . . . It can correspond to an instantaneous shock ("give five dollars of credit at time t") or an action that takes place over an extended period of time ("promote this podcast from time t to time t + 5"): we just assume the meaning to be directly encoded within arbitrary category labels 0, 1, 2, 3, . . . . 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 a second time" and so on).In Section 2, we provide an 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 4 before we perform experiments in Section 5, highlighting the shortcomings of black-box alternatives. 
 
 Appendix: Much of traditional causal inference aims at estimating the outcome of a single-shot treatment that is not intertwined within a sequence with uncontrolled intermediate outcomes [24] . However, many causal questions involve sequential treatments. For example, we might be interested in 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 to as time-varying treatments [22] .In many designs, it is assumed that a time-fixed treatment applied to individuals in the population happens at the same time (i.e. there is a 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 at different time steps, potentially resulting in no 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/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) [35, 36] of time-varying treatments is to have d0 = (0, 0, • • • , 0) ("never treated") and d1 = (1, 1, • • • , 1) ("always treated") as the two options in the design space. This allows us to answer what the average expected outcome contrasted between these two options. A fine-grained estimation at every time step will require further levels of contrast for other regimes, such as d = (0, 1, • • • , 1) and d = (1, 0, • • • , 0). 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 a 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.Proposition 1 Let (i) f (n, d 1:T , x 1:T , z) := f T d (d 1:T )f nxz (n, x 1:T , z), where function sequences f d (d 1:T ) and f nxz (n, x 1:T , z) are defined for all T ∈ N + and have codomain R rg ; (ii) f dl (d 1:T ) be given as in Eq. ( 6) with r h l = 1; (iii) the l-th entry of f nxz (n, x 1:T , z) be given by f nxzl (n, x 1:T , z) := u l (n) T v l (x 1:T , z), where both u l (•) and v l (•, •) have codomain R r s ′ ; (iv) v l (x 1:t , z) be defined analogously to f dl (d 1:t ) as in Eq. ( 6), carrying out the fixed-size z as an extra argument. Then there exists some integer r and three functions a, b and c with codomain R r so that f (n,d 1:T , x 1:T , z) = r l=1 a l (x 1:T , z) × b l (n) × c l (d 1:T , T ). Moreover c l (d 1:T , T ) = T t=1 m l (d t , t, T ) for some function m l : N 3 → R.Proof. The proof is tedious and the result intuitive, but they help to formalize the main motivation for Eqs. ( 1) and (3).f T d (d 1:T )f nxz (n, x 1:T , z) = rs l ′ =1 f dl ′ (d 1:T )f nxzl ′ (n, x 1:T , z) = rs l ′ =1 f dl ′ (d 1:T ) r s ′ l ′′ u l ′ l ′′ (n)v l ′ l ′′ (x 1:T , z) = r s ′ l ′′ =1 rs l ′ =1 u l ′ l ′′ (n)v l ′ l ′′ (x 1:T , z)f dl ′ (d 1:T ) = r l=1 a l (x 1:T , z) × b l (n) × c l (d 1:T ), where r = r s × r s ′ and a l (x 1:T , z) := v l ′ l ′′ (x 1:T , z) for l = l ′ × l ′′ ; b l (n) := u l ′ l ′′ (n) for l = l ′ × l ′′ ; and c l (d 1:T ) = f dl ′ (d 1:T ) for l ′ = ⌈l/r s ′ ⌉, where ⌈•⌉ is the ceiling function.Finally, each c l (d 1:T ) is by definition some f dl ′ (d 1:T ). From Eq. ( 6) and the assumptions, f dl ′ (d 1:T ) = g l ′ (d 1:T , T, T ), and g l ′ (d 1:t , t ′ , t) := h l ′ (g l (d 1:t ′ -1 , t ′ -1, t), d t ′ , t ′ , t). For t ′ > 0 and d > 0, we have that h ′ l is given by hl ′ 1 (g l ′ (d 1:t , t ′ -1, t))h l ′ 2 (d t , t ′ , t). For t ′ = 0, h l ′ (v, d, t ′ , t) ≡ 1, and for t ′ > 0 and d = 0, h l ′ (v, d, t ′ , t) ≡ g l ′ (d 1:t ′ -1 , t ′ -1, t). By recursive application, this implies f dl ′ (d 1:T ) = T t ′ =1:d t ′ >0 h l ′ 2 (d t ′ , t ′ , T) and the result follows. □ Proposition 2 Given a fixed value of T , suppose we have a measurable function f :T t=2 X 1:t-1 × T t=1 D t × Z → R d , t ≤ T , where X 1:t-1 and Z are compact, D t is finite. Let x ∈ T t=2 X 1:t-1 and let d ∈ T t=1 D t . For d = 1, • • • , d max , let m d : T t=1 D t → {1, • • • , T } with m d (d) = t if ∃t s.t. d t = d else 0.Then, for any ϵ > 0, there exist r = O(T dmax+1 ) and measurable functions ϕ l : X 1:t-1 × Z → R d and ψ ld ′ : D t → R, and real numbers β l , such thatsup x∈ T t=2 X 1:t-1 ,d∈ T t=1 D t ,z∈Z f (x, d, z) - r l=1 ϕ l (x, z) × β l × dmax d=1 ψ ld ′ (m d ′ (d)) ≤ ϵ. (13) Proof. First, let f (x, d, z) := r l=1 ϕ l (x, z) × β l × dmax d=1 ψ ld ′ (m d ′ (d)).Notice that we can construct an equivalent definition of D t for some fixed t as follows. We for now call the new object Dt :EQUATION)Note that there is a one-to-one correspondence between Dt and D t , and denote the bijection connecting them by g : Dt → D t . Since both sets are finite, we can approximate f with f , if and only if we can approximatef • i × g × i by f • i × g × i.So, without loss of generality, just consider Dt as D t .The set of all sequences of treatments is the union of D t over t = 1, ..., T , for some known T :EQUATION)Note thatEQUATION)EQUATIONtherefore the size of D T is polynomial in T .Fixing the ground truth function f at D ∈ D T gives us a new function T t=1 X 1:t-1 × Z → X . Analogously, fixing f at some D ∈ D T also gives a new function in the same function space. Since D T is finite, we can enumerate its elements. Thus, for everyEQUATION)We can then express the vector field as:EQUATION1 ) . . . . . . . . .EQUATIONBy choosing r = |D T | and appropriate values for β l and ψ ld (D (i) ), we can make W invertible. Now, consider the vector field given by the ground truth function fixed at all values of D ∈ D T :EQUATIONand consider its transformation under W -1 :EQUATION)Since recurrent neural networks are universal approximators [40] , for every δ > 0 we can chooseEQUATION)Since W is a finite matrix, its operator norm is well-defined. Then 1) , z) . . .EQUATIONEQUATION)EQUATIONEQUATION)EQUATIONHowever, Proof. Let b n be a (T 0n -1)×1 vector where each entry t is given by E[X t+1 n | x 1:t n , z n , do(0 1:t+1 )], where 0 1:t+1 denotes all actions kept at the default level of 0 from time 1 to time t + 1. We know from standard least-squares theory, and the given assumptions, that as long as the number T 0 ≥ r of observations is larger than or equal to the dimension r and the matrix is full rank, that the system is (over-)determined andEQUATIONβ n = A + n b. □Theorem 1 Assume we have a dataset of N units and T time points (d 1:T n , x 1:T n , z n ) generated by a model partially specified by Eq. ( 1). Assume also knowledge of the conditional expectations E[X t n | x t-1 n , z n , do(d 1:t n )] and basis functions ϕ l (•, •) for all l = 1, 2, . . . , r and all 1 ≤ n ≤ N and 1 ≤ t ≤ T . Then, under Assumptions 1-4 applied to all individuals and all intervention levels d that appear in our dataset, we have that all β n and all ψ l (d, •, •) are identifiable.Proof. By Assumptions 1 and 2, and Proposition 3, all β n are identifiable.For the next step, letEQUATIONα ∆t dl := ∆t t ′ =1 ψ l (d, 1, t ′ ).Consider now the first time t where an intervention level d > 0 is assigned to someone in the entire dataset. By Assumption 3, there exists a dataset N d with the given properties where, for all n ∈ N d , we haveη t nl = β nl × t t ′ =1 ψ l (d t ′ n , t ′ , t) = β nl , since all interventions prior to t have been at level 0. Let b dt ′ be a N d × 1 vector where each entry i is given by E[X t+t ′ ni | x 1:t+t ′ -1 ni , z ni , do(d 1:t+t ′ ni )], for t ′ = 1, 2, . . . , k d -1.Consider the time-bounded case with free parameters w dl,1 , l = 1, 2, . . . , r. Given Assumption 4, we can solve for η t+1 n = A + d1 b d1 , the vector with entries η t+1 nl . As η t+1 nl = η t nl × ψ l (d, t + 1, t + 1) = η t+1 nl × w dl,1 , with non-zero η nl by Assumption 3, this also identifies w dl,1 . Also by Assumption 3, at time point t+2 we have that η t+1 nl = η t nl ×ψ l (d, t+1, t+2) for all units in N d , which by analogy to the previous paragraph, identifies ψ l (d, t + 1, t + 2) and consequently w dl,2 . As this goes all the way to t + k d -1, this identifies all of ψ l (d,• • • , •).If intervention d level is the time-unbounded model of Eq. (3), by Assumption 3, we also identify α 1 d , α 2 d and α 3 d , as they are defined by functions ψ l (d, t + 1, t + t ′ ) = ψ l (d, 1, t ′ ) for t ′ = 1, 2, 3 that we established as identified by the reasoning in the previous paragraph. Asσ(w 1dl ) × w 2dl - σ(w 1dl ) 3 × w 2dl = α 1 dl -α 3 dl and σ(w 1dl ) × w 2dl -σ(w 1dl ) 2 × w 2dl σ(w 1dl ) 2 × w 2dl -σ(w 1dl ) 3 × w 2dl = α 1 dl -α 2 dl α 2 dl -α 3 dl ,solving for the above results inw 1dl = σ -1 α 2 dl -α 3 dl α 1 dl -α 2 dl , w 2dl = α 1 dl -α 3 dl σ(w 1dl ) -σ(w 1dl ) 3 , w 3dl = α 1 dl -σ(w 1dl ) × w 2dl .So far we have established identification of the ψ(d, •, •) 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 our dataset after time t.Assume the induction hypothesis that t s is the s-th unique time an intervention level d > 0 is assigned to some unit in the dataset, and all intervention levels d ⋆ that appeared up prior to t d have theirTo sample data from this simulator, we take the initial state of X 0 n and F 0 n from the pre-processed data. We use the Z n vector which comes from the real-world user embedding. Then, we iteratively generate the next F n via categorical draws from the softmax output based on user embedding Z n all previous F n and X n . Once a categorical value of F n is drawn, we replace it with its embedding form ϕ(F n ), which is the centroid of its cluster. We use this ϕ(F n ) to generate X n 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 0:T n }, intervention sequence {D 1:T n } and covariates Z n .In this Section, we present more experimental results for our method. This includes: (1) the effect of choosing r;(2) further examples of visualization and empirical results.The effect of different choices of r is presented in Fig. 3 , with r = 5 being the ground truth parameter from which the data is generated. We have also plotted the performance of our method under r = 3 and r = 10. We observed that when r is smaller than the ground truth, the model tend to under-fit the data with a slightly higher MSE error. When r 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 as applying L 1 and L 2 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 r without much concern about what an exact value should be. Using a large r leads to higher variance but in general performs reasonably well compare to knowing the true r (here, r = 5). We also notice that applying regularization techniques, even when we know the right dimension of r, can be beneficial. This leads to a more stable performance (that is, with lower variance), but does not change much of the mean estimate as shown in Fig. 3 . This maybe partially caused by the gradient optimization process when training the deep neural network. Hence, we claim that the main price to be paid for a large choice of r is that this makes data requirements stricter (such as the number of steps prior to the first intervention) in order to guarantee identification. 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: r = 5, T = 20, T 0 = 10, K = 5, N = 50, 000, z dim = 5, = 3. We also generated two further datasets of M = 2, 000, one for calibration and one for testing.We first present the model-based uncertainty quantification at level of α = 0.95, in Fig. 8 . Given estimates of the conditional mean and homoscedastic error variance σ 2 learned from data, we can calculate predictive intervals as (µ -1.96 × σ, µ + 1.96 × σ). For purely model based uncertainty quantification with plug-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. NeurIPS Paper Checklist• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. (d) We recognize that reproducibility may be tricky in some cases, in which case are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Justification: The introduction lines up why the problem is important, while the Conclusions mentions briefly possible misuses to be avoided.• The answer NA means that there is no societal impact of the performed.• If the authors answer NA or No, they should explain why their work has no societal impact or 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 considerations. • The conference expects that many papers will be foundational research and not 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 correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or 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 a system learns from feedback over time, improving the efficiency and accessibility of ML).We include here models which explicitly parameterize causal contrasts in longitudinal interventions, such as the blip models of[36].The code is available at https://github.com/jialin-yu/CSI-VAEhttps://open.spotify.com/https://research.atspotify.com/datasets/https://open.spotify.com/https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.htmlhttps://github.com/joshlk/k-means-constrained 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: WEAKNESSES

* 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?
questions: 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: limitations: see above
paper_summary: paper_summary: This paper studies treatment effect estimation under sequential, discrete interventions.
It is assumed that all treatments are fully independent and impact all future outcomes.
It studies the problem of "causal extrapolation" where some combinations of sequential interventions may not have been observed during training.
Parametric assumptions on the causal mechanisms are made and used for generalising to unseen interventions.
A matrix factorisation approach is developed to fill in the gaps for the unseen interventions.
The approach is empirically tested on synthetic and semi-synthetic data.
strength_weakness: strength_weakness: - 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.
- 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).

**Minor:**

- L21 and L 86: Typo "an unit"
questions: 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?
- L136: What are ANOVA models?
- L123: What does the last sentence in this paragraph mean?
- 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?
- Prop. 3: Does this assume that equation 1 holds?
limitations: 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.
- 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?
paper_summary: paper_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.
strength_weakness: strength_weakness: 1. Section 2.1 is really hard to parse. 
2. The assumptions in section 2.2 is not explained intuitively and I don’t know how necessary they are. 
3. 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.
questions: questions: 1. In equation (1), do you need the lower case x to be a potential outcome? 
2. Is the notation changed in 2.1? f does not appear in equation (1).
3. For CSI-VAE, why use GRU instead of more modern arch like attention or even LSTM?
limitations: limitations: The only limitation is how practical these methods are.
paper_summary: paper_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.
strength_weakness: strength_weakness: 1) IMO the main drawback of the work is in clarity.

1a) 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.

1b) 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.

1c) 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?)

2) While experiments generally suggest the relevance of the method, I'd have a few concerns :

2a) 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?

2b) 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?

Note : 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
questions: questions: 1) I struggle to understand what these mean, can you develop a bit?

a) l.38-39 : "including the bursts of non-stationarity that the application of an intervention brings in"

b) 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$".

2) 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?

3) How stringent are Assumptions 3 and 4?

4) Section 3.1 :
a) Is identifiability preserved with variational inference in CSI-VAE?
b) What is $\mu$ in the equation of l.250, and how does this incorporate the functional form of Equation 1?
limitations: limitations: Limitations have a dedicated paragraph, and social impact is no different than any other work on treatment effect estimation more broadly.

**EDIT (2024/08/11)** : increased my score following authors' rebuttal and discussion.

Comments: title: Thank you: closing comment
Closing comments
Closing comments
Thank you
Answer to Aug 9 answer
Thank you
Answer to Aug 7 rebuttal
Re: Official comment
comment: Thank you for your assessment. As the discussion started from our paper, we will close it, showing our disagreement,

* "Eq. 1 is a reasonable approach" Eq. 1 follows a matrix-factorisation approach that is well-received by the community. We are still not sure where the misunderstanding is coming from, or how the reviewer initially missed the discussion in Section 2.1, or the source behind statements such as "Eq. 1 as an estimator" following our rebuttal. Perhaps in the same way the reviewer is unfamiliar with ANOVA, they are also unfamiliar with the long history behind structures such as Eq. 1. 

*  Predictions under control of a full sequence of actions is a fundamental problem of control theory, dynamic treatment regimes, reinforcement learning, among other fields, including the multi-million user company which gave us feedback on the writing of this paper. It's a literature which we assume our readers feel comfortable with. Once again, we repeat that we do not require the data itself to follow actual randomized trials.

Unfortunately, besides the itemized questions in the review which were responded to in detail above and cross-referenced with the literature and passages in the paper, we do not believe we got any piece of concrete advice that would allow us to address the subjective judgement call of "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". As such, we would like to leave our own judgement that this statement is unsubstantiated.

To end on a positive note, we are mindful that reviewing is a volunteer, time-consuming job, we are still genuinely thankful for the time and engagement of the reviewer.
Thank you for providing additional details and clarifications.

However, I'm afraid that my main concerns aren't addressed:

- After reading the re-reading the relevant sections in the paper and the rebuttal, it is still unclear to me why Eq. 1 is a reasonable approach to model generic sequences. This also mirrors the overall lack of clarity in the submission.
- Now, even if we suppose that Eq. 1 was a good ansatz, I think that the assumption of controlling all treatments is overly restrictive for the approach to have an impact for real problems.

Therefore, I will keep my score.
Thanks again for your time and engagement. We hope we have addressed all questions about clarify and experiments.
Thanks again for your review and we hope we have addressed the question about benchmarks - including the easiness by which further comparisons can be added to our codebase.
Thanks for you reply! I don't have any further question at this point.
Thank you very much for engaging with us!
Many thanks for your new answer. All is clear now. I will increase my score and more firmly recommend admission.
Thank you, this follow-up is super-useful to understand where you (and other potential readers) are coming from! Yes, we used $x_n^{1:t - 1}(d_n^{1:t - 1})$ as notation for a realization of $X_n^{1:t - 1}(d_n^{1:t - 1})$. Although not unsound, we can see how this may invite unnecessary confusion. Even though using $v$ as short-hand notation for the event $V = v$ when following a conditioning bar is standard, the extra piece of notation about the potential outcome index gets in the way and it's not very common. We wouldn't want to use $\mathbb E[X_n^t(d_n^{1:t}) | X_n^{1:t-1}(d_n^{1:t-1}), Z_n]$ (capital letters) as later on in e.g. Proposition 3 we mention realizations, but  $\mathbb E[X_n^t(d_n^{1:t}) | X_n^{1:t-1}(d_n^{1:t-1}) = x_n^{1:t-1}, Z_n = z_n]$ will help to establish a convention about the implicit regime invoked when referring to $x_n^{1:t-1}$, including our references for $\phi_l$ (where we do need to instantiate as $\phi_l(x_n^{1:t-1}, d_n^{1:t-1}, z_n)$ when defining $\phi_{nl}^t$ in line 109). Thank you!

Concerning your point 2, yes those are posterior parameters. They are not independent parameters of the encoder per se, strictly speaking the parameters are the ones in the MLP/GRU pipeline (lines 250-251), which then define the mean and variance of the variational Gaussian distribution in the usual amortized variational inference sense. Line 233 is the likelihood function indeed, where $f(\cdot)$ is the conditional mean as in Eq. 1 and the error variance is just another parameter we can optimize with respect to. Finally, regarding your point c), our theory can be respected by basically freezing $\phi$ after the initial $T_0$ period (in this case, $\beta$ is considered identified with respect to the learnt $\phi$ -- it's not fundamental that a different $\phi$ would imply a different $\beta$, the choice of basis $\phi$ is problem-dependent anyway). In practice though, as we say in line 258, we allow our implementation to just backprop through $\phi$ even after $T_0$ (based on preliminary experiments with different simulations), although this can be easily switched off.

Many thanks again, very useful to get suggestions on how to tweak the notation to increase accessibility. Good stuff!
Many thanks to authors for the extensive rebuttal. My questions have been answered, except the following points :

1) I still do not understand exactly how $E[X_n^t(d_n^{1:t})|x_n^{1:t−1}(d_n^{1:t−1}),z_n]$ in Equation 1 and $\phi_l(x_n^{1:t−1}(d_n^{1:t−1}), z_n)$ in Equation 109 depend on $x_n^{1:t−1}, d_n^{1:t−1}, z_n$, or more precisely what the notation $x_n^{1:t−1}(d_n^{1:t−1})$ means. It suggests that $x_n^{1:t−1}$ is a function of $d_n^{1:t−1}$, while of course also depending on $n$ and $t$. Further, from this interpretation, it is at first unclear whether this function is known a priori (before observing the data) or pre-specified, or it is observed in the data, or it has to be learnt, etc... This is what I mean in the statement raised by your own question "we have one question for the reviewer: we didn’t quite get what the reviewer meant by" : if $x_n^{1:t−1}$ is a function that is pre-specified or known a priori, then indeed $x_n^{1:t−1}(d_n^{1:t−1})$ only depends on $d_n^{1:t−1}$ in the argument, and $n,t$ in the index.

It is later in the paper that I understood that these $x_n^{1:t−1}$ are actually values taken by the (potential) covariates, thus these $x_n^{1:t−1}$ are vectors or scalars but not functions (if I am not mistaken). Thus, it now seems to me that the $x_n^{1:t−1}(d_n^{1:t−1})$ notation for small $x$ is incorrect and should be scraped altogether and replaced with either $X_n^{1:t−1}(d_n^{1:t−1})$ with capita $X$ in statistical quantities, and $x_n^{1:t−1}, d_n^{1:t−1}$ in functions. To be more precise, I understand that one should write

a) $E[X_n^t(d_n^{1:t})|X_n^{1:t−1}(d_n^{1:t−1}) = x_n^{1:t−1}, z_n]$ in Equation 1, which makes it clear that $X_n^{1:t−1}(d_n^{1:t−1})$ is a potential covariate that is indexed at the observed intervention $d_n^{1:t−1}$ and is here taking the observed value $x_n^{1:t−1}$ ;

b)  $\phi_l(x_n^{1:t−1}, d_n^{1:t−1}, z_n)$ in Equation 109, as this is a function as in l.173 and it does not formally rely on random variables, including potential covariates.

Is all of this correct? If not, is $x_n^{1:t−1}$ for a small $x$ actually a function of $d_n^{1:t−1}$?

2) My mistake, I had misread Section 3.1 and missed the point that the mean and conditional variance of l.250-251 are actually the posterior parameters of $\beta_n$ (feel free to correct me if I am wrong here). Here are some questions, mostly for clarity:
a) Are these parameters for the encoder then?
b) By the conditional mean and variance of the decoder, do you refer to Equation 233?
c) Do you also confirm that the VAE departs from the setup of Theorem 1, as notably the basis function is learnt?
Thank you very much for the reply. We hope the following addresses the remaining follow-up questions.

**Eq. 1: So are you saying that the RHS of Eq. 1 is an estimator for the LHS? ...**

Apologies, but we are genuinely confused where this conclusion is coming from. We don't see anything in our rebuttal that suggests that, and it would be helpful to have the quote of the passage leading to it. Eq. 1 is indeed a "a statement about the true data generating process". 

**Fig. 1: ...in your rebuttal you say that confounding between action variables is allowed, but not preceding action variables causing future actions...**

Again, we don't see anywhere in our rebuttal a statement that $D$ variables cannot be causing others. In fact, we explicitly say the opposite ("*...nothing in our results change at all if ... past $D$ variables pointing to all future $D$ variables.*"). Perhaps it's useful to recall the difference between interventions and random variables within our context.

In a fully connected DAG according to ordering $(D_1, X_1, D_2, X_2)$ where $D_1$ and $D_2$ are random instead of controlled, we have the edges 

$D_1 \rightarrow X_1, D_1 \rightarrow X_2, D_1 \rightarrow D_2$, $X_1 \rightarrow D_2$, $X_1 \rightarrow X_2$, $D_2 \rightarrow X_2$.

When $D_1$ is controlled to $d_1$ and $D_2$ to $d_2$, one graphical characterization is

$d_1 \rightarrow X_1, d_1 \rightarrow X_2$,  $X_1 \rightarrow X_2$, $d_2 \rightarrow X_2$,

where lower case here indicates that we are talking about exogenous fixed indices (squares in Fig. 1). 

There isn't anything else to be said if $d_2$ is functionally independent of the past, but even then this is not at all an issue. This is because whether $d_2$ is a function or not of $(d_1, X_1)$, it won't affect identification strategies such as the sequential back-door/g-formula see e.g. chapter 4 of [32]. We can still carry $d_2$ symbolically into any averaging over the past, even if $d_2$ is functionally related to $(D_1, X_1)$ and we treat $D_1$ as uncontrolled and average over it (since we don't average over the past, this point is redundant anyway). Other representations such as SWIGs suggest adding edges between fixed indices with functional dependencies see e.g. Chapter 19 of [22], but if we were to adopt SWIGs in Fig. 1 the diagram would become an incomprehensible mess. Technically speaking, even edges like $d_1 \rightarrow X_1$ are "unnecessary", as the graphical model is meant to represent the independence structure of a distribution over random variables, and $d_1$ isn't one (preserving "$d_1 \rightarrow \dots$") saves us from having to label the nodes as $X_1(d_1)$ etc.). 

To summarize, our Fig. 1 is a cosmetic choice among other plausible choices and exemplifies already the case of non-dynamic regimes with no ambiguity. There wasn't really much of a deep point we were trying to make (other than not requiring any Markovian assumptions connecting past and future), and we are earnestly surprised that this is raising a discussion...

**Why did you make the assumption of randomization in the first place?**

This was just a way of saying that we assume to have access to the distribution of (single-world) potential outcomes, where our predictions lie. Whether we obtain it by (say) controlled experiments, sequential ignorability, proxies, instrumental variables etc. is orthogonal to our main results. We thought readers would appreciate if we focused on the main novel aspects of our contribution. We are happy to make this point more explicitly in the introduction.

We really appreciate your engagement, and we hope the above has been helpful.
Thank you for giving a detailed rebuttal. I have read the other reviews, rebuttals and had another look at the paper and will try to answer to this rebuttal below.

## Eq. 1

So are you saying that the RHS of Eq. 1 is an estimator for the LHS?
The equal sign would suggest that this is a statement about the true data generating process rather than a way to estimate the expectation.


## Fig. 1

I'm afraid the explanation did not resolve the confusion.

In L39 you say “We consider the case where each $D_n^t$ behaves as if it was randomized…”, a similar quote in L81, which you provide.
But in your rebuttal you say that confounding between action variables is allowed, but not preceding action variables causing future actions (and I suppose no confounding between actions and outcomes).

Why did you make the assumption of randomization in the first place? In other words, which results in the paper do not hold when you do not make the assumption of randomization?