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.
46,736 characters · 17 sections · 53 citation commands
Lasso under multi-way clustering: Estimation and Post-selection Inference
\address[Harold Chiang]{Department of Economics, Vanderbilt University, United States}\address[Yuya Sasaki]{Department of Economics, Vanderbilt University, United States}
This paper studies a method of estimation and post-selection inference for regression parameters in high-dimensional linear models by lasso under multi-way clustering. The objective is motivated by recently increasing demands from applied economic research. On one hand, economists often use multi-way cluster sampled data. Examples include, but are not limited to, network data, matched employer-employee data, matched student-teacher data, scanner data where observations are double-indexed by stores and products, market share data where observations are double-indexed by market and products, and growth/development data where observations are double-indexed by ethnicity and geographical units -- see Section (ref) for specific applications of the last example. On the other hand, researchers also often use machine learning methods of estimation and inference for high-dimensional models in today's big data environments. There are a number of useful methods in the literature that deal with each of these two issues (multi-way clustering and high dimensionality) separately, but the existing methods do not seem to provide a solution to dealing with both of these practically relevant issues simultaneously. In this light, we present lasso under multi-way clustering, and propose a post-selection inference method for regression parameters under this sampling assumption.
In the important branch of the literature following the seminal work by BCCH12, post-selection inference with lasso has been widely studied under various settings by BCH14\footnote{See also BCH11.}, JM14, vdG14, ZZ14, BCK15, and many others. For empirical researchers, lasso has become a powerful machine learning tool under data-rich environments. Most of the papers in this literature assume i.i.d. or independent sampling. In many empirical applications, it is sometimes more plausible to assume multi-way cluster sampling (e.g., network data, matched employer-employee data, and matched student-teacher data). Building upon BCCH12, BC13 and BCH14, this paper generalizes lasso and post-double-selection procedure by allowing for multi-way cluster sampling. To our best knowledge, the present paper is the first in the literature of high-dimensional models to consider lasso under multi-way cluster sampling.
The influential work by CGM11 proposes multi-way cluster-robust inference methods for linear and nonlinear regression models -- also see CM15 for a survey. Formal analysis of asymptotic properties and bootstrap validity under multi-way clustering is studied by menzel17 using the Aldous-Hoover representation -- see kallenberg05 for example. Under the assumptions of separable exchangeability, the method of menzel17 covers both degenerate and non-degenerate cases. Using the same representation, while focusing on the non-degenerate cases, DDG18 develop empirical process theory under multi-way cluster sampling which applies to a large class of econometric models. Building upon the asymptotic framework of these two papers, MNW19 propose several wild bootstrap procedures for linear regression models, and examine their finite-sample performances under several different cluster sampling scenarios. In this paper, we take advantage of the innovations by these preceding papers to develop a multi-way cluster-robust inference method for high-dimensional models. To our best knowledge, the present paper is the first in this literature on multi-way clustering to consider high-dimensional models.
The rest of this paper is organized as follows. Section (ref) introduces the model. Section (ref) presents an overview of the proposed methodology. Section (ref) discusses a formal asymptotic theory. Section (ref) presents an extension of the baseline results to cases of heterogeneous cluster sizes. Section (ref) presents simulation studies. Section (ref) presents an empirical illustration with development and growth economics. Section (ref) concludes. The appendix contains mathematical proofs and auxiliary lemmas.
Consider the high-dimensional regression model
where $Y_{ij}$ is an observed outcome variable, $(D_{ij},X_{ij}')'$ is an observed vector of regressors, and $R_{ij}^Y$ is an approximation error for the unit of observation with the double index $(i,j)$. We set $\alpha$ as a scalar parameter of interest. The dimension $p$ of the nuisance parameter vector $\beta \in \mathbb R^p$ is potentially increasing in the sample size. Following the literature on high-dimensional post-selection inference BCH14, we also consider the auxiliary projection
where $R^D$ is an approximation error. The dimension $p$ of the nuisance parameter vector $\gamma\in \mathbb R^p$ is the same as that of $\beta$, and is potentially increasing in the sample size.
In the absence of two-way clustering, the system ((ref))--((ref)) would be the same as the model considered in BCH14. We first consider two-way clustering where each cell contains one observation. Section (ref) presents an extension to the case of heterogeneous cluster sizes.
In this section, we present an overview of the proposed method, namely estimation and post-selection inference. Formal theoretical justifications are discussed in Section (ref).
A researcher observes a sample $\left\{\left. (Y_{ij},D_{ij},X_{ij}') \right\vert i \in \{1,...,N\}, j \in \{1,...,M\}\right\}$ of size $NM$. The estimation procedure consists of two steps. First, define the lasso estimates for ((ref)) and ((ref)) by
respectively, for some regularization parameters $\lambda_1$ and $\lambda_2$, valid choices of which are discussed in the statement of Theorem (ref) ahead. Denote the supports of the lasso estimates by $\widehat I_1={\rm support}(\widehat \beta)$ and $\widehat I_2={\rm support}(\widehat \gamma)$, and let $\widehat I= \widehat I_1 \cup \widehat I_2$. In the second step, define the post-double-selection lasso estimate $\widetilde \alpha$ by
Let $\underline C =N\wedge M$, $\mu_N=\underline C/N$ and $\mu_M=\underline C/M$. Under suitable conditions to be formally stated in Section (ref), we have the asymptotic normality
where the asymptotic variance is given by $\sigma^2= Q^{-1}\Gamma Q^{-1}$ with
and $\bar\mu_N$ and $\bar\mu_M$ denoting the limits of $\mu_N$ and $\mu_M$, respectively.
The asymptotic variance is estimated by the sample counterpart $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$, where
$\widehat v_{ij}= D_{ij}- X_{ij}'\widehat \gamma$, and $\widehat \varepsilon_{ij}=Y_{ij}- \widehat \alpha D_{ij} - X_{ij}'\widehat\beta$.
In summary, we propose to report the post-double-selection lasso estimate $\widetilde\alpha$ as an estimate of $\alpha$ with its standard error given by $\widehat\sigma / \sqrt{\underline C}$. The $\alpha^\ast$-level confidence interval can be constructed as $\left[\widetilde\alpha + \Phi^{-1}(\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C},\widetilde\alpha + \Phi^{-1}(1-\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C}\right]$, where $\Phi^{-1}$ denotes the quantile function of the standard normal distribution.
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$. For simplicity, each size of intersection $n_{ij}$ is assumed to be uniformly bounded by a positive integer $\bar n$ that is independent of $n$. With this said, we will suppress the index notation and write $(N,M)$ for simplicity. We fix a number of notations. For each $n$, let ${\mathrm{P}}_{n}$ denote the law with respect to sample size $(N,M)$ -- note that we allow the dimension $p$ of $X_{ij}$ to grow with $n$. Let $a := p\vee (NM)$. Also recall the notations $\underline C =N\wedge M$, $\mu_N=\underline C/N$, and $\mu_M=\underline C/M$ from Section (ref). We use the short-hand notation $[k] = \{1,...,k\}$ and $[k]^c=\mathbb{N} \backslash [k]$ for any $k \in \mathbb{N}$. For a sequence $(t_{ij})_{i\in[N],j\in [M]}$, denote $\|t_{ij}\|_n = \sqrt{\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} t_{ij}^2 }$. Thus, $\|X'_{ij}\delta\|_n = \sqrt{\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} \delta' X_{ij}X_{ij}'\delta}$ is the prediction norm of $\delta$. Let $\|A\|_\infty=\max_{k,l}|A_{k,l}|$ denote the max norm of matrix $A$. We write $a \lesssim b$ to mean $a \leqslant cb$ for some $c > 0$ that does not depend on $n$. We also write $a \lesssim_{\mathrm{P}} b$ to mean $a = O_P(b)$. We write $Z_{ij}=(Y_{ij},D_{ij},X'_{ij})'$ for the $(p+2)$-dimensional random vector in data. Throughout, we assume that this random vector $Z_{ij}$ is Borel measurable -- see kallenberg05. With these notations, we state the following four assumptions.
In this section, we show that an empirical process in multi-way clustered samples can be represented as a sum of independent variables via H\'ajek projection. Furthermore, its variance can be shown to be approximated by covariances of observed variables.
For any $f:{\rm support}(Z)\to \mathbb R$, we let
denote its empirical process.
A proof of this lemma can be found in Appendix (ref). The first part of the lemma shows that an empirical process $\mathbb{G}_C f$ under multi-way cluster sampling can be represented as a sum of independent unobserved variables via H\'ajek projection $H_n f$. While $U_{i0}$ and $U_{0j}$ are unobserved, the second part of this lemma in turn shows that the variance of the H\'ajek projection can be approximated by covariances of observed variables. Note that, since $H_n f$ is a H\'ajek projection, the lemma implies $\frac{\mathbb{G}_C f}{\sqrt{V(\mathbb{G}_C f)}}=\frac{H_n f}{\sqrt{V(H_n f)}}+o_{\mathrm{P}}(1)$ if $\bar\mu_N Cov(f(Z_{11}),f(Z_{12})) + \bar\mu_M Cov(f(Z_{11}),f(Z_{21}))$ is bounded and bounded away from zero uniformly in $n$.
Our Lemma (ref) can be seen as an extension to Lemma D.2 in DDG18. Specifically, while DDG18 consider a fixed data generating process over the sample size $n$, our Lemma (ref) allows the data generating process to vary with $n$ in particular for the sake of accommodating the increasing of dimensionality $p$ for high-dimensional models. The lemma serves as a main building block for all the asymptotic results to be presented ahead.
We next show the convergence rates of the lasso estimator $(\widehat\alpha,\widehat\beta',\widehat\gamma')'$ and the post-lasso estimator $(\widetilde\alpha,\widetilde\beta',\widetilde\gamma')'$ under multi-way clustering.
A proof can be found in Appendix (ref), and is based on the previous result (Lemma (ref)). In the multi-way sampling, this lemma can be viewed as a counterpart of Lemma 6 and Lemma 7 in BCCH12.
In this section, we present the main result of this paper. The limit normal distribution of the post-double-selection lasso estimate $\widetilde\alpha$ is established based on the previous two results (Lemma (ref) and Theorem (ref)).
A proof can be found in Appendix (ref). This result provides a theoretical justification for the asymptotic variance proposed in the overview in Section (ref). In practice, we do not know the components, $Q$ and $\Gamma$, of the asymptotic variance. The following subsection proposes estimators of them.
In this section, we propose an analog variance estimator. The components, $Q$ and $\Gamma$, of the asymptotic variance can be estimated by
respectively, where $\widehat v_{ij}= D_{ij}- X_{ij}'\widehat \gamma$ and $\widehat \varepsilon_{ij}=Y_{ij}- \widehat \alpha D_{ij} - X_{ij}'\widehat\beta$ are the residuals. With these component estimators, we propose that the asymptotic variance $\sigma^2 = Q^{-1} \Gamma Q^{-1}$ be estimated by $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$. The following theorem provides a theoretical support for this variance estimator.
A proof is found in Appendix (ref). In light of this result, we propose to compute the standard error by $\widehat\sigma / \sqrt{\underline C}$. Similarly, in light of this result together with Theorem (ref), we propose to construct the $\alpha^\ast$-level confidence interval by $\left[\widetilde\alpha + \Phi^{-1}(\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C},\widetilde\alpha + \Phi^{-1}(1-\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C}\right]$, where $\Phi^{-1}$ denotes the quantile function of the standard normal distribution.
Thus far, we focus on the case where each cluster contains one observation. In this section, we presented an extension of the baseline results to situations where the numbers of observations are heterogeneous across clusters. Suppose that we have $n_{ij}$ observations for each cell $(i,j) \in [N] \times [M]$, where $n_{ij}$ is a random variable that is allowed to depend on $(X_{ij,\ell})_{\ell\ge 1}$. To deal with the situation of $n_{ij}=0$, for any sequence $(t_\ell)_{\ell \ge 1}$, define $\sum_{\ell=1}^0 t_\ell=0$. Consider the model
where $Y_{ij,\ell}$ is an observed outcome variable, $(D_{ij,\ell},X_{ij,\ell}')'$ is an observed vector of regressors, and $R_{ij,\ell}^Y$ is an approximation error for the unit $\ell\in [n_{ij}]$ with the double index $(i,j)$ indicating $i$-th cluster in the first clustering dimension and $j$-th cluster in the second clustering dimension. Using matrix notations, we can rewrite the model as
where each of $Y_{ij}=(Y_{ij,\ell})_{\ell\in [n_{ij}]}$, $D_{ij}=(D_{ij,\ell})_{\ell\in [n_{ij}]}$, $R^Y_{ij}=(R^Y_{ij,\ell})_{\ell\in [n_{ij}]}$, and $\varepsilon_{ij}=(\varepsilon_{ij,\ell})_{\ell\in [n_{ij}]}$ is of dimension $n_{ij}\times 1$, and $X_{ij}=(X_{ij,\ell}')_{\ell\in [n_{ij}]}$ is of dimension $n_{ij}\times p$. We similarly write the accompanying auxiliary projection as
where $R^D$ is of dimension $n_{ij}\times 1$ representing approximation errors, and $v_{ij}$ is of dimension $n_{ij}\times 1$ representing projection errors.
Under this setting, the first step of estimation procedure consists of
In turn, the second-step estimates are obtained by
The asymptotic variance estimator for $\widetilde\alpha$ is given by $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$, where
$\widehat v_{ij}= D_{ij}- X_{ij}\widehat \gamma$, and $\widehat \varepsilon_{ij}=Y_{ij}- \widehat \alpha D_{ij} - X_{ij}\widehat\beta$.
We now formally state assumptions for the extended theory to support the asymptotic validity of this procedure. Define $W_{ij}=(n_{ij},(Z_{ij,\ell})_{\ell \ge 1})$ and $\ddot M=\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M} X_{ij}'X_{ij}$.
The following statement provides a theoretical guarantee for the estimation and inference procedure for the extended model outlined above.
A proof of Corollary (ref) closely follows that of the results in Section (ref), and are therefore omitted. The key difference is that we now apply Aldous-Hoover representation on $W_{ij}$ rather than on $Z_{ij}$.\footnote{For more insights on this extension, see Section 3.1 of DDG19.}
In this section, we present simulation studies of finite-sample performance of the proposed method of estimation and post-selection inference. We compare the performance of our method against existing alternatives from the lasso literature that do not account for multi-way clustering.
We consider the linear model
The parameter values are fixed at $(\alpha,\beta')' = \left(0.5, 0.5^2, \cdots, 0.5^{\text{dim}(X)+1}\right)'$. The random vector $(D_{ij}, X_{ij}', \varepsilon_{ij})$ is constructed by
with two-way clustering weights $(\omega^{x}_1, \omega^{x}_2)$ and $(\omega^{\varepsilon}_1,\omega^{\varepsilon}_2)$, where $\upsilon^{x}_{ij}$, $\upsilon^{x}_{i}$, and $\upsilon^{x}_{j}$ are independently generated according to
and $\upsilon^{\varepsilon}_{ij}$, $\upsilon^{\varepsilon}_{i}$, and $\upsilon^{\varepsilon}_{j}$ are independently generated according to
Note that the weights $(\omega^{x}_1, \omega^{x}_2)$ and $(\omega^{\varepsilon}_1,\omega^{\varepsilon}_2)$ specify the extent of dependence in two-way clustering in $(D_{ij},X_{ij}')$ and $\varepsilon_{ij}$, respectively. Also, the parameter $\rho$ specifies the extent of collinearity among the high-dimensional covariates $(D_{ij},X_{ij}')$. We set $(\omega^{x}_1, \omega^{x}_2) = (0.25,0.25)$, $(\omega^{\varepsilon}_1,\omega^{\varepsilon}_2) = (0.25,0.25)$, and $\rho = 0.50$.
We compare the performance of our multi-way cluster-robust variance estimator with two existing alternative benchmarks. One is the heteroskedasticity robust variance estimator (such as the one in BCH14) without accounting for cluster sampling, i.e., $\Gamma$ is estimated by
We will refer to this variance estimator $\widehat Q^{-1} \widehat \Gamma_{HC} \widehat Q^{-1}$ as the `0-Way' estimator. The other is the one-way cluster-robust variance estimator (similar to those of BCHK16 and kock16) clustered at one (e.g., second) dimension, i.e., $\Gamma$ is estimated by
We will refer to this variance estimator $\widehat Q^{-1} \widehat \Gamma_{CR} \widehat Q^{-1}$ as the `1-Way' estimator.
Table (ref) summarizes simulation results. The first two columns indicate the two-way sample sizes $(N,M)$. The third column indicates the dimension (Dim) of $(\alpha,\beta')'$. The next four columns report simulation statistics for $\widetilde\alpha$. These statistics include the average (Avg), bias (Bias), standard deviation (SD), and root mean square error (RMSE). The last three columns report 95% coverage frequencies of $\alpha$ based on three variance estimators. The first is the heteroskedasticity robust variance estimator (0-Way). The second is the one-way cluster-robust variance estimator (1-Way). The third is our multi-way cluster-robust variance estimator (2-Way). The results are based on 25,000 Monte Carlo iterations for each row in the table.
In view of the statistics columns, observe that the post-double-selection lasso estimate $\widetilde\alpha$ behaves well in larger sample sizes (e.g., $N, M \geqslant 20$) both in terms of bias and variance. Next, observe the 95% coverage frequencies by the three alternative variance estimators. Both the 0-Way and 1-Way variance estimators significantly underestimate the variances of the post-double-selection lasso estimate $\widetilde\alpha$. On the other hand, the coverage frequency based on our 2-Way variance estimator approaches the nominal probability (95%) as the sample size increases. These results demonstrate that, when the true sampling process entails multi-way clustering, traditional variance estimators may bias the inference and our multi-way cluster-robust variance estimator performs robustly well.
In this section, we illustrate our proposed method with applications to a couple of empirical studies. There is a sequence of recent growth and development economic studies using empirical data that are clustered at ethnic and geographical levels NunnWantchekon11,Michalopoulos_Papaioannou2013,MP14,MP16,gershman16,anderson18,dickens18. The next two subsections present how our method can enrich the model flexibility and robustness of such studies, focusing on the cases of NunnWantchekon11 and Michalopoulos_Papaioannou2013.
NunnWantchekon11 analyze the effects of slave trade on mistrust in Africa, controlling for various demongraphic and geographical covariates including age, age squared, ethnic fractionalization, gender, urban residence, occupation, religion, and living conditions as well as country fixed effects in their baseline model. Estimates of these effects are obtained by running regressions with a sample that pools $n_{ij}$ individuals $\ell \in [n_{ij}]$ in ethnic group $i$ and districts $j$ across the cells $(i,j) \in [N] \times [M]$ of $N (=185)$ ethnic groups and $M (=1257)$ districts. Standard errors are computed by the two-way cluster-robust method of CGM11 for the ethnic group and district as two ways of clustering.
With our proposed method that is applicable to both high-dimensional models and multi-way clustering, they could consider even more flexible model specifications, for example, allowing for higher orders of age rather than just the quadratic specification and interactions of the age polynomials with various other dummy variables. We present estimates with standard errors under such extended models with flexible specifications, demonstrate that qualitatively similar results continue to be obtained without substantial loss of statistical significance, and thus confirm further robustness of the main empirical findings by NunnWantchekon11.
Consider the model
where $Y_{ij,\ell}$ denotes a measure of trust, $D_{ij,\ell}$ denotes an intensity measure of slave trade, $X_{ij,\ell}$ contains polynomial basis elements of age up to degree 10, ethnic fractionalization, gender, urban residence, occupation, religion, living conditions, the interactions of the polynomial basis of age with all the dummy variables, and country fixed effects, consisting of 597 dimensions of covariates in total. Note that the total number of regressors ($p+1=598$) is much larger than the effective sample size ($\underline C = N \wedge M = 185$) of two-way clustering in this extended setting.
Table (ref) summarize the estimates of the effects of slave trade on mistrust as measured by the “trust of neighbors,” corresponding to Table 1 of NunnWantchekon11. The last two columns in the table show the original estimates obtained under the prototypical model by NunnWantchekon11 and corresponding lasso estimates obtained under more flexible model specification by our method. Across all the measures of slave exports, the original estimates and our lasso estimates are similar with similar levels of statistical significance. These results demonstrate that, even for flexible model specifications entailing high-dimensional covariates, the proposed method allows to produce qualitatively similar results without extensive loss of significance, and we can thus confirm further robustness of the main empirical findings by NunnWantchekon11.
Michalopoulos_Papaioannou2013 analyze the effects of pre-colonial institutions on contemporary regional developments in Africa, controlling for various population, locational and geographic covariates including population density, distance to capital, distance to sea coast, distance to border, water area, land area, elevation, land suitable for agriculture, ecological suitability, petrolium, and diamond mine as well as country fixed effects in their baseline model. Estimates of these effects are obtained by running regressions with a sample that pools $n_{ij}$ populated pixels $\ell \in [n_{ij}]$ in ethnic group $i$ and country $j$ across the cells $(i,j) \in [N] \times [M]$ of $N (=93)$ ethnic groups and $M (=48)$ countries. Standard errors are computed by the two-way cluster-robust method of CGM11 for the ethnic group and district as two ways of clustering.
With our proposed method that is applicable to both high-dimensional models and multi-way clustering, they could consider even more flexible model specifications, for example, allowing for interactions of all combinations of geographical covariates and locational covariates. We present estimates with standard errors under such extended models with flexible specifications, demonstrate that qualitatively similar results continue to be obtained without substantial loss of statistical significance, and thus confirm further robustness of the main empirical findings by Michalopoulos_Papaioannou2013.
Consider the model
where $Y_{ij,\ell}$ denotes a regional development measured by night light intensity, $D_{ij,\ell}$ denotes an intensity measure of pre-colonial ethinic institutions, $X_{ij,\ell}$ contains population density, interactions of all combinations of locational controls (distance to capital, distance to sea coast, and distance to border), interactions of all combinations of geographical controls (water area, land area, elevation, land suitable for agriculture, ecological suitability, petrolium, and diamond mine), and country fixed effects, consisting of 82 or 83 dimensions of covariates in total. Note that the total number of regressors ($p+1=83$ or 84) is much larger than the effective sample size ($\underline C = N \wedge M = 48$) of two-way clustering in this extended setting.
Table (ref) summarize the estimates of the effects of pre-colonial institutions on regional development as measured by the “light density,” corresponding to parts of Table 3 of Michalopoulos_Papaioannou2013. The last two columns in the table show the original estimates obtained under the prototypical model by Michalopoulos_Papaioannou2013 and corresponding lasso estimates obtained under more flexible model specification by our method. Across all the measures of pre-colonial institutions and all specifications, the original estimates and our lasso estimates are similar with similar levels of statistical significance. These results demonstrate that, even for flexible model specifications entailing high-dimensional covariates, the proposed method allows to produce qualitatively similar results without extensive loss of significance, and we can thus confirm further robustness of the main empirical findings by Michalopoulos_Papaioannou2013.
In this paper, we investigate high-dimensional regression models when data is sampled under multi-way clustering. We establish the convergence rates for the lasso and post-lasso estimators under multi-way clustering. We then propose an inference method based on a post-double-selection procedure and show that it is asymptotically valid under multi-way clustering. Simulation studies demonstrate that the proposed procedure works well in finite sample under multi-way clustering in comparison with existing alternatives. We demonstrate that our method can enrich the flexibility of regression models and robustness of empirical results through a couple of empirical applications in growth and development economics.
Indeed, both multi-way clustering and high dimensionality are two important issues which concern applied research. The existing literature provide solutions to each of multi-way clustering and high-dimensionality separately. To our best knowledge, the literature does not seem to provide a solution to both of these issues simultaneously. In this paper, we filled this void in the literature.