We consider experimentation in the presence of non-stationarity, inter-unit (spatial) interference, and carry-over effects (temporal interference), where we wish to estimate the global average treatment effect (GATE), the difference between average outcomes having exposed all units at all times to treatment or to control. We suppose spatial interference is described by a graph, where a unit's outcome depends on its neighborhood's treatments, and that temporal interference is described by an MDP, where the transition kernel under either treatment (action) satisfies a rapid mixing condition. We propose a clustered switchback design, where units are grouped into clusters and time steps are grouped into blocks, and each whole cluster-block combination is assigned a single random treatment. Under this design, we show that for graphs that admit good clustering, a truncated Horvitz-Thompson estimator achieves a O~(1/NT) mean squared error (MSE), matching the lower bound up to logarithmic terms for sparse graphs. Our results simultaneously generalize the results from \citet{hu2022switchback,ugander2013graph} and \citet{leung2022rate}. Simulation studies validate the favorable performance of our approach.
Figures & tables
Figure 1: Clustered Switchback Experiments . The image illustrates clustered switchback on DoorDash [ 34 ] . The time and geographical locations are grouped into blocks. Each spatio-temporal cluster (i.e., Cartesian product) is independently assigned a treatment arm. The goal is to estimate the difference in the average (counterfactual) “outcomes” (e.g., revenues) between the all-treatment and all-control policy.
Interference Graph
(Spatio-) Clustering
Timeblocks
MSE
No edges (no interference)
one node per cluster
r=O~(tmix)
Corollary 2 : O(TNtmix)
Arbitrary
all nodes in single cluster (pure switchback)
r=O~(tmix)
Corollary 3 : O(Ttmix) (also Hu and Wager [16] )
Max Degree D
max cluster degree dΠmax max cluster size Nc
r=O~(tmix)
Corollary 4 : O(TNtmixDdΠmaxNc22dΠmax)
2D-grid spatial graph edges w/in distance h
uniform 2D grid
r=O~(tmix)
Corollary 6 : O(TNtmixh2)
2D-grid spatial graph edges w/in distance
uniform 2D grid
r=T (pure GCR)
Leung [23] : O(Nh2)
κ -Restricted Growth Max Degree D
1 -hop-max randomized clustering
r=O~(tmix)
Corollary 7 : O(TNtmixD3κ5)
Table 1: Implications of Our Main Result alongside Known Results: Our main theorem implies the following MSE bounds in several fundamental special cases. We omit polylogarithmic terms in N,T in these MSE bounds. We also show our results alongside several known results for “pure” switchback and “pure” A/B testing under interference.
Figure 2: Positioning of this work. By and large, a model for experimentation involves a subset of four key features, as illustrated above. Prior work has addressed spatial and temporal interference separately, often incorporating heterogeneity across individuals as well. However, in real-world applications, all four features often arise simultaneously. This work aims to lay the theoretical foundation for cluster switchback experiments, which are increasingly used in practice to navigate such complex interference.
Figure 3: Dependence Graph: The regions correspond to the clusters in a partition Π . Units i,j intersect a common cluster C , so (i,j)∈EΠ (or i⊥⊥j ).
Design
When to Use
MSE bound
Pure Switchback
High temporal mixing; environment stable over time
O~(Ttmix)
Pure Spatial Clustering
Strong spatial interference; sparse or weakly connected graph
O(2−2dΠmaxN2∣{(i,i′):i⊥⊥i′}∣)
Clustered Switchback
Both significant spatial interference and fast temporal dynamics
O~(Ttmix⋅2−2dΠmaxN2∣{(i,i′):i⊥⊥i′}∣)
Table 2: When to use each experimental design.
Figure 4: Locations of 169 Hotels in Washington DC provided in dataset [ 8 ] .
Statistic
Value
Nodes
169
Edges
4,332
Mean degree
51.3
Median degree
56
Std. degree
30.3
Min / Max degree
0 / 95
Table 3: Summary statistics of the DC hotel interference graph ( κ=0.1 ).
Figure 8
Figure 7: MSE vs. total (spatio-temporal) clusters K , with HT estimator, δ=0.2 . Joint spatio-temporal design has lower MSE than both pure designs.
Figure 8: Heatmaps of bias (left), standard deviation (centre), and MSE (right) for the HT estimator as a function of time-block length ℓ (rows) and fractional recency fr=r/ℓ (columns), fixing 1/s=4 and δ=0.2 (Uniform network, n=100 , T=10,000 ). MSE is U-shaped in fr for every ℓ , with the optimum near fr∗≈0.5 .
Figure 9: Heatmaps of bias (left), standard deviation (center), and MSE (right) for the HT estimator as a function of spatial cells per dimension 1/s (rows) and fractional recency fr (columns), fixing ℓ=100 and δ=0.2 (Uniform network). MSE is minimized at an interior value 1/s∈{3,4} ; beyond this, variance rises as each spatial cluster contains fewer units.
Figure 10: Heatmaps of bias, standard deviation, and MSE for the HT estimator over (ℓ,1/s) ( δ=0.2 , best r per cell). Rows: 1/s∈{1,…,8} ; columns: ℓ∈{40,…,1000} . This is expected, since when ℓ is small, “contaminated” data (rounds) make up a higher proportion of the entire data.
Figure 11: MSE vs. time-block length ℓ for the Hajék (left), HT (center), and DM (right) estimators, fixing 1/s=3 and δ=0.2 (Uniform network, x -axis truncated to ℓ≤400 , y -axis to [0,0.5] ). Each curve corresponds to a different fractional parameter: fr∈{0.25,0.5,1,1.5,2} for Hajék/HT and fb∈{0,0.1,0.3,0.5} for DM. Hajék uniformly dominates HT; DM with fb=0 is competitive at large ℓ but biased at small ℓ .
Figure 12: Bias and standard deviation of the Hajék, HT, and DM estimators as a function of time-block length ℓ , fixing 1/s=10 and δ=0.2 . Each curve corresponds to a different fractional parameter: fr∈{0.25,0.5,1,1.5,2} for Hajék/HT and fb∈{0,0.1,0.3,0.5} for DM. Shaded bands indicate ±1 standard deviation.
Figure 13: MSE of the Hajék, HT, and DM estimators as a function of time-block length ℓ , fixing 1/s=10 and δ=0.2 . Each curve corresponds to a different fractional parameter: fr∈{0.25,0.5,1,1.5,2} for Hajék/HT and fb∈{0,0.1,0.3,0.5} for DM.
Figure 14: Bias and MSE of the DM estimator on the Hotel network as a function of ℓ , for burn-in fractions b/ℓ∈{0,0.1,0.3,0.5} . Each MSE curve has an interior minimum, and a small positive burn-in fraction substantially reduces bias with modest variance cost.
Figure 15: Bias and MSE of the HT estimator as a function of ℓ , for recency values r∈{5,10,20,40,80,…} with δ=0.2 . The HT estimator exhibits large MSE relative to Hajék (note the y -axis scale), and the optimal r depends on ℓ .
Figure 16: Bias and MSE of the Hajék estimator as a function of ℓ , for recency values r with δ=0.2 . The MSE scale ( ≤0.02 ) is far smaller than for HT ( Fig. 15 ), highlighting the efficiency advantage of the Hajék correction.
This paper studies experimental designs for estimation and inference on policies with spillover effects. Units are organized into a finite number of large clusters and interact in unknown ways within each cluster. First, we introduce a single-wave experiment that estimates the marginal effect of a change in treatment probabilities, taking spillover effects into account. Using the marginal effect, we propose a policy optimality test and improve the population's average outcome through marginal changes. Second, we design a multi-wave experiment to estimate optimal treatment rules, with strong regret guarantees. We provide an implementation in a large-scale field experiment studying climate adaptation.
Davide Viviano, Jess Rudder
Department of Economics, Harvard University, 1805 Cambridge St, Cambridge, MA 02138. · Department of Economics, Oregon State University, 1500 SW Jefferson Way, Corvallis, OR 97331.
Adaptive experimentation under unknown network interference requires solving two coupled problems: (i) learning the underlying dynamics of interference among units and (ii) using these dynamics to inform treatment allocation in order to maximize a cumulative outcome of interest (e.g. revenue). Existing adaptive experimentation methods either assume the interference network is fully known or bypass the network by operating on coarse cluster-level randomizations. We develop a Thompson sampling algorithm that jointly learns the interference network and adaptively optimizes individual-level treatment allocations via a Gibbs sampler. The algorithm returns both an optimized treatment policy and an estimate of the interference network; the latter supports downstream causal analyses such as estimation of direct, indirect, and total treatment effects. For additive spillover models, we show that total reward is linear in the treatment vector with coefficients given by an n-dimensional latent score. We prove a Bayesian regret bound of order nT⋅Blog(en/B) for exact posterior sampling; empirically, our Gibbs-based approximate sampler achieves regret consistent with this rate and remains sublinear when the additive spillovers assumption is violated. For general Neighborhood Interference, where this reduction is unavailable, we analyze an explore-then-commit variant with O(n2logT) graph-discovery cost. An information-theoretic Ω(nlogT) lower bound complements both results. Empirically, our method achieves more than an order-of-magnitude reduction in regret in head-to-head comparisons. On two real-world networks, the algorithm achieves sublinear regret and yields downstream effect estimates with small RMSE relative to the truth.
Estimating causal effects from real-world spatiotemporal data is challenging due to hidden confounders and interference. Standard causal identification methods assume conditional exchangeability given observed covariates, which fails whenever hidden confounders affect both treatment and outcomes - a common setting in domains such as climate, environmental policy, epidemiology, and regional economics. In this paper, we propose a novel spatiotemporal proximal causal inference framework that extends proximal identification theory to spatiotemporal settings. The proposed method jointly captures local and neighborhood-level confounding information by introducing treatment- and outcome-inducing proxies, and we derive a spatiotemporal outcome confounding bridge function that identifies the potential outcome without requiring direct recovery of the hidden confounder. We establish the identifiability of this bridge function under proxy exclusion restrictions and a spatiotemporal completeness condition, and show that the resulting estimator recovers the outcome through a proximal generalization of the g-computation formula. To operationalize this identification result, we propose a neural architecture that learns proxies via transformer-based spatiotemporal encoders - coupled with a conditional mutual information critic to enforce exclusion restrictions and a moment-matching network to guarantee that the learned bridge function satisfies the underlying identifying equation. We further introduce a stabilized weighting scheme to address treatment support imbalance. Experiments on synthetic datasets demonstrate that our approach achieves comparable performance to baseline causal inference methods, while providing, to our knowledge, the first theoretically grounded outcomes for the hidden confounding in the presence of spatiotemporal interference through a proximal causal inference framework.