Multiway Cluster Robust Double/Debiased Machine Learning

Multiway Cluster Robust Double/Debiased Machine Learning

Introduction

We propose a novel multiway cross fitting algorithm and a double/debiased machine learning (DML) estimator based on the proposed algorithm. This objective is motivated by recently growing interest in use of dependent cross sectional data and recently increasing demand for DML methods in empirical research. On one hand, researchers frequently use multiway cluster sampled data in empirical studies, such as network data, matched employer-employee data, matched student-teacher data, scanner data where observations are double-indexed by stores and products, and market share data where observations are double-indexed by market and products. On the other hand, we have witnessed rapidly increasing popularity of machine learning methods in empirical studies, such as random forests, lasso, post-lasso, elastic nets, ridge, deep neural networks, and boosted trees among others. To date, available DML methods focus on i.i.d. sampled data. In light of the aforementioned research environments today, a new method of DML that is applicable to multiway cluster sampled data may well be of interest by empirical researchers.

The DML was proposed by the recent influential paper by CCDDHNR18. They provide a general DML toolbox for estimation and inference for structural parameters with high-dimensional and/or infinite-dimensional nuisance parameters. In that paper, the estimation method and properties of the estimator are presented under the typical microeconometric assumption of i.i.d. sampling. We advance this frontier literature of DML by proposing a modified DML estimation procedure with multiway cross fitting, which accommodates multiway cluster sampled data. Even for multiway cluster sampled data, we show that the proposed DML procedure works under nearly identical set of assumptions to that of CCDDHNR (CCDDHNR18). To our best knowledge, the present paper is the first to consider generic DML methods under multiway cluster sampling.

Another branch of the literature following the seminal work by CGM11 proposes multiway cluster robust inference methods. Menzel17 conducts formal analyses of bootstrap validity under multiway cluster sampling robustly accounting for non-degenerate and degenerate cases. DDG18 develop empirical process theory under multiway cluster sampling which applies to a large class of models. We advance this practically important literature by developing a multiway cluster robust inference method based on DML. In deriving theoretical properties of the proposed estimator, we take advantage of the Aldous-Hoover representation employed by the preceding papers. To our knowledge, the present paper is the first in this literature on multiway clustering to develop generic DML methods.

Relations to the Literature

The past few years have seen a fast growing literature in machine learning based econometric methods. For general overviews of the field, see, e.g., AtheyImbens19 or MullainathanSpiess17. For a review of estimation and inference methods for high-dimensional data, see BCH14review. For an overview of data sketching methods tackling computationally impractically large number of observations, see LeeNg19. The DML of CCDDHNR (CCDDHNR18) is built upon BCK15, which proposes to use Neyman orthogonal moments for a general class of Z-estimation statistical problems in the presence of high-dimensional nuisance parameters. This framework is further generalized in different directions by BCFH17 and BCCW18. CCDDHNR (CCDDHNR18) combine the use of Neyman orthogonality condition with cross fitting to provide a simple yet widely applicable framework that covers a large class of models under i.i.d. settings. The DML is also compatible with various types of machine learning based methods for nuisance parameter estimation.

Driven by the need from empiricists, the literature on cluster robust inference has a long history in econometrics. For recent review of the literature, see, e.g., CM15 and MacKinnon2019. On the other hand, coping with cross-sectional dependence using a multiway cluster robust variance estimator is a relatively recent phenomenon. CGM11 first provide a multiway cluster robust variance estimator for linear regression models without imposing additional parametric assumptions on the intra-cluster correlation structure. This variance estimator has significantly reshaped the landscape of econometric practices in applied microeconomics in the past decade.\footnote{As of December 31, 2019, CGM11 has received over 2,500 citations. The majority of such citations came from applied economic papers.} In contrast to the popularity among empirical researchers, theoretical justification of the validity of this type of procedures was lagging behind. The first rigorous treatment of asymptotic properties of multiway cluster robust estimators are established by Menzel17 using the Aldous-Hoover representation under the assumptions of separable exchangeability and dissociation. The asymptotic theory of Menzel17 covers both non-degenerate and degenerate cases. Focusing on non-degenerate situations, DDG18 further extend this approach to a general empirical process theory.\footnote{See also DDG19 for further generalization of the empirical process theory for dyadic data under joint exchangeability assumption.} Using this asymptotic framework, MacKinnonNielsenWebb2019 study linear regression models under the non-degenerate case and examine the validity of several types of wild bootstrap procedures and the robustness of multiway cluster robust variance estimators under different cluster sampling settings.

