Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,653 characters · 14 sections · 36 citation commands
Estimating the Long-Term Effects of Novel Treatments: The Dynamically Adjusted Surrogate Index
\icmltitle{Estimating the Long-Term Effects of Novel Treatments:\\ The Dynamically Adjusted Surrogate Index}
\icmlsetsymbol{equal}{*}
\icmlaffiliation{msr}{Microsoft Research, New England}
\icmlcorrespondingauthor{Vasilis Syrgkanis}{[email removed]}
\icmlkeywords{long-term effects, dynamic treatment effects, surrogates, high-dimensional, double machine learning}
\vskip 0.3in
\printAffiliationsAndNotice
Businesses frequently invent new ways of interacting with their customers. Marketing departments frequently devise new marketing campaigns. Pharmaceutical companies typically roll out trials of new drugs. In a multitude of domains, policy makers want to understand the long-term effects of novel treatments, recently deployed, while only having access to short-term data following their deployment. Such policy makers only have access to long-term historical data where other treatments had been deployed.
We propose an estimation methodology that leverages historical data to derive estimates of the long-term treatment effects of novel treatments. The seminal paper of athey2020estimating proposes a surrogate index methodology as one solution to this problem. The fundamental assumption of the surrogate index is that there exist short-term proxies that are observed in the short-term dataset and that causal effects on long-term outcomes have to be primarily channeled through these short-term signals. In other words, the treatment has a long-term causal effect, if and only if it has an effect on a variety of short terms signals. This allows one to use the historical long-term data set to learn a mapping from short-term signals to a projected long-term reward --- referred to as the surrogate index --- and subsequently, estimate the causal effect of novel treatments on the surrogate index.
But often historical treatment policies are dynamic: treatments are assigned repeatedly, and their assignments depend on past treatments and short-term outcomes. This can easily break the assumptions needed for the surrogate index to work and introduce bias. For example, suppose that a firm offers multiple investments to a particular customer in the historical data and these investments are auto-correlated, i.e. if a customer receives an investment this month, then they will receive an investment with high probability in one of the subsequent months. These future investments can substantially increase the long-term outcome of interest, and this increase will be attributed to the short-term proxies. The surrogate index thus formed will tend to over-predict long-run outcomes, and so when it is used to measure the treatment effect of some new treatment in the short-run data, the estimated treatment effect will be bigger (in absolute magnitude) than the truth.
The main methodological innovation of this paper is to suggest using the dynamic treatment effect analysis of lewis2020doubledebiased on the historical data in order to create an unbiased dynamically adjusted surrogate index. Our dynamically adjusted surrogate index takes the interpretation of the projected long-term reward in the absence of any future treatments. Applying this dynamically adjusted surrogate model to the short-term dataset leads to unbiased causal effect estimates of the long-term effects of the novel treatments.
A second contribution of the paper is to generalize the causal analysis step that uses the experimental sample to allow for multiple continuous treatments rather than a single binary treatment. We do this by proposing a new estimator for this expanded surrogate approach which allows for the construction of valid confidence intervals, even when using flexible machine learning models at both stages of the estimation to deal with high-dimensional data. In short, we show how one can combine three recently developed techniques, i.e. i) the surrogate index approach of athey2020estimating, ii) the double machine learning approach of chernozhukov2018 and iii) the dynamic treatment effect estimation approach of lewis2020doubledebiased, in a single data analysis pipeline to estimate treatment effects in the presence of dynamic treatment policies.
Our work lies in the broader field of estimating causal effects with machine learning and Neyman orthogonality Neyman:1979,robinson:88,Ai2003,Chernozhukov2016locally,chernozhukov2018double. Moreover, it relates to the work on machine learning estimation of treatment effects in the dynamic treatment regime nie2019learning,Thomas2016,petersen2014targeted,kallus2019efficiently,kallus2019double,lewis2020doubledebiased,bodoryevaluating,singh2020kernel and on structural nested models in biostatistics robins1986new,robins1992g,robins1994correcting,robins1997toward,robins2000marginal,lok2012impact,vansteelandt2014structural,Vansteelandt2016. Finally, it relates to the surrogacy literature in causal inference prentice1989surrogate,begg2000use,frangakis2002principal,freedman1992statistical. Our work builds on insights in these works and proposes the first complete method that combines all three lines of work, so as to estimate the long-term effect of novel treatments from observational data that stem from a dynamic observational policy and in a manner that allows for high-dimensionality.
Though our methodology applies in many domains, for concreteness we consider the running example of a firm making investments in its customers in order to increase subsequent purchases by those customers.
A firm has a number of distinct investments/treatments $T_1, T_2 \ldots T_k$ it offers to its customers. At each period $t$ (e.g. the period of a month), a vector of treatments $T_{i,t}=(T_{i, t,1}, \ldots, T_{i,t,k})\in \ensuremath{\mathbb{R}}^k$, is applied to each customer $i$. We also observe a vector of $p$ characteristics $X_{i, t}=(X_{i, t, 1}, \ldots, X_{i,t,p})\in \ensuremath{\mathbb{R}}^p$, some of which are constant within customer (e.g. industry) and some of which vary over time and customer (e.g. last month's revenue). We observe the outcome of interest $Y_{i, t}\in \ensuremath{\mathbb{R}}$ (e.g. monthly revenue).
We are interested in identifying the average effect of each treatment at some period $t$, on the cumulative outcome in the subsequent $M$ periods, i.e. $\bar{Y}_{i, t} = \sum_{\kappa=0}^{M} Y_{i,t + \kappa}$. We assume that for most customers and periods we do not observe the subsequent $M$ periods following a treatment. This makes direct inference from the historical data impractical. Instead, we assume that for every customer $i$ and period $t$, we have access to a vector of $d$ short-term proxies/surrogates, $S_{i, t}=(S_{i, t, 1}, \ldots, S_{i, t, d})$. In practice, this can for instance include the next few months of revenue and other measures that are indicative of a customer's trajectory.
Given these surrogate measures, athey2020estimating propose the following estimation strategy:
This standard approach will be consistent for the treatment effects under two assumptions. The first is that the relationship between the surrogate variables and the target outcome remains unchanged over our historical sample period and the periods into which we project future outcomes. We follow athey2020estimating in maintaining this assumption. The model is able to account for steady growth in a natural way: as the surrogates grow, the predicted outcomes grow too.\footnote{It is not robust to changes in the mapping between the surrogates and outcomes over time, and so we would recommend refitting the model periodically to account for this.}
The second main assumption is that the only causal path from the treatment to the outcome goes through the surrogates. In other words, the treatment has an effect on $M$-period outcomes if and only if it affects these surrogates. This critical “sufficiency” assumption typically requires that a few periods of post-treatment data are available to build a relatively good understanding of the $M$-period path.
Serially correlated treatment policies will often violate this second assumption. A customer who receives a treatment in period $t$ may be more likely to receive a second treatment over the next year, even controlling for observed features $X_{it}$. This positive correlation could be either because a customer that is treated subsequently becomes more salient to the firm, or because the treatment successfully drives surrogate metrics higher, thereby improving the perceived return on subsequent treatments of this customer (an example could be a firm that targets its fast-growing customers).
In either case, the standard estimated surrogate index model $\hat{g}$ will now be biased; we would wrongly estimate that even small increases in our proxy measures forecast strong growth in $M$-period outcomes on average. If we subsequently estimate the causal effect of novel treatments on the surrogate index, we will over-estimate the causal effect on $M$-period outcomes. Essentially, causal effect that should be attributed to the future treatment gets doubly attributed to the current treatment. To remove this bias we include a dynamic adjustment step based on lewis2020doubledebiased that alters our observational data-set to remove the effect of future treatments from the target outcome before estimating our surrogate index model. To show this why this method is both necessary and sufficient for removing bias, we frame the problem and solution in approach in a two-period example before formally presenting the general case in section (ref).
Assume that $\bar{Y}_{i,t} = Y_{i,t} + Y_{i, t+1}$ (i.e. the long-term outcome contains two periods). We assume a linear Markovian model of how the random variables evolve over time. We collapse $X_{i,t}$ and $S_{i,t-1}$ so that the surrogates at each period and the controls in the next period are all denoted by $S_{i,t-1}$. In other words, all customer observable characteristics in the next period are candidate surrogates and become controls for the period after next. This is without loss of generality. Finally, we focus on a single scalar investment.
With these simplifications, the structural model that describes how all the random variables evolve can be written in three equations:
where $A$ is an $p\times 1$ matrix, $B$ is a $p\times p$ matrix, $\gamma, \lambda$ are $p$-dimensional vectors and $\kappa$ is a scalar. The terms $\epsilon_{i,t}, \zeta_{i,t}, \eta_{i,t}$ are exogenous mean-zero independent noise terms. For conciseness we drop the customer index $i$ in the remaining equations.
\paragraph{Target effect estimand} The effect of treatment $T_{t}$ on the long-term outcome $\bar{Y}_{t}$ can be derived as follows:
Thus we see that the effect of $T_t$ on $Y_{t+1} + Y_{t}$, keeping all other variables fixed, is:
A one-unit increase in investment $T_t$, leads to a total of $\theta_0$ units more revenue in the next two periods assuming future treatments are held constant.
\paragraph{Surrogate index without dynamic adjustment} If we train a surrogate index by regressing $\bar{Y}_{t}=Y_{t} + Y_{t+1}$ on $S_{t}$, as is the method proposed in athey2020estimating, then this would result in the following surrogate index:
Subsequently, if we estimate the causal effect of $T_t$ on $g_0(S_{t})$, then this effect would be:
The standard estimate contains a bias stemming from the fact that investment today can lead to higher investment tomorrow, either directly (i.e. that $\kappa \geq 0$) or indirectly through the surrogates (i.e. that $\lambda'A \geq 0$). The standard surrogate approach is valid only when the investment policy is not adaptive, i.e. $\kappa=0$ and $\lambda=0$, since then $T_{t+1}$ would be independent of $S_{t}$. We see here that the bias is larger if the investment policy is highly auto-correlated, which is typical in practice. The two possible channels for bias, through direct auto-correlation or surrogate-dependent treatment policies, can result in substantially biased estimates.
\paragraph{Dynamic adjustment} The goal of the dynamic adjustment is to remove the effect of the next-period treatment from the long-term outcome $\bar{Y}_{i, t}$. We achieve this by estimating a separate causal effect of $T_{i,t+1}$ on $Y_{i,t+1}$, controlling for $X_{i,t+1}, S_{i,t}$. Observe that this conditional expectation is equal to:
We then subtract the causal effect, $\alpha_{t+1}$, from $Y_{i, t+1}$ to get a dynamically adjusted outcome, $Y_{i, t+1}^{\ensuremath{\text{adj}}} := Y_{i, t+1} - \alpha_{t+1}' T_{t+1}$. Finally, we can create an adjusted long-term outcome, $\bar{Y}_{i, t}^{\ensuremath{\text{adj}}} := Y_{i, t} + Y_{i, t+1}^{\ensuremath{\text{adj}}}$.
We then build a new dynamically adjusted surrogate index, which is the projected adjusted long-term outcome, conditional on the observed surrogates:
This new index captures the projected $M$-period outcomes as if the customer was offered no future treatments.
When we estimate the effect of the treatment based on this adjusted surrogate index we recover:
With the dynamic adjustment, the coefficient in front of $T_t$ that we recover is the true causal effect $\theta_0=\gamma'(I+B)A$. lewis2020doubledebiased show how this adjustment approach can be extended to many periods, via a recursive peeling process, and also to high-dimensional surrogates and controls via a dynamic double machine learning approach.
\paragraph{Estimating causal effects of new treatments} With this adjusted surrogate index in hand, we can also estimate the long-term effect of any other treatment that was introduced more recently. In this example, the stationarity assumption of both proposed surrogate approaches requires that $B$, which governs how the surrogate evolves, and $\gamma$, which governs how surrogates translate to per-period outcomes, do not change between the observational and experimental data-sets. These two parameters govern how surrogates today relate to future outcomes in the absence of any treatment.
Under such a condition, if we introduce a new treatment $T_t^{\ensuremath{\text{new}}}$, which has a different effect $A^{\ensuremath{\text{new}}}$ on the surrogates (and hence on the long-term outcome), then the effect of this treatment on the long-term outcome is $\theta_0^{\ensuremath{\text{new}}}=\gamma'(I+B)A^{\ensuremath{\text{new}}}$. This $\theta_0^{\ensuremath{\text{new}}}$ is exactly the outcome of estimating the causal effect of $T_t^{\ensuremath{\text{new}}}$ on $g_0^{\ensuremath{\text{adj}}}(S_{t+1})$ controlling for $S_t$.
In this section we present the general problem formulation and in the subsequent sections we present our formal main results. There are a number of innovations beyond the basic strategy presented in the two period example above. First, we develop a generalization of the doubly robust estimation method in athey2020estimating to the case of multiple continuous treatments, under a semi-parametric assumption (c.f. Section (ref)) and in the presence of a dynamic treatment policy in the observational data (c.f. Section (ref)). Second, we make use of orthogonal machine learning techniques chernozhukov2018 throughout, to allow for a rich set of potential confounders and valid analytic confidence intervals.
We assume that we have access to two sample populations: an experimental population, denoted as $e$, and an observational population, denoted as $o$. A sample from each population consists of a finite horizon time-series $(S_0, T_1, S_1, Y_1, T_2, S_2, Y_2, \ldots, T_M, S_M, Y_M)$. We observe the full $M$-period time series for each sample in the observational population, but we only observe $(S_0, T_1, S_1)$ for each sample from the experimental population. \vsedit{Moreover, the random variables in the two populations could be distributed differently and even have different support (e.g. treatments in the population $e$ could be different from treatments in population $o$).} As in the two period example, we simplify without loss of generality by merging the control variables in period $t$ and the surrogates for period $t-1$. In other words, all next-period control variables serve as potential surrogates and vice versa. We assume that the data obey the Markovian assumptions depicted in the causal graph in Figure (ref).
Our goal is to estimate the causal effect of treatment vector $T_1$ on the long term outcome:
in the experimental/short-term sample, \vsedit{assuming all future treatments take some baseline value. For simplicity we use the $0$ value as the baseline, but this could be replaced by any baseline treatment vector value. In other words, if we set the future treatments $T_{>1}:=(T_1,\ldots, T_M)$ that each sample receives to the baseline level of $0$ and we change the treatment $T_1$ from some value $t_0$ to some other value $t_1$, then what is the change on the long term outcome $\bar{Y}$, i.e.:
}
We present our theoretical results in two steps. In the first setting, we assume that treatments happen only at period $1$ (Section (ref)). This is the setting analyzed in athey2020estimating, albeit only for the case of a single binary treatment $T_1$.\footnote{We note that the work of athey2020estimating also allowed for estimation of average treatment effects, even in the case when there is arbitrary treatment effect heterogeneity. In this work, we assume that treatment effects are constant. A generalization to the case of arbitrary treatment effect heterogeneity is feasible, but would require the estimation of conditional covariance matrices, which would make the estimation algorithm more brittle and the exposition much more complex.} We then show how this approach can be modified to incorporate a dynamic treatment policy in the observational and experimental sample (Section (ref)).
\paragraph{Notation} Throughout, we will denote with $\ensuremath{\mathbb{E}}_e[\cdot]$ the expectation conditional on the experimental population and $\ensuremath{\mathbb{E}}_o[\cdot]$ the expectation conditional on the observational population. Moreover, for any vector-valued function $f$ that takes as input a random variable $Z$, we denote with:
and analogously $\|f\|_{2,e}$. We denote with $\ensuremath{\mathbb{E}}_n[\cdot]$, the empirical expectation over all the samples, i.e. for any random variable $Z$, $\ensuremath{\mathbb{E}}_n[Z]:=\frac{1}{n}\sum_i Z_i$, and with $\ensuremath{\mathbb{E}}_{e,n}$ and $\ensuremath{\mathbb{E}}_{o,n}$ the empirical expectation over the experimental and observational samples correspondingly, i.e. $\ensuremath{\mathbb{E}}_{e,n}[Z] = \frac{1}{n_e} \sum_{i\in e} Z_i$ and $\ensuremath{\mathbb{E}}_{o,n}[Z]=\frac{1}{n_o} \sum_{i\in o} Z_i$.
For expository purposes, we begin by analyzing the setting where \vsedit{$T_{>1}=0$ almost surely, i.e. treatments occur only in period $1$, but are multi-dimensional and potentially continuous.} We will further assume a partially linear relationship between the treatment $T_1$ and the long-term outcome in the experimental sample:
for some known feature map $\phi(\cdot, \cdot)$, but arbitrary function $b_0(\cdot)$.
Formally, the invariance of the surrogate-outcome relationship requires that that the mean-relationship between the surrogates $S_1$ and the long-term outcome does not change between the observational and the experimental sample:
We denote with $g_0$ the surrogate index model and with $g_0(S_1)$ the surrogate index.
\vsedit{Finally, for the surrogate approach to be valid we need that the long-term outcome $\bar{Y}$ is independent of $S_0, T_1$, conditional on $S_1$. In fact we simply need conditional mean independence:
i.e. $\ensuremath{\mathbb{E}}[\bar{Y} \mid S_1, T_1, S_0] = \ensuremath{\mathbb{E}}[\bar{Y}\mid S_1]$. The latter is satisfied under the causal graph assumption of Figure (ref), when $T_{>0}=0$ a.s..}
\vsedit{Under the PLR assumption and the causal graph governing our data, we have by a standard $g$-formula (see e.g. hernan2010causal) that:
} To estimate our treatment effect of interest, it suffices to find an estimate $\hat{\theta}$ of the parameter vector $\theta_0$. Subsequently, we can also estimate:
We establish three valid estimands for $\theta_0$ that follow similar intuition as athey2020estimating, adapted to consider linear effects of continuous treatments rather than a single binary treatment. A graphical depiction of the different identification strategies is depicted in Figure (ref) \vsedit{
}
The core estimation challenge that the surrogate approach resolves is that the treatments and outcome of interest are not observed in a single dataset. Intuitively, the first surrogate index representation approaches this challenge by using realized treatments from the experimental sample and, in place of realized outcomes, substitutes the expected outcome conditional on the surrogates, $g_0\left(S\right)$, which can be identified from the observational sample and then constructed in the experimental sample.
The second surrogate score representation reverses this substitution. The second term pairs an expectation of the featurized treatment conditional on surrogates, $\ensuremath{\mathbb{E}}_e[\tilde{\Phi}_1\mid S_{1}]$, with the realized outcomes from the observational sample. This representation requires an added ratio of probabilities of appearing in each sample, $\Pr(e\mid S_1)$ and $\Pr(o\mid S_1)$, to adjust for variation of the marginal surrogate distribution across the two datasets.
The third orthogonal representation blends the first two representations and satisfies Neyman orthogonality, which allows the construction of confidence intervals and double robustness.\footnote{One practical difficulty with this third doubly robust approach is that it is less transparent and requires access to the raw historical dataset whenever estimating a new treatment option. In contrast, the surrogate index representation allows for segmentation: one can estimate the surrogate index in the observational data once and store only the $g_0$ model. Treatment effects can then be estimated using only these stored parameters and the experimental dataset, or multiple experimental datasets. This explicit construction of expected outcomes in the experimental data also makes the first approach particularly easy to interpret. However, one must then be careful to account for the additional uncertainty stemming from estimating outcomes in the first step, as standard confidence intervals in the second step will not account for this pre-estimated component.}
We note that the parameter identified in the equations in Theorem (ref) is interpretable even if the partially linear assumption is violated. In this case, the equations are identifying the best linear projection of the variation in the long-term outcome that is not explained by the initial state, i.e. $\bar{Y} - \ensuremath{\mathbb{E}}_e[\bar{Y}\mid S_0]$, on the variation in the feature map $\phi(T_1, S_0)$, that is also un-explained by the initial state, i.e. $\tilde{T}_1$. That is the quantity: \vsedit{
} We formulate the estimation of $\theta_0$ based on the orthogonal representation as a $Z$-estimator based on a vector of moment equations that depends on a vector of nuisance functions $f_0$, i.e.:
and such that it satisfies the Neyman orthogonality condition:
Subsequently, this will allow us to invoke the results in chernozhukov2018double, to derive an asymptotic normal estimator, even when high-dimensional, regularized approaches are used to estimate the nuisance functions $f_0$.
\vsedit{
}
\vsedit{ Given the latter orthogonal moment formulation of the target parameter of interest, one can achieve a root-$n$ asymptotically normal estimate and accompanied asymptotically valid confidence intervals by invoking the results in chernozhukov2018double and verifying that the general conditions required by the main theorems in chernozhukov2018double are satisfied. Given that the estimate that we present in the next section (see Theorem (ref)) is a generalization of the setting presented in this section, we omit this result and refer the reader to the more general theorem of the next section. } \vsdelete{
Asymptotically valid confidence intervals can be constructed using empirical analogues of the co-variance matrix. }
In this section we deal with the case where the treatment policy in the observational and the experimental data is dynamic and we want to estimate only the effects of the treatment at period $1$, under zero future treatments, i.e. the part of the effect that does not go through future treatments but solely through the surrogates/control variables.
\paragraph{Preliminary definitions.} To present the identification and estimation strategy we will need to introduce some notation from the dynamic treatment regime literature. Consider an arbitrary time-series process $\{S_{t-1}, T_t, Y_t\}_{t=1}^{M}$, with $S_t\in \ensuremath{{\cal S}}_t$ and $T_t\in \ensuremath{{\cal T}}_t$. For any time $t$, let $\bar{S}_t=\{S_1,\ldots, S_t\}$ and $\bar{T}_t=\{T_1,\ldots, T_t\}$ denote the sequence of the variables up until time $t$ and similarly, let $\underline{S}_t = \{S_t, \ldots, S_M\}$ and $\underline{T}_t=\{T_t,\ldots, T_M\}$. We will also denote with $\bar{s}_t, \bar{\tau}_t, \bar{y}_t, \underline{s}_t, \underline{\tau}_t, \underline{y}_t$, corresponding realizations of the latter random sequences. Moreover, we will be denoting with $(\bar{\tau}'_t, \underline{\tau}_{t+1})$, the sequences of treatments that follows $\tau'$ up until time $t$ and then continues with $\tau$. We let $0\in \ensuremath{{\cal T}}_t$ denote a baseline policy value, which could be appropriately instantiated based on the context.
\paragraph{Target quantity.} For any sequence of treatment $\tau=(\tau_1,\ldots, \tau_M)$, let $Y_{t}^{(\tau)}$ denote the counterfactual outcome at period $t$ under such a sequence of interventions, equivalently in do-calculus notation $Y_{t} \mid do(\bar{T}_M=\bar{\tau}_M)$. Note that $Y_t^{(\tau)}$ is only a function of $\bar{\tau}_t$, i.e. $Y_t^{(\tau)} \equiv Y_t^{(\bar{\tau}_t)}$, but for simplicity of notation we use the overall vector of treatments. We will also denote with $\bar{Y}_t^{(\tau)}:=\sum_{j=t}^M Y_j^{(\tau)}$, the counterfactual cumulative outcome from period $t$ and onwards, and with $\bar{Y}^{(\tau)} = \sum_{j=1}^M Y_j^{(\tau)}$, the total counterfactual cumulative outcome. Under this counterfactual notation, we can re-write our target quantity of interest from Equation (ref) as:
We show that the target quantity of interest is non-parametrically identified if the data generating processes adhere to the causal graph depicted in Figure (ref) and satisfy a regularity condition on overlap, as well a a dynamic analogue of an invariance relationship between the observational and experimental setting. We first present a set of high-level conditions that lead to non-parametric identification and then present the main identification result.
We assume that the data generating process satisfies the following sequential conditional randomization assumption:
This condition is for instance satisfied if the data generating process adheres to the causal graph presented in Figure (ref), as can be easily verified from the single-world-intervention graph (SWIG) in Figure (ref). Moreover, we will assume a surrogacy assumption, that under a zero future treatment policy, the effect of $T_t$ on future outcomes only goes through $S_t$. This is again satisfied if the data generating process adheres to the causal graph presented in Figure (ref), as can be easily verified from the single-world-intervention graph (SWIG) in Figure (ref). In fact, we will only require a conditional mean-independency assumption.
Since we do not observe long-term outcomes from the experimental setting, we will need to assume a dynamically adjusted analogue of the invariance property, so that we can use long-term outcomes from the observational dataset to “impute” long-term outcomes in the experimental dataset.
Observe that the dynamic invariance Assumption (ref) is much more permissive in practice than the standard invariance assumption as we no longer require that the dynamic treatment policy in the observational data be the same as in the experimental data, but simply that the adjusted outcomes under baseline treatment levels retain the same relationship with the surrogates. Moreover, for conveniency to reader's more familiar with do-calculus notation, we can equivalently express this assumption as:
Finally, we also require a regularity condition of sequential positivity (aka overlap), which essentially states that the density of treatment is bounded away from zero a.s.. To define sequential positivity, we will denote with $\pi_d(\tau_t, s_{t-1})$ the marginal densities of the random variables $(T_t, S_{t-1})$, for any setting $d\in \{e, o\}$ and period $t\in [1,M]$. Then sequential positivity is defined as:
Under these high-level assumptions, we can show that the target outcome of interest is non-parametrically identified using a variant of the $g$-formula, based on a recursively defined estimand.
The non-parametric identification argument of Theorem (ref) requires the estimation of quantities of the form $\ensuremath{\mathbb{E}}[f_{t+1,j}(S_t)\mid S_{t-1}, T_t=0]$. When treatment $T$ is binary or discrete, then such quantities can be estimated in a relatively accurate manner by fitting nested regression models on the sub-population for which $T_t=0$. Moreover, we can also employ the great variety of doubly robust estimators for the quantity $\ensuremath{\mathbb{E}}_o[\bar{Y}^{(\underline{0}_2)}\mid S_1]$ (see e.g. tran2019double) combined with a doubly robust estimator for the surrogate part, to arrive at an overall doubly robust estimator. For instance, we can adapt the efficient influence function (EIF) representation of dynamic treatment effects proposed in scharfstein1999adjusting,van2011targeted,robins2000comment, to the case of a surrogate index setting as follows:
Since, we have that $\tau(t_1,t_0)=\sum_{m=1}^M \ensuremath{\mathbb{E}}_e\left[Y_m^{(t_1, \underline{0})}\right] - \ensuremath{\mathbb{E}}_e\left[Y_m^{(t_0, \underline{0})}\right]$, we can combine the doubly robust representations prescribed by Equation (ref), for each $m\in [1, M]$ and $\tau\in \{t_1, t_0\}$, to get an overall doubly robust representation of the target quantity. Similar to existing augmented inverse propensity methods in the dynamic treatment regime, the latter representation will lead to a consistent estimation if either all the models that go into the inverse propensity weights $\{W_{t}\}_{t=0}^M$ are consistent, or if all the nested regression functions $\{f_{t,j}\}_{1\leq t\leq m\leq M}$ are consistent. Moreover, this variant of the double robustness property also implies Neyman orthogonality (local robustness) of the moment implicitly defined by Equation (ref). Thus using the general results in chernozhukov2018double, we can devise an estimation strategy that enables valid inference while using machine learning, adaptive and regularized estimators for the auxiliary regression and classification models required by the above identification strategy. One could also adapt and apply alternative adaptive estimation frameworks, that also allow for the use of machine learning, adaptive estimators for the auxiliary models, such as the longitudinal targeted minimum loss estimation approach rotnitzky2012improved,van2011targeted, based on the latter representation of the target quantity.\footnote{We omit these details for succinctness and since the main estimation algorithm we propose in this work, which applies to both discrete and continuous treatments, appears in Section (ref) under a semi-parametric assumption.}
However, when treatment $T$ is continuous and potentially multi-dimensional, then non-parametric estimation rates for the quantities described in Theorem (ref) are required, without further assumptions, and can be potentially very slow and prohibitive. Moreover, finite sample performance will heavily depend on the number of samples observed in a region around the baseline treatment level at each period $T_t=0$, which could be very small and impact statistical power. Since our main application of interest (return-on-investments) involves multiple continuous treatments, being able to handle this setting is of primary practical importance.
To achieve parametric estimation rates, with valid confidence intervals, and more stable finite sample performance for the quantities of interest, even in the case of multiple continuous treatments, we will make further semi-parametric assumptions on the data-generating processes, i.e. that some parts of the data-generating process adhere to a known parametric form. One option for instance, would be to assume that the regression functions $\ensuremath{\mathbb{E}}[f_{t+1,j}(S_t)\mid S_{t-1}, T_t]$ adhere to some known parametric form, e.g. $\theta^\top\phi(T_t, S_{t-1})$, for a known feature map $\phi$. However, this essentially assumes a fully parametric model: even in the absence of any treatment, the world behaves in a simple manner. Unlike, for instance, in the classic partially linear model, where the baseline behavior under no-treatment is left non-parametric and only the effect of the treatment on the baseline behavior is modeled in a parametric manner. Instead, we could only model how these regression functions behave as the treatment $T_t$ deviates from the baseline, i.e.
Hence, analogous to the partially linear model, we are leaving un-modeled, the baseline behavior at each period, conditional on the past (the nested conditional mean). This is exactly the approach taken in the line of work on structural nested mean models (SNMMs), which we explore in the subsequent sections. As it will be shown below, the structural parameters $\theta$ of these nested means, can be identified without the need to estimate local non-parametric regression quantities of the form $\ensuremath{\mathbb{E}}[f_{t+1,j}(S_t)\mid S_{t-1}, T_t=0]$ and hence wont suffer from low sample sizes near the baseline treatment. Moreover, the target quantity of interest can be expressed in terms of these structural parameters $\theta$ of the SNMM.
The aforementioned semi-parametric assumption, can be expressed in terms of primitive counterfactual quantities, using the notion of a blip function.
These functions are a variant of what are known as the blip functions Chakraborty2013,Robins2004 and can be shown to be non-parametrically identifiable, assuming sequential conditional exogeneity and a sequential analogue of the positivity (aka overlap) assumption Robins2004. Theorem 3.1 of Robins2004 combines a telescoping sum argument and the sequential randomization condition to express counterfactual outcomes in terms of blip functions. We restate this result here, adapting it to our notation and our variant of sequential conditional exogeneity and blip function definition and providing a proof for completeness:
Intuitively, each term $\gamma_j$, removes from the outcome the blip effect of the observed action $T_j$. Consider any target outcome $Y_j$. When we remove $\gamma(T_j, S_{j-1})$ from $Y_j$, then what remains is, in-expectation (and crucially, even conditional on $S_{j-1}, T_{j}$), equal to the counterfactual outcome, where the sample received zero-treatment at period $j$. Subsequently, removing $\gamma(T_{j-1}, S_{j-2})$ from this remnant, then what remains is in-expectation (and crucially, even conditional on $S_{j-2}, T_{j-1}$), equal to the counterfactual outcome, where the sample received zero-treatment at periods $\{j-1,j\}$, and so on and so forth.
Note that if we denote with $\gamma_{o, t,j}$ the blip functions of the observational setting, then the latter lemma immediately gives an alternative identification strategy to the one presented in Theorem (ref), since we can write:
Thus if we can identify the blip functions, then the target quantity is also immediately identified without further assumptions.
One strategy for identifying the blip functions is to assume that they obey some known parametric form and then identify the parameters via a set of moment restrictions that the blip functions need to satisfy. In particular, by Lemma (ref), we know that the quantity $H_{t,j}(\theta^*)$ is equal in-expectation, and conditional on $S_{t-1}, T_t$ to the counterfactual outcome $Y_j^{(\bar{T}_{t-1}, \underline{0}_t)}$. However, this counterfactual outcome, by the sequential conditional exogeneity implied by the causal graph assumption, is independent of the treatment $T_t$, conditional on $S_{t-1}$, i.e. $Y_j^{(\bar{T}_{t-1}, \underline{0}_t)}\perp \!\!\! \perp T_t\mid S_{t-1}$. Thus for any function $f$ of $T_{t}, S_{t-1}$:
Moreover, by the conditional mean equivalence of this counterfactual outcome and the “remnant of the blip effects” $H_{t,j}(\theta^*)$, the same conditional mean independence property needs to hold for $H_{t,j}(\theta^*)$.
This leads to the following lemma:
Hence, if we have found the right $\theta^*$, then the infinite set of conditional moment restrictions in Equation (ref) need to be satisfied. Lemma (ref) is an adaptation of Theorem 3.2 of Robins2004 to our notation and we include its proof for completeness. Methods that estimate the structural parameters by utilizing such conditional mean independence moment restrictions are typically referred to in the literature on dynamic treatment effects as $g$-estimation methods.\footnote{$g$-estimation is a different term than $g$-computation, which typically refers to using the $g$-formula for dynamic treatment effects and estimating effects in a plug-in manner by estimation all conditional densities, and conditional means.}
One approach to operationalize Lemma (ref) would be to perform a grid search over some discretization of the parameter space and check that this set of conditional moment restrictions holds. In the full generality of structural nested mean models, without any further assumptions on the blip functions, such an exhaustive grid search could be inevitable, and renders the method impractical from a computational perspective.
For this reason, a typical approach in structural nested mean models, to render the methodology practical, is to assume a linear parametric form for the blip functions, leading to the class of linear structural nested mean models.\footnote{We note that the literature on $g$-estimation has also analyzed other forms of generalized linear parametric forms and provided practical methods (see e.g. Robins2004,Chakraborty2013,vansteelandt2014structural).}
Assuming that the expected conditional covariance matrix $\ensuremath{\mathbb{E}}[\ensuremath{\mathtt{Cov}}(\phi_t(T_t, S_{t-1})\mid S_{t-1})]$ of the feature map $\phi_t(T_t, S_{t-1})$ conditional on $S_{t-1}$, is full rank, then we can uniquely identify $\theta^*$ by finding a parameter vector $\theta$ that satisfies a small subset of the moment restrictions of the form:
What is most appealing about linear SNMMs is that the latter system of moment equations has a recursive closed form solution. In particular, we can express parameter $\theta_{t,j}$ as a function of parameters $\theta_{\tau,j}$ for $\tau>t$, in a closed form manner:
This immediately portrays the practicality of the method and the sufficiency of this subset of moment restrictions.
We will assume that both the data generating processes that generated the observational dataset and the experimental dataset obey a SNNM model with linear blip functions. Albeit, we allow both the treatments to change in between the two environments, as well as the blip function parameterizations to be different. We will denote with $\gamma_{e,t,j}, \gamma_{o,t,j}$, the blip functions in the two settings, with $\theta_{e,t,j},\theta_{o,t,j}$ the structural parameters of the blip functions in the two settings and with $\phi_{e,t},\phi_{o,t}$ the corresponding feature maps.
We start by presenting an identification argument for the target quantity of interest, as a function of the blip functions $\gamma_{o,t,j}$ in the observational dataset. Subsequently, in Theorem (ref), we combine it with a separate identification argument for the structural parameters $\theta_{o,t,j}$ of the blip functions, to arrive at a complete identification strategy.
One caveat of Theorem (ref) is that $\bar{Y}^{o,\ensuremath{\text{adj}}}$ and $g_{\ensuremath{\text{adj}}}^*$ are defined in terms of the dynamic effects $\{\theta_{o,t,j}\}_{2\leq t\leq j\leq M}$ of the observational setting, which are parameters that also need to be estimated. In particular, if we denote with $\Phi_{o,t} := \phi_{o,t}(T_t, S_{t-1})$, then we can write:
However, we can combine the Neyman orthogonal moment equations developed in lewis2020doubledebiased (which are an orthogonal variant of the moment equations in Equation (ref) and a variant of the doubly robust version of this equation introduces by Robins2004), with the Neyman orthogonal moment equation from Theorem (ref) to arrive at an overall Neyman orthogonal strategy for simultaneously identifying $\theta_0$ and these auxiliary dynamic effects. In particular, the parameters $\theta_{o,t}$ are identified recursively by the moment restrictions:
where $\bar{Y}_t := \sum_{j=t}^M Y_j$.
Collecting all the aforementioned discussion, we find that in order to identify the structural parameters of interest we need to estimate the following auxiliary nuisance models:
Note that all the nuisance functions $f$ are estimable from the observed data. All nuisance functions except $q$ correspond to a regression problem and $q$ can be decomposed into a classification problem for estimating the odds ratio $\frac{\Pr(e\mid S_1)}{1 - \Pr(e\mid S_1)}$ and a regression problem for estimate $\ensuremath{\mathbb{E}}_e\left[\tilde{\Phi}_1\mid S_1\right]$. Given these nuisance models we can define the parameter $\theta_0$ of interest as the solution to a set of moment restrictions that are Neyman orthogonal with respect to all the nuisance functions. To state our theorem we first define the vector of orthogonal scores.
We are now ready to state our main semi-parametric identification theorem via Neyman orthogonal moment restrictions:
Given that we have formulated the target structural parameters of interest as the solution to a vector of Neyman orthogonal moment equations, we can now easily transfer this identification argument to an estimation strategy, by invoking standard approaches. In particular, our estimation strategy will first estimate and apply the nuisance functions in a cross-fitting manner and subsquently solve a plug-in empirical analogue of the moment equations. Algorithm (ref) provides a formal description of the process.
To guarantee that our estimator is root-$n$ consistent and asymptotically normal, we need to assume that our first stage estimates of the nuisance functions are sufficiently accurate. In particular, we need to make the following nuisance rate assumptions:
Note that these nuisance rate assumptions possess almost a doubly robust flavor. With the exception of the nuisance quantities $\hat{p}_{o,t,t}$ and $\hat{p}_{e,1}$, which need to admit $o_p(n^{-1/4})$ root-mean-squared-error (RMSE) rates, for the remainder of the nuisance functions it suffices that the product of their RMSE rates with some other nuisance function be $o_p(n^{-1/2})$ and not that they individually satisfy $o_p(n^{-1/4})$ rates. For instance, if we knew the treatment policy in the experimental sample (captured by the propensity $p_{e,1}$) and the dynamic treatment policy in the observational sample (captured by the dynamic porpensity $p_{o, t,t}$), then we don't need any rates for $\hat{h}, \hat{p}_{e,t}, \hat{b}_{o,t}, \{\hat{p}_{o,j,t}\}_{t<j}$. Moreover, it suffices that the product of the surrogate score $\hat{q}$ error and the surrogate indices $\hat{g}, \hat{g}_t$ error, be small. Subject to these nuisance rate conditions we can show asymptotic normality of our estimate and provide asymptotically valid confidence intervals.
\vscomment{Add potential discussion on variance and jacobian.} \vsdelete{denote a $(m\,d)\times (m\,d)$ upper triangular matrix consisting of $d\times d$ blocks, such that the $(t,j)$ block, for $t\geq 2$ is defined as:
be an estimate of $J$ such that the $(t,j)$ block, for $t\geq 2$ is defined as:
}
The asymptotic linearity of our estimate also allows for alternative computationally convenient resampling methods for the construction of confidence intervals, with potentially better finite sample properties. For instance, constructing intervals by running the Bootstrap on the final stage estimation (keeping the nuisance estimates fixed), will be asymptotically valid. Moreover, the computationally even more convenient multiplier Bootstrap can also be used Chatterjee2005,Chernozhukov2013,Chernozhukov2014,Spokoiny2015,Zhilova2020, which can also be used for joint inference on multiple parameters, such as for constructing uniform confidence bands on dose response curves, i.e. the curve of the form $t \to \tau(t, 0)$, for $t$ in some bounded range $[U, L]$.
As a simple example where the linear SNMM assumption holds, consider the following linear Markovian (albeit high-dimensional) data generating process:
where $\epsilon_t, \eta_t, \zeta_t$ are i.i.d. random shocks. Our assumptions are satisfied if the quantities $B, C$ remain unchanged between the experimental and the observational setting, while the quantities $A, D, G$, as well as the distributions of mean-zero random shocks, can change arbitrarily, in the two settings, denoted as $A_d, D_d, G_d$ for $d\in\{o, e\}$.
In this case, the blip functions take the simple form: $\psi(\tau_t, s_{t-1})=\tau_t$ and $\theta_{t,j}=C B^{j-t} A$. Moreover, note that in this case, for any non-adaptive sequence of treatments $\tau_{>1}$, we have that:
Thus the quantity that our algorithm estimates is valid, irrespective of the baseline future policy that one considers and is a universal effect quantity that holds under any non-adaptive future sequence of treatments. This is practically convenient, as the causal effect derived is not heavily dependent on the future treatment sequence that a sample will receive in the short-term data set. Finally, note, that even though the surrogates/controls can be high-dimensional objects and hence the matrices $B, C, G$ are high-dimensional objects, our estimation strategy allows to estimate the target parameter $\theta_0 = \sum_{j=1}^m \theta_{1,j}$, which is low-dimensional at parametric root-$n$ rates and with asymptotically normal distributional limits. The intuition is that our analysis and estimation strategy, never really identifies or argues about estimation errors of these intermediate high-dimensional quantities.
We evaluate the performance of our proposed estimation strategy on a semi-synthetic dataset. The semi-synthetic data retain qualitative characteristics of data on real-world incentive investments in customers at a major corporation, although all data series and relationships have been perturbed to retain confidentiality.
The semi-synthetic dataset, like the real-world dataset on which it is based, displays several patterns that are common across many potential applications. The treatments, in this case incentive investments, are lumpy: in most periods most customers get no investments. Proxies, which include single period values of the outcome of interest, are highly auto-correlated over time. Treatments are also auto-correlated, and correlated with past values of proxies. Finally, we include a set of time-invariant controls that affect both proxies/outcomes and treatments.
To build the semi-synthetic data we estimate a series of moments from a real-world dataset: a full covariance matrix of all proxies, treatments, and controls in one period and a series of linear prediction models (lassoCV) of each proxy and treatment on a set of 6 lags of each treatment, 6 lags of each proxy, and time-invariant controls. Using these values, we draw new parameters from distributions matching the key characteristics of each family of parameters. Finally, we use these new parameters to simulate proxies, treatments, and controls by drawing a set of initial values from the covariance matrix and forward simulating to match intertemporal relationships from the transformed prediction models. For further details on the data generation process, see Appendix Section (ref).
We now compare multiple possible approaches for estimating the effects of our three synthetic treatments on a long-term outcome. To construct this outcome we select one proxy to be the outcome of interest. We consider the effect of each treatment in period $t$ on the cumulative sum of the outcome from period $t$ to $t+3$, four periods, or $t$ to $t+7$, eight periods. We can calculate the true treatment effects in the synthetic data as a function of parameters from the linear prediction models.
Because we construct a single, long synthetic dataset for this exercise it is possible to estimate the treatment effects on realized long-term outcomes directly, unlike the typical use case for a surrogate approach. Following a likely approach, we estimate the effect of each treatment at time $t$ on outcomes over the next 4 or 8 periods using double machine learning and controlling for invariant customer characteristics and contemporaneous and lagged values of all proxies and other treatments. The blue, “total" bars in each panel of Figure (ref) show the distribution of the estimation error\footnote{We use the $\ell_2$ error $\|\hat{\theta}-\theta_0\|_2$.} in the estimated treatment effects from this method across 100 simulated datasets. The top row plots the estimation error when estimating the effect on four periods of outcomes, increasing the sample size of each simulation from left to right, while the bottom row shows the same for the effect on eight periods of outcome. As predicted, the auto-correlation in treatments causes this method, which does not control for future treatments, to substantially overestimate treatment effects relative to their true values.
We then estimate the same set of treatment effects using the unadjusted surrogate approach described in Section (ref). The distribution of estimation errors from this approach is represented in the orange “surrogate" bars in each panel of Figure (ref). Since this approach still fails to control for future treatments when estimating the surrogate index, the estimated treatment effects are still substantially larger than the true effects on average. Note that the surrogate model exhibits slightly less bias than the direct “total" approach. Intuitively, because the surrogate approach is only capturing the relationship between treatment and outcome that passes through the surrogates it picks up less of the bias resulting from future correlated treatments than the direct approach.
The third set of green “adj. total" bars plot the distribution of estimation errors when estimating treatment effects on adjusted realized outcomes using the method of lewis2020doubledebiased. When a dataset containing both all treatments of interest and realized long-run outcomes is available, this should be the preferred approach. This third methodology, which removes the effects of future treatments from the long-run outcome in a first step, exhibits significantly less bias than the first two methods, particularly for reasonably large samples in the right two columns.
The final two bars in each panel of Figure (ref) illustrate the success of the adjusted surrogate approach described in Section (ref). We recommend this approach in the case when treatments are serially-correlated, as in the synthetic data, and it is not possible to collect a single dataset that contains both long-term outcomes and all the treatments of interest. As illustrated by the red “adj. surrogate" bars, this adjusted surrogate approach is highly accurate in predicting long-term effects with a performance comparable to that of having access to the raw long-term outcome itself. The final purple “new treat." bars show that the approach works equally well when considering the effect of a novel treatment that appears only in the experimental sample and was not part of the dynamic adjustment. Overall, this methodology overcomes a common data limitation when considering long-term effects of novel treatments and expands the surrogate approach to consider a common, and previously problematic, pattern of serially correlated treatments.