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.
80,901 characters · 17 sections · 45 citation commands
Estimation and Inference for Causal Functions with Multiway Clustered Data
Multiway cluster dependence is ubiquitous in empirical data. Units in such data exhibit strong dependence within geographical locations, industrial sectors, and other clusters. For instance, data for airline and automobile industries used for demand estimation berry1992estimation,berry1995automobile are indexed by markets and products, which represent two cluster sources of strong dependence through demand and supply shocks, respectively. Datasets used to analyze the causal effects of trade and technology shocks on employment autor2013growth,autor2015untangling,acemoglu2016import exhibit two-way cluster dependence by commuting zones and job occupations. Financial data used to analyze causal effects of technological innovation hsu2014financial exhibit two-way cluster dependence by countries and industries.
Accounting for this type of cross-sectional dependence is essential for conducting accurate statistical inference as such dependence may invalidate the conventional asymptotic theory developed for iid sampling and can produce substantial size distortions for causal inference cameron2011robust. The challenges that arise from multiway clustering necessitate special care to address complicated cross-sectional dependence structures.
The recent literature develops methods of inference about various estimands under multiway clustering, but they do not directly apply to inference for causal functions. This paper addresses this gap. Under multiway clustering, we propose novel methods for estimating and conducting uniform inference for a general class of non-parametric causal functions $\tau_0(\cdot)$, which can be indentified as a conditional expectation of the form
where the signal $\psi(\eta_0)$ may depend on a high-dimensional nuisance parameter $\eta_0$. Examples of such functions $\tau_0(\cdot)$ include the conditional average treatment effect (CATE) and continuous treatment effects (CTE), among others, as discussed in more detail later. Typical examples of nuisance parameters $\eta_0$ include the propensity score, conditional density, etc.
There are extensive discussions in the literature concerning causal functions with high-dimensional nuisance parameters. However, nearly all the results are provided under the iid sampling scheme. For instance, semenova2021debiased and fan2022estimation consider estimation and inference for the CATE function under iid sampling while kennedy2017non consider estimation and inference for the CTE function under the framework of iid sampling. To our knowledge, this literature is silent about sampling schemes with strong cross-sectional dependence, such as multiway clustering, which is relevant to many empirical data used in economics. The current paper aims to address this gap in the literature.
We propose a two-step procedure to estimate the causal function $\tau_0(\cdot)$. The first step involves estimating the high-dimensional nuisance parameters $\eta_0$ by a machine learner (ML; e.g., lasso, neural network, random forest). The second step involves a sieve estimation chen2007large of the causal function. We consider two approaches to implement this two-step procedure, namely, the full-sample and multiway cross-fitting approaches. In the full-sample approach, the nuisance parameters $\eta_0$ and causal function $\tau_0(\cdot)$ are estimated based on all the observations. On the other hand, the multiway cross-fitting approach splits the data into multiway folds according to the clustering scheme, so the estimation of causal functions relies on one multiway fold independent of the multiway folds in which the ML estimates nuisance parameters. Our full-sample approach is similar to belloni2017program under the iid case. On the other hand, our cross-fitting estimator differs from the iid counterpart chernozhukov2018double, and we apply a variant of the multiway cross-fitting method chiang2022multiway.
Since our sieve estimator entails a non-Donsker issue andrews1994empirical, a functional central limit theorem (CLT) is not applicable. Instead, we apply the high-dimensional CLT for the separately exchangeable arrays chiang2021inference and approximate the standardized process by an intermediate Gaussian process of an increasing dimension. As the limiting null distribution of the sup-test statistic is non-Gaussian, we also develop a resampling method that consistently approximates its critical values. Related bootstrap methods that work under multiway clustering include the pigeonhole bootstrap mccullagh2000resampling, polyadic bootstrap davezies2021empirical, and wild bootstrap mackinnon2021wild,menzel2021bootstrap, but none of these methods applies to a high-dimensional score vector. Therefore, we extend the sieve score bootstrap for iid data chen2018optimal to our novel setting of multiway clustering, introducing the multiway cluster-robust sieve score bootstrap. Finally, we prove the uniform probability coverage of the uniform confidence bands (UCBs) achieved by this novel bootstrap method.
In addition to the theoretical developments, we conduct simulations to examine the finite-sample performance of our proposed methods. The results reveal that both our full-sample and multiway cross-fitting UCBs perform well with reasonable size controls. Overall, the multiway cross-fitting UCBs achieve more accurate probability coverage in finite samples than the full-sample UCBs. Therefore, we suggest using the multiway cross-fitting UCBs over the full-sample ones in practice. Furthermore, with multiway clustered data, our proposed sieve score bootstrap method demonstrates superior probability coverage compared to the conventional method developed for the iid case, further corroborating the necessity of incorporating cross-sectional dependence into inference. Finally, we illustrate an empirical application to the analysis of the causal relationship between the mistrust levels in Africa and the history of slave trade. Our findings reveal heterogeneous treatment effects, with significant impacts at higher trade volumes. These findings significantly enhance our understanding of this critical empirical question.
{\bf Structure:} The rest of this paper proceeds as follows. Section (ref) introduces the setup. Section (ref) explores the estimation procedures. Section (ref) outlines the assumptions and provides the limit theory. Section (ref) details the uniform inference method based on our multiway cluster-robust sieve score bootstrap. Section (ref) evaluates the finite-sample performance of the proposed method. Section (ref) demonstrates an empirical application. Finally, Section (ref) provides the conclusion. Technical and additional details are relegated to the appendix and the online supplementary appendix.
{\bf Notations:} The quantities, $N$ and $M$, denote the cluster sizes. We write $\left[N\right] =\{1,2,...,N\}$ and $\left[M\right] = \{1,2,...,M\}$. The two-way sample sizes $\left(N,M\right)$ will be index by $n \in \mathbb{N}$ as $(N,M)=(N(n),M(n))$, where $N(n)$ and $M(n)$ are non-decreasing in $n$ and $N(n)M(n)$ is increasing in $n$; for instance $n=\min\{M,N\}$. For simplicity, we write the cluster sizes as $(N,M)$ omitting its dependence on $n$. Let $\{\mathcal{P}_n\}$ be a sequence of spaces of empirical probability laws of $\{Z_{ij}\}$ where we allow for increasing dimensions of $Z_{ij}$ in the sample size $n$. Let $\mathbb{P}=\mathbb{P}_n\in \mathcal{P}_n$ denote the law under the sample size $(N,M)$. The $\ell^2$-norm of a vector $v$ is denoted by $\left\Vert v\right\Vert$, and the operator norm of a matrix $Q$ is denoted by $\left\Vert Q\right\Vert$. For a given matrix-form quantity $Q(x)$ with $x\in\mathcal{X}$ {and $\mathcal{X}$ being a compact support}, $\left\Vert Q\right\Vert_{\mathbb{P},q}$ denotes $(\int_{x\in\mathcal{X}}\left\Vert Q(x)\right\Vert^qd\mathbb{P}\left(x\right))^{1/q}$. The notation $a \lesssim b$ indicates $a\leq cb$ for some constant $c$ that does not depend on the sample size $\left(N,M\right)$. The notation $a\lesssim_{p} b$ indicates $a=O_{p}\left(b\right)$. The notation $a\asymp b$ indicates (i) $a\lesssim b$ and $b\lesssim a$ or (ii) $a\lesssim_p b$ and $b\lesssim_p a$. The notations $\rightarrow_{p}$ and $\Rightarrow$ denote convergence in probability and weak convergence, respectively, in Euclidean or functional space. The vector $\ell_{1\times d}$ refers to a $d$-dimensional row vector of ones and the notation $\mathbf{0}_{d\times 1}$ denotes a $d$-dimensional (column) vector of zeros. We also use standard notations from the empirical process literature:
for any subsets $I_k\subseteq \{1,...,N\}$ and $J_\ell\subseteq \{1,...,M\}$.
We consider a class of causal functions $\tau_0(\cdot)$ identified in the form of (ref) with the Neyman-orthogonal signals, $\psi(\cdot)$, that satisfy the Neyman-orthogonality condition
for all $x$ and $\eta$. If the signal $\psi(\cdot)$ satisfies the Neyman-orthogonality (ref), then $\psi(\eta)$ is immune to the first-order bias of $\eta$, helping deliver a more efficient estimation of the causal function $\tau_0(\cdot)$. For the sake of generality, $\psi(\cdot)$ and $\tau_0(\cdot)$ are left unspecified for the moment, but we provide two concrete examples shortly in Sections (ref)--(ref).
Suppose that an econometrician observes copies of $Z$ including $X$ as a subvector, where the copies may exhibit cross-sectional dependence in the form of multiway clustering. Without loss of generality and for brevity, we consider the two-way clustering -- see Appendix (ref) in the supplementary material for a general multiway cluster sampling framework. That is, $\{Z_{ij}\}_{i\in[N], j\in[M]}$ is a double-indexed sample of sizes $N$ and $M$ satisfying the following condition.
Assumption (ref)(ref) imposes the condition of separate exchangeability davezies2021empirical, which is comparable with the identical distribution assumption of the iid case, that is, perturbing labels $i$ and $j$ separately does not affect the joint distribution of $\{Z_{ij}\}$. Assumption (ref)(ref) requires that $\{Z_{ij}\}$ is independent if they do not share a common index. However, the random vectors $\{Z_{ij}\}$ are allowed to be dependent in any form if they share the same index, $i$ or $j$.
Throughout this paper, we consider scenarios in which researchers use cross-sectionally dependent data $\{Z_{ij}\}_{i\in[N], j\in[M]}$ satisfying Assumption (ref) to estimate $\eta_0$ satisfying (ref), and thence $\tau_0$ via (ref). We propose two estimation strategies -- one based on the full sample and the other based on the multiway cross fitting -- and develop their asymptotic theory. We further propose a novel method of bootstrap resampling to construct uniform confidence bands of $\tau_0(\cdot)$.
The next two subsections overview two leading examples of our setup: one is the conditional average treatment effect (CATE); and the other is the continuous treatment effect (CTE). Each subsection introduces a causal function $\tau_0$ and an associated Neyman-orthogonal score $\psi$ satisfying (ref), that in turn identifies $\tau_0$ through (ref).
For each $i\in[N]$ and $j\in[M]$, let $D_{ij}$ indicate whether the individual is treated, $Y_{ij}(1)$ and $Y_{ij}(0)$ denote the potential outcomes under treatment and no treatment, respectively, and $W_{ij}$ denote a vector of covariates. The observed outcome is constructed by $Y_{ij}=D_{ij}Y_{ij}(1)+(1-D_{ij})Y_{ij}(0)$. The observed data $Z_{ij}$ consist of $Z_{ij} = (D_{ij},W_{ij}^{\prime}, Y_{ij})^{\prime}$ for $i\in[N]$, $j\in[M]$. The quantity of our interest is the conditional average treatment effect (CATE) function $\tau_0(\cdot)$ defined by
where $X_{ij}$ is a $d_X$-dimensional subvector of $W_{ij}$ of a fixed dimension $d_X$. We allow for high-dimensional data, where the dimension of the covariates $W_{ij}$ grows with the sample size.
Assume the observational unconfoundedness rosenbaum1983central. That is, the treatment status $D_{ij}$ is independent of the potential outcomes $(Y_{ij}(1)$, $Y_{ij}(0))$ conditional on observed controls $W_{ij}$:
Define $\mu_{0}(l,w)=\mathbb{E}[Y_{ij}|D_{ij}=l,W_{ij}=w]$ for $l=0,1$. The law of iterated expectations and unconfoundedness yield $\mathbb{E}[Y_{ij}(l)|X_{ij}=x]=\mathbb{E}[\mu_0(l,w)|X_{ij}=x]$. Thus, the CATE function $\tau_0(\cdot)$ can be identified by $\tau_0(x)=\mathbb{E}[\mu_0(1,w)-\mu_0(0,w)|X_{ij}=x]$. Moreover, we consider a more robust result based on the Neyman-orthogonal signal fan2022estimation adapted to our multiway clustered setup:
where $\pi_{0}(w)=\mathbb{P}(D_{ij}=1|W_{ij}=w)$ denotes the propensity score. Lemma (ref) in the supplementary appendix shows that the signal (ref) identifies (ref) via (ref) and is both locally and doubly robust with respect to the nuisance parameter $\eta_0(w):=(\pi_0(w),\mu_0(1,w),\mu_0(0,w))$.
Whenever the CATE is discussed throughout the rest of this paper, we shall assume its identification condition, formally stated in the supplementary appendix as Assumption (ref).
For each $i\in[N]$ and $j\in[M]$, let $X_{ij}$ be a continuous treatment intensity that takes values in $\mathcal{X} \subset \mathbb{R}$, $\{Y_{ij}(x)\}_{x \in \mathcal{X}}$ be the potential outcomes indexed by $x \in \mathcal{X}$, and $W_{ij}$ be a covariate vector. An econometrician observes the actual outcome $Y_{ij}=Y_{ij}(X_{ij})$ instead of the potential outcome $Y_{ij}(x)$ per se. Let the observed data $Z_{ij}=(X_{ij}, W_{ij}^{\prime}, Y_{ij})^{\prime}$ be a multiway clustered sample supported on $\mathcal{Z}=\mathcal{X} \times \mathcal{W}\times\mathcal{Y} \subseteq \mathcal{R}\times\mathcal{R}^{d_W}\times \mathcal{R}$. We consider the setup in which the dimension $d_W$ of the covariates, $W_{ij}$, grows with sample size, while the continuous treatment $X_{ij}$ of interest is a scalar random variable. The main target function is the continuous treatment response, $\tau_0(\cdot)$, defined by
In causal inference, we are often interested in its linear transformations that measure continuous treatment effects (CTE), such as the derivative function $d\tau_0(\cdot)/dx$ or a difference function $\tau_0(\cdot) - \tau_0(\bar x)$ from a fixed benchmark value $\bar x$ of $x$.
Assume the observational unconfoundedness. That is, the potential outcomes $\{Y_{ij}(x), x\in \mathbb{R}\}$ are independent of the continuous treatment $X_{ij}$ conditional on the observed controls $W_{ij}$:
Similarly to the notations in Section (ref), we define $\mu_{0}(x,w)=\mathbb{E}[Y_{ij}|X_{ij}=x,W_{ij}=w]$. The law of iterated expectations and the unconfoundedness yield $\tau_0(x):=\mathbb{E}[Y_{ij}(x)]=\mathbb{E}[\mu_0(x,W_{ij})]$ for any given $x$. We consider to extend the doubly robust signal from kennedy2017non to the multiway clustering setting:
that relies on the marginal distribution ${P}_{W}(\cdot)$ of $W_{ij}$, the generalized propensity score
and the marginal treatment density
In particular, Lemma (ref) in the supplementary appendix shows that the signal (ref) can identify the causal function via (ref) in both locally robust and doubly robust manners with the nuisance parameter
The corresponding estimation and inference procedures for the CTE can be further established based on the doubly robust signal ${\psi}(Z_{ij};\eta)$ defined in (ref).
Whenever the CTE is discussed throughout the rest of this paper, we shall assume its identification conditions, formally stated in the supplementary appendix as Assumptions (ref)--(ref).
This section presents estimation approaches for the general causal functions $\tau_0(\cdot)$ (encompassing the CATE and the CTE in particular) identified via (ref) based on the Neyman-orthogonal signal $\psi(\cdot,\cdot)$ satisfying (ref). We propose two versions of non-parametric estimators, depending on whether the first-step ML estimation and the second-step sieve regression use the same set of empirical data. These two approaches to estimation, referred to as the full-sample and multiway cross-fitting approaches, are separately introduced below.
This section provides a uniform limit theory for the sieve estimator of causal functions. To facilitate our discussions, we first introduce some additional notations.
The sup-norm of the sieve basis functions will be denoted by $\xi_{p}:=\sup_{x\in\mathcal{X}}\Vert p(x)\Vert$. Let $\alpha(x):=p(x)/\Vert p(x)\Vert$ denote the normalized basis $p(x)$, and let
denote its Lipschitz constant. With these notations, we now introduce some conditions for the second-step non-parametric estimation.
Assumption (ref) collects several standard regularity conditions used in the non-parametric estimation literature. Assumption (ref)(ref), as an identification condition that eliminates collinearity in the population signal matrix, is standard in the literature of non-parametric estimation andrews1991heteroskedasticity,newey1997convergence. Assumptions (ref)(ref)--(ref) bound the non-parametric approximation error. Specifically, the sequence of constants $r_{p}$ bound the $\ell^2$-norm while $\ell_{p}r_p$ bounds the $\ell^{\infty}$-norm. Assumption (ref)(ref) bounds the conditional variance of $\left\{u_{ij}\right\}$ uniformly, facilitating the uniform limit theory. In particular, Assumption (ref)(ref) implies that $2<m\leq q$.
Define the lower and upper bounds of the conditional second moment of $u_{ij}$ given $X_{ij}$ by
We impose the following condition on the upper bound.
In what follows, we state two conditions that bound the estimation errors of the Neyman-orthogonal moments for each of the two cases of full-sample and multiway cross-fitting estimators. First, we give high-level conditions for the full-sample estimator, lower-level sufficient conditions for which will be presented later in the contexts of the CATE and the CTE.
Assumption (ref) assumes the robustness condition for the unobserved orthogonal signal $\psi(Z_{ij};\eta_0)$ regarding the biased estimation of the nuisance parameter. The above condition is comparable to the upper bound for estimation errors in fan2022estimation and semenova2021debiased.
Next, we also bound the estimation errors of the Neyman-orthogonal moments in the cross-fitting case. For simplicity, we assume that all the observations are partitioned into blocks of equal size; that is, $\vert I_k\vert=\vert I\vert\asymp N$ and $\vert J_{\ell}\vert=\vert J\vert\asymp M$ for all $k$, $\ell\in[K]$.
We will show the plausibility of Assumptions (ref) and (ref) in the context of CATE and CTE with lower-level sufficient conditions in Sections (ref)--(ref). Given these conditions, Theorem (ref) below establishes the uniform linear representation and uniform convergence rate for the estimator of general causal functions.
See Appendix (ref) for a proof.
Theorem (ref) provides the linear representations of the full-sample and cross-fitting estimators with uniform upper bounds of their remainder terms. These linear representation forms show the leading terms of the empirical process of interest and serve as foundations for strong Gaussian approximation and uniform inference procedures for the causal functions.
Before introducing the high-dimensional central limit theorem (CLT), we remark that two-way clustered data satisfying Assumption (ref) are known to accommodate the so-called Aldous-Hoover-Kallenberg representation kallenberg1989representation:
where $U_{(i,0)}$, $U_{(0,j)}$ and $U_{(i,j)}$ follow the iid $U[0,1]$ distribution and are mutually independent. The above factor representation implies that we may focus on the three iid leading components in orthogonal directions, ($U_{(i,0)},U_{(0,j)},U_{(i,j)}$).
Now, define $g_{i0}(U_{(i,0)})=\mathbb{E}[p_{ij}u_{ij}\vert U_{(i,0)}]$ and $g_{0j}(U_{(0j)})=\mathbb{E}[p_{ij}u_{ij}\vert U_{(0,j)}]$, where $U_{(i,0)}$ and $U_{(0,j)}$ follow $\text{iid }U(0,1)$. For simplicity, we write $g_{i0}(U_{(i,0)})$ and $g_{i0}(U_{(i,0)})$ as $g_{i0}$ and $g_{0j}$ when there is little risk of ambiguity. With these definitions and notations, we now introduce the following necessary conditions for the high-dimensional CLT.
Assumption (ref) imposes conditions that ensure the applicability of the high-dimensional CLT for the separately exchangeable array chiang2021inference. Assumption (ref)(ref) requires that each coordinate of $p_{ij}u_{ij}$ is bounded. Assumption (ref)(ref) requires the maximum of the third moment across coordinates to be increasing at a speed no faster than the polynomials of $D_{n}$. Assumption (ref)(ref) further ensures that the H\'{a}jek projection of the multiway clustering empirical process of interest is nondegenerate.
Following the H\'{a}jek projection given in (ref)--(ref) in the supplementary appendix, we define the asymptotic variance:
with $\overline{\mu}_{N}=\lim(N\wedge M)/N$ and $\overline{\mu}_{M}=\lim(N\wedge M)/M$.
Define the approximating Gaussian variate, $\gamma_{\Sigma}\overset{d}{=} \mathcal{N}(\mathbf{0}_{p\times1},\Sigma)$. In what follows, Theorem (ref) establishes a strong Gaussian approximation to the non-parametric sieve estimator by a sequence of zero-mean Gaussian processes, $\gamma_{\Sigma}$.
In the proof of this theorem, provided in Appendix (ref), we apply the high-dimensional CLT for separately exchangeable arrays chiang2021inference and check all its prerequisite conditions.
To further facilitate the uniform inference procedure for general causal functions, we need to estimate the covariance matrix $\gamma_{\Sigma}$ by considering the iid components of the H\'{a}jek projection in Proposition (ref) in the supplementary appendix. Specifically, when we have the full sample, the covariance matrix of the Gaussian approximating term, $\Sigma$, can be estimated by
When we consider the multiway cross-fitting method, the covariance matrix of the Gaussian approximating term of the $(k,\ell)$-block can be estimated by
Averaging across partitions, we have the covariance estimator given by
We also define the population and sample signal matrices by
with $\widehat{Q}_{k\ell}:=\mathbb{E}_{n,k\ell}[p_{ij}p^{\prime}_{ij}]$. We then define the population variance $\sigma^{2}_{\tau}(x)=p(x)^{\prime}Q^{-1}\Sigma Q^{-1}p(x)$ with the full-sample and multiway cross-fitting estimators:
See Appendix (ref) for a proof.
With consistent variance estimates, we can stabilize the test statistic for various causal functions and further asymptotically approximate this test statistic by a Gaussian process whose distribution depends on data only by the variance function $\sigma_{\tau}(\cdot)$. Enhanced by this theory of the strong Gaussian approximation, we can conduct uniform inference accompanied by the bootstrapped critical values, which will be provided in Section (ref).
Before proceeding with the inference, we discuss lower-level sufficient conditions for the high-level statements of Assumptions (ref)--(ref) in the context of the CATE and the CTE in the following two subsections.
This subsection provides low-level sufficient conditions on the regression functions $\mu_0(1, \cdot)$, $\mu_0(0,\cdot)$ and propensity score $\pi_0(\cdot)$ such that the uniform convergence rate of Theorem (ref) and uniform Gaussian approximation of Theorem (ref) hold for the CATE.
Throughout, we assume the identifying condition for the CATE, stated as Assumption (ref) in the supplementary appendix. See Lemma (ref) in the supplementary appendix for the formal identification result.
We are going to restrict the estimation errors of the first-step nuisance parameter estimations. The upper bounds of estimation errors are shown separately for the cases of full-sample and multiway cross-fitting estimators. We first impose the following conditions for the full-sample estimator $\widehat{\eta}$.
Assumptions (ref)(ref) and (ref)(ref) control the errors of first-step nuisance parameter estimation $\widehat{\eta}$ in various norms. Assumption (ref)(ref) imposes conditions on the complexity of the functional space for nuisance parameters by entropy conditions, which can be satisfied by random forest wager2018estimation, deep neural network farrell2021deep and other machine learners.
Next, we provide similar conditions for the multiway cross-fitting estimator.
We impose a finite value $K$ and slightly abuse the notations for the multiway cross-fitting case and use $\delta_{1n,k\ell}$, $\delta_{2n,k\ell}$ and $\delta_{3n,k\ell}$ to bound the first-step estimation errors. It is evident that Assumption (ref) imposes no conditions on the entropy conditions of the space of nuisance functions since the first-step and second-step estimations rely on independent subsamples. To be specific, in the theoretical proofs of cross-fitting, we can treat the estimators of the nuisance parameters as fixed by conditioning on the subsample that estimates causal functions. This idea helps simplify our theoretical discussions substantially. Since we assume that all the observations are split into blocks of equal sizes, then $\vert I_k\vert=\vert I\vert\asymp N$, $\vert J_{\ell}\vert=\vert J\vert\asymp M$ for each $k$, $\ell\in[K]$.
Under the above conditions, the following theorem establishes the consistency of estimating the Neyman-orthogonal signal, $\psi(\cdot,\cdot)$.
As a result, the statements of Theorems (ref)--(ref) hold for the conditional average treatment effects (CATE) estimation and inference, further enabling the uniform convergence rate in estimating the non-parametric causal function and the associated coupling principle.
This subsection provides low-level sufficient conditions on the nuisance parameters
for the CTE to ensure the high-level condition of Assumptions (ref) and (ref), so that the uniform rate of Theorem (ref) and Gaussian approximation of Theorem (ref) continue to hold for the case of the CTE estimation and inference.
Throughout, we assume the identifying condition for the CTE, stated as Assumptions (ref)--(ref) in the supplementary appendix. See Lemmas (ref)-(ref) in the supplementary appendix for the formal identification result.
Similarly to the discrete case of the CATE, we are going to restrict the errors of the first-step nuisance parameter estimations. The upper bounds of estimation errors will be shown for the full-sample and cross-fitting estimators separately. Then, we proceed to show that the Neyman-orthogonal signal (ref) in the supplementary appendix satisfies the error bound given in Assumptions (ref)--(ref) for both the full-sample and cross-fitting estimators.
In what follows, we first discuss the conditions for the full-sample estimator, $\widehat{\eta}$.
Similarly to the case of the CATE, Assumptions (ref)(ref)--(ref)(ref) bound the estimation errors for the nuisance parameters regarding the full-sample estimator. Assumption (ref)(ref) restricts the complexity of the space of nuisance parameters.
Next, we provide low-level conditions for the nuisance parameters in the case of the multiway cross-fitting estimator.
As the cross-fitting method estimates the causal functions by the subsample that stays independent from the first-step subsample, imposing conditions for the space of nuisance parameters is unnecessary, similar to the case of the CATE. To simplify our discussions, we also assume that all the observations used for the CTE estimation and inference are partitioned into blocks of equal sizes. That is, $\vert I_k\vert=\vert I\vert\asymp N$, $\vert J_{\ell}\vert=\vert J\vert\asymp M$ for all $k$, $\ell\in[K]$.
The following theorem shows that the Neyman-orthogonal signal, $\psi(\cdot,\cdot)$ can be estimated accurately for the CTE inference.
Once the conditions for Lemma (ref) hold, the conclusions of Theorems (ref)--(ref) hold for the CTE estimation, further enabling the uniform rate of convergence in estimating the CTE function and the associated strong Gaussian approximation. The proof of Lemma (ref) is provided in Appendix (ref) in the supplementary appendix.
This section provides a uniform inference method by extending the sieve score bootstrap approach chen2018optimal to the cases of multiway clustered data. We call this extended method the “multiway cluster-robust sieve score bootstrap”. Our testing approach facilitated by this bootstrapping method uses the estimated nuisance parameters of the first step and randomizes the empirical process of the sieve estimator in the second step, avoiding the heavy computation burdens for high-dimensional covariance matrices. Our multiway cluster-robust sieve score bootstrap method also differs from its iid version chen2018optimal in what types of empirical processes are involved in re-sampling. In the following discussion, we will detail the challenges the conventional sieve score bootstrap faces under multiway clustered data and introduce the uniform inference procedure that relies on our multiway cluster-robust sieve score bootstrap method.
Under random sampling, chen2018optimal propose to use the sieve score bootstrap method to calculate critical values $cv_{n}^{b}(1-\alpha)$ so that the confidence bands can attain the uniform probability coverage. Adapted to our two-index notations, they define $\{\omega_{ij}\}_{i\in[N],j\in[M]}$ as the iid standard Gaussian variates which are independent of data $\{Z_{ij}\}_{i\in[N],j\in[M]}$. Then, chen2018optimal define the sieve score bootstrapping $t$-statistic empirical process:
with the residuals $\widehat{u}_{ij}$ from Section (ref). The test statistic process can capture the distribution of $\sup_{x\in\mathcal{X}}t^{b}_{\tau}(x)$ and calculate its critical values.
The above sieve score bootstrapping method cannot be applied to data with strong cross-sectional dependence. The iid perturbation process $\{\omega_{ij}\}$ utilizes degrees of independence (or weak dependence) in $Z_{ij}$ across both $i\in[N]$ and $j\in[M]$, so the validity of resampling relies crucially on the condition of independence (or weak dependence). However, the presence of strong within-cluster dependence violates this prerequisite condition when we consider the multiway clustered data. Alternatively, we rely on the independence across clusters and further modify the bootstrapping method of chen2018optimal to ensure its applicability in the present context.
In the following subsection, we first revisit the key machinery, the H\'{a}jek projection, that facilitates the establishment of a novel sieve score bootstrap method that stays robust to multiway clustering. Further, we will apply this bootstrapping method and discuss how it helps construct the bootstrapping uniform confidence bands (UCBs) for causal functions.
Though two-way clustered data introduces strong cross-sectional dependence in two dimensions, Proposition (ref) in the supplementary appendix shows that this type of data can be projected onto two orthogonal directions expanded by iid random variables $U_{(i,0)}$ and $U_{(0,j)}$, which follow uniform distributions, as shown by (ref) in the supplementary appendix. (See the Aldous-Hoover-Kallenberg representation (ref) for $U_{(i,0)}$ and $U_{(0,j)}$.) This intuition motivates us to randomize the projected score vector rather than the original score as in (ref).
Even though Proposition (ref) in the supplementary appendix illustrates that the score vector of interest can be projected on two orthogonal spaces, it is infeasible to directly approximate the unknown projections $g_{i0}(U_{(i,0)})$ and $g_{0j}(U_{(0,j)})$ since $U_{(0,i)}$ and $U_{(0,j)}$ are latent. Using the idea of menzel2021bootstrap, we approximate $g_{i0}(U_{(i,0)})$ and $g_{0j}(U_{(0,j)})$ by the following averages
respectively, where $\widehat{u}_{ij}$ is given in Section (ref). For simplicity, we write $\widehat{g}_{i0}(U_{(i,0)})$ as $\widehat{g}_{i0}$ and $\widehat{g}_{0j}(U_{(0,j)})$ as $\widehat{g}_{0j}$.
The multiway cluster-robust sieve score bootstrap method can be summarized as follows. Let $\{\omega_{1,i}\}_{i\in[N]}$ and $\{\omega_{2,j}\}_{j\in[M]}$ be independent $\mathcal{N}(0,1)$ random variables independent of the data. Then, we obtain the multiway cluster-robust sieve score bootstrap empirical process:
To compute the $1-\alpha$ critical value, $cv_{n}^{b}(1-\alpha)$, one can calculates $\sup_{x\in\mathcal{X}}\vert t^{b}_{\tau}(x)\vert$ based on two sequences of independent draws of $\{\omega_{1,i}\}$ and $\{\omega_{2,j}\}$, which are independent between each other. Specifically, $cv_{n}^{b}(1-\alpha)$ is given by the $(1-\alpha)$-quantile of $\sup_{x\in\mathcal{X}}\vert t^{b}_{\tau}(x)\vert$ over the multiple random draws of $\{\omega_{1,i}\}$ and $\{\omega_{2,j}\}$:
Then, we compute the bootstrapping UCBs for the causal function $\tau_0(\cdot)$ by
where $cv^{b}_{n}(1-\alpha)$ ensures that $\tau(x)\in[\widehat{\tau}^{b}_{l}(x),\widehat{\tau}^{b}_{u}(x)]$ for all $x\in\mathcal{X}$ with the level of confidence $100(1-\alpha)\%$ asymptotically, as formally shown in the following theorem.
See Appendix (ref) for a proof. Theorem (ref) extends the conventional sieve score bootstrap method chen2018optimal that suits iid data to the case of multiway clustering. In such a case, though the cross-sectional correlation that can result in additional technical challenges exists, we can still apply the high-dimensional central limit theorem of chiang2021inference to recover the Gaussian approximation and prove the uniform probability coverage of the multiway cluster-robust sieve score bootstrapping UCBs. Also, the asymptotic validity of the UCBs holds regardless of whether the estimator relies on the full sample. Regardless of whether the full-sample or cross-fitting estimator is considered, the uniform size control brought by our bootstrapping method always holds.
In this section, we examine the finite-sample performance for our bootstrap uniform inference procedures on the CATE and the CTE through numerical simulations.
The current subsection focuses on the CATE function. The observed outcome is constructed by
for each $i\in[N]$ and $j\in[M]$, where the potential outcomes are generated by
With $x$ denoting the first coordinate of $w$, the structure $\mu_1(\cdot)$ can be specified via (i) the polynomial function $\mu_1(w)=x$; (ii) the exponential functions $\mu_1(w)=e^x/(1+e^x)$ and $e^{3x}/(1+e^{3x})$; and (iii) the trigonometric functions $\mu_1(w)=\cos(x)$, $\sin(x)$ and $\sin(x)+\cos(x)$. On the other hand, we set $\mu_0(w)=x$ throughout. The treatment allocation is determined by
with the logistic link function $\Lambda(\cdot)$ and the $d$-dimensional parameters $\zeta=\left(0.7, 0.7^2, \ldots, 0.7^d\right)^{\prime}$. Our target estimand is the CATE function $\tau_0(\cdot)$ defined by $\tau_0(x)=\mathbb{E}[Y(1)|X=x]-\mathbb{E}[Y(0)|X=x]$, or equivalently, $\tau_0(x)=\mu_1(x)-\mu_0(x)$.
The primitive two-way clustered random variables $W_{ij}$, $\varepsilon_{ij}$ and $v_{ij}$ are constructed by
where $( r_1^W, r_2^W)$, $( r_1^{\varepsilon}, r_2^{\varepsilon})$ and $( r_1^v, r_2^v)$ denote the two-way clustering weights; the $d$-dimensional latent variables $\alpha_{ij}^W$, $\alpha_{i}^W$, and $\alpha_{j}^W$ are generated jointly according to
and $\alpha_{ij}^{\varepsilon}$, $\alpha_{ij}^{\varepsilon}$, $\alpha_{ij}^{\varepsilon}$, $\alpha_{ij}^v$, $\alpha_{i}^v$ and $\alpha_{j}^v$ are independently generated according to
The three pairs of weights $\left( r_1^W, r_2^W\right)$, $\left( r_1^{\varepsilon}, r_2^{\varepsilon}\right)$, and $\left( r_1^v, r_2^v\right)$ dictate the strength of the cross-sectional dependence. The dependence parameter $\rho$ specifies the magnitude of collinearity within the covariates $W_{ij}$. The values of these parameters are set to $\left( r_1^W, r_2^W\right)=( r_1^{\varepsilon}, r_2^{\varepsilon})=( r_1^v, r_2^v)=(0.4,0.4)$ and $\rho=0.25$ throughout. We compute empirical probability coverage of our test statistics under the null hypothesis based on $1,000$ Monte Carlo iterations for each scenario, with the number of observations fixed at $N=M=25$. The uniform coverage of $\tau_0(\cdot)$ is tested across $100$ grid points over the 1%--99% quantiles of the sample distribution of $X_{ij}$.
Table (ref) reports the results when the dimension of the covariates is $d=4$. Several findings are given in order. First, compared to the pointwise and iid sieve score bootstrap approaches, the UCBs based on the multiway cluster-robust sieve score bootstrap method attain the nominal level more accurately for all the data-generating processes. The size distortion induced by pointwise critical values is widely shown in the literature, and we refer interested readers to belloni2015some for more details. Also, the iid bootstrap method fails to capture the cluster dependence and generates spurious statistical significance. Second, the multiway robust cross-fitting technique can help reduce over-rejection in the finite sample. Table (ref) shows that the full-sample uniform inference procedure still produces substantial size distortions when the function of interest is highly nonlinear (e.g., $\sin(x)+\cos(x)$). Comparatively, the multiway cross-fitting method can generate empirical probability coverage nearly identical to the nominal level. Such observations echo the motivation of the cross-fitting estimator to improve the finite sample performance in the iid setting chernozhukov2018double.
This subsection investigates the finite sample performance of our multiway cluster-robust sieve score bootstrapping inference for the CTE. The outcome is constructed by
for $i\in[N]$ and $j\in[M]$, where $X_{ij}$ that represents the one-dimensional continuous treatment. The function $g(x)$ can be specified as (i) the polynomial function $g(x)=x$; (ii) the exponential functions $g(x)=e^x/(1+e^x)$ and $e^{3x}/(1+e^{3x})$; and (iii) the trigonometric functions $g(x)=\cos(x)$, $\sin(x)$ and $\sin(x)+\cos(x)$. The continuous treatment $X_{ij}$ follows a beta distribution:
The $d$-dimensional covariates $W_{ij}$ and a scalar error $\varepsilon_{ij}$ are generated in exactly the same ways as in the CATE setup presented in Section (ref). The parameters are set to $\gamma=(0.5,0.5^2,...,0.5^d)^{\prime}$, $\zeta=(0.7,0.7^2,...,0.7^d)^{\prime}$, and we let $d=4$. Our targeted function is the continuous treatment response function $\tau_0(\cdot)$ given by $\tau_0(x):=E[Y(x)]=g(x)$.
We compute the empirical probability coverage of our bootstrapping test statistics under the null hypothesis with the nominal level of $95\%$ based on 1,000 Monte Carlo iterations for each scenario, with the number of observations fixed at $N=M=25$. The uniform coverage of $\tau_0(\cdot)$ is tested across $100$ grid points over the 20%--80% quantiles of the sample distribution of $X_{ij}$.
Table (ref) evidently shows that the conventional inference approach that relies on the standard Gaussian critical values and the iid sieve score bootstrapping method can lead to severe size distortions for both the full sample and cross-fitting cases. Our multiway cluster-robust sieve score bootstrapping method also incurs size distortions with the full sample case. The reason might lie in the complexity of the nuisance parameter associated with the marginal treatment density in the first step. On the other hand, the multiway cluster-robust sieve score bootstrapping method relying on the multiway cross-fitting procedure stays valid in the setting of multiway clustered data, though some extent of conservativeness can be observed for this novel bootstrapping method under the cross-fitting case. This evidence again corroborates the robust behavior of multiway cross-fitting methods in improving finite-sample performance in the setting of cross-sectionally correlated data.
In summary, we propose employing the multiway cluster-robust sieve score bootstrap method based on the multiway cluster-robust cross-fitting estimation for practical applications.
In this section, we present an empirical application of our multiway cluster-robust inference procedures. Revisiting the empirical debate on the causal relationship between the level of mistrust in Africa and the history of slave shipments, we provide non-parametric estimates of continuous treatment effects (CTEs), accompanied by both pointwise confidence intervals (CIs) and uniform confidence bands (UCBs).
We use the empirical data from nunn2011slave, defining Trust of Neighbors as the observed outcome and slave Exports as the treatment variable. Specifically, we examine the CTEs of slave Exports on Trust of Neighbors while controlling for a high-dimensional set of covariates. Our empirical analysis considers eleven covariates, including the respondent’s age, age squared, a gender indicator variable, and an indicator variable for whether the respondent resides in an urban area, among others. The dataset consists of 20,027 observations clustered along two dimensions: ethnicities ($i$) and regions ($j$). The number of clusters is $N = 185$ for the ethnicity index ($i$) and $M = 171$ for the region index ($j$). The effective number of clusters, given by $N \wedge M = 171$, is sufficiently large to justify the use of our asymptotic theory.
To investigate the causal relationship of interest, we estimate the CTEs, defined as the derivative function $d\tau_0(\cdot)/dx$ of the continuous treatment response function $\tau_0(\cdot)$. Our estimation procedure leverages the Neyman-orthogonal signal for $\tau_0(\cdot)$ from (ref). We construct pointwise CIs based on standard Gaussian critical values and multiway cluster-robust sieve score bootstrap UCBs using the method outlined in Section (ref). Our study may be considered to revisit the findings presented in Columns 1 and 2 of Table 2 in nunn2011slave by relaxing their parametric assumptions on continuous treatment effects. We conduct uniform inference across 100 grid points, covering the full range of the treatment variable from its minimum to maximum values. The resulting estimates, along with both pointwise CIs and UCBs, are displayed in Figures (ref) and (ref).
Figure (ref) illustrates the causal relationship between Trust of Neighbors and the nonlinear transformation of slave Exports, specifically $\ln(1 + Exports/Area)$. The vertical axis measures the CTE measured by $d\tau_0(\cdot)/dx$. The figure highlights heterogeneous treatment effects that would not be captured by linear parametric models. The treatment effects are not significant at lower levels of the slave export but become significantly negative at higher levels. This significance is indicated not only by the pointwise CIs but also by the UCBs. Notably, these results suggest that the null hypothesis, $H_0: d\tau_0(x)/dx = 0 \ \forall x$, which posits uniformly zero treatment effects, is rejected. Figure (ref), which illustrates the causal relationship between Trust of Relatives and the nonlinear transformation of slave Exports, yields qualitatively the same conclusion.
These new empirical findings, enabled by our proposed method of uniform inference, advance our understanding of the causal relationship between mistrust and the history of the slave trade in Africa. Specifically, we reject the joint hypothesis of uniformly zero treatment effects, with the heterogeneous effects being particularly strong at high levels of the slave trade.
This paper introduces novel methods for the estimation and uniform inference of causal functions under the conditions of multiway clustering. Our focus is on a general class of causal functions, such as the Conditional Average Treatment Effect (CATE) and Continuous Treatment Effect (CTE), which are characterized as conditional expectations of a Neyman-orthogonal signal dependent on high-dimensional nuisance parameters. We propose a two-step procedure that leverages machine learning techniques to estimate nuisance parameters and then applies a sieve estimation method for the parameter of interest.
Our primary methodological contribution is the development of the multiway cluster-robust sieve score bootstrap, an extension of the conventional sieve score bootstrap that is specifically tailored to account for multiway clustering in data. This resampling procedure is crucial for accurate inference in contexts where strong dependencies exist within clusters, which can otherwise result in significant size distortions if conventional asymptotic methods are applied. The proposed method achieves uniform confidence bands with desirable finite-sample properties, as evidenced by both our theoretical developments and extensive Monte Carlo simulations.
Overall, our findings contribute to the growing literature on causal inference under complex dependence structures. By developing tools that accommodate multiway clustering, we provide researchers with more reliable methods for drawing inferences from data characterized by such dependencies. Future research could extend these methods to other types of clustering structures or explore their applicability in broader empirical settings. The methodological innovations presented here open new avenues for robust causal analysis, particularly in fields where data are inherently clustered across multiple dimensions.
Our work underscores the importance of accounting for cross-sectional dependencies in empirical research. The multiway cluster-robust sieve score bootstrap not only improves the accuracy of statistical inference but also prevents misleading conclusions that can arise from ignoring the complex dependence structures often present in real-world data. We encourage further exploration and application of these methods in diverse empirical contexts, aiming to enhance the robustness and reliability of causal inferences in social science research.