Despite of the popularity of both machine learning and cluster robust inference among empirical researchers, relatively limited cluster robust inference results exist for machine learning based methods. Inference for machine learning based methods with one-way clustering is studied by BCHK16, Kock2016, KockTang2018, SGCT18 and HansenLiao19 for different variations of regularized regression estimators and AtheyWager19 for random forests. ChiangSasaki2019 investigate the performance of lasso and post-lasso in the partially linear model setting of BCH14 under multiway cluster sampling. To our best knowledge, there is no general machine learning based procedures with known validity under multiway cluster sampling environments.

Overview

Setup

Suppose that the researcher observes a sample $\left\{\left. W_{ij} \right\vert i \in \{1,...,N\}, j \in \{1,...,M\}\right\}$ of double-indexed observations of size $NM$. Let $P$ denote the probability law of $\{W_{ij}\}_{ij}$, and let ${\rm E}_{P}$ denote the expectation with respect to $P$. Let $\underline C= N \wedge M$ denote the sample size in the smaller dimension. We consider two-way clustering where each cell contains one observation for simplicity of notations, but results for higher cluster dimensions and random cluster sizes can be obtained at the expense of involved notations -- see Appendix (ref) for a general case.

The structural model is assumed to entail the moment restriction

for some score $\psi$ that depends on a low-dimensional parameter vector $\theta \in \Theta \subset \mathbbm R^{d_\theta}$ and a nuisance parameter $\eta \in T$ for a convex subset $T$ of a normed linear space. The nuisance parameter $\eta$ may be finite-, high-, or infinite-dimensional, and its true value is denoted by $\eta_0 \in T$. In this setup, the true value of the low-dimensional target parameter, denoted by $\theta_0 \in \Theta$, is the object of interest.

Let $\widetilde T=\{\eta - \eta_0 : \eta \in T\}$, and define the Gateaux derivative map $D_r: \widetilde T \rightarrow \mathbbm R^{d_\theta}$ by

for all $r\in[0,1)$. Also denote its limit by

We say that the Neyman orthogonality condition holds at $(\theta_0,\eta_0)$ with respect to a nuisance realization set $\mathcal T_n \subset T$ if the score $\psi$ satisfies ((ref)), the pathwise derivative $D_r[\eta-\eta_0]$ exists for all $r\in[0,1)$ and $\eta\in \mathcal T_n$, and the orthogonality equation

holds for all $\eta\in \mathcal T_n$. Furthermore, we also say that the $\lambda_n$ Neyman near-orthogonality condition holds at $(\theta_0,\eta_0)$ with respect to a nuisance realization set $\mathcal T_n\subset T$ if the score $\psi$ satisfies ((ref)), the pathwise derivative $D_r[\eta-\eta_0]$ exists for all $r\in[0,1)$ and $\eta\in \mathcal T_n$, and the orthogonality equation

holds for all $\eta\in \mathcal T_n$ for some positive sequence $\{\lambda_n\}_n$ such that $\lambda_n=o(\underline C^{-1/2})$.

Throughout, we will consider structural models satisfying the moment restriction ((ref)) and either form of the Neyman orthogonality conditions, ((ref)) or ((ref)). Consider linear Neyman orthogonal scores $\psi$ of the form

A generalization to nonlinear score follows from linearization with Gateaux differentiability as in Section 3.3 of CCDDHNR (CCDDHNR18). We focus on linear scores as they cover a wide range of applications.

The Multiway Double/Debiased Machine Learning

