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.
52,545 characters · 14 sections · 25 citation commands
Non-Robustness of the Cluster-Robust Inference: with a Proposal of a New Robust Method
\onehalfspacing
{7.45mm}
Cluster-robust (CR) standard errors account for within-cluster correlations. Such correlations often arise by construction within an industry He1998 or within a state BeDuMu2004, to list a couple of the earliest examples. Today, even if a model may not induce cluster dependence by construction, applying CR methods by observable group identifiers is quite common in practice.
While CR methods have been used for decades, “it is only recently that theoretical foundations for the use of these methods in many empirically relevant situations have been developed” MaNiWe2022. The initial theory fixes the (maximum) cluster size $N_g$ to a small number and assumes some large number $G$ of clusters Wh84,LiZe86,Ar87. More recent theory DjMaNi19,HaLe19,Ha22 allows for large $N_g$, but still imposes the restriction $\sup_g N^2_g / N \rightarrow 0$ of vanishing maximum cluster size relative to the whole sample size $N = \sum_{g=1}^G N_g$ as $G \rightarrow \infty$, for establishing the convergence rate and asymptotic distribution of the estimator.
These conventional assumptions may not be always satisfied. They do not hold for certain important situations, such as the case in which a researcher uses the 51 states in the U.S. as clusters. As documented later in this article, the distribution of such cluster sizes $N_g$ appears to follow the power law. Consequently, $\sup_g N^2_g / N \rightarrow 0$ will not hold. To fix ideas, think of $\sup_g N_g^2 / N \gg \sup_g N_g / N \approx 0.1 \gg 0$ for the state of California as an outlier in terms of the cluster size in typical survey data.
We show that the assumption of bounded second moments for the CR score fails if the distribution of $N_g$ follows the power law with the tail exponent less than two. For the example of using the 51 states as clusters, unfortunately, we cannot rule out the possibility of the tail exponent being less than two. Furthermore, we show that the conventional assumption $\sup_g N_g^2/N \rightarrow 0$ fails too. Therefore, we cannot recommend using the conventional CR methods for such data. We highlight several examples from a list of original research articles recently published in a top journal.
We propose a simple fix in light of these adverse findings concerning the conventional CR methods. Specifically, we propose a weighted CR (WCR) method to suppress heavy tails in the distributions of cluster sums of scores. We provide supporting theories to guarantee that this proposed method works. Simulation studies demonstrate that our proposed WCR method achieves significantly more accurate coverage than the conventional CR methods, thus supporting our arguments that the WCR method is robust while the conventional CR methods are not.
The initial theory \citep*[][]{Wh84,LiZe86,Ar87} for CR methods assumes small cluster sizes $N_g$ as $G \rightarrow \infty$. It is implemented by the `cluster()' option or the `vce(cluster)' option by Stata, and is used by almost all, if not all, empirical papers that report CR standard errors. It has been known that a large cluster size $N_g$ in data could lead to a large CR standard error \citep*[e.g.,][p. 324]{CaMi15}.
More recently, \citet*{Ha07}, \citet*{IbMu10,IbMu16}, \citet*{BeCoHa11}, and \citet*{CaSaSh21} propose alternative theory to accommodate large $N_g$ while $G$ is fixed. Their frameworks differ from ours in a few aspects. First, they exploit within-cluster central limit theorem (CLT) which is suitable for panel data, while we focus on the cross-sectional settings which accommodate arbitrary within-cluster dependence as opposed to weak dependence. Second, they consider sequences of increasing cluster sizes $N_g$ for all $g$ which is suitable for long panels, while we treat $N_g$ as random variables drawn from a widely supported distribution to accommodate both small clusters (like the state of Wyoming) and large clusters (like the state of California) in cross-sectional settings. Third, they consider small $G$ while we consider large $G$.
As such, more closely related to this article are the recent developments by \citet*[][]{DjMaNi19}, \citet*{HaLe19} and \citet*{Ha22}. Both their frameworks and our framework consider large $G$ with arbitrary within-cluster dependence as opposed to weak dependence. The asymptotic theories in these papers essentially require $\sup_g N^2_g/N \rightarrow 0$. As mentioned in Section (ref), however, we have $\sup_g N^2_g / N \gg \sup_g N_g / N \approx 0.1 \gg 0$ for the state clusters in typical survey data in the U.S. More generally, we find that the assumption of $\sup_g N^2_g/N \rightarrow 0$ fails if the distribution of $N_g$ follows the power law, as is the case with the state clusters among others. Consequently, we propose the alternative WCR method, accommodating cases with $\sup_g N^2_g/N \not\rightarrow 0$. With this said, we want to stress that our framework with the power law does not necessarily nest those of these preceding papers.
Furthermore, it has been known that the conventional CR standard errors are biased downward, and the jackknife estimator exhibits a better performance Ha2022,MaNiWe2022fast,MaNiWe2022leverage. When the cluster size $N_g$ follows a power law, however, the self-normalized means may converge to a non-Gaussian limit in general, and hence even the jackknife is not guaranteed to work in such cases. That being said, we suggest to combine our proposed WCR method with the jackknife variance estimation. Specifically, we suggest the HC3 estimator MaWh85 with the weighted cluster sum as an effective unit of observation.
Finally, this paper is also closely related to the recent literature on cluster randomized experiments \citep*[e.g.][]{BaLiShTa2022,BuCaShTa2022,CrToVa2022} in a couple of ways, even though the main objectives in this literature are different from ours. First, while the existing literature discussed thus far treats $N_g$ as deterministic sizes, this new literature treats $N_g$ as stochastic sizes similarly to our paper. While this recent literature assumes $\sup_g N_g^2/N \stackrel{p}{\rightarrow} 0$ (and equivalently $\mathbb{E}[N_g^2]<\infty$), however, we allow $N_g$ to be drawn from heavy-tailed distributions so $\sup_g N_g^2/N \stackrel{p}{\rightarrow} 0$ need not hold. Second, our proposed WCR estimator encompasses one of the estimators proposed by \citet*{BuCaShTa2022}. This implies that their estimator in fact works under our assumption as well, and consequently, their approach does not require $\sup_g N_g^2/N \stackrel{p}{\rightarrow} 0$ under the set of our assumptions.
Besides those discussed above, there are many important papers in the extensive literature of CR methods. We refer readers to the comprehensive surveys by \citet*{CaMi15,CaMi22} and \citet*{MaNiWe2022}.
Often available for empirical research is a sample of $G$ clusters, where the $g$-th cluster consists of $N_g \in \mathbb{N}$ observations for each $g \in \{1,\cdots,G\}$. A common assumption in this setting is that observations within a cluster may be arbitrarily correlated, but they are assumed to be independent across clusters. For instance, the fifty-one states are often treated as $G=51$ clusters based on a random sample drawn from the population in the U.S., allowing state-specific factors to induce dependence among observations within each state possibly.
To fix ideas, suppose that a researcher uses a clustered sample $\{\{(Y_{gi},X_{gi}')'\}_{i=1}^{N_g}\}_{g=1}^G$ to study the linear model
where $U_g = (U_{g1},\cdots,U_{gN_g})'$ and $X_g = (X_{g1},\cdots,X_{gN_g})'$ for each $g \in \{1,\cdots,G\}$. Also write $Y_g = (Y_{g1},\cdots,Y_{gN_g})'$ for each $g \in \{1,\cdots,G\}$. Then, the ordinary least squares (OLS) estimator for the parameter vector $\theta$ under the cluster sampling takes the form of
Furthermore, commonly employed cluster-robust (CR) variance estimators take the form of
where $a_n \rightarrow 1$ is a suitable finite-sample adjustment and $ \widehat S_g = \sum_{i = 1}^{N_g} X_{gi} \widehat U_{gi} $ with $ \widehat U_{gi} = Y_{gi} - X_{gi}'\widehat\theta. $ In particular, a majority of empirical economics papers is based on the `cluster()' option or the `vce(cluster)' option in Stata, which uses (ref) and (ref) with the finite-sample adjustment
Note that this adjustment factor satisfies $a_n \rightarrow 1$ as $G \rightarrow \infty$.
Also used is the jackknife variance estimator defined by
where $\widehat\theta_{-g}$ denotes the leave-one-cluster-out estimator defined by $$ \widehat\theta_{-g} = \left(\sum_{h \neq g} X_h'X_h\right)^{-1} \left(\sum_{h \neq g} X_h'Y_h\right). $$
While they provide robustness to cluster dependence, the CR methods are vulnerable to common cross-sectional situations in which there are a small number of large clusters.
The standard econometric theory to guarantee that the OLS estimator (ref) and the CR variance estimator (ref) behave well under the cluster sampling are based on the asymptotic property
as $G \rightarrow \infty$, where $\stackrel{p}{\rightarrow}$ stands for convergence in probability, $\stackrel{d}{\rightarrow}$ stands for convergence in distribution, and
Indeed, the cluster-robust approach allows for robustness against within-cluster dependence. However, this robustness is not cost-free. Namely, suitable moments of $\Xi_g$ and $S_g$ need to be finite for the weak law of large numbers and the central limit theorem to be invoked in (ref). Specifically, the asymptotic normality (ref) and the validity of the CR variance estimator (ref) under the cluster sampling require $S_g$ to have finite second moments. If the second moment of $S_g$ does not exist, then $\sqrt{G}\left(\widehat\theta-\theta\right)$ diverges and hence would not converge in distribution as in (ref). However, we are going to argue that this bounded second moment condition for the cluster score $S_g$ can fail under the cluster sampling even if the second moment of the individual's score $X_{gi}U_{gi}$ were finite.
For ease of illustration, we first introduce a few definitions and notations. For a given $j \in \{1,\cdots,\dim\{X_{gi}\}\}$, let $\Sigma_g$ and $Z_{gi}$ be short-hand notations for the $j$-th coordinate of $S_g$ and the $j$-th coordinate of $X_{gi}U_{gi}$, respectively -- we omit the dependence on $j$ in these notations for brevity. For any distribution function $F$, we say that $F$ is regularly varying (RV) at infinity if it satisfies the following property that
for any $x>0$ and some constant $\alpha>0$. The constant $\alpha$ is referred to as the tail exponent, which measures the tail heaviness of $F$. Let $F$ denote the marginal distribution of $Z_{gi}$. For each $n \in \{2,3,\cdots\}$, let $C^n$ denote the copula that characterizes the joint distribution of $(Z_{g1},\dots,Z_{gn})$ such that $$ \mathbb{P}\left( Z_{g1}\leq z_1,\dots,Z_{gn} \leq z_n \right) = C^n(F(z_1),\dots,F(z_n)). $$ With these definitions and notations, we make the following set of assumptions.
We provide some discussions about Assumption (ref). Assumption (ref).(ref) requires that the CDF of $Z_{gi}$ implies a finite mean and its tail is regularly varying. The regularly varying tail is mild and satisfied by many commonly used heavy-tailed distributions such as Pareto, Student-t, Cauchy, F distributions. Without loss of generality, we focus on the right tail and consider that $Z_{gi}$ is non-negative. Otherwise, our conclusion still holds provided ZhShWe09
Given the regularly varying tail, we can elegantly characterize the moment conditions mikosch1999regular as
When $Z_{gi}$ has bounded infinite-order moments (such as Gaussian distribution) or even has a compact support, the corresponding tail exponent can be considered to be arbitrarily large. See, for example, de07 for a more precise statement in terms of $1/\alpha$. These distributions are commonly referred to as thin-tailed distributions, which are still covered by our conclusion. For conciseness, we suppress them in this assumption but still implement them in the simulation studies in Section (ref). Assumption (ref).(ref) allows for dependence among $Z_{g1},\dots,Z_{gN_g}$ within each cluster $g$. This condition allows for very general formats of correlations within each cluster.
Assumption (ref).(ref) requires that the distribution of the cluster size $N_g$ is also regularly varying. This condition is key to our results and is plausibly satisfied as shown in the following section. Since $N_g$ takes integer values, we may consider it as the integer part of some continuous random variable, say $N^*_g$ whose tail is regularly varying. Given that $N^*_g-1 \leq N_g\leq N^*_g$, we have that $\mathbb{P}(N_g>y)\sim \mathbb{P}(N^*_g>y)$ as $y\rightarrow\infty$ and hence continue with $N_g$ for conciseness. In some applications, the requirement of independent cluster sizes in Assumption (ref).(ref) may be too strong. \citet*{BuCaShTa2022} for instance consider non-independent cluster sizes for cluster randomized experiments. For later results, we can relax this requirement at the cost of strengthening Assumption (ref).(ref) -- see Remarks (ref) and (ref) ahead for details.
Suppose that $\beta < \alpha$ under Assumption (ref).(ref) and (ref).(ref). In other words, the tail of the distribution of $N_g$ dominates that of $Z_{gi}$. In this case, the tail heaviness of the summation $\Sigma_g=\sum_{i=1}^{N_g}Z_{gi}$ is dictated by that of $N_g$. Therefore, even if $Z_{gi}$ has a finite $r$-th moment for $r<\alpha $, the $r$-th moment of $\Sigma_g$ might still be infinite if $r>\beta $. The following theorem formalizes this argument.
We provide a proof in Appendix (ref).
Under i.i.d. sampling, the `robust' standard errors (e.g., those based on Eicker-Huber-White variance estimator or the robust option in Stata) exist and behave well if $Z_{gi}$ has a bounded second moment. Under cluster sampling, however, the cluster-`robust' standard errors (i.e., those based on (ref) or the cluster()/vce(cluster) options in Stata) may fail to exist even if $Z_{gi}$ has a bounded second moment. The first part of Theorem (ref) shows that the assumption of $\sup_g N_g^2/N \rightarrow 0$ which has been imposed by even the recent literature (cf. Section (ref)) on cluster-robust inference is also implausible under the power law.\footnote{DjMaNi19 discuss the scenario in which $\sup_g N_g^2/N \rightarrow 0$ is relaxed to $\sup_g N_g/N \rightarrow 0$ under some additional assumptions. We find that such relaxation could be violated under our framework with a random $N_g$. See Appendix (ref) for details. } The second part of Theorem (ref) further implies that $\Sigma_g$ may have an infinite second moment even if $Z_{gi}$ has a bounded second moment. This occurs when $\beta < 2 < \alpha$.
A natural question now is whether the distribution of cluster sizes $N_g$ has a heavy tail with tail index $\beta < 2$. An answer will of course depend on data and the context of empirical research. Let us focus on one of the most common situations of empirical economic analyses with clustering. Specifically, suppose that a researcher clusters a sample of individuals in the U.S. Panel Study of Income Dynamics (PSID) by the 51 states. We extract the sample of all the male individuals from the most recent wave of 2019. This sample consists of $4808$ observations across $G=51$ states.
In this data set, the five largest clusters have their sizes of $N_{4} = 629$ (California), $N_{42} = 430$ (Texas), $N_{32} = 366$ (North Carolina), $N_{21} = 310$ (Michigan), and $N_{39} = 304$ (South Carolina). Figure (ref) plots the logarithm of the rank of $N_g$ (in descending order) against the logarithm of $N_g$ for all the 51 clusters. The gray line indicates the fitted line using the largest 25 observations, i.e., 50% of the 51 states. Note that these largest 25 states follow the line of slope $-1.76$, while the remaining 26 states follow another line of a less steep slope. This shape of the log-log plot implies that the distribution of the cluster sizes $N_g$ is approximately asymmetric double Pareto. The slope of $-1.76$ can be interpreted as the negative value of the Pareto exponent $\beta$ in the right tail of the distribution $H$ of $N_g$. Hence, the first part of Theorem (ref) implies that the conventional assumption $\sup_g N_g^2/N \rightarrow 0$ is likely to fail.
From this graphical analysis, it seems quite possible that $\beta < 2$. We now resort to a more formal and widely practiced test to reinforce this heuristic observation. Figure (ref) displays the so-called Hill plot DrReDe2000, consisting of 95% confidence intervals for the Parteto exponent $\beta$ of the distribution of $N_g$ using the largest $k$ order statistics for $k \in \{2,\cdots,25\}$. Observe that we cannot rule out the possibility of $\beta < 2$ for any $k$.
Recall from Theorem (ref) that the Pareto exponent $\beta$ of the distribution $H$ of $N_g$ dictates the heaviness of the distribution of the coordinates $\Sigma_g$ of the score $S_g$. If $\beta < 2$, then the asymptotic convergence (ref) fails even if $\alpha > 2$ holds for each coordinate $Z_{gi}$ of the score $X_{gi}U_{gi}$. In light of the conclusion from the previous two paragraphs, we cannot rule out $\beta < 2$. Therefore, it follows that the widely employed CR variance estimator (ref) may lead to misleading inference. In summary, we do not recommend the common practice of using the CR standard error with the 51 states as clusters for this data set from the U.S. PSID.
We highlight several examples based on recent original research articles published in 2020 in Econometrica that use cluster-robust methods of inference for which we cannot rule out the possibility of $\beta < 2$. Namely we showcase in point with state clusters in the U.S. BuHaTiVo20, region clusters in Russia EnMaPe20, and NGO branch clusters AlBaBaBuRaSuVi20.
We would like to emphasize that the main point of this section is not to criticize a specific list of papers -- the authors of these papers are in fact great for allowing us to replicate their analyses more easily than many others. Our selection from one journal from only one recent year reflects only the tip of the iceberg -- furthermore, there are even more articles published in Econometrica and other journals for which data are unavailable to us or replication was difficult for us. We suspect that an enormous number of empirical research papers from a wide variety of journals and from a long history of the literature are in fact subject to this problem of non-robustness.
BuHaTiVo20 use 56 state clusters in the U.S.\footnote{The cluster variable is the state as in Section (ref). We selected this example, however, because the data are different and the numbers of state categories are different too.} for their analysis summarized in their Tables 1--2. We extracted their cluster variable and draw its Hill plot for the largest 28 (50%) clusters in Figure (ref).
For every number $k$ of the top order statistics except for $k=20$ and 21, the 95% confidence interval does not exclude the possibility of $\beta < 2$. Therefore, we cannot plausibly assume $\beta > 2$ for this data set to use the conventional CR methods. Furthermore, using the formal statistical test presented in Appendix (ref), we reject the hypothesis of $\mathbb{E}[S_g^2]<\infty$ for a number of regressions whose results are reported in their Tables 1--2.
EnMaPe20 use 78 region clusters in Russia for their analysis summarized in their Tables 1--3. We extract their cluster variable and draw its Hill plot for the largest 39 (50%) clusters in Figure (ref).
For every number $k$ of the top order statistics, the 95% confidence interval does not exclude the possibility of $\beta < 2$. Therefore, we cannot plausibly assume $\beta > 2$ for this data set to use the conventional CR methods. Furthermore, using the formal statistical test presented in Appendix (ref), we reject the hypothesis of $\mathbb{E}[S_g^2]<\infty$ for a number of regressions whose results are reported in their Tables 1--3.
AlBaBaBuRaSuVi20 use 108 NGO (BRAC) branch clusters for for their analysis summarized in their Table 7. We extract their cluster variable and draw its Hill plot for the largest 54 (50%) clusters in Figure (ref).
For every number $k$ of the top order statistics, the 95% confidence interval does not exclude the possibility of $\beta < 2$. Therefore, we cannot plausibly assume $\beta > 2$ for this data set to use the conventional CR methods. Furthermore, using the formal statistical test presented in Appendix (ref), we reject the hypothesis of $\mathbb{E}[S_g^2]<\infty$ for a number of regressions whose results are reported in their Table 7.
Once again, we want to stress that the above list of examples is just the tip of the iceberg. Our selection of the three papers is due only to the relative ease of replication thanks to the efforts by the authors of these papers. There were many other papers for which replication was difficult or infeasible, sometimes due to limited data availability. It is worthy to note that there are at least nine papers that cluster standard errors among those published in American Economic Review in year 2020, and at least eight papers that cluster standard errors among those published in Econometrica in year 2020, among others.\footnote{The results of the test of the null hypothesis $\mathbb{E}[S_g^2]<\infty$ for all these papers, as well as those in other top five journals, based on the method presented in Appendix (ref), will become available in spreadsheet upon request from the authors under certain conditions. We thank our graduate research assistants for their great efforts to replicate the analyses by the large number of papers.}
When the score does not have a finite variance, one might still want to resort to the self-normalized central limit theorem (CLT) to establish the asymptotic normality under a possibly slower convergence rate of convergence than the standard $\sqrt{G}$ rate. Under our framework in which the distribution of $\{N_g\}$ has a heavy tail, however, such a self-normalized CLT is not guaranteed to work. It may or may not hold depending on the structure of within-cluster dependence.
For simplicity of illustration, suppose that one is interested in obtaining the standard error of the estimator $\hat\theta = N^{-1} \sum_{g=1}^G \sum_{i=1}^{N_g} Y_{gi}$ for the mean $\theta = \mathbb{E}[Y_{gi}]$. In this setting, one may hope to establish the self-normalized asymptotic normality $\mathbb{E}[(\hat\theta - \theta)^2]^{-1/2} (\hat\theta - \theta) \stackrel{d}{\rightarrow} \mathcal{N}(0,1)$, or more practically, $\widehat{\mathbb{E}}[(\hat\theta - \theta)^2]^{-1/2} (\hat\theta - \theta) \stackrel{d}{\rightarrow} \mathcal{N}(0,1)$. Such a self-normalized CLT may fail to hold, however, if the distribution of $N_g$ has a heavy tail and within-cluster correlation is strong.
For instance, consider the data generating process $ Y_{gi} = \rho_G R_g + e_{gi}, $ where $R_g$ is a cluster-specific effect which induces within-cluster dependence, and $e_{gi}$ is a random noise which is i.i.d. across $i$ and $g$. In this case, the limit distribution of the self-normalized mean $\mathbb{E}[(\hat\theta - \theta)^2]^{-1/2} (\hat\theta - \theta)$ in general fails to be asymtotically Gaussian if $\rho_G \gg G^{1/2-1/\beta} \rightarrow 0$ as $G\rightarrow\infty$ when $\beta<2$. See Appendix (ref) for details about this argument. Hence, the rate-adaptive CR standard errors and even the jackknife standard errors may fail to work in our framework of cluster sampling involving cluster-specific effects and heavy-tailed distributions of cluster sizes.
Recall that cluster-specific variables were the initial motivations for applied researchers to start using cluster-robust standard errors in cross-section data He1998. Cluster-specific variables are ubiquitous in economic data set. This is also true by design for cluster randomized experiments.
This section proposes a new approach to cluster-robust inference. It does not suffer from the aforementioned problem with the standard CR methods characterized by Theorem (ref). Our proposal is simple, but we provide a complete theoretical rationale for why it works unlike the standard CR approach.
One simple way to fix the problem of the non-robustness of the standard CR methods is to modify (ref) and (ref) by
respectively, where $a_n \overset{p}\rightarrow 1$.
Likewise, we may also modify the jackknife estimator (ref) by
where $$ \widehat\theta_{-g}^{\text{WCR}} = \left(\sum_{h \neq g} N_h^{-1}X_h'X_h\right)^{-1} \left(\sum_{h \neq g} N_h^{-1}X_h'Y_h\right) $$ denotes the weighted leave-one-cluster-out estimator.
We are not the first to suggest a weighted estimator. In the different context of cluster randomized experiments, \citet*{BuCaShTa2022} suggest a weighted difference-in-mean estimator similarly to ours. While they propose such a weighting for the purpose of identifying and estimating certain parameter of interest, our estimator $\widehat\theta^{\text{WCR}}$ in fact encompasses their estimator $\hat\theta_{1,G}$ as a special case.\footnote{They propose two kinds of estimators. While their $\hat\theta_{1,G}$ is analogous to our estimator $\widehat\theta^{\text{WCR}}$, their $\hat\theta_{2,G}$ is analogous to the conventional CR estimator $\widehat\theta$.} This implies that their estimator works under our assumptions as well. See Remark (ref) ahead.
As in Section (ref), for given $j \in \{1,\cdots,\dim\{X_{gi}\}\}$, let $\Sigma_g$ and $Z_{gi}$ be short-hand notations for the $j$-the coordinate of $S_g$ and the $j$-th coordinate of $X_{gi}U_{gi}$, respectively. Define $\widetilde \Sigma_g := \Sigma_g/N_g = \sum_{i=1}^{N_g} Z_{gi}$. In contrast to Theorem (ref), the following theorem establishes that $\widetilde \Sigma_g$ has the same tail index as the original score $Z_{gi}$.
We present a proof in Appendix (ref).
Recall that the non-robust methods (e.g., those based on Eicker-Huber-White standard errors) requires only bounded second moments of the score $X_{gi} U_{gi}$. Also recall from Section (ref) that the conventional CR methods can fail when the second moment of $N_g$ is infinite even if the second moments of the score $X_{gi} U_{gi}$ were finite. Now, Theorem (ref) states that our proposed robust approach works as far as the second moments of the score $X_{gi} U_{gi}$ is bounded regardless of whether the second moment of $N_g$ if finite or not. Accordingly, we allow $\sup_g N_g^2/N$ to diverge.
To derive the asymptotic properties of our WCR estimator, we make the following additional assumption.
Part 1 of this assumption requires the i.i.d. sampling across clusters, as is assumed in general under cluster sampling environments. It does, however, allow for arbitrary dependence within each cluster as in the general cluster sampling environments. Part 2 of the above assumption rules out multi-collinearity. In summary, Assumption (ref) is a standard requirement.
Recall that our main claim from Section (ref) is that the asymptotic normality (ref) for the conventional CR inference may not work if $\beta < 2$ even if $\alpha > 2$. Section (ref) demonstrates that we cannot rule out such a pathetic case with $\beta < 2$ even for the most common setting of cluster sampling, namely using 51 states as clusters. Furthermore, Section (ref) argues that even the self-normalized CLT generally fails for the conventional CR inference. On the other hand, the following theorem shows that our proposed CR approach based on $\widehat\theta^{\text{WCR}}$ and $\widehat{V}_{\hat\theta}^{\text{WCR}}$ works as far as $\alpha > 2$ is true, regardless of whether $\beta < 2$ is true or not.
A proof is presented in Appendix (ref). A few remarks are in order.
As a consequence of Theorem (ref), one can use $$ \widehat\theta^{\text{WCR}}_j \pm 1.96 \sqrt{\widehat V_{\hat\theta, jj}^{\text{WCR}}} \qquad\text{ or }\qquad \widehat\theta^{\text{WCR}}_j \pm 1.96 \sqrt{\widehat V_{\hat\theta, jj}^{\text{WCR,JACK}}} $$ as a robust 95% confidence interval for the $j$-th coordinate of $\theta$ even if $\beta < 2$ as is likely the case with 51 states used as clusters.
To justify the use of the jackknife WCR variance estimator $\widehat V_{\hat\theta, kk}^{\text{WCR,JACK}}$, we close this section with the following proposition.
A proof is presented in Appendix (ref).
In this section, we present simulation studies to evaluate the finite sample performance of our proposed WCR methods in comparison with the conventional CR methods. In light of the discussions in Section (ref), we focus on data that contain cluster-specific variables.
The data generating design is defined as follows. We consider the cluster treatment model with individual covariates $$ Y_{gi} = \theta_0 + \theta_1 T_g + \sum_{j=1}^K \theta_j X_{g,j+1} + U_{gi} $$ following \citet*[][Equation (40)]{MaNiWe2022fast} among others. The binary treatment variable $T_g$ takes the value of one for $\lceil 0.2G \rceil$ clusters and zero for the remaining clusters $G - \lceil 0.2G \rceil$, where $\lceil a \rceil$ denotes the smallest integer greater than or equal to $a$. We draw cluster sizes $ N_g \sim \lceil 10 \cdot \text{Pareto}(1,\beta) \rceil $ independently for $g \in \{1,\cdots,G\}$. For each $g \in \{1,\cdots,G\}$, we independently draw $N_g$-variate random vectors, $ (\tilde X_{g1j},\cdots, \tilde X_{gN_gj})' \sim \mathcal{N}(0,\Omega) $ for $j \in \{1,\cdots,K\}$ and $ (\tilde U_{g1},\cdots,\tilde U_{gN_g})' \sim \mathcal{N}(0,\Omega), $ where $\Omega$ is an $N_g \times N_g$ variance-covariance matrix such that $\Omega_{ii}=1$ for all $i \in \{1,\cdots,N_g\}$ and $\Omega_{ij}=1/2$ whenever $i \neq j$. The controls are constructed by $ X_{gik} = 0.2 F_{\text{Beta}(2,2)}^{-1} \circ \Phi (\tilde X_{gik}), $ where $F_{\text{Beta}(2,2)}$ and $\Phi$ denote the CDFs of the $\text{Beta}(2,2)$ and standard normal distributions, respectively. The errors are heteroskedastically constructed by $ U_{gi} = 0.2 \tilde U_{gi} $ if $T_g=0$ and $ U_{gi} = \tilde U_{gi} $ if $T_g=1$.
We vary values of the exponent parameter $\beta \in \{1,2,4\}$ across sets of simulations. The regression coefficients are fixed to $(\theta_0,\theta_1,\theta_2,\cdots,\theta_{K+1})'=(1,1,1,\cdots,1)'$ throughout, whereas the dimension $K$ of covariates vary as $K \in \{0,1,5\}$. We set the sample size (i.e., the number of clusters) to $G = 50$ across sets of simulations, which is close to the number of states in the U.S. we discussed as an example earlier. Each set of simulations consists of 10,000 Monte Carlo iterations.
Figures (ref)--(ref) draw Q-Q plots under $\beta=2$ and $\beta=1$, respectively, of the self-normalized statistics:
For these figures, we focus on the case with $K=0$. In each figure, the dashed line indicates the $45^{\circ}$ line, and the solid line indicates the fitted line.
Observe that the self-normalized statistics based on the conventional CR methods suffer from farther deviation away from the theoretical quantiles, whereas those based on our proposed WCR methods more precisely follow the theoretical quantiles. This observation is true for both the analytic standard error estimator and the jackknife estimator. The deviations for the conventional CR methods further exacerbate as the distribution of $N_g$ becomes heavier, as in the transition from $\beta=2$ (in Figure (ref)) to $\beta=1$ (in Figure (ref)). These results are consistent with the general non-Gaussianity of the conventional CR methods as discussed in Section (ref), as well as the guaranteed Gaussianity of our proposed WCR methods as discussed in Section (ref).
Table (ref) summarizes simulation results. Displayed are the mean square error (MSE), the rejection frequencies based on the analytic standard error estimation in round brackets, and the rejection frequencies based on the jackknife standard error estimation in square brackets. The nominal probability of rejection is set to $p=0.05$ throughout. The first three columns show the results for the conventional CR methods. The last three columns show the results for our proposed WCR methods.
We find the following three observations in these simulation results. First, focus on the MSE. While the MSE of the OLS estimator $\widehat\theta_1$ exponentially blows up as $\beta$ decreases, the MSE of the WCR estimator $\widehat\theta_1^{\text{WCR}}$ remains stable as $\beta$ varies. Second, consider the rejection frequencies reported in the round brackets based on the analytic standard error estimators, $\widehat V_{\hat\theta,11}^{\text{CR}}$ and $\widehat V_{\hat\theta,11}^{\text{WCR}}$. While the rejection frequencies for the conventional CR method based on $\widehat V_{\hat\theta,11}^{\text{CR}}$ blows up as $\beta$ decreases, the rejection frequencies for our proposed WCR method based on $\widehat V_{\hat\theta,11}^{\text{WCR}}$ remain stable and closer to the nominal rejection probability of $p=0.050$ as $\beta$ varies. Third, consider the rejection frequencies reported in the square brackets based on the jackknife standard error estimators, $\widehat V_{\hat\theta,11}^{\text{CR,JACK}}$ and $\widehat V_{\hat\theta,11}^{\text{WCR,JACK}}$. While these jackknife standard error estimators deliver more desirable rejection frequencies than the analytic standard error estimators for each of the conventional CR method and the new WCR method, we continue to observe the same qualitative pattern as in the case of the analytic estimators. Namely, while the rejection frequencies for the conventional CR method based on $\widehat V_{\hat\theta,11}^{\text{CR,JACK}}$ blows up as $\beta$ decreases, the rejection frequencies for our proposed WCR method based on $\widehat V_{\hat\theta,11}^{\text{WCR,JACK}}$ remain stable and closer to the nominal rejection probability of $p=0.050$ as $\beta$ varies.
From these observations, it seems more desired to use the new WCR methods over the conventional CR methods for estimation accuracy as well as the property of robust inference.
CR standard errors are used extensively in the practice of empirical economic analyses. Despite its wide practice, the literature has been silent about a potential pitfall of the conventional CR methods. The theory for the conventional CR methods assumes that (i) $N_g$ is fixed while $G \rightarrow \infty$, or at best (ii) $\sup_g N_g^2/N \rightarrow 0$ as $G \rightarrow \infty$. Even the most common empirical settings, such as the case of using the 51 states in the U.S. as clusters, violate these existing assumptions.
We establish three theoretical results in this paper. First, if the distribution of $N_g$ follows the power low with exponent less than two, then the second moment of the CR score fails to exist even if the second moment of the individual score were finite. Second, the second moment of our proposed WCR score exists as far as the second moment of the individual score exists regardless of the exponent of the distribution of $N_g$. Third, consequently, our proposed WCR methods enjoy robust inference unlike the conventional CR methods. It is also worthy to note that, in the special case where the tail exponent $\beta$ is large, the WCR methods behave similarly to the conventional CR methods.
One may want to resort to self-normalized CLT to validate the conventional CR standard errors and the jackknife estimators. However, the limit distribution of self-normalized sums is not guaranteed to be Gaussian under the power law of cluster sizes. Hence, we cannot recommend using the conventional CR methods if one is interested in robustness in inference under a small number of large clusters, such as the case of using the 51 states in the U.S. as clusters. A main disadvantage of our proposed methods is that the weighting alters the estimator from our familiar OLS estimator, which is a cost that a researcher needs to pay to enjoy the robustness.
{7.5mm} {7.9mm}