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.
42,549 characters · 9 sections · 62 citation commands
Neighborhood Stability in Double/Debiased Machine Learning with Dependent Data
\onehalfspacing
Double/debiased machine learning (DML) methods enable $\sqrt{n}$-consistent estimation of low-dimensional parameters, such as average treatment effects, using modern machine learning methods to approximate complex nuisance parameters. Large-sample theory for DML is mostly limited to i.i.d.\ data, but there is recent interest in applications to dependent data. brown2024inference studies DML trained on time series data. gilbert2024causal employ DML to control for spatial confounding. emmenegger2025treatment, leung2025graph, and wang2024graph apply DML to network data.
A key theoretical challenge is verifying the high-level condition
where $\theta_0$ is the low-dimensional parameter of interest, $M_n$ is the empirical moment, $M$ the population moment, $g_0$ the nuisance parameter, and $\hat{g}$ the machine learning estimate of $g_0$. The goal of this paper is to provide lower-level conditions for (ref), without cross-fitting, when the data is a $\beta$-mixing process. We allow observations to be situated in a generic metric space that can accommodate temporal, spatial, and network data.
There are three main approaches to verifying (ref). First, classical semiparametric estimation theory takes the approach of verifying stochastic equicontinuity, which implies (ref) andrews1994asymptotics. This condition replaces $\hat{g}$ on the left-hand side with a deterministic function $g$ and takes the supremum over all $g$ in the nuisance space $\mathcal{G}$.\footnote{Several papers verify stochastic equicontinuity for neural networks without cross-fitting. See farrell2021deep for multilayer perceptrons and i.i.d.\ data, wang2024graph for graph neural networks and clustered network data, and brown2024inference for multilayer perceptrons and time series data.} As discussed in chernozhukov2018double, this requires restrictions on the complexity of $\mathcal{G}$ that can be too strong for modern machine learning applications. chen2022debiased note that these restrictions can be conservative since they only concern properties of $\mathcal{G}$ rather than the machine learner $\hat{g}$.
Second, we can employ cross-fitting, which essentially entails training $\hat{g}$ on a random subsample of the data (training fold) and evaluating $M_n$ on the remaining subsample (evaluation fold). When the data is i.i.d., this ensures $M_n(\theta_0, \cdot)$ and $\hat{g}$ are independent, making (ref) straightforward to verify. With dependent data, however, it is an open question how cross-fitting should be modified in general to establish (ref).
emmenegger2025treatment and gilbert2024causal consider settings with local dependence where observations are independent when separated by some minimum distance in space or a network. They naturally propose to cross-fit by excluding from the training fold observations that are sufficiently far from the evaluation fold. This modification, however, can substantially reduce the sample size of the training fold in practice, particularly with network data since the number of units at distance $k$ can grow rapidly with $k$. In the simulation study in (ref), the training fold size can be less than {\em ten} for reasonable parameterizations of the network formation model.
We therefore opt to avoid cross-fitting and follow a third and more recent approach due to chen2022debiased that builds on chernozhukov2024adversarial. They establish (ref) using a natural “stability” condition on the machine learner $\hat{g}$, rather than conditions on $\mathcal{G}$. Stability essentially requires $\hat{g}$ to be insensitive to perturbations of the training data that replace any observation with an independent copy. It has been commonly employed in the statistical learning literature austern2020asymptotics,bayle2020cross,bousquet2002stability,celisse2016stability,elisseeff2003leave and is advantageous due to its intuitive appeal and relative ease of verification for certain machine learners.
We show that the analysis of chen2022debiased can be extended to dependent data. A key piece of the argument is strengthening stability to “neighborhood” stability. This requires insensitivity to perturbations that replace the observations in a slowly growing neighborhood of any given observation with an independent copy. Resampling neighborhoods is somewhat analogous to modifications of the bootstrap and subsampling for dependent data that resample blocks of contiguous observations. Our conditions intuitively connect the degree of dependence in the data to the strength of the restriction on the machine learner. Specifically, when the dependence between observations decays faster with distance, the neighborhood stability condition is weaker in that the machine learner needs only to be invariant to resampling smaller neighborhoods of observations.
Stability has been verified for a variety of machine learners, including decision trees arsov2019stability, regularized M-estimators bousquet2002stability, M-estimators trained via stochastic gradient descent hardt2016train, and bagged learners chen2022debiased. We show that several of these results can be adapted to verify neighborhood stability. The sufficient conditions on the machine learner are satisfied when dependence between observations at distance $r$ is small relative to the sizes of $r$-neighborhoods. This is a familiar theme for central limit theorems (CLTs) for dependent data jenish2009central,kojevnikov2021limit. When $r$-balls in Euclidean space contain more observations, a given observation is proximate to many others, so a spatial CLT requires dependence to decay more rapidly with distance to compensate.
The paper is organized as follows. The next section defines the model and estimator. Section (ref) states our assumptions and the main result. We provide primitive sufficient conditions for neighborhood stability in (ref). Section (ref) presents results from a simulation study of the emmenegger2025treatment network interference model. Finally, (ref) concludes. All proofs are relegated to the appendix.
Let $[n] = \{1, \ldots, n\}$ be the set of observations and $\rho_n\colon [n] \times [n] \rightarrow \mathbb{R}_+$ be a semi-metric, so that $\rho_n(i,j)$ represents the distance between $i,j \in [n]$. For example, observations can be located in Euclidean space with $\rho_n(i,j)$ denoting the Euclidean distance between the locations of $i$ and $j$. Alternatively, observations may be connected by a social network with $\rho_n(i,j)$ denoting the shortest path distance between $i$ and $j$. Abusing notation, for $S_1, S_2 \subseteq [n]$, let $\rho_n(S_1,S_2)=\min \{\rho_n(i,j)\colon i\in S_1, j\in S_2 \}$. For $r \in \mathbb{R}_+$, let $\mathcal{N}(i,r) = \{j \in [n]\colon \rho_n(i,j) \leq r\}$ denote the {\em $r$-neighborhood} of observation $i$.
The data consists of the triangular array $\bm{Z} = \{Z_i\}_{i=1}^n$ where $Z_i \in \mathcal{Z} \subseteq \mathbb{R}^{d_z}$ for all $i \in [n]$. Let $\mathcal{G}$ be a set of functions $g\colon \mathcal{Z} \rightarrow \mathbb{R}^{d_g}$. We consider the moment model
studied, for example, by chen2022debiased and chernozhukov2018double. Here $\theta \in \mathbb{R}^p$ is the parameter of interest, $g$ a nuisance function, $\psi\colon \mathcal{Z} \times \mathcal{G} \rightarrow \mathbb{R}^{p\times p}$, and $\nu\colon \mathcal{Z} \times \mathcal{G} \rightarrow \mathbb{R}^p$. We suppose there exist “true values” $g_0 \in \mathcal{G}$ and $\theta_0 \in \mathbb{R}^p$ such that
Let $\hat{g}$ denote the machine learner, an estimate of $g_0$ trained on $\bm{Z}$. Formally, $\hat{g}$ is a mapping from $\mathcal{Z}^n$ to $\mathcal{G}$, but we suppress the dependence of $\hat{g}$ on the training set $\bm{Z}$ and simply treat $\hat{g}$ as an element of $\mathcal{G}$, writing $\hat{g}(Z_i)$. Let $\hat\theta$ be the two-step method of moments estimator that satisfies
Our objective is to provide lower-level conditions for the following version of (ref), which we will refer to as “SE.” For any $f \in \{\psi, \nu\}$, denote by $f_{k,l}$ the $kl$-th component of $f$. For the case of $f=\nu$, which is vector- rather than matrix-valued, we define $f_{k,l}=0$ whenever $k\neq l$. Let $F_{k,l}(i, g) = {\bf E}[f_{k,l}(Z_i, g)]$. SE corresponds to
For i.i.d.\ data, Theorem 1 of chen2022debiased or Theorem 3.1 of chernozhukov2018double demonstrate that $\hat\theta$ is approximately normal under SE ((A.13) of the latter reference is analogous to SE). These results can be extended to dependent data given an appropriate CLT brown2024inference,gilbert2024causal,leung2025graph, so our focus is solely on verifying SE.
Let $\lVert\cdot\rVert$ denote the Euclidean norm, and for any function $g$, let $\lVertg\rVert_{2,i} = {\bf E}[\lVertg(Z_i)\rVert^2]^{1/2}$.
Let $\tilde{\bm{Z}} = \{\tilde{Z}_i\}_{i=1}^n$ denote an independent copy of $\bm{Z}$ and $\hat{g}^{(-i),0}$ be the machine learner trained on $\{Z_1, \ldots, Z_{i-1}, \tilde{Z}_i, Z_{i+1}, \ldots, Z_n\}$ rather than $\bm{Z}$. The latter is a perturbed version of the data that replaces $Z_i$ with an independent copy. Let $\bm{Z}^* = \{Z_i^*\}_{i=1}^n$ denote another copy of $\bm{Z}$ independent of both $\tilde{\bm{Z}}$ and $\bm{Z}$. chen2022debiased impose the following stability condition on the moments and $\hat{g}$ (see their Lemma 2):
When the machine learner is trained on i.i.d.\ data using cross-fitting, $Z_i$ and $\hat{g}$ are independent, and SE is straightforward to verify. Stability essentially weakens full independence to requiring that $Z_i$ has a negligible influence on $\hat{g}$.
To adapt (ref) to dependent data, we replace a {\em set} of observations in a slowly growing neighborhood of $Z_i$ with an independent copy. Formally, let $\{r_n\}_{n\in\mathbb{N}} \subseteq \mathbb{R}_+$ be a sequence of {\bf neighborhood radii}. Construct from $\bm{Z}$ a new dataset by replacing $\{Z_k\colon k \in \mathcal{N}(i,r_n) \cup \mathcal{N}(j,r_n)\}$ with $\{\tilde{Z}_k\colon k \in \mathcal{N}(i,r_n) \cup \mathcal{N}(j,r_n)\}$, and let $\hat{g}^{(-i,-j)}$ be the machine learner trained on this data, leaving its dependence on $r_n$ implicit for economy of notation. For the case $i=j$ we abbreviate $\hat{g}^{(-i)} \equiv \hat{g}^{(-i,-j)}$.
This strengthens (ref) to demand invariance of $\hat{g}$ to resampling larger blocks of observations. The faster $r_n$ diverges, the stronger the requirement. How fast we require $r_n$ to diverge will depend on the degree of dependence in $\bm{Z}$, as we will see shortly. We verify (ref) for various machine learners in (ref).
Let $\bar{N}_n(r) = n^{-1} \sum_{i=1}^n \lvert\mathcal{N}(i,r)\rvert$ be the average $r$-neighborhood size.
This requires the machine learner to converge at a sufficiently fast rate. For the case of $q=2$, with standard machine learners and weakly dependent data, we expect
which corresponds to the usual $n^{-1/4}$ rate of convergence chernozhukov2018double.\footnote{Note that (ref) contains $\hat{g}$ rather than $\hat{g}^{(-i)}$. The analogous statement for $\hat{g}^{(-i)}$ should be easier to verify since the data on which it is trained exhibits a greater degree of independence relative to $\bm{Z}$.} Then in order for (ref) to not be unduly restrictive, we require
which is a restriction on the metric space and $r_n$. Choosing $r_n = O(1)$ would satisfy (ref) for most spaces of interest, but our final assumption below may require $r_n$ to diverge. As we will discuss in the examples below, this together with (ref) amounts to requiring sufficiently fast decay of a certain measure of dependence.
For any $\sigma$-fields $\mathcal{A},\mathcal{B}$, define the $\beta$-mixing coefficient
where the supremum is taken over all $J,K\geq 1$ and partitions $\{A_j\}_{j=1}^J$ and $\{B_k\}_{k=1}^K$ of the sample space such that $A_j \in \mathcal{A}$ and $B_j \in \mathcal{B}$. For any $S \subseteq [n]$, let $\sigma(S)=\sigma(\{Z_i\}_{i\in S})$ be the $\sigma$-field generated by $\{Z_i\}_{i\in S}$. Define
This measures the dependence between one observation and sets of observations at least distance $r$ away. We require sufficiently fast decay with respect to $r$ in the following sense.
We now state the main result.
To prove the theorem, we adapt the proof of the analogous i.i.d.\ result, Lemma 2 of chen2022debiased. We make use of a powerful $\beta$-mixing coupling lemma due to berbee1979random ((ref)), which allows us to reduce certain expressions to the independent case at the cost of an approximation error involving $\beta_n$.
This section shows that existing results on the stability of machine learners can be adapted to verify neighborhood stability. We first note that (ref) is a restriction on both the moment function and machine learner. The next proposition shows that, under smoothness conditions on the moment functions, specifically a strengthened version of (ref)(b), (ref) is a consequence of a more primitive neighborhood stability condition imposed directly on the machine learner.
This is analogous to Corollary 4 of chen2022debiased, and the proof is immediate. The remainder of this section provides primitive sufficient conditions for (ref).
Suppose $\mathcal{G}$ is a reproducing kernel Hilbert space with kernel $\kappa$ and norm $\lVert\cdot\rVert_\kappa$ berlinet2011reproducing. Let
where $\ell\colon \mathcal{G} \times \mathcal{Z} \rightarrow \mathbb{R}_+$ is the loss function and $\lambda \geq 0$ the penalty parameter. This setup includes bounded SVM regression, soft margin SVM classification, and ridge regression.
Suppose all $g \in \mathcal{G}$ share a common codomain $\text{Co}(\mathcal{G})$. The following result is a simple application of Theorem 22 of bousquet2002stability.
The result implies (ref) when
In the standard regime where $\lambda = O(n^{-1/2})$, which results in a nondegenerate limiting distribution knight2000asymptotics, this is similar to, though stronger than, (ref). See Examples (ref) and (ref).
Suppose the machine learner is parameterized by $\theta \in \Theta \subseteq \mathbb{R}^{d_\theta}$ so that $g(z) \equiv g(z; \theta)$, and the estimated learner is $\hat{g}(z) \equiv g(z; \hat\theta)$. To obtain $\hat\theta$, we choose an objective $\psi\colon \Theta \times \mathcal{Z} \rightarrow \mathbb{R}$ and seek to solve
To compute the solution, we consider the use of stochastic gradient descent (SGD) where we initialize $\hat\theta_0 = 0$ and then apply the updating rule
sequentially for all $t \in [n-1]$. Here $\alpha_{t+1}$ is the step size and $\psi'$ is the derivative of $\psi$ with respect to its first argument. The estimator is $\hat\theta = \hat\theta_{n-1}$.
hardt2016train and kissel2023black assume $\psi$ satisifes the following smoothness and convexity conditions.
Define $\hat\theta^{(-i,-j)}$ such that $\hat{g}^{(-i,-j)}(\cdot) = g(\cdot,\hat\theta^{(-i,-j)})$. The following result is a simple application of Proposition 5.1 of kissel2023black.
Suppose $L/\beta$ is nondegenerate. Then the right-hand side of (ref) is $o(n^{-1/2})$ if
Since the left-hand side is typically non-degenerate, the learning rate $a$ must exceed 0.5. With stronger dependence ($r_n$ diverges faster) or larger neighborhood sizes ($\mathcal{N}(i,r)$ grows faster with $r$), the condition requires a faster learning rate $a$. In the case of $a$ close to one, the requirement is similar to, though stronger than, (ref).
Neighborhood stability (ref) follows under appropriate smoothness conditions on $g(\cdot; z)$, as illustrated in the next example.
We lastly consider machine learners that average the predictions of “base” learners trained on resampled data, for example random forests. Following the setup of chen2022debiased, we draw $m$ observations from $[n]$, either with or without replacement, independently $B$ times to obtain $\{\bm{Z}_{[m]}^b\}_{b=1}^B$, which we will to refer to as “subsamples.” Let $\hat{h}(\cdot; \bm{Z}_{[m]}^b)\colon \mathcal{Z} \rightarrow \mathcal{G}$ be a base machine learner trained on $\bm{Z}_{[m]}^b$. The bagged learner is
Let $\bm{Z}_{[m],(-i,-j)}^b$ be defined by replacing each element of $\{Z_k\colon k \in \mathcal{N}(i,r_n) \cup \mathcal{N}(j,r_n)\}$ in $\bm{Z}_{[m]}^b$ with the corresponding element of $\{\tilde{Z}_k\colon k \in \mathcal{N}(i,r_n) \cup \mathcal{N}(j,r_n)\}$. The following result is a simple extension of Theorem 5 of chen2022debiased.
The second part of the parameters condition holds with a sufficiently large number of subsamples $B$. The first part restricts the subsample size $m$ in a manner that depends on how fast dependence decays relative to neighborhood growth rates. The condition allows $m$ to grow polynomially with $n$ under the following strengthened version of (ref):
for some $c > 0$. Then the condition holds with $m = o(n^c)$. The strengthened condition is satisfied if the mixing rate decays fast enough relative to the neighborhood growth rate. In the case of spatial data ((ref)), this holds if $\alpha^{-1} a \delta < 0.5$. In the case of network data ((ref)), this holds if $\alpha^{-1} a d < 0.5$.
We replicate the design used in emmenegger2025treatment, which features a potential outcomes model with network interference. We find that computing their estimator without cross-fitting improves the bias and standard deviation. We also find that increasing the density of the network used in their design can result in extremely small training fold sizes when cross-fitting.
In this design, observations are units connected through a network which is considered non-random or conditioned upon. Then $\rho_n$ is shortest path distance, and $\mathcal{N}(i,r)$ is the set of units is the set of units no more than path distance $r$ from $i$. For each unit we observe a scalar outcome $Y_i$, binary treatment $W_i$, covariate $C_i$, and network feature $X_i$. Outcomes are given by
where ${\bf E}[\varepsilon_i \mid W_i, X_i, C_i] = 0$. The network features are functions of treatments and covariates of units connected to $i$, generating what the authors refer to as “$X$-spillovers.” Treatments are unconfounded in that $W_i \perp\!\!\!\perp \varepsilon_i \mid X_i, C_i$, and the estimand of interest is the average treatment effect
We replicate the data-generating process in section 3.1 of emmenegger2025treatment that uses an Erd\H{o}s-R\'{e}nyi network. See (ref) for specifics of the model. The network formation model draws links between pairs in an i.i.d.\ fashion with linking probability $\Delta/n$ with $\Delta=3$. This means $\Delta$ is the limiting expected degree (number of connections involving a given unit), which measures network density.
Due to the network features, the outcomes of units that share a common network neighbor are correlated since both depend upon the common neighbor's treatment. As discussed in (ref), this induces local dependence, which satisfies $\beta$-mixing. We can take $r_n=3$ to satisfy (ref) since units are independent if they are at least path distance 3 apart and therefore lack a common neighbor.
As in emmenegger2025treatment we estimate $\theta_0$ using the doubly robust estimator in (ref) using random forests to approximate the propensity score and outcome regressions. We replicate the parameters of the learner specified in their section 3.1, which uses standard bootstrapped trees; see (ref) for details. We also compare with subsampled trees satisfying the setup of (ref). Since $r_n = O(1)$, we treat $N^*(r_n)$ in (ref) as fixed and take the subsample size to be $m = \min\{10\tilde n^{1/3}, \tilde n\}$ to satisfy assumption (a) of (ref) where $\tilde n$ is the size of the training set.
We report results with and without cross-fitting. In the latter case, the random forests are trained on the entire dataset with the outcome regressions ${\bf E}[Y_i \mid W_i=w, X_i, C_i]$ estimated separately for $w \in \{0,1\}$. Hence $\sum_{i=1}^n W_i$ is the sample size used in the outcome regression for $w=1$. The data-generating process is calibrated such that this is about $n/2$.
Due to local dependence, emmenegger2025treatment propose the following modification of cross-fitting to ensure independence between evaluation and training folds. First, they randomly partition $[n]$ into $K=5$ folds. For a given fold $k$ used to compute the estimator, the training fold for the random forest is the subset of the units not in fold $k$ that are not connected to a unit in $k$ and do not share a neighbor in common with any unit in $k$ (see their eq.\ (6)). The estimates are averaged across the $K$ folds to obtain $\hat\theta$. emmenegger2025treatment further reduce the randomness of this procedure by repeating the procedure $B$ times and taking $\hat\theta$ to be the median result, but we only consider the standard case of $B=1$.
(ref) presents the results. Whether using bootstrapped or subsampled trees, the emmenegger2025treatment estimator sees a reduction in bias and variance when using the full sample instead of cross-fitting. We also find that bootstrapped trees perform better in terms of bias and variance, particularly when cross-fitting. Such trees correspond to a subsample size $m$ (in the notation of (ref)) equal to $n$, which violates assumption (a) of the proposition. This indicates that, while the condition is sufficient, it may not be necessary.
(ref) shows the average size of the training fold when cross-fitting for different values of $\Delta$ and $K$. We see that for the base design in (ref), the fold size is small, only about 10 percent of the full sample, which explains why cross-fitting results in higher variance. When we increase the expected degree to 8, the fold size becomes negligible. This is because the network is denser, and almost all units in the training fold are less than distance 3 from the evaluation fold and end up excluded from training. When the expected degree is 5 but we reduce the number of folds to two, we again obtain close to negligible sample sizes. Despite the network being sparser, half of all units are in the evaluation fold, so the vast majority of units in the training fold are again less than distance 3 from the evaluation fold.
When cross-fitting with dependent data, existing proposals reduce correlation between training and evaluation folds by excluding from training observations too close to the evaluation set. This can result in a large number of excluded observations, particularly with network data. We establish the validity of DML methods using dependent data without cross-fitting when the machine learner is neighborhood-stable. This requires the learner to be robust to resampling data within a slowly growing neighborhood of any observation.
When dependence is weaker in the sense that the $\beta$-mixing coefficient decays faster with distance, neighborhood stability imposes weaker restrictions on the learner in that the neighborhood can be smaller relative to the sample size. CLTs for dependent data typically require mixing coefficients to decay with distance $r$ sufficiently fast relative to the sizes of $r$-neighborhoods. We show that, under similar requirements, neighborhood stability holds for several classes of learners, including regularized $M$-estimators and bagged learners.