For the class of models introduced in Section (ref), we propose a novel $K^2$-fold multiway cross fitting procedure for estimation of $\theta_0$. For any $r \in \mathbb N$, we use the notation $[r]=\{1,...,r\}$. With a fixed positive integer $K$, randomly partition $[N]$ into $K$ parts $\{I_1,...,I_K\}$ and $[M]$ into $K$ parts $\{J_1,...,J_K\}$. For each $(k,\ell) \in [K]^2$, obtain an estimate $$\widehat \eta_{k\ell}=\widehat \eta\left((W_{ij})_{(i,j)\in ([N]\setminus I_k )\times ([M]\setminus J_\ell)}\right)$$ of the nuisance parameter $\eta$ by some machine learning method (e.g., lasso, post-lasso, elastic nets, ridge, deep neural networks, and boosted trees) using only the subsample of those observations with multiway indices $(i,j)$ in $([N]\setminus I_k ) \times ([M]\setminus J_\ell)$. In turn, we define $\widetilde \theta$, the multiway double/debiased machine learning (multiway DML) estimator for $\theta_0$, as the solution to

where $\mathbbm E_{n,k\ell} [f(W)] = \frac{1}{|I_k||J_\ell|}\sum_{(i,j)\in I_k\times J_\ell} f(W_{ij})$ denotes the subsample empirical expectation using only the those observations with multiway indices $(i,j)$ in $I_k \times J_\ell$.

We call this procedure the $K^2$-fold multiway cross fitting. Note that, for each $(k,\ell)\in [K]^2$, the nuisance parameter estimate $\widehat\eta_{k\ell}$ is computed using the subsample of those observations with multiway indices $(i,j) \in ([N]\setminus I_k ) \times ([M]\setminus J_\ell)$, and in turn the score term $\mathbbm E_{n,k\ell}[\psi(W; \cdot,\widehat\eta_{k\ell})]$ is computed using the subsample of those observations with multiway indices $(i,j) \in I_k \times J_\ell$. This two-step computation is repeated $K^2$ times for every partitioning pair $(k,\ell)\in [K]^2$. Figure (ref) illustrates this $K^2$-fold cross fitting for the case of $K=2$ and $N=M=4$, where the cross fitting repeats for $K^2 (= 2^2 = 4)$ times.

We propose to estimate the asymptotic variance of $\sqrt{\underline C}(\widetilde\theta-\theta_0)$ by

where $\widehat \Gamma$ and $\widehat J$ are given by

accounting for multiway cluster dependence. For a $d_\theta$-dimensional vector $r$, the $(1-a)$ confidence interval for the linear functional $r'\theta_0$ can be constructed by

Example: Partially Linear IV Model with Multiway Cluster Sample

For an illustration, consider as a concrete example the partially linear IV model (cf. Okui, Small, Tan and Robins, OkuiSmallTanRobins2012 ; CCDDHNR, CCDDHNR18, Section 4.2) adapted to the multiway cluster sample data:

A researcher observes the random variables $Y_{ij}$, $D_{ij}$, $X_{ij}$, and $Z_{ij}$, which are typically interpreted as the outcome, endogenous regressor, exogenous regressors, and instrumental variable, respectively. The low-dimensional parameter vector $\theta_0$ is an object of interest.

A Neyman orthogonal score $\psi$ for such model is given by

as in OkuiSmallTanRobins2012 and CCDDHNR (CCDDHNR18), where $w=(y,d,x,z)$, $\eta=(g_1,g_2,m)$ and $g_1$, $g_2$, $m\in L^2(P)$. It is straightforward to verify that this score satisfies both the moment restriction ((ref)), ${\rm E}_{P}[\psi(W_{11};\theta_0,\eta_0)]=0$, and the Neyman orthogonality condition ((ref)), $\partial_\eta {\rm E}_{P} \psi(W_{11};\theta_0,\eta_0)[\eta - \eta_0]=0$ for all $\eta \in \mathcal{T}_n$ at $\eta_0=(g_{10},g_{20},m_0)$, where $g_{10}(X)={\rm E}_{P}[Y|X]$, $g_{20}(X)={\rm E}_{P}[D|X]$, and $m_0(X)={\rm E}_{P}[Z|X]$.

The following algorithm is our proposed multiway DML procedure introduced in Section (ref), specifically applied to this partially linear IV model.

