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.
54,180 characters · 17 sections · 32 citation commands
Modeling Covariate Transition for Efficient Estimation of Longitudinal Treatment Effects in Randomized Experiments
\etocdepthtag.toc{main}
\icmltitlerunning{Modeling Covariate Transition for Efficient Estimation of Longitudinal Treatment Effects in Randomized Experiments}
\twocolumn[ \icmltitle{Modeling Covariate Transition for Efficient Estimation of \\ Longitudinal Treatment Effects in Randomized Experiments}
\icmlsetsymbol{intern}{\textdagger}
\icmlaffiliation{sanken}{SANKEN, The University of Osaka, Osaka, Japan.} \icmlaffiliation{cyber}{CyberAgent, Inc., Tokyo, Japan} \icmlaffiliation{keio}{Keio University, Tokyo, Japan.}
\icmlcorrespondingauthor{Naoki Chihara}{[email removed]}
\icmlkeywords{Machine Learning, ICML}
\vskip 0.3in ]
\printAffiliationsAndNotice{\textsuperscript{\textdagger}Work done during research internship at CyberAgent AI Lab.}
Randomized experiments rct have been the gold standard for inferring causal effects and are widely employed across various domains, including medicine Rubin1997-ph and economics Duflo2007-yx,Athey2017-zx. However, the statistical power is limited, and large sample sizes are needed to obtain reliable results, as reported in online experiments Lewis2015-ke. To address this challenge, regression adjustment has been widely used to reduce the variance of treatment effect estimators by utilizing auxiliary covariates collected in randomized experiments Lin2013-xv,Deng2013-ig,Imbens2015-mr,Bloniarz2016-bi,Oka2026-eh.
Although regression adjustment under randomized experiments is widely employed for variance reduction in treatment effect estimation, standard methods typically utilize only pre-treatment covariates to estimate the average treatment effect (ATE). As a result, these methods overlook important insights into the temporal evolution of effects induced by early treatments. This limitation leads to a significant missed opportunity, especially in modern applications, where randomized experiments are often conducted over extended periods under static regimes. For example, the ATE does not capture the essential temporal dynamics of the response, such as whether the treatment yields an immediate impact, exhibits a delayed effect, or maintains persistence over time. Capturing these temporal dynamics requires incorporating information from intermediate outcomes and evolving post-treatment covariates. However, naively conditioning on such post-treatment variables allows for variance reduction but introduces post-treatment bias, distorting the target estimand even in randomized settings Freedman2008-qx,Montgomery2018-tn. For instance, in streaming services, user preferences dynamically evolve in response to recommendation algorithms. Adjusting for these evolving preferences as covariates blocks the causal pathway, thereby failing to capture the time-evolving impact of treatments on long-term engagement.
In this paper, we propose a regression adjustment method for the estimation of longitudinal treatment effects under randomized experiments. A key advantage of our framework is its ability to incorporate dynamic transitions in covariates induced by treatments into the estimation procedure. To account for temporal dynamics, we exploit transition kernels to represent how past treatments shape future histories and propagate outcome predictions forward along this kernel. This forward integration aggregates information contained in post-treatment covariates while preserving the marginal target estimand, thereby avoiding the post-treatment bias inherent in naive conditioning. As a by-product, this explicit modeling allows us to directly incorporate domain knowledge regarding system dynamics, such as seasonality or known exogenous shocks, into the transition mechanism.
In addition, the proposed method builds upon the Neyman-orthogonal moment conditions Chernozhukov2018-jn,Chernozhukov2022-ar, which provide robustness against first-order estimation errors in high-dimensional or complex nuisance components. These nuisance functions are estimated using flexible machine learning methods, making our method model-agnostic. Incorporating cross-fitting further strengthens robustness against estimation errors. With this design, we derive the asymptotic distribution of the proposed estimator and establish its semiparametric efficiency bound under mild regularity conditions. These results enable statistical inference, including standard error estimation and the construction of confidence intervals.
\noindentContributions. The main contributions of this paper are: (i) we present the dynamic regression-adjustment method using transition kernels and forward integration under randomized experiments; (ii) we establish asymptotic normality and provide valid inference under cross-fitted machine learning nuisance estimation, generalizing beyond traditional settings in causal inference; (iii) we derive the semiparametric efficiency bound for our estimator and show that it attains this bound; (iv) simulation studies and empirical analyses demonstrated the practical advantages of our approach.
\noindentOutline. The rest of our paper is organized in a conventional format. After the introduction, we review related works in (ref) and provide our problem definition in (ref). Next, we introduce the dynamic regression-adjusted estimator and its estimation procedure in (ref) and the asymptotic properties of our results in (ref). We then provide our experimental results and discussion in (ref), followed by a conclusion in (ref). All of the proofs for our theoretical results are provided in (ref).
\noindentLongitudinal treatment effects. Early seminal works in epidemiology and biostatistics established principled frameworks for estimating longitudinal treatment effects g-computation,msms,snmms. Recently, many researchers in the machine learning community have been studying sequential decision making in adaptive experimental designs, such as multi-armed bandits bandit. The adaptive nature of these approaches poses significant challenges for statistical inference, and it has been reported that they may introduce bias and make it harder to estimate causal effects Shin2021-li,Hadad2021-nr. In this context, policy evaluation in contextual bandits and reinforcement learning Dudik2011-pb,Jiang2016-hi,Kallus2020-ei,Hadad2021-nr, as well as dynamic treatment regimes dynamic-treatment,Robins2004-gi,Zhang2013-zg,Lewis2021-et,Bradic2024-ss, primarily aims to estimate the value of adaptive decision rules. Deep learning approaches have also been proposed for counterfactual prediction from temporal trajectories r-msn,crn,g-net,causal-transformer,Frauen2025-id, while our focus is on estimators with formal statistical inference and efficiency analysis under randomized experiments. From a different perspective, surrogate index methods Prentice1989-eh,surroage-index infer long-term treatment effects from short-term experimental results under strong surrogacy assumptions. Despite their respective strengths, the majority of these existing methods focus on estimating a single cumulative outcome, thereby obscuring the rich temporal dynamics of how treatment effects evolve at each specific time point. While the practical importance of tracking how treatment effects evolve over time has been widely recognized in large-scale online experiments Hohnhold2015-so, formal efficiency theory for the full trajectory of marginal treatment effects in static randomized experiments remains limited.
\noindentRegression adjustment. Randomized experiments have been widely used for the unbiased identification of the parameter of interest, but they often suffer from relatively low statistical power, as reported in online experiments Lewis2015-ke. To tackle this challenge, regression adjustment using auxiliary covariates to improve precision in treatment effect estimation has been extensively studied Fisher1925-sa,Cochran1977-ex,Yang2001-gj,Rosenbaum2002-ry,Freedman2008-yf,Tsiatis2008-ki,Rosenblum2010-sv,Lin2013-xv,Deng2013-ig,Ghadiri2023-lr. Recent work extends covariate adjustment to high-dimensional settings and modern machine learning, enabling variance reduction with many auxiliary variables Wager2016-ym,Bloniarz2016-bi,Poyarkov2016-mf,Guo2021-ug,Lei2021-ab,Masoero2023-lm,Chiang2025-bb. However, post-treatment covariates dynamically evolve in response to past treatments, and naive conditioning on them can distort the target estimand, even in randomized experiments Freedman2008-qx,Montgomery2018-tn. Our research tackles this difficulty using recursive forward integration based on transition kernels to exploit post-treatment histories for variance reduction while preserving the marginal estimand under longitudinal randomized experiments.
\noindentSemiparametric estimation. Our approach is grounded in semiparametric efficiency theory, which addresses the challenge of estimating low-dimensional parameters of interest in the presence of possibly infinite-dimensional nuisance components Klaassen1987-bq,Robinson1988-xc,Bickel1993-zx,Andrews1994-al,Newey1994-ni,Robins1995-ua. Building on this line of work, recent developments in double/debiased machine learning (DML) adapt semiparametric methods to accommodate flexible machine learning estimators for nuisance functions Ichimura2022-qp,Byambadalai2024-fa,Byambadalai2025-zz,Byambadalai2025-xo,Ahrens2025-ow. Specifically, our work leverages Neyman-orthogonal moment conditions Neyman1959-po,Chernozhukov2022-ar combined with cross-fitting Chernozhukov2018-jn. This strategy mitigates the impact of slow convergence rates and overfitting of complex nuisance estimators, yielding valid asymptotic distributions.
To the best of our knowledge, this is the first work to develop a regression-adjustment estimator for longitudinal treatment effects in randomized experiments that explicitly leverages covariate transition modeling using transition kernels and to establish that our algorithm attains the semiparametric efficiency bound for our estimator under this formulation.
First, we introduce key notations used in this paper. Please see (ref) for details. We consider randomized controlled trials (RCTs) with longitudinal treatments $\{W_{t}\}_{t=1}^{T}$ under static regimes, where $T$ is the time horizon and $W_{t}\in\mathcal{W}_{t}\coloneqq\{1,\ldots,\abs{\mathcal{W}_{t}}\}$ is the treatment assignment at time $t$. Let $Y_{t}\in\mathcal{Y}_{t}\subseteq\mathbb{R}$ denote the scalar-valued observed outcome of interest at time $t$ and $X_{t}\in\mathcal{X}_{t}\subseteq\mathbb{R}^{d}$ denote the covariates at time $t$. For each time $t$, we observe $n$ random samples $\{Z_{i,t}\}_{i=1}^n=\{(X_{i,t},W_{i,t},Y_{i,t})\}_{i=1}^n$ from a distribution on the product space $\mathcal{Z}_{t}\coloneqq\mathcal{X}_{t}\times\mathcal{W}_{t}\times\mathcal{Y}_{t}$. We denote $\bar{X}_{i,t}=\{X_{i,1},\ldots,X_{i,t}\}$ with $\bar{W}_{i,t},\bar{Y}_{i,t}$ defined analogously, and $Z_{t}\in\mathcal{Z}_{t}$ denotes the generic data point at time $t$. The probability of treatment assignment $\bar{w}_{t}$ is denoted as $\pi_{\bar{w}_{t}}=P(\bar{W}_{i,t}=\bar{w}_{t})$ satisfying $\sum_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}\pi_{\bar{w}_{t}}=1$, while $n_{\bar{w}_{t}}$ indicates the number of observations in treatment history $\bar{w}_{t}$ satisfying $\sum_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}n_{\bar{w}_{t}}=n$. Our work follows the potential outcome framework Rubin1974-gx,Imbens2015-mr in accordance with conventional settings and we let $Y_{t}(\bar{w}_{t})$ and $X_{t}(\bar{w}_{t-1})$ denote the potential outcome and covariate under longitudinal treatments $\bar{w}_{t}$.
Here, we provide the assumptions used in this paper.
Our parameters of interest are the expected potential outcomes at time $t$, as follows:
where $\bar{w}_{t}\in\bar{\mathcal{W}}_{t}$. Based on the parameters, we define the longitudinal average treatment effect (ATE) between treatments $\bar{w}_{t},\bar{w}_{t}'$ as
Under \crefrange{assume:sutva}{assume:positivity}, the expected outcomes $\mu_{\bar{w}_{t}}(t)$ are identifiable without any covariates while potential outcomes $\{Y(\bar{w}_{t})\}_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}$ are unobserved variables. Specifically, a simple estimator for $\mu_{\bar{w}_{t}}(t)$ is given by
While this estimator is unbiased and consistent, we aim to enhance its precision using observed historical data.
In this section, we present a new estimator that leverages observed historical data through covariate transition modeling. We also provide cross-fitted estimation and statistical inference procedures.
We begin by introducing the historical data $\bar{H}_{t}(\bar{w}_{t-1})=(\bar{X}_{t}(\bar{w}_{t-1}),\bar{Y}_{t-1}(\bar{w}_{t-1}))\in\bar{\mathcal{H}}_{t}$. Note that the historical data at time $t=1$ only has pre-treatment covariates, i.e., $H_1=X_1$ and $Y_0=0$. We adopt a regression adjustment framework to incorporate historical data $\bar{H}_{t}(\bar{w}_{t-1})$ into an estimation procedure. For each time $t\in[T]$ and longitudinal treatment $\bar{w}_{t}\in\bar{\mathcal{W}}_{t}$, we define the mean regression function $m_{\bar{w}_{t}}^{(t)}(\bar{h}_{t})$ as
The conditional mean function can be estimated using various supervised learning methods, such as LASSO, random forests, boosted trees, or deep neural networks. Under (ref), we rewrite $\mu_{\bar{w}_{t}}(t)$ using the mean regression function $m_{\bar{w}_{t}}^{(t)}(\bar{h}_{t})$ as
where $P_{\bar{H}_{t}}^{\bar{w}_{t-1}}$ is the probability measure on $\bar{\mathcal{H}}_{t}$ under $\bar{w}_{t-1}$. A natural estimator of $\mu_{\bar{w}_{t}}(t)$ is given by
where $\est m_{\bar{w}_{t}}^{(t)}$ is an estimator for $m_{\bar{w}_{t}}^{(t)}$. The estimator takes the form of the well-known augmented inverse-propensity weighting estimator. However, the empirical support for $\bar{H}_{t}$ under the longitudinal treatments $\bar{w}_{t-1}$ can be sparse, which leads to a high-variance estimate of (ref). In particular, when the units in the group $\{i:\bar{W}_{i,t-1}=\bar{w}_{t-1}\}$ coincide with those in $\{i:\bar{W}_{i,t}=\bar{w}_{t}\}$, we have $\est\mu_{\bar{w}_{t}}^{emp}(t)=\est\mu_{\bar{w}_{t}}^{adj}(t)$, which implies that no variance reduction is achieved.
To mitigate the aforementioned problem, we complement (ref) with a forward integration scheme over the dynamic trajectories in the historical data. Let $p^{(\tau)}_{\bar{w}_{\tau}}(dh_{\tau+1}\,|\,\bar{h}_\tau)$ be the transition kernel of $H_{\tau+1}$ given $\bar{h}_{\tau}$ under $\bar{w}_{\tau}$ and let $\estp^{(\tau)}_{\bar{w}_{\tau}}$ be an estimator for $p^{(\tau)}_{\bar{w}_{\tau}}$, where $\bar{w}_{\tau}\coloneqq(w_1,\ldots,w_{\tau})$ is the length-$\tau$ prefix of $\bar{w}_{t}$. The transition kernel can be learned using models that support conditional sampling, such as vector autoregressions, Gaussian processes, or deep neural networks. We decompose the probability measure using the transition kernels as follows:
Substituting the above decomposition into (ref), we obtain
We define the iterated conditional expectation $\Gamma_{\bar{w}_{t}}^{(1)}(h_1)$ under treatment history $\bar{w}_{t}$ as
and, recursively for $\tau=t - 1,\ldots,1$,
and $\Gamma_{\bar{w}_{t}}^{(t)}(\bar{h}_t)=m_{\bar{w}_{t}}^{(t)}(\bar{h}_{t})$. Intuitively, $\Gamma_{\bar{w}_{t}}^{(\tau)}(\bar{h}_\tau)$ is the conditional expectation of the terminal regression $m_{\bar{w}_{t}}^{(t)}(\bar{h}_{t})$ after propagating the remaining future history from time $\tau+1$ to $t$ according to the transition kernels under $\bar{w}_{t}$. Consequently, the target estimand satisfies $\mu_{\bar{w}_{t}}(t)=\mathbb{E}[\Gamma_{\bar{w}_{t}}^{(1)}(X_1)]$. Then, the resulting estimator augments the residual average with a fully marginalized iterated conditional expectation:
where $\est\Gamma^{(1)}_{\bar{w}_{t}}$ is an estimator for $\Gamma^{(1)}_{\bar{w}_{t}}$ and $\widehat{A}^{(t)}_{\bar{w}_{t}}$ is an estimator for the auxiliary correction term $A^{(t)}_{\bar{w}_{t}}$ defined as follows:
The second term in (ref) compensates for discrepancies between realized and one-step transitions via the auxiliary term $A^{(t)}_{\bar{w}_{t}}$, and the third integrates $\est m_{\bar{w}_{t}}$ along the transition under $\bar{w}_{t}$. The transition $\estp^{(\tau)}_{\bar{w}_{\tau}}$ supplies a forward model for the trajectory of covariates, enabling a stable approximation to (ref). The conditional regression functions $m_{\bar{w}_{t}}$ and the transition kernels $p^{(\tau)}_{\bar{w}_{\tau}}$ are treated as nuisance functions.
We rewrite our estimation problem as a moment condition problem. This formulation is crucial for establishing the asymptotic properties of our estimator when nuisance functions are estimated via machine learning (ML). We define $m_{t}\coloneqq(m_{\bar{w}_{t}}^{(t)}(\cdot))_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}$ and $p_{t}\coloneqq(p_{\bar{w}_{t}})_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}$, where $p_{\bar{w}_{t}}\coloneqq\{p^{(\tau)}_{\bar{w}_{\tau}}\}_{\tau=1}^{t-1}$. In addition, let $\theta_{t}=(\mu_{\bar{w}_{t}}(t))_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}$ be the target estimand. We define the moment functions as
In this equation, for each $\bar{w}_{t}\in\bar{W}_{t}$,
where $1\!\mathrm{l}\{\cdot\}$ represents the indicator function. The following lemma states what moment conditions are implied by our setup with a randomized experiment.
The key to valid inference with ML estimators lies in the robustness of this moment condition to local perturbations in the nuisance functions.
Neyman orthogonality implies that the moment condition is first-order robustness of the moment condition to errors in the nuisance estimates. Leveraging this property, coupled with cross-fitting, allows us to obtain asymptotic normality of the dynamic regression-adjusted estimator under mild conditions, despite using ML models for the nuisance functions.
This lemma shows that the solution to (ref) coincides with the vector whose components are the dynamic regression-adjusted estimators $\est\mu_{\bar{w}_{t}}^{dy\text{-}adj}(t)$ defined in (ref). (ref) ensure that the estimator $\est\mu_{\bar{w}_{t}}^{dy\text{-}adj}(t)$ is insensitive to the estimation errors in the nuisance functions.
Here, we explain our algorithm for computing the dynamic regression-adjusted estimator $\est\mu_{\bar{w}_{t}}^{dy\text{-}adj}(t)$ defined in (ref). (ref) shows an overview of our estimation procedure. To ensure the validity of our asymptotic results and avoid overfitting biases, we employ an $L$-fold cross-fitting strategy. This procedure decouples nuisance function estimation from moment condition evaluation. The estimation proceeds in two main stages:
\noindentStep 1. Nuisance estimation. We randomly partition the observation indices into $L$ disjoint folds. For each fold $\ell$, we utilize the units in the remaining $L-1$ folds as the training set. On this training set, we fit the supervised learning model $\mathcal{M}$ to estimate the conditional mean function $\est m_{\bar{w}_{t}}^{(t)}$ and the transition model $\mathcal{T}$ to learn the transition kernel $\estp_{\bar{w}_{\tau}}^{(\tau)}$. These trained models are then applied to the units in fold $\ell$ to generate out-of-sample predictions and transition estimates.
\noindentStep 2. Recursive integration. We calculate the integral over the high-dimensional historical covariate space defined in (ref). However, since analytical integration is often intractable, we approximate it using Monte Carlo integration. Specifically, for each unit in fold $\ell$, we draw $S$ samples from the estimated transition kernel $p_{\bar{w}_{\tau}}^{(\tau)}$ to recursively compute $\est\Gamma^{(\tau)}_{\bar{w}_{t}}$ backwards from $\tau=t-1$ down to $1$.
Lastly, these components are combined to compute the auxiliary correction term $\widehat{A}_{\bar{w}_{t}}^{(t)}$ and the target estimator $\mu_{\bar{w}_{t}}^{dy\text{-}adj}(t)$.
Our estimator is easy to use for statistical inference, such as standard error estimation and the construction of confidence intervals, because it is based on Neyman-orthogonal moment conditions in (ref) and has an explicit influence function in (ref). In this section, we briefly explain two practical inference methods for the target parameter $\theta_{t}=(\mu_{\bar{w}_{t}}(t))_{\bar{w}_{t}\in\bar{\mathcal{W}}_{t}}$ and for contrasts such as longitudinal ATEs. The asymptotic validity of these procedures is established in the next section.
\noindentAnalytical variance estimation. The asymptotic normality result in (ref) allows us to estimate standard errors using the sample variance of the estimated influence functions. Since the moment function $\psi_{t}^\pi$ derived in (ref) coincides with the efficient influence function ((ref)), the asymptotic covariance matrix $\Sigma_{t}$ is consistently estimated by the empirical second moment of the scores:
Using this consistent variance estimator, a $(1-\alpha)$ Wald-type confidence interval for a linear contrast $c^\top\theta_{t}$ can be constructed as
where $z_{1-\alpha/2}$ denotes the $(1-\alpha/2)$-quantile of the standard normal distribution.
\noindentMultiplier bootstrap. While the analytical approach yields closed-form standard errors, it relies on asymptotic approximations that may be inaccurate in finite-sample settings. The multiplier bootstrap multiplier-bootstrap provides an alternative approximation that can improve finite-sample accuracy relative to the normal approximation. Importantly, in our semiparametric setting, this method avoids the computationally expensive step of refitting nuisance models for each bootstrap iteration Chernozhukov2018-jn. Specifically, we first draw multipliers $\{\xi^{(b)}_i\}_{i=1}^n$ for $b\in[B]$ independently of data, from a distribution with mean zero and unit variance (e.g., standard normal or Rademacher). Then, for each $t\in[T]$, we compute
Conditional on the observed data, the empirical distribution of $\{G_{t}^{(b)}\}_{b=1}^B$ serves as an approximation to the law of $\sqrtn(\est\theta_{t}-\theta_{t})$. A two-sided $(1-\alpha)$ confidence interval for a scalar contrast is then given by
where $\widehat{q}_{\alpha}$ is the $\alpha$-th empirical quantile of $\{c^\top G_{t}^{(b)}\}_{b=1}^{B}$. This bootstrap is computationally attractive in our setting because it operates only on the already computed influence function evaluations and therefore avoids repeating the training of the nuisance learners.
In this section, we derive the asymptotic distribution of the proposed estimator, which enables statistical inference and the construction of confidence intervals. Additionally, we establish the semiparametric efficiency bound for our dynamic regression-adjusted estimator and demonstrate that our estimator achieves this bound under the specified assumptions. We begin by introducing additional assumptions to formalize our results.
(ref) rules out extremely heavy-tailed outcomes and allows us to apply standard limit theorems. (ref) requires that the cross-fitted nuisance functions become sufficiently accurate as the sample size grows. Note that transition kernels are viewed as maps from the conditioning history to the vector space of finite signed measures, equipped with the variation norm. (ref) states that small estimation errors in the nuisance functions do not get amplified when constructing $\est\Gamma^{(\tau)}_{\bar{w}_{t}}$ and $\widehat{A}^{(t)}_{\bar{w}_{t}}$.
We now establish the weak convergence of our proposed estimator in the following theorem, which serves as the theoretical foundation for statistical inference.
(ref) and \crefrange{theorem:moment_condition}{theorem:sample_moment_condition} show that the moment function $\psi^{\pi}_{t}(\bar{Z}_{t})$ is the influence function. The next theorem reveals that $\psi^{\pi}_{t}(\bar{Z}_{t})$ is also the efficient influence function.
(ref) demonstrate that the dynamic regression-adjusted estimator is asymptotically normal and semiparametrically efficient, reaching the efficiency bound.
We conducted simulation study to validate our theoretical results and demonstrated the effectiveness of our estimator using real-world A/B test data collected from a large-scale streaming platform in Japan.
\noindentSetup. We explain our synthetic data generating process. We generate the outcome $Y_{i,t}$ and covariate $X_{i,t}$ for the $i$-th unit at time point $t$ according to the following process:
where the initial states $X_{i,1}\sim\mathcal{N}(F,\sigma_0^2I_d)$ and $Y_{i,0}=0$, the error terms $\varepsilon_{i,t}\sim\mathcal{N}(0,\sigma_Y^2)$, $B_t$ is a standard Wiener process representing system noise, and $Q$ is a random orthogonal matrix sampled uniformly from the Haar measure. Here, let $g(\cdot)$ be a nonlinear outcome function depending only on a subset of covariates, and let $\tau(t)$ be a time-varying treatment effect. The covariate transition $f(\cdot)$ is constructed based on the rotated forced Lorenz-96 model lorenz96, which is commonly used as a benchmark. This function includes external time-dependent forcing $F_{t}$ coupled with both historical outcomes and treatment assignments, which together shape the evolution of the covariates over time. In addition, we drew a subject-level treatment indicator $W_{i}\sim\text{Bernoulli}(0.5)$ independently across $i$ and set $W_{i,t}=W_i$ for all $t\in\![T]$, which follows a standard A/B testing setup. Detailed descriptions of the data generating process are provided in (ref). This design incorporates nonlinear dependencies, irrelevant covariates, influence from past treatments, and time-varying treatment effects, so the synthetic data generated by this process reflects the complexities of real-world scenarios.
We used synthetic data of size $n\in\{500,1000,2000\}$ and length $T=7$ with $10$-dimensional covariates and estimated the longitudinal ATE defined by (ref) with $1000$ simulations. We approximated the ground-truth values using Monte Carlo sampling with $10^6$ draws. We used an OLS model (linear) and an RF model (nonlinear) as supervised learning models $\mathcal{M}$, and a VAR model (linear) and a Gaussian MLP model (nonlinear) as transition learning models $\mathcal{T}$. For instance, the estimator combining a RF model and Gaussian MLP model is referred to as RF-MLP adjustment. All adjusted estimators used $5$-fold cross-fitting. Our source code for the simulation study is publicly available at: \url{https://github.com/C-Naoki/dynamic-adjustment}.
\noindentResult. (ref) shows the statistical properties (RMSE, average $95\%$ confidence interval (CI) length, and its coverage probability) of different estimators using synthetic data with a sample size of $n=1000$, based on $1000$ simulation runs. Additional results for other sample sizes are found in (ref). The pointwise confidence intervals are calculated using sample estimates of the asymptotic variance. All adjusted estimators exhibit lower RMSE and shorter average 95% confidence interval lengths compared to the empirical estimator while maintaining the nominal coverage probability. In particular, RF-MLP adjustment achieves the best performance, demonstrating the benefits of incorporating machine learning techniques. OLS-MLP adjustment yields performance comparable to RF-MLP adjustment because the lagged data itself contains a certain amount of information, allowing even a simple linear baseline to capture the underlying trend to a large extent. In contrast, RF-VAR adjustment offers only a subtle improvement because a VAR model fails to capture the complex nonlinear dynamics, resulting in sample trajectories that deviate significantly from the true data distribution. (ref) shows RMSE reduction in % relative to the empirical estimator for different sample sizes $n\in\{500,1000,2000\}$. Across all sample sizes, dynamic regression adjustment consistently yields substantial precision gains with performance improving as sample size $n$ increases. These results validate the theoretical discussion presented in (ref), and emphasize the value of flexible dynamic regression adjustment in improving finite-sample efficiency for longitudinal treatment effect estimation.
\noindentSetup. In this section, the proposed dynamic regression adjustment framework is applied to evaluate the results of a longitudinal A/B test at a streaming platform in Japan, using non-public proprietary data provided by the company. The A/B test was conducted to investigate how different types of up-next content recommendations influenced user viewing time over a 10-day period. The logged data on the platform is organized as panel data, with each observation showing the daily viewing behavior during the experiment. Users were randomly assigned to one of two groups. The control group received recommendations based on user interaction data, such as past clicks, while the treatment group received recommendations based on content-level similarity derived from item features. The outcome variable is the daily viewing time for a specific series, and the covariate variables are the daily viewing times for the five series, including lagged outcome variables. To ensure reliable analysis, we selected a total of 6,802 users with sufficient observation coverage during the experimental period. For regression adjustment, we used an RF model and a zero-inflated log-normal (ZILN) MLP model based on the domain knowledge that users did not watch the target series daily, leading to a large number of zero observations. We used 5-fold cross-fitting as well as the simulation studies.
\noindentResult. (ref) (a) visualizes the estimated trajectories of the longitudinal treatment effects on daily viewing time over the 10-day experiment. The left panel displays the results from the empirical estimator, while the right panel presents those from our proposed dynamic regression-adjusted estimator. The shaded regions represent the $95\%$ pointwise confidence intervals computed using multiplier bootstrap multiplier-bootstrap with $1000$ repetitions. We can see that the treatment effect initially showed slightly positive values but subsequently decreased. This pattern may reflect that users initially engaged with the new recommendations but did not sustain their viewing time, possibly due to a mismatch with their longitudinal preferences. (ref) (b) quantifies the corresponding precision gains in terms of the standard errors (SE) and the reduction rates over time. Specifically, the regression adjustment reduces the standard errors by approximately $0.4\%$ to $20.4\%$ across the experimental period. These empirical findings demonstrate the practical utility of incorporating covariate transition dynamics to enhance statistical power in large-scale online experiments.
We present a regression adjustment framework for the estimation of longitudinal treatment effects under randomized experiments. By explicitly modeling covariate dynamics via transition kernels and employing recursive forward integration, our method successfully extracts information from evolving post-treatment histories without distorting our target estimand. Moreover, we derive several theoretical foundations of our proposed estimator, including asymptotic normality ((ref)) and a semiparametric efficiency bound ((ref)). Simulation studies and empirical validations using A/B test data from a large-scale streaming platform demonstrated the practical advantages of our approach.
\noindentLimitations. Our approach has some limitations, which are interesting to tackle in the future. First, this paper considers experimental data under perfect compliance and no interference. Such conditions are realistic in standard A/B tests, but these assumptions may be restrictive in other contexts, such as covariate-adaptive randomization (CAR) Tsiatis2008-ki and observational settings. In these scenarios, time-varying confounding and various latent regimes Saggioro2020-rp,rahmani2025causal,modeplait,Mameche2025-hx may further complicate the estimation of nuisance components. In addition, any learned components would be treated as additional nuisances, whose uncertainty is propagated to inference. Second, adapting our framework for distributional treatment effects (DTE) Doksum1974-ns,Lehmann2006-wm provides richer insights into heterogeneous impacts than solely focusing on overall average effects. Lastly, addressing irregular sampling intervals and missing data through continuous-time modeling Chen2018-xe may broaden empirical applicability.
The goal of this paper is to advance the field of machine learning. There are many potential societal consequences of our work. Here, we highlight key points to ensure our method works well. At first, our framework builds on the Neyman-orthogonal moment conditions, which provide first-order robustness against small errors in nuisance estimation. However, severe model misspecification can still yield unexpected or unreliable results, making the rigorous validation of nuisance functions essential. We recommend a step-by-step approach to model selection and validation. When modeling the nuisance functions $m_t$ and $p_t$, it is best to establish a stable foundation by starting with simple linear baselines. If the underlying system dynamics turn out to be complex or the linear models underperform, practitioners should then scale up to more flexible, nonlinear models. Furthermore, while the theoretical assumptions cannot be directly evaluated on real-world data, the cross-fitting procedure offers an empirical diagnostic. Practitioners should closely monitor out-of-fold predictive performance (e.g., out-of-fold RMSE). Achieving high predictive accuracy without severe overfitting on the validation folds serves as a highly reliable, practical proxy for confirming that the necessary structural assumptions are reasonably satisfied. Also, we emphasize that our work focuses on static regimes (i.e., pre-specified treatment sequences). Therefore, adaptive experiments, such as contextual bandits, fall outside the scope of our method. Applying our framework to such settings could lead to substantial bias. The recursive integration and Monte Carlo sampling used in our algorithm are often computationally intensive, particularly for long time horizons, large sample sizes, and high-dimensional covariates. If the data is too large, our method may be impractical due to limited computational budgets. Lastly, although the proposed estimator is asymptotically normal, finite-sample bias can cause the coverage of confidence intervals to deviate from nominal levels. To improve finite-sample accuracy, the multiplier bootstrap is often an effective approach, but it still relies on empirical evaluations of the influence functions. Consequently, if the nuisance components are poorly estimated due to insufficient sample sizes, the bootstrap cannot fully correct for the resulting inference errors.
We would like to thank the anonymous referees for their valuable and helpful comments. This work was partly supported by “Program for Leading Graduate Schools” of the Osaka University, Japan, JSPS KAKENHI Grant-in-Aid for Scientific Research Number JP25KJ1729, JP26H02499, JST CREST JPMJCR23M3, JST K Program JPMJKP25Y6, JST COI-NEXT JPMJPF2009, JST COI-NEXT JPMJPF2115, Future Social Value Co-Creation Project - Osaka University.