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.
50,385 characters · 12 sections · 24 citation commands
-2.5cm A Variance-Based Test for Heterogeneous Treatment Effects
\bgroup \let\footnoterule\relax
\thispagestyle{empty}
\egroup \setcounter{page}{1}
The analysis of causal effects has traditionally centered on the Average Treatment Effect (ATE), which summarizes the mean impact of a policy or intervention across an entire population. While nonparametric estimators for the ATE with valid statistical inference are now well-established under the unconfoundedness assumption robins1994estimation,chernozhukov2018double, the ATE often masks substantial heterogeneity in individual responses. Recognizing this heterogeneity is crucial for understanding the underlying mechanisms of a treatment and for designing optimal policies that tailor interventions to specific subpopulations heckman1997making, athey2017state.
Recent advances in causal machine learning have facilitated flexible estimation of the Conditional Average Treatment Effect (CATE) function, \(\tau(x) = \mathbb{E}[Y(1) - Y(0) | X=x]\), even in high-dimensional settings wager2018estimation, nie2021quasi. However, obtaining valid statistical inference for the full CATE function remains a formidable challenge. The complexity of modern machine learning algorithms often precludes the use of classical empirical process theory, and the regularization bias required for estimation makes formally testing hypotheses about the shape of \(\tau(x)\) difficult. Consequently, researchers often face a trade-off between the robust inference available for the ATE and the granular, yet often unstable, characterization of the full CATE curve.
To bridge this gap, we propose a robust hypothesis test for the existence of heterogeneous treatment effects. Rather than attempting to estimate the shape of the heterogeneity immediately, we ask a preliminary question: Is the treatment effect constant across subpopulations defined by covariates? We answer this by conducting inference on a single scalar parameter, the variance of the CATE, \(\theta_0 = \operatorname{Var}(\tau(X))\). This parameter serves as a natural omnibus measure. If \(\theta_0 = 0\), the effects are homogeneous almost surely, and if \(\theta_0 > 0\), relevant heterogeneity exists, justifying further granular investigation.
Developing a valid test for \(\theta_0\) presents a theoretical challenge at the intersection of causal inference and machine learning. To robustly estimate this variance, one must rely on the doubly robust pseudo-outcome, which is the influence function of the ATE. However, testing the null hypothesis of homogeneity using these pseudo-outcomes introduces an impasse characterized by two problems, null degeneracy and the breakdown of Neyman orthogonality.
First, under the null hypothesis, the true parameter lies on the boundary of the parameter space, and the CATE collapses to the ATE. In this boundary case, the true variance components of the pseudo-outcome become identical, and the influence function of the standard variance difference degenerates to zero almost surely. If the test statistic is computed on a single full sample, the asymptotic variance collapses and standard Gaussian inference breaks down. The modern semiparametric resolution to such degeneracy is to decouple the empirical processes of the variance components via sample-splitting williamson2023general.
Second, resolving null degeneracy via standard sample-splitting inadvertently causes a breakdown of Neyman orthogonality, which is a key ingredient for developing \(\sqrt{n}\)-consistent semiparametric estimators chernozhukov2018double. By the Law of Total Variance, \(\theta_0\) is identified as the difference between the total and residual variances of the pseudo-outcome. While the pseudo-outcome itself is doubly robust, its squared loss is not, introducing \(O_P(n^{-1/4})\) first-order regularization biases into both variance components. These biases can cancel each other when both components are constructed using the exact same nuisance estimators. Standard splitting destroys this symmetry. By evaluating the components on separate folds to decouple their empirical processes, it forces the use of independently trained machine learning models, and thus, their biases fail to cancel. When the test statistic is scaled by \(\sqrt{n}\), this uncancelled residual error diverges to infinity, invalidating asymptotic inference.
To resolve this methodological impasse, a valid test must evaluate the variance components on disjoint observations to resolve degeneracy, while applying identical nuisance estimators to preserve bias cancellation. We achieve this via a novel Intra-Fold Sample-Splitting algorithm. We partition the data into \(K\) main folds and train a single set of nuisance functions out-of-fold. We randomly bisect each in-fold dataset into two mutually disjoint halves and compute the total variance exclusively on the first half and the residual variance on the second, applying the exact same out-of-fold nuisance estimators to both. This paired structure ensures a strictly positive asymptotic variance under the null and cancels out the non-orthogonal squared biases. We formally prove this procedure yields consistent and asymptotically normal estimators, and guarantees valid Type I error control on the boundary.
Our work contributes to the growing literature on testing for treatment effect heterogeneity. Existing methods largely fall into two categories, projection-based tests and distributional tests. A prominent strand of literature focuses on testing whether the projection of the CATE onto a specific set of basis functions of covariates is zero. crump2008nonparametric propose a nonparametric test based on sieve estimation, while semenova2021debiased develop a Double/Debiased Machine Learning (DML) inference framework for the coefficients of a linear projection of the CATE. Theoretically, these projection-based methods are consistent against general nonlinear alternatives provided the number of basis functions grows sufficiently with the sample size. However, in practice, this approach faces a fundamental trade-off between approximation error and statistical power. Testing the joint significance of a high-dimensional vector of coefficients consumes degrees of freedom, diluting statistical power. Conversely, specifying a parsimonious basis to maximize power risks inconsistency if the true heterogeneity is orthogonal to the chosen subspace. In contrast, our CATE Variance Test targets a single scalar parameter. Because \(\theta_0 = 0\) is a necessary and sufficient condition for a constant CATE, our test remains consistent against any deviation from the null without incurring the power penalty associated with high-dimensional coefficient testing.
A second strand of literature focuses on distributional effects, testing for differences in the marginal distributions or variances of potential outcomes ding2016randomization, chung2021permutation. While observing a difference in marginal distributions of potential outcomes implies the existence of individual treatment effect heterogeneity, it is not a direct test of moderation by observables. It is possible for individual effects to vary while the conditional average effect \(\tau(x)\) remains constant. Our test specifically isolates the heterogeneity explained by covariates, making it directly relevant for policy evaluation and design. Beyond testing for treatment effect heterogeneity, our algorithm offers a generalizable framework for conducting valid hypothesis testing on nonlinear transformations of doubly robust scores.
The remainder of the paper is organized as follows. Section (ref) establishes the econometric framework and the identification of the target parameter via pseudo-outcomes. Section (ref) formalizes the theoretical tension between null degeneracy and Neyman orthogonality, introduces our Intra-Fold Sample-Split algorithm and establishes its asymptotic properties. Section (ref) presents Monte Carlo simulation results comparing our test to existing projection-based alternatives. Section (ref) applies the test to empirical data from the NSW job training program, and Section (ref) concludes. All proofs are collected in the Appendix.
In this section, we define the causal parameters of interest, state the assumptions required for identification and inference, and derive the variance decomposition that forms the basis of our test statistic.
We follow the potential outcomes framework rubin1974estimating. We observe a random sample of \(n\) independent and identically distributed units \(O_i = (Y_i, D_i, X_i)\) for \(i=1, \dots, n\), drawn from an unknown distribution \(P_0\). Here, \(D_i \in \{0, 1\}\) is a binary treatment indicator, \(X_i \in \mathcal{X} \subset \mathbb{R}^p\) is a vector of covariates, and \(Y_i \in \mathbb{R}\) is the observed outcome. Let \(Y_i(1)\) and \(Y_i(0)\) denote the potential outcomes under treatment and control, respectively. The observed outcome relates to the potential outcomes via the consistency condition \(Y_i = D_i Y_i(1) + (1-D_i)Y_i(0)\). The fundamental problem of causal inference is that for any unit \(i\), we observe only one of the two potential outcomes. Our primary focus is the CATE, defined as \[ \tau_0(x) = \mathbb{E}[Y_i(1) - Y_i(0) | X_i = x]. \] We also define the ATE, denoted by \(\tau_{\mathrm{ATE}} = \mathbb{E}[\tau_0(X_i)]\). To facilitate identification, we define the nuisance functions \(\mu_0(d, x) = \mathbb{E}[Y_i | D_i=d, X_i=x]\) for \(d \in \{0, 1\}\) representing the conditional outcome means, and \(e_0(x) = P(D_i=1 | X_i=x)\) representing the propensity score.
We invoke the standard assumptions for causal identification in observational studies rosenbaum1983central, alongside regularity conditions required for valid asymptotic inference.
Assumptions (ref) and (ref) allow for the identification of the CATE function \(\tau_0(x) = \mu_0(1, x) - \mu_0(0, x)\). Assumption (ref)(i) ensures finite moments needed for the central limit theorem and for the empirical-process arguments underlying double machine learning. Condition (ii) directly guarantees that the asymptotic variance of our test statistic is bounded away from zero on the boundary of the parameter space, eliminating the pathological degeneracy that would otherwise arise under the null hypothesis of homogeneity.
We investigate whether the treatment effect is constant across the population defined by \(X\). Formally, we define the CATE Variance parameter: \[ \theta_0 = \operatorname{Var}(\tau_0(X_i)). \] The variance serves as an omnibus measure of heterogeneity. If \(\theta_0 = 0\), the CATE is constant almost surely (i.e., \(\tau_0(x) = \tau_{\mathrm{ATE}}\) for all \(x\)). If \(\theta_0 > 0\), there exists variation in the treatment effect explained by the covariates. Accordingly, we test the null hypothesis of homogeneity against the one-sided alternative of heterogeneity:
Since \(\theta_0\) is non-negative, the null hypothesis lies on the boundary of the parameter space. We address the inferential implications of this boundary condition in Section (ref).
A direct estimator of \(\operatorname{Var}(\tau_0(X))\) based on a plug-in estimate of the function \(\hat{\tau}(x)\) would suffer from first-order regularization bias, particularly when \(X\) is high-dimensional. To address this, we utilize a doubly robust pseudo-outcome, also known as the uncentered influence function for the ATE. Define the pseudo-outcome \(\psi(O_i)\) as
This pseudo-outcome possesses two critical properties. First, it is an unbiased signal of the CATE \[ \mathbb{E}[\psi(O_i) | X_i] = \tau_0(X_i), \] which also implies \(\mathbb{E}[\psi(O_i)] = \tau_{\mathrm{ATE}}\). Second, it allows us to identify \(\theta_0\) through a variance decomposition. By the Law of Total Variance applied to \(\psi(O_i)\), we have \[ \operatorname{Var}(\psi(O_i)) = \operatorname{Var}(\mathbb{E}[\psi(O_i)|X_i]) + \mathbb{E}[\operatorname{Var}(\psi(O_i)|X_i)]. \] Substituting the conditional expectation with \(\tau_0(X_i)\), we can rearrange this to identify the CATE variance \[ \theta_0 = \operatorname{Var}(\psi(O_i)) - \mathbb{E}[(\psi(O_i) - \tau_0(X_i))^2]. \] Or, expressed in terms of Mean Squared Error which facilitates our estimation strategy
Equation (ref) provides the identification result for our test. It expresses the CATE variance as the difference between \(V_{\mathrm{tot}}\), the MSE of the best constant predictor of the pseudo-outcome (\(\tau_{\mathrm{ATE}}\)), and \(V_{\mathrm{res}}\), the MSE of the best conditional predictor (\(\tau_0(X)\)).
In this section, we develop a formal hypothesis test for the presence of heterogeneous treatment effects. Having identified the CATE variance, \(\theta_0 = \operatorname{Var}(\tau_0(X))\), as our target parameter in Section (ref), we test the hypotheses in Equation (ref). This test leverages the identification result derived in Equation (ref). While estimating \(\theta_0\) fits within the general framework of semiparametric inference, the null hypothesis poses a unique theoretical challenge known as null-degeneracy. Below, we derive the influence function for \(\theta_0\), analyze its properties, and detail the algorithm for the hypothesis test.
We first derive the influence function for \(\theta_0\).
Proof. See Appendix.
Based on Proposition (ref), a standard "one-step" efficient estimator can be constructed by solving the empirical equation \(n^{-1} \sum \hat\phi_{\theta}(O_i) = 0\), where \(\hat\phi_{\theta}(O_i)\) is obtained by plugging in the estimated nuisance parameters. Under the alternative hypothesis, \(H_1: \theta_0 > 0\), standard semiparametric theory guarantees that such an estimator is \(\sqrt{n}\)-consistent and asymptotically normal \[ \sqrt{n}(\hat\theta - \theta_0) \xrightarrow{d} \mathcal{N}(0, \operatorname{Var}(\phi_\theta)), \] provided that the nuisance parameters converge at sufficiently fast rates, typically \(n^{-1/4}\) chernozhukov2018double. This allows us to employ modern machine learning methods to estimate \(\mu_0\) and \(e_0\), and plug in the pseudo-outcome \(\psi(O_i)\). Constructing \(\hat\theta\) also requires feasible estimators for \(\tau_{\mathrm{ATE}}\) and \(\tau_0(X)\), which we review in the next section.
Existing strategies for estimating \(\tau_0(x)\) and \(\tau_{\mathrm{ATE}}\) largely fall into two categories, the T-learner and the DR-learner. The T-learner estimates the conditional means \(\mu_0(1, x)\) and \(\mu_0(0, x)\) and computes their differences to obtain an estimator for the treatment effect. For the ATE, a T-learner is \(n^{-1}\sum_{i=1}^n (\hat\mu(1, X_i) - \hat\mu(0, X_i))\), and for the CATE, a T-learner is \(\hat\mu(1, X_i) - \hat\mu(0, X_i)\). By the triangle inequality, the \(L_2\) error of the T-learner is bounded by the errors of the baseline outcome models. Therefore, provided the nuisance estimators \(\hat{\mu}(1, \cdot)\) and \(\hat{\mu}(0, \cdot)\) satisfy the \(o_P(n^{-1/4})\) rate, we can show that the T-learner also satisfies this rate and can be applied in our algorithm. However, as noted by kunzel2019metalearners, T-learners can suffer from regularization bias, particularly when the CATE function is sparser than the baseline outcome functions or when there is poor overlap between treatment groups.
The DR-learner, on the other hand, treats ATE and CATE estimation as a direct regression of the pseudo-outcome \(\psi(O_i)\) on the covariates. The DR-learner for the ATE is \(n^{-1}\sum_{i=1}^n \hat\psi(O_i)\), which is also known as the Augmented Inverse Propensity Weighting (AIPW) estimator. Because the pseudo-outcome is Neyman orthogonal, \(\hat{\tau}_{\mathrm{ATE}}\) is \(\sqrt{n}\)-consistent provided the product of the \(L_2\) estimation errors for the propensity score and outcome mean vanishes at an \(o_P(n^{-1/2})\) rate robins1994estimation, chernozhukov2018double. For the CATE, it is the minimizer of the mean squared error \(\sum_{i=1}^n (\hat{\psi}(O_i) - f(X_i))^2\). Its estimation error is bounded by the oracle smoothing error of the CATE plus the product of the nuisance errors \(\|\hat{e} - e_0\|_{P,2} \times \|\hat{\mu} - \mu_0\|_{P,2}\) kennedy2023towards, where \(\|\cdot\|_{P,2}\) denotes the \(L_2(P_0)\) norm. This imparts a "double robustness of rates." Even if the baseline outcome model \(\hat{\mu}\) converges at a rate slower than \(n^{-1/4}\) due to complex confounding, the DR-learner can still achieve the requisite \(o_P(n^{-1/4})\) rate, provided the propensity score converges sufficiently fast and the true CATE is sufficiently smooth.
In this paper, we adopt the DR-learner for both ATE and CATE estimation. It often yields more stable estimates than differencing two regression functions, and the convergence rate depends on the product of nuisance errors, making it robust to misspecification of the nuisance models. However, simply plugging the nuisance estimators into a standard full-sample or cross-fitting empirical analogue of \(\theta_0\) fails to yield valid inference. We formalize this fundamental breakdown of Neyman orthogonality in the next section.
To develop a valid semiparametric test based on the variance of CATE, the first hurdle is the problem of null degeneracy. Under the null hypothesis of homogeneity \(H_0: \theta_0 = 0\), the true CATE is constant almost surely, i.e., \(\tau_0(X_i) = \tau_{\mathrm{ATE}}\). Consequently, the true total and residual losses are identical, and their corresponding influence functions coincide perfectly. If one computes the empirical analogues of \(V_{\mathrm{tot}}\) and \(V_{\mathrm{res}}\) using the same sample of observations, the empirical processes become perfectly correlated, and the asymptotic variance of their difference collapses to zero. This degeneracy violates the regularity conditions required for standard Gaussian approximations and destroys the size calibration of the test: rather than attaining its nominal level, the same-sample statistic degenerates and becomes severely conservative (Appendix (ref)). williamson2023general suggest that this degeneracy can be resolved by evaluating the components on disjoint subsets of the data via sample-splitting.
However, resolving null degeneracy via standard sample-splitting breaks down the Neyman orthogonality of the influence function in Proposition (ref). To formalize this, consider the pathwise Gâteaux derivative of the expected squared residual loss, \(\mathbb{E}[(\psi - \tau_0(X))^2]\), with respect to the propensity score \(e(x)\). The expected first-order bias depends on the cross-term conditional on \(X_i\), \[ \mathbb{E}\left[ 2(\psi(O_i) - \tau_0(X_i)) \frac{\partial \psi}{\partial e}(O_i) \bigg| X_i \right]. \] Substituting the residual error \[ \psi(O_i) - \tau_0(X_i) = \frac{D_i(Y_i-\mu_0(1, X_i))}{e_0(X_i)} - \frac{(1-D_i)(Y_i-\mu_0(0, X_i))}{1-e_0(X_i)} \] and its partial derivative \[ \frac{\partial \psi}{\partial e}(O_i) = - \frac{D_i(Y_i-\mu_0(1, X_i))}{e_0(X_i)^2} - \frac{(1-D_i)(Y_i-\mu_0(0, X_i))}{(1-e_0(X_i))^2}, \] the cross-products strictly vanish since the treatment indicator satisfies \(D_i(1-D_i) = 0\). Using the unconfoundedness assumption to replace the expected squared residual outcomes with the true conditional variances \(\sigma_1^2(X_i)\) and \(\sigma_0^2(X_i)\), \[ \mathbb{E}\left[ 2(\psi(O_i) - \tau_0(X_i)) \frac{\partial \psi}{\partial e}(O_i) \bigg| X_i \right] = 2 \left( - \frac{\sigma_1^2(X_i)}{e_0(X_i)^2} + \frac{\sigma_0^2(X_i)}{(1-e_0(X_i))^2} \right) \equiv g(X_i). \] Crucially, this derivative \(g(X_i)\) is generally non-zero. Because this derivative does not vanish, the squared pseudo-outcome is not Neyman orthogonal. Any plug-in estimator for the residual variance \(V_{\mathrm{res}}\) is therefore contaminated by a first-order regularization bias of order \(O_P(\|\hat{e} - e_0\|_{P,2})\). The estimator for the total variance, \(V_{\mathrm{tot}} = \mathbb{E}[(\psi - \tau_{\mathrm{ATE}})^2]\), suffers from the same non-orthogonal bias \(g(X_i)\).
The target parameter \(\theta_0\) remains \(\sqrt{n}\)-consistent only because of an exact algebraic cancellation. If both variance components are evaluated using the exact same nuisance estimators, their respective first-order biases \(g(X_i)\) are mathematically identical and cancel one another when taking the difference \(V_{\mathrm{tot}} - V_{\mathrm{res}}\). Standard sample-splitting, which evaluates the two variance components on different data folds, structurally destroys this delicate symmetry by forcing the use of independently trained machine learning nuisance estimators (e.g., evaluating \(V_{\mathrm{tot}}\) with an out-of-fold propensity score \(\hat{e}_{\mathrm{odd}}\) and \(V_{\mathrm{res}}\) with \(\hat{e}_{\mathrm{even}}\)). Because these independently trained nuisance estimators differ in finite samples, their induced non-orthogonal biases no longer match. The uncancelled first-order bias in the split-sample estimator becomes approximately \[ \text{Bias}(\hat{\theta}_{\mathrm{split}}) \approx \int g(X) \big( \hat{e}_{\mathrm{odd}}(X) - \hat{e}_{\mathrm{even}}(X) \big) dP_0(X). \] Standard rates only guarantee that independently trained nuisance estimators differ by \(o_P(n^{-1/4})\), so the \(\sqrt{n}\)-scaled test statistic inherits a residual bias of order \(o_P(n^{1/4})\) — a quantity that need not converge to zero, invalidating asymptotic inference.
To resolve the methodological impasse formalized in Section (ref), a valid testing procedure must simultaneously evaluate the total and residual variance components on strictly disjoint sets of observations and construct these components using the same nuisance estimators to preserve the algebraic cancellation of the squared pseudo-outcomes, thereby restoring Neyman orthogonality.
We achieve these requirements via a novel Intra-Fold Sample-Split (IF-SS) algorithm, detailed in Algorithm (ref). Instead of splitting the evaluation of the variance components across entirely different main folds, our algorithm introduces an internal data partition. We first partition the data into \(K\) main folds and train a single set of nuisance functions on the out-of-fold data. We randomly bisect each in-fold evaluation dataset into two mutually disjoint halves. We compute the total variance exclusively on the first half and the residual variance exclusively on the second half, applying the identically trained out-of-fold nuisance estimators to both.
To establish the asymptotic validity of Algorithm (ref), we impose regularity conditions on the estimators used. We maintain Assumptions (ref)--(ref) from Section (ref) and further introduce the following regularity conditions on nuisance estimators.
Because Algorithm (ref) algebraically cancels the non-orthogonal bias, the only remaining estimation errors depend strictly on the doubly robust linear pseudo-outcome terms and the Mean Squared Error of the CATE estimator itself (\(\|\hat{\tau}_k - \tau_0\|_{P,2}^2\)). Provided Assumption (ref) holds, these remaining errors rigorously vanish at an \(o_P(n^{-1/2})\) rate. We formalize the asymptotic validity of this test in Theorem (ref).
Proof. See Appendix.
We evaluate the performance of the proposed test using Monte Carlo simulations. In all designs, we generate \(n \in \{250, 500, 1000, 2000\}\) independent and identically distributed observations \(O_i = (Y_i, D_i, X_i)\) where \(X_i \in \mathbb{R}^p\). The outcome follows a common structural model \[ Y_i = \mu_0(X_i) + D_i \cdot \tau(X_i) + \varepsilon_i, \] where \(\mu_0(x)\) is the baseline outcome function, \(\tau(x)\) is the CATE, and \(\varepsilon_i \sim N(0, 1)\). The covariates are drawn from a multivariate normal distribution \(X_i \sim N(0, \Sigma)\). The treatment assignment \(D_i\) follows a Bernoulli distribution conditional on \(X_i\) with propensity score \(e(x) = (1 + \exp(-x'\alpha))^{-1}\).
We adopt a sparse setting with \(p = 50\) and uncorrelated covariates, \(\Sigma = I_p\). The propensity score depends on the first three covariates, with \(\alpha = (0.2, 0.2, 0.2, 0, \dots, 0)' \in \mathbb{R}^p\). The baseline outcome is a sparse linear function of the first five covariates, \[ \mu_0(x) = x_1 + 0.5 x_2 + 0.5 x_3 + 0.3 x_4 + 0.3 x_5. \]
We examine four specifications of the CATE function \(\tau(x) = \mu(1, x) - \mu(0, x)\):
Figure (ref) provides visualizations of the data generating processes through the scatter plots of \(Y_i\) against \(X_{1i}\), alongside the true conditional mean functions \(\mu(1, x_1)\) and \(\mu(0, x_1)\) evaluated at the mean of all other covariates. The models are designed to reflect qualitatively different patterns of treatment effect heterogeneity. The constant CATE model falls under the null of Equation (ref), while the other three models fall under the alternative. The linear and nonlinear models have a zero ATE by construction, so conventional ATE-targeted approaches such as OLS or IPW would fail to detect the existence of treatment effects.
We implement our proposed Algorithm (ref) using \(K=5\) folds. We estimate the nuisance parameters and the DR-learner for the CATE function using two machine learning algorithms: Lasso and XGBoost\footnote{We use the glmnet R package for Lasso and the xgboost package for XGBoost.}. To demonstrate the theoretical necessity of our IF-SS structure, we introduce a Naive DML benchmark. This benchmark utilizes standard DML cross-fitting but omits our internal sample-splitting step. Specifically, for each fold \(k\), it computes both the total variance \(\hat{V}_{\mathrm{tot},k}\) and the residual variance \(\hat{V}_{\mathrm{res},k}\) on the entire evaluation fold \(\mathcal{I}_k\) using the identically trained nuisance estimators \(\hat{\eta}_k\). The Naive DML benchmark uses XGBoost for nuisance estimation. While this naive approach preserves Neyman orthogonality, it fails to solve the null degeneracy problem. Under the null hypothesis, the influence function \(\phi_i = (\hat{\psi}_i - \hat{\tau}_{\mathrm{ATE}})^2 - (\hat{\psi}_i - \hat{\tau}(X_i))^2\) converges to zero for all \(i\), so that both the point estimate and the estimated standard error degenerate. Because the two variance components are evaluated on the same observations, the flexible CATE learner contributes a spurious dispersion that biases \(\hat{\theta}\) downward, and dividing this negative bias by a standard error of even smaller order drives the standardized statistic to \(-\infty\); the rejection probability of the one-sided test converges to zero (Proposition (ref)). The naive test is therefore severely undersized rather than unreliable in an unpredictable direction. We characterize this conservative degeneracy formally in Appendix (ref).
We also compare the performance with two existing tests in the literature, the nonparametric test proposed by crump2008nonparametric (hereinafter CHIM) and the debiased machine learning test proposed by semenova2021debiased (hereinafter SC).
We implement the sieve-based nonparametric test proposed by CHIM to evaluate the null hypothesis of a constant conditional average treatment effect. This method approaches the problem by comparing the shapes of the conditional outcome mean functions for the treated and control groups. We approximate these functions, \(\mu_1(x)\) and \(\mu_0(x)\), using a sieve basis expansion \(\mathbf{P}(x) = (1, p_1(x), \dots, p_K(x))'\), where the basis terms are constructed as a linear function of the covariates \(X_i\). This vector includes an intercept and \(K\) covariate-dependent basis terms. We estimate the coefficients by running two separate OLS regressions of the observed outcome \(Y_i\) on \(\mathbf{P}(X_i)\) for the treated and control subsamples, yielding the coefficient vectors \(\hat{\boldsymbol{\xi}}_1 = (\hat{\alpha}_1, \hat{\boldsymbol{\beta}}_1')'\) and \(\hat{\boldsymbol{\xi}}_0 = (\hat{\alpha}_0, \hat{\boldsymbol{\beta}}_0')'\). Under the null hypothesis, the treatment effect is constant, implying that the outcome functions are parallel and their slope coefficients are identical (\(\boldsymbol{\beta}_1 = \boldsymbol{\beta}_0\)). The test statistic evaluates the quadratic distance between these estimated slopes \[ T_{\mathrm{Crump}} = (\hat{\boldsymbol{\beta}}_1 - \hat{\boldsymbol{\beta}}_0)' \widehat{\mathbf{V}}_{\beta}^{-1} (\hat{\boldsymbol{\beta}}_1 - \hat{\boldsymbol{\beta}}_0), \] where \(\widehat{\mathbf{V}}_{\beta}\) is the robust covariance matrix for the difference in slope estimates.
As a benchmark for high-dimensional settings, we implement the Best Linear Predictor (BLP) test proposed by SC, following Example 2.2 in their paper. This framework approximates the CATE by projecting it onto a linear dictionary of covariates. The core of the method is the construction of a Neyman-orthogonal signal, which is the pseudo-outcome \(\psi(O_i)\) in Equation (ref), which serves as an unbiased proxy for the latent individual treatment effect. We employ the cross-fitting procedure proposed in their Definition 2.1. The sample is split into \(K\) folds, and for each observation \(i\) in fold \(k\), the signal \(\psi(O_i)\) is constructed using nuisance parameters estimated on the complementary folds. In the second stage, we project this cross-fitted signal onto a vector of covariates \(Z_i\) constructed as a second-order polynomial expansion of \(X_i\) (including interaction terms) to estimate the BLP coefficients. We solve the Lasso optimization problem \[ (\hat{\beta}_{0}, \hat{\boldsymbol{\beta}}_{\mathrm{Lasso}}) = \arg\min_{\beta_0, \boldsymbol{\beta}} \frac{1}{n} \sum_{i=1}^n (\hat{\psi}(O_i) - \beta_0 - Z_i'\boldsymbol{\beta})^2 + \lambda \|\boldsymbol{\beta}\|_1. \] The null hypothesis of a constant treatment effect implies that the best linear predictor is constant, or equivalently, that the slope coefficients are zero (\(\boldsymbol{\beta} = \mathbf{0}\)). We test this hypothesis using the debiased Lasso estimator to account for regularization bias, constructing a Wald statistic for the joint significance of the slope coefficients.
The empirical rejection proportions over \(1,000\) Monte Carlo replications at the nominal \(\alpha = 0.05\) level are presented in Table (ref). Under the constant CATE model, the results demonstrate the impasse detailed in Section (ref). The Naive DML estimator is severely undersized under the null: reusing the same evaluation fold biases its point estimate downward while its standard error degenerates even faster, so the standardized statistic drifts to \(-\infty\) and the one-sided test almost never rejects. The drift is slow for regularized learners, which is why the rejection rates remain small but non-zero and essentially flat across the sample sizes considered (see Appendix (ref)). Furthermore, CHIM and SC fail severely, with rejection rates far exceeding the nominal level even at large sample sizes. In contrast, our proposed IF-SS-CVT maintains excellent size control across all sample sizes.
Under the alternative hypotheses, all tests show consistent high power when \(n\) is large, while our IF-SS-CVT has lower power than the other tests when \(n\) is small. This is expected, as by randomly bisecting each evaluation fold to decouple empirical processes, the IF-SS-CVT operates on an effective sample size of \(n/2\). Despite this inherent finite-sample penalty, the IF-SS-CVT remains remarkably powerful when the sample size is large.
In this section, we demonstrate the application of the proposed test to the NSW job training program data. In this program, participants were randomly assigned to either a job training program or a control group, and the treatment effect on future earnings can be estimated by directly comparing outcomes of the treated and control groups. In order to evaluate the validity of econometric estimators of treatment effects, lalonde1986evaluating compared the treated individuals from the experiment to control groups drawn from two survey datasets: the Panel Study of Income Dynamics (PSID) and the Current Population Survey (CPS). The resulting datasets have been extensively analyzed in the influential works by dehejia1999causal, smith2005does, angrist2009mostly, sloczynski2022interpreting, among others. In the context of CATE hypothesis testing, the dataset was analyzed by hsu2017consistent and dai2023nonparametric, who focused specifically on heterogeneity with respect to age. Using the proposed test, we examine heterogeneity with respect to all available covariates.
The dataset we use is NSW-CPS, which contains 185 treated units from the experiment and 15992 control units from the CPS. The outcome \(Y_i\) is the earnings in 1978, and the treatment \(D_i\) is a binary indicator of whether the individual received the job training. We consider the same set of covariates as those in column 4 of Table 3.3.3 in angrist2009mostly, which includes age, age squared, education, dummy variables for black and Hispanic, marital status, a dummy indicator for high school degree, and pre-treatment earnings in 1974 and 1975. For this set of covariates \(X_i\), we test \(H_0: \tau(x) = c\) for some constant \(c\) and all covariate values \(x\). For nuisance parameter estimation in the IF-SS-CVT, we employ XGBoost. We compare the results with the CHIM and SC tests introduced in the simulation section, maintaining the same specifications for the basis functions (linear basis for CHIM and second-order polynomials with interactions for SC).
The test results are presented in Table (ref). The CHIM test fails to reject the null hypothesis of constant treatment effects at the 5% significance level (\(p = 0.40\)). This lack of rejection might be attributed to the test's lower power in finite samples with moderate-dimensional covariates, as observed in our simulations. In contrast, both the SC test and our proposed IF-SS-CVT with Lasso or XGBoost strongly reject the null hypothesis (\(p < 0.01\)), providing robust evidence for the presence of heterogeneous treatment effects. The rejection by the SC test suggests that some of the heterogeneity is linear in the covariates, while the consistent rejection by both Lasso- and XGBoost-based IF-SS-CVT confirms that this finding is not an artifact of a specific machine learning method. These findings complement the conventional ATE-focused analyses by highlighting that the treatment effect of job training likely varies across individuals with different characteristics.
This paper develops a hypothesis test for the presence of heterogeneous treatment effects by targeting a single omnibus parameter, the variance of the CATE, \(\theta_0 = \operatorname{Var}(\tau_0(X))\). In developing this test, we identify a fundamental theoretical impasse in semiparametric inference at the boundary of the parameter space. On one hand, evaluating variance components on the identical sample leads to null degeneracy, where the asymptotic variance collapses to zero and invalidates standard Gaussian approximations. On the other hand, decoupling the empirical processes via standard sample-splitting destroys the Neyman orthogonality of the squared pseudo-outcomes.
To resolve this impasse, we develop a novel Intra-Fold Sample-Split algorithm. By randomly bisecting each evaluation fold and computing the total and residual variance components on mutually disjoint halves, our procedure guarantees a positive asymptotic variance under the null. By strictly coupling both evaluation halves to identically trained nuisance estimators, the non-orthogonal squared biases cancel out. We formally prove that this algorithm restores Neyman orthogonality, yields asymptotic normality, and guarantees valid Type I error control under the null hypothesis.
Monte Carlo simulations and an empirical application to the NSW job training program confirm the robust finite-sample performance of the proposed test. Our simulations provide empirical proof that, on the boundary, standard cross-fitted DML statistics degenerate and become severely conservative while projection-based HTE tests severely over-reject, whereas our IF-SS-CVT attains the nominal size. Beyond testing for treatment effect moderation, our algorithm provides a general framework for conducting robust hypothesis testing on nonlinear transformations of doubly robust scores.