For the sake of concreteness, we present this algorithm specifically based on lasso (in the three sub-steps under step 2), but another machine learning method (e.g., post-lasso, elastic nets, ridge, deep neural networks, and boosted trees) may be substituted for lasso.

Theory of the Multiway DML

In this section, we present formal theories to guarantee that the multiway DML method proposed in Section (ref) works. We first fix some notations for convenience. The two-way sample sizes $(N,M) \in \mathbb{N}^2$ will be index by a single index $n \in \mathbb{N}$ as $(N,M) = (N(n),M(n))$ where $M(n)$ and $N(n)$ are non-decreasing in $n$ and $M(n)N(n)$ is increasing in $n$. With this said, we will suppress the index notation and write $(N,M)$ for simplicity. Let $\{\mathcal P_n\}_n$ be a sequence of sets of probability laws of $\{W_{ij}\}_{ij}$ -- note that we allow for increasing dimensionality of $W_{ij}$ in the sample size $n$. Let $P=P_{n}\in \mathcal P_n$ denote the law with respect to sample size $(N,M)$. Throughout, we assume that this random vector $W_{ij}$ is Borel measurable. Recall the notations $\underline C =N\wedge M$, $\mu_N=\underline C/N$, and $\mu_M=\underline C/M$, and suppose that $\mu_N\to \bar \mu_N$, $\mu_M\to \bar \mu_M$. We write $a \lesssim b$ to mean $a \leq cb$ for some $c > 0$ that does not depend on $n$. We also write $a \lesssim_P b$ to mean $a = O_P(b)$. For any finite dimensional vector $v$, $\|v\|$ denotes the $\ell_2$ or Euclidean norm of $v$. For any matrix $A$, $\|A\|$ denotes the induced $\ell_2$-norm of the matrix. For any set $B$, $|B|$ denotes the cardinality of the set.

We state the following assumption on multiway clustered sampling.

Recall that we focus on the linear Neyman orthogonal score of the form

Let $c_0>0$, $c_1>0$, $s>0$, $q\ge 4$ be some finite constants with $c_0\le c_1$. Let $\{\delta_n\}_{n\ge 1}$ (estimation errors) and $\{\Delta_n\}_{n\ge 1}$ (probability bounds) be sequences of positive constants that converge to zero such that $\delta_n \ge \underline C^{-1/2}$. Let $K\ge 2$ be a fixed integer. Let $W_{00}$ denote a copy of $W_{11}$ that is independent from the data and the random set $\mathcal T_n$ of nuisance realization. With these notations, we consider the following assumptions.

The following result presents the main theorem of this paper, establishing the linear representation and asymptotic normality of the multiway DML estimator. It corresponds to Theorem 3.1 of CCDDHNR (CCDDHNR18), and is an extension of it to the case of multiway cluster sampling.

As is commonly the case in practice, we need to estimate the unknown asymptotic variance. The following theorem shows the validity of our proposed multiway DML variance estimator.

Theorems (ref) and (ref) can be used for constructing confidence intervals.

As in Section 3.4 of CCDDHNR (CCDDHNR18), we can also repeatedly compute multiway DML estimates and variance estimates $S$-times for some fixed $S\in \mathbbm N$ and consider the average or median of the estimates as the new estimate. This does not have an asymptotic impact, yet it can reduce the impact of a random sample splitting on the estimate.

Simulation Studies

Simulation Setup

Consider the partially linear IV model introduced in Section (ref). We specifically focus on the following high-dimensional linear representations

where the parameter values are set to $\theta_0 = \pi_{10} = 1.0$ and $\zeta_0 = \pi_{20} = \xi_0 = (0.5,.0.5^2,\cdots,0.5^{\text{dim}(X)})'$ for some large $\text{dim}(X)$. The primitive random vector $(X_{ij}',\epsilon_{ij},\upsilon_{ij},V_{ij})'$ is constructed by

with two-way clustering weights $(\omega_1^X,\omega_2^X)$, $(\omega_1^\epsilon,\omega_2^\epsilon)$, $(\omega_1^\upsilon,\omega_2^\upsilon)$, and $(\omega_1^V,\omega_2^V)$, where $\alpha_{ij}^X$, $\alpha_{i}^X$, and $\alpha_{j}^X$ are independently generated according to

$(\alpha_{ij}^\epsilon,\alpha_{ij}^\upsilon)'$, $(\alpha_{i}^\epsilon,\alpha_{i}^\upsilon)'$, and $(\alpha_{j}^\epsilon,\alpha_{j}^\upsilon)'$ are independently generated according to

and $\alpha_{ij}^V$, $\alpha_{i}^V$, and $\alpha_{j}^V$ are independently generated according to

The weights $(\omega_1^X,\omega_2^X)$, $(\omega_1^\epsilon,\omega_2^\epsilon)$, $(\omega_1^\upsilon,\omega_2^\upsilon)$, and $(\omega_1^V,\omega_2^V)$ specify the extent of dependence in two-way clustering in $X_{ij}$, $\epsilon_{ij}$, $\upsilon_{ij}$, and $V_{ij}$, respepctively. The parameter $s_X$ specifies the extent of collinearity among the high-dimensional regressors $X_{ij}$. The parameter $s_{\epsilon\upsilon}$ specifies the extent of endogeneity. We set the values of these parameters to $(\omega_1^X,\omega_2^X) = (\omega_1^\epsilon,\omega_2^\epsilon) = (\omega_1^\upsilon,\omega_2^\upsilon) = (\omega_1^V,\omega_2^V) = (0.25, 0.25)$ and $s_X = s_{\epsilon\upsilon} = 0.25$.

Results

Monte Carlo simulations are conducted with 2,500 iterations for each set. Table (ref) reports simulation results. The first four columns in the table indicate the data generating process ($N$, $M$, $\underline C$, and dim$(X)$). The next column indicates the integer $K$ for our $K^2$-fold cross fitting method. We use $K=2$ and $3$ in the simulations for the displayed results, since $2^2 (\approx 5)$ and $3^2 (\approx 10)$ are close to the common numbers of folds used in cross fitting in practice. The next column indicates the machine learning method for estimation of $\widehat\eta_{k\ell}$. We use the ridge, elastic net, and lasso. The last four columns of the table report Monte Carlo simulation statistics, including the bias (Bias), standard deviation (SD), root mean square error (RMSE), and coverage frequency for the nominal probability of 95% (Cover).

For each covariate dimension $\text{dim}(X) \in \{100,200\}$, for each choice $K \in \{2,3\}$ for the number $K^2$ of multiway cross fitting, and for each of the three machine learning methods, we observe the following patterns as the effective sample size $\underline C=N \wedge M$ increases: 1) the bias tends to zero; 2) the standard deviation decreases approximately at the $\sqrt{\underline C}$ rate; and 3) the coverage frequency converges to the nominal probability. These results confirm the theoretical properties of the proposed method. We ran several other sets of simulations besides those displayed in the table, and this pattern remains the same across different sets.

Comparing the results across the three machine learning methods, we observe that the ridge entails larger bias and smaller variance relative to the elastic net and lasso in finite sample. This makes the coverage frequency of the ridge less accurate compared with the elastic net and lasso. This result is perhaps specific to the data generating process used for our simulations. On one hand, the choice $K=3$ (i.e., $9$-fold) of the multiway cross fitting contributes to mitigating the large bias of the ridge relative to the choice $K=2$, and hence $K=3$ produces more preferred results for the ridge. On the other hand, the choice $K=2$ tends to yield preferred results in terms of coverage accuracy for the elastic net and lasso. In light of these results, we recommend the elastic net or lasso along with the use of $2^2$- fold (i.e., $4$-fold) cross fitting. This number of folds in cross fitting is in fact similar to that recommended by CCDDHNR (CCDDHNR18) for i.i.d. sampling -- see their Remark 3.1 where they recommend 4- or 5-fold cross fitting.

Empirical Illustration: Demand Analysis with Market Share Data

Let us revisit the demand model of Example (ref) in Section (ref). Recall that, for the consumer demand model of Berry94 introduced in Example (ref), LuShiTao19 derive the partial-linear equation

for estimation of $\theta_0$, where $Y_{ij} = \log( S_{ij} ) - \log( S_{0j} )$ denotes the observed log share of product $i$ relative to the log of the outside share in market $j$, $D_{ij}$ denotes the log price of product $i$ in market $j$, and $X_{ij}$ denotes a vector of observed attributes of product $i$ in market $j$. To deal with the likely endogeneity of $D_{ij}$, researchers often use instruments $Z_{ij}$ such that ${\rm E}_{P}[\epsilon_{ij}|X_{ij},Z_{ij}]=0$. Such instruments often consist of observed attributes of other products in the market.

The implied equation ((ref)) together with this mean independence assumption yields the reduced-form model ((ref)). Furthermore, we write the innocuous nonparametric projection equation ((ref)). Therefore, we apply Algorithm (ref) in Section (ref) for the two-way cluster robust DML estimation of $\theta_0$ with a robust standard error.

We present an application of the proposed algorithm to the U.S. automobile data of BLP95. The sample consists of unbalanced two-way clustered observations with $N=557$ models of automobiles and $M=20$ markets. The observed attributes $X_{ij}$ consist of horsepower per weight, miles per dollar, miles per gallon, and size. The instrument $Z_{ij}$ is defined as the sum of the values of these attributes of other products.

For the purpose of highlighting the effect of clustering assumptions, we report estimates and standard errors under the zero-way cluster robust DML (based on the i.i.d. assumption) and the one-way cluster robust DML (based on clustering along each of the product and market dimensions), as well as the two-way cluster robust DML (along both of the product and market dimensions). The number $K=4$ of folds of cross fitting is used for the zero- and one-way cluster robust DML, while the number $K^2=4$ of folds of two-way cross fitting is used for the two-way cluster robust DML following the recommendations from Section (ref) and those by CCDDHNR (CCDDHNR18, Remark 3.1). To mitigate the uncertainty induced by sample splitting, we compute estimates based on the average of ten rerandomized DML following CCDDHNR (CCDDHNR18, Section 3.4) with variance estimation according to CCDDHNR (CCDDHNR18, Equation 3.13) adapted to our two-way cluster-robustness.

Table (ref) summarizes the results. For each of the zero-, one-, and two-way cluster robust DML, both the point estimates and standard errors are similar across all the choices of instrument. Furthermore, the point estimates are also similar across all of the zero-, one-, and two-way cluster robust DML. On the other hand, the standard errors tend to increase as the assumed number of ways of clustering increases. In other words, the zero-way cluster robust DML reports the smallest standard error while the two-way cluster robust DML reports the largest standard error. To robustly account for possible cross-sectional dependence of observations in such two-way cluster sampled data as this market share data, we recommend that researchers use the two-way cluster robust DML although it may incur larger standard errors as is the case with this application.

Conclusion

In this paper, we propose a multiway DML procedure based on a new multiway cross fitting algorithm. This multiway DML procedure is valid in the presence of multiway cluster sampled data, which is frequently used in empirical research. We present an asymptotic theory showing that multiway DML is valid under nearly identical reguarity conditions to those of CCDDHNR (CCDDHNR18). The proposed method covers a large class of econometric models as is the case with CCDDHNR (CCDDHNR18), and is compatible with various machine learning based estimation methods. Simulation studies indicate that the proposed procedure has attractive finite sample performance under various multiway cluster sampling environments for various machine learning methods. To accompany the theoretical findings, we provide easy-to-implement algorithms for multiway DML. Such algorithms are readily implementable using existing statistical packages.

There are a couple of possible directions for future research. First, whereas we focused on linear orthogonal scores that cover a wide range of applications, it may be possible to develop a method and theories for non-linear orthogonal scores as in CCDDHNR (CCDDHNR18; Section 3.3). Second, whereas we focused on unconditional moment restrictions, it may be possible and will be important to develop a method and theories for conditional moment restrictions AiChen2003,AiChen2007,ChenLintonKeilegom2003,ChenPouzo2015. We leave these and other extensions for future research.