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.
99,051 characters · 16 sections · 61 citation commands
Clustering with Potential Multidimensionality: Inference and Practice
In our survey of articles published in American Economic Review in the years 2021 and 2022, 70% of 133 articles containing empirical specifications reported some cluster-robust standard errors. Among these papers, 12 reported two-way clustered standard errors. Despite the common use of cluster-robust inference in empirical work, little guidance has been provided on the motivation for clustered standard errors and the level of clustering for linear regressions and nonlinear estimators. Even less is known about the formal reasoning underlying two-way clustered standard errors. In this paper, we establish asymptotic properties of M-estimators under finite populations with potentially multiway cluster dependence, allowing for unbalanced and unbounded cluster sizes in the limit.
We use the finite population framework for the choice of appropriate inference because we can combine sampling-based uncertainty that arises from possibly not observing the entire population with design-based uncertainty caused by the stochastic assignment of treatment or policy variables. Following abadie2023should, we distinguish between two situations that justify computing clustered standard errors: \romannumeral 1) cluster sampling induced by random sampling of groups of units, and \romannumeral 2) cluster assignment caused by correlated assignment of “treatment” within the same group. While abadie2023should justify cluster-robust standard errors for the difference-in-means estimator with one-way clustering, we generalize their setup both in considering general M-estimators and in allowing for multiway clustering.
We show that for conducting inference with general M-estimators, one-way clustering is only necessary when there is either cluster sampling or cluster assignment, or both on nested or same dimensions. Multiway clustering can be justified when clustered assignment and clustered sampling occur on different dimensions, or when either sampling or assignment is multiway clustered. The same results are also shown for functions of M-estimators with the estimator of the average partial effect (APE) as a leading example. In the special case of linear regression on a binary treatment variable, one-way clustered standard errors on the assignment dimension is sufficient under homogeneous treatment effects even if the sampling and assignment dimensions are non-nested.
Until recently, accounting for the large sample behavior of design-based settings with multiway clustering has been a difficult problem. Asymptotic theory for variables that have multi-dimensional dependence has thus far relied on separate exchangeability (e.g., davezies2018asymptotic). Separate exchangeability implies that the marginal distributions of clusters are exchangeable mackinnon2021wild. However, by construction, separate exchangeability is violated within a design-based framework. Consider a binary assignment variable $X_i$. Since the error term $u_i = X_i u_i(1) + (1 - X_i) u_i(0)$ depends on the nonstochastic potential error $u_i(x)$, even if treatments $X_i$ are identically distributed across clusters, the marginal distribution of $u_i$ differs because $X_i$ is weighted differently. Hence, a limit theory that accommodates heterogeneity of clusters and observations is required. By building on the central limit theorem in yap2023general, we derive results on large-sample behavior of standard estimators in this environment.
When estimating the variance-covariance matrix, the usual variance estimators are typically overly conservative for the finite population variance-covariance matrix in one-way clustering. However, we find that the two-way clustered variance estimator as proposed in cameron2011robust, henceforth CGM, can be anticonservative. In response to the anticonservativeness of the CGM estimator in design-based settings, there are two approaches that empirical researchers may take. The first approach is to make an assumption on how the individual treatment effects are correlated within the same cluster: CGM is conservative when the correlation is positive, which is reasonable in most applications. The second approach is to remain agnostic and to use CGM2, a more conservative version of the CGM variance estimator proposed by davezies2018asymptotic. Since simulations show CGM2 is often unnecessarily conservative, we propose a simple shrinkage variance estimator relying on adjustments using covariates. The probability limit of our adjusted variance estimator is guaranteed to be no smaller than the finite population variance matrix and provides a smaller upper bound than CGM2.
We discuss several practical settings involving clustering and justify the validity of cluster-robust inference within a design-based framework. These examples include a standard difference-in-means estimator, fixed effects regressions, a triple difference estimator, and a linear regression on two assignment variables clustered on different dimensions. We also show the performance of our proposed shrinkage variance estimators in simulations and an empirical illustration. In particular, our adjusted standard errors can be substantially smaller than CGM2 while still maintaining correct coverage.
Within the growing literature on design-based inference (e.g., abadie2020sampling, xu2020potential, athey2022design, abadie2023should, de2023level), there are at least four nuances unique to our environment that we find from our theory and from analyzing applications of cluster-robust inference. First, there are special cases where two-way clustering reduces to one-way clustering in a way that one-way clustering does not reduce to a heteroskedasticity-robust variance-covariance matrix. Second, CGM can be anticonservative while standard one-way robust variances are conservative. Third, beyond the special case of abadie2023should, estimands from a fixed effects regression cannot be interpreted as the average treatment effect in general. A broader lesson here is that cluster dependence not only affects variance estimation in inference; it can also affect the interpretation of estimands. Fourth, with multiple assignment variables on different clustering dimensions, we find that in certain cases it suffices to use one-way cluster variances on the respective dimensions. In contrast, this setting is not allowed or cannot be discussed in one-way clustering.
Our paper contributes to at least three strands of literature that have gained recent attention. First, we contribute to the literature on design-based inference (cited and summarized above). Second, we contribute to the literature on multiway clustering (e.g., davezies2018asymptotic, mackinnon2019cluster, menzel2021bootstrap, chiang2023using, yap2023general, chiang2024standard). To the best of our knowledge, we are the first to consider any design-based environment with multiway clustering, which is a difficult problem as separate exchangeability used in most multiway clustering limit theorems does not hold in the general design-based setting. Third, we contribute to the literature on causal panel data (e.g., de2020two, callaway2021difference, sun2021estimating, wooldridge2021two, athey2022design, borusyak2024revisiting, dechaisemartin2024difference, gardner2024two,arkhangelsky2024causal). By using a design-based setup, we make an assumption on assignment instead of the potential outcome as most of these papers do with parallel trends. Considering multiway clustered assignment is new relative to the design-based setting of athey2022design, and our setting provides new insights into the interpretation of the fixed effect estimands in these designs.
Consider a sequence of finite populations indexed by population size $M$, where $M$ diverges to infinity in deriving the asymptotic properties. Suppose there are $G$ mutually exclusive clusters in population $M$ defined as either the primary sampling units in the sampling scheme or the partition in the assignment design, where each cluster has $M^G_g$ units, $g=1,2,\dots,G$. Further, suppose the population can also be partitioned into $H$ mutually exclusive clusters according to either the sampling scheme or assignment design on dimensions possibly different from that of $G$ clusters. Each cluster $H$ contains $M^H_h$ units, $h=1,2,\dots,H$. We use $\mathcal{N}^C_c$ to denote the set of observations in the $c$th cluster on the $C \in \{ G, H, G\cap H \}$ dimension, and $c(i)$ to denote the cluster that $i$ belongs to on the relevant dimension. If there is only one way of partitioning, $H$ and $G$ clusters coincide with each other.
For each unit $i$ within cluster $g$ and cluster $h$, we observe $(X_{iM}, z_{iM}, Y_{iM})$. The vector $X_{iM}$ is the vector of stochastic assignment variables, $z_{iM}$ is a set of non-stochastic attributes, and $Y_{iM}$ is the realized outcome. The categorization of assignments and attributes depends on the empirical question. Typically, the key variables of interest in an empirical study could be viewed as assignment variables, and the rest covariates as attribute variables. With the potential outcome framework, there exists a mapping, denoted by the potential outcome function $y_{iM}(x)$, from the assignment variables to the potential outcomes. For example, $y_{iM}(x)=x\theta_{01}+z_{iM}\theta_{02}+e_{iM}$ for continuous outcomes, and $y_{iM}(x)=\mathbbm{1}[x\theta_{01}+z_{iM}\theta_{02}+e_{iM}>0]$ for binary outcomes, where $z_{iM}$ and $e_{iM}$ are observed and unobserved attributes respectively.\footnote{We emphasize $x$ as the argument of the potential outcome function because it is the only stochastic variables in the function.} The potential outcome function $y_{iM}(\cdot)$ is non-stochastic.\footnote{This implies that the unobserved attributes are non-stochastic.} Nevertheless, the realized potential outcome, $Y_{iM}=y_{iM}(X_{iM})$, is random. Hence, the finite population setting can be understood as a setting that conditions on the potential outcomes and attributes of the $M$ units in the population. We use $W_{iM}:=(X_{iM}, Y_{iM})$ to denote the random vector for brevity.
We study solutions to a population minimization problem, where the estimand of interest is a $k\times 1$ vector denoted by $\theta^*_M$:
The expectation $\mathbb{E}$ in ((ref)) is taken over the distribution of $X$ since $X$ is the source of randomness here. The function $q_{iM}(\cdot, \cdot)$ is the objective function for a single unit. Examples include (nonlinear) least squares, weighted least squares, maximum likelihood, and instrumental variables estimation in the just identified case (based on the first-order condition of the minimization problem). The subscripts of the objective function imply its dependence on the non-stochastic attribute variables $z_{iM}$, so covariates are allowed in the model. Interpreting this estimand is context-dependent, so we abstract from this discussion in our general framework.
Let $R_{iM}$ denote the binary sampling indicator, which is equal to one if unit $i$ is sampled. Hence, the sample size is $N=\sum^G_{g=1}\sum^{H}_{h=1}\sum_{i\in \mathcal{N}^{G\cap H}_{(g,h)}}R_{iM}=\sum^M_{i=1}R_{iM}$. The sample size is random unless the sample is population. The estimator of $\theta^*_M$ is denoted by $\hat{\theta}_N$, which solves the minimization problem in the sample:
Our random variables are two-way clustered in that random variables for indices $i$ and $j$ are independent if $g(i)\ne g(j)$ and $h(i) \ne h(j)$, formalized in the following assumptions.
Various sampling schemes are allowed under Assumption (ref). When all sampling probabilities are equal to one, we observe the entire population; when $\rho_{gM}=\rho_{hM}=1$ but $\rho_{uM}<1$, we have independent sampling; if the dimensions of $H$ and $G$ coincide or $H$ is nested within $G$ without loss of generality (e.g., zip code areas nested within counties), $\rho_{gM}<1$ or $\rho_{hM}<1$ implies one-way cluster sampling; lastly, if the dimensions of $H$ and $G$ differ, $\rho_{hM}<1$ and $\rho_{gM}<1$ implies two-way cluster sampling. For instance, one can first sample according to occupations and industries and then sample individuals from the chosen intersections of occupations and industries.
Assumption (ref) allows assignments to be correlated within clusters on the dimensions of $G,H$, or both. Cluster assignment is another source of within-cluster correlation in addition to cluster sampling. The assignment variables $X_{iM}$ are not necessarily identically distributed, which allows the assignments to depend on the fixed attributes $z_{iM}$. Independent sampling and assignment processes imply Assumption (ref).
Cluster sizes in our theorems are allowed to be unbalanced and unbounded in the limit. Nevertheless, in order to apply asymptotic theory, we have assumptions that restrict cluster heterogeneity and the growth rate of the cluster sizes relative to the population size and variances.
This assumption implies $G, H \to \infty$, and rules out the case where a particular subset of clusters dominates the population. In the results that follow, for matrices $A$ and $B$, when we say $A \geq B$, we mean that $A- B$ is positive semidefinite.
To fix ideas, we start with the simpler case of one-way clustering. Namely, there is only one way of partitioning so that $G$ and $H$ clusters coincide. Without loss of generality, let $\rho_{hM}=1$ in this case. Let $m_{iM}(W_{iM},\theta)$ denote the score function of $q_{iM}(W_{iM}, \theta)$. The variance matrix of M-estimators is defined as
where
and
It can be shown that:
where
and
account for heteroskedasticity (Eicker-Huber-White (EHW)) and within-cluster correlation respectively. The terms
and
are the finite population counterparts of $\Delta_{ehw,M}(\theta^*_M)$ and $\Delta_{cluster,M}(\theta^*_M)$.
The conventional superpopulation variance matrix is denoted by
Notice that the middle part of the sandwich form of $V_M$ is different from that of $V_{SM}$ due to two “extra" (E) terms $\Delta_{E,M}$ and $\Delta_{EC,M}$ scaled by the composite sampling probability.
The usual cluster-robust variance estimator (CRVE) that uses the estimator from liang1986longitudinal for $V_{\Delta M}$ is given by
where
and
For one-way clustering, we use a stronger version of Assumptions (ref) to enhance interpretability. This assumption ensures convergence at rate $N^{-1/2}$ and follows hansen2019asymptotic. The results are extended to an arbitrary convergence rate and two-way clustering in the next subsection.
Theorem (ref) shows asymptotic normality with the finite population cluster-robust asymptotic variance (CRAV). In the variance-covariance matrices, the term $\Delta_{cluster,M}(\theta^*_M)$ is scaled by the sampling probability $\rho_{uM}$ because of the two-stage sampling scheme. Nevertheless, the usual CRVE, $\hat{V}_{SN}$, converges to $V_{SM}$, in which the estimation of $\rho_{uM}$ has been accounted for.
The term related to clustering in the variance formula, $\Delta_{cluster,M}(\theta^*_M) - \rho_{gM} \Delta_{EC,M}$, is zero if we have both independent sampling and independent assignment. Otherwise, these components in the variance must be accounted for. Hence, Remark (ref) suggests that we should adjust standard errors of M-estimators for clustering at the level of cluster sampling or cluster assignment. It generalizes the results in abadie2023should: they prove the case for the difference-in-means estimator, while the remark above holds for all M-estimators with either continuous or discrete assignment variables.
Since the sum of the two additional terms, $\Delta_{E,M}+\Delta_{EC,M}$, is positive semidefinite, we reach the conclusion in Remark (ref). Remark (ref) together with Theorem (ref)(2) imply that the usual CRVE is often too conservative. There are exceptions where using the usual CRVE for inference is approximately correct, with the leading scenario summarized in the remark below.
Another special case for the usual CRVE to be correct for inference is when $\Delta_{E,M}+\Delta_{EC,M}=\textbf{0}$, which is true if either $\mathbb{E}\big[m_{iM}(W_{iM},\theta^*_M)\big]=\textbf{0},\ \forall\ i=1,\dots,M_g^G,\ g=1,\dots,G$ or $\sum_{i \in \mathcal{N}^G_g}\mathbb{E}\big[m_{iM}(W_{iM},\theta^*_M)\big]=\textbf{0},\ \forall\ g=1,\dots,G$. The former is true for the variance of the coefficient estimator on the assignment variables under the sufficient conditions provided by abadie2020sampling, including constant treatment effects and other linearity conditions. The latter holds if the finite population is composed of repetitions of the smallest cluster. With this kind of data structure, $\theta^*_M$ that solves $\mathbb{E}\Big[\sum^G_{g=1}\sum_{i \in \mathcal{N}^G_g}m_{iM}(W_{iM},\theta^*_M)\Big]= \textbf{0}$ is also the solution to $\mathbb{E}\Big[\sum_{i \in \mathcal{N}^G_g}m_{iM}(W_{iM},\theta^*_M)\Big]=\textbf{0}$ for each cluster $g$. However, these kinds of special cases rarely hold in practice.
Sometimes, we are interested in the functions of M-estimators rather than M-estimators themselves. Let $f_{iM}(W_{iM},\theta^*_M)$ be a $q \times 1$ function of $W_{iM}$ and $\theta^*_M$. Suppose we wish to estimate $\gamma^*_M=\frac{1}{M}\sum^M_{i=1}\mathbb{E}\big[f_{iM}(W_{iM},\theta^*_M)\big]$. As an example, $\gamma^*_M$ could be the APE from nonlinear models, where $f(\cdot, \cdot)$ is some partial derivative for continuous variables or some difference function for discrete variables.
Let $\hat{\gamma}_N=\frac{1}{N}\sum^M_{i=1}R_{iM}f_{iM}(W_{iM},\hat{\theta}_N)$ be the estimator of $\gamma^*_M$. Denote the finite population variance matrix by
The superpopulation variance matrix is then $V_{f,SM}=\Delta^f_{ehw,M}+\rho_{uM}\Delta^f_{cluster,M}$. And the usual CRVE is denoted by $\hat{V}_{f,SN}=\hat{\Delta}^f_{ehw,N}+\hat{\Delta}^f_{cluster,N}$. The detailed definition of each term can be found in Appendix A.
Theorem (ref) shows that the conservative property of the usual CRVE of M-estimators also applies to the usual CRVE of any functions of M-estimators.
Now, suppose the $H$ and $G$ clusters are partitioned on different dimensions. The variance matrix of M-estimators is defined as
where
By grouping the cross products of score functions into different cases: individual units ($\Delta_{ehw,M}(\theta^*_M)$), units belonging to the intersection of $G$ and $H$ ($\Delta_{(G\cap H),M}(\theta^*_M)$), units belonging to $G$ but in different $H$'s ($\Delta_{G,M}(\theta^*_M)$), units belonging to $H$ but in different $G$'s ($\Delta_{H,M}(\theta^*_M)$), it can be shown that:
where
and
As before, $\Delta_{E(G\cap H),M}$, $\Delta_{EG,M}$, and $\Delta_{EH,M}$ are the finite population counterparts of $\Delta_{(G\cap H),M}(\theta^*_M)$, $\Delta_{G,M}(\theta^*_M)$, and $\Delta_{H,M}(\theta^*_M)$, respectively.
Assumption (ref) ensures that the overall variance is not driven by a few large clusters. A stronger way of stating Assumptions (ref) and (ref) is that $\frac{1}{M} \sum\limits^G_{g=1}\left(M^G_g\right)^2\leq C<\infty$ and $\max\limits_{g\leq G}\frac{\left(M^G_g\right)^2}{M}\to 0$, as $M\to \infty$ with an analogous condition in the $H$ dimension where $\lambda_M \geq \underline{c} M$ for some $\underline{c} >0$. This assumption is more similar to the setting of hansen2019asymptotic, but rules out two-way balanced clusters where there is one unit in every intersection: if there are $G$ clusters on both the $G$ and $H$ dimensions, then $M=G^2$ so $\frac{1}{M} \sum\limits^G_{g=1}\left(M^G_g\right)^2= G^3/ G^2 = G \rightarrow \infty$. Assumption (ref) as stated makes no such restriction as $\frac{1}{M^2} \sum\limits^G_{g=1}\left(M^G_g\right)^2 = G^3/ G^4 =1/G \rightarrow 0$. With more flexible cluster sizes, the convergence rate depends on the variance of the sum, which motivates Assumption (ref). The stronger version of the assumption only allows a convergence rate of $M^{-1/2}$, which is not necessarily true in the weaker version. For instance, the weaker version can allow a slower convergence rate of $(\sum_{g=1}^G (M^G_g)^2/M^2)^{-1/2} = G^{-1/2}$ instead of $M^{-1/2}=G^{-1}$. Since we can allow for different convergence rates, we use the scale $\lambda_M$ to restrict cluster heterogeneity in Assumption (ref) to derive the asymptotic distribution.
The proof of this theorem is largely analogous to Theorem (ref), just that we apply the central limit theorem (CLT) from yap2023general instead of hansen2019asymptotic. chiang2023using (Table 1) pointed out that the two-way cluster-robust standard errors are usually valid. A notable exception is when the additive components are degenerate, such that the random variable can be written as $D_{it}=\alpha_i \gamma_t$, where $\alpha_i,\gamma_t$ are cluster-specific random variables on the respective dimensions. As noted in Remark 1 of yap2023general, this data generating process (earlier pointed out by menzel2021bootstrap) is ruled out by our summability condition in Assumption (ref).
Remark (ref) coupled with Remark (ref) imply that two-way clustering is only justified if there is (\romannumeral 1) two-way clustered sampling (i.e., $\rho_{gM}<1$ and $\rho_{hM}<1$); (\romannumeral 2) two-way clustered assignments (i.e., $\Delta_{G,M}(\theta^*_M)\neq \Delta_{EG,M}$ and $\Delta_{H,M}(\theta^*_M)\neq \Delta_{EH,M}$); or (\romannumeral 3) clustered sampling and clustered assignments on different dimensions. For the third case, there could be many combinations of sampling schemes and assignment designs, possibly combining two-way sampling and two-way assignments at the same time.
In a special case, $G$ and $H$ clusters could be partitioned at different but nested levels. Without loss of generality, suppose $H$ clusters are nested in $G$ clusters. The variance-covariance matrix can be simplified to be
where (\romannumeral 1) $\rho_{1M}=\rho_{uM}\rho_{hM}$ and $\rho_{2M}=\rho_{uM}\rho_{gM}\rho_{hM}$ for nested sampling; (\romannumeral 2) $\rho_{1M}=\rho_{2M}=\rho_{uM}$ for nested assignment; (\romannumeral 3) $\rho_{1M}=\rho_{uM}$ and $\rho_{2M}=\rho_{uM}\rho_{gM}$ for sampling at the $G$ level and assignment at the $H$ level; (\romannumeral 4) $\rho_{1M}=\rho_{2M}=\rho_{uM}\rho_{hM}$ for assignment at the $G$ level and sampling at the $H$ level. As a consequence, one-way clustering at the higher level $G$ is sufficient.
While the usual variance estimators are generally conservative for the finite population variance-covariance matrix in one-way clustering, this is not necessarily true with multiway clustering. We denote the usual two-way cluster-robust variance estimator initially proposed in cameron2011robust by
where
for $C\in \{G,H,G\cap H\}$. Define the superpopulation two-way CRAV to be
For variance estimators to converge to the variance-covariance matrices in the general environment, we impose an additional assumption.
The condition in Assumption (ref) is required in the following propositions so that the asymptotic error incurred by using the matrix estimator $\hat{V}$ relative to the true matrix $V$ converges to zero. The difference between Assumption (ref) and the existing assumptions is that the previous assumptions defined $\lambda_M$ as the variance of the sum, which includes all two-way clustered terms, but here, $\lambda^C_M$ only includes terms from one of the two dimensions. Since the strategy for showing such convergence is similar to yap2023general, an analogous summability condition and a condition on the largest cluster having a negligible contribution to the variance are required. Since Assumption (ref) only accounts for one-way clustering in the denominator, Assumption (ref) is stronger than Assumption (ref).\footnote{Assumption (ref) in its present form requires that the cluster correlation on each clustering dimension be of comparable scale. However, this can be moderated in situations where one clustering dimension dominates the other dimension in the estimation of variance matrix.}
The anticonservativeness results from the subtraction of the correlation terms within intersection clusters $\hat{\Delta}_{{G\cap H},N}(\hat{\theta}_N)$ so that the difference in the meat of the variance sandwich is $\rho_{uM} \rho_{gM} \rho_{hM} (\Delta_{E,M} + \Delta_{E(G\cap H),M} + \Delta_{EG,M} + \Delta_{EH,M})$, which can be positive or negative in general. In Example B.1 in Appendix B, we give an example where $\hat{V}_{CGM}$ is anticonservative. Nevertheless, if all within-cluster correlation of $E[m_{iM}(W_{iM},\theta^*_M)]$ is positive, then $\hat{V}_{CGM}$ is still a conservative variance estimator. davezies2018asymptotic propose an alternative variance estimator that does not adjust for double counting the intersection clusters. Let
and
Hence, in contrast to CGM, CGM2 is asymptotically conservative.
Even though we restore conservativeness of the usual variance estimator by using $\hat{V}_{CGM2}$, it can be too conservative. The terms in the usual CRAV can be estimated in the standard way. Taking one-way clustering as an example, it is more challenging to estimate the two extra terms, $\Delta_{E,M}$ and $\Delta_{EC,M}$, because $\mathbb{E}\big[m_{iM}(W_{iM},\theta^*_M)\big]$ is generally non-identifiable due to the missing data problem of the potential outcome framework. For instance, with a binary assignment variable, $\mathbb{E}\big[m_{iM}(W_{iM},\theta^*_M)\big]=P(X_{iM}=1)\cdot m_{iM}\big((1,y_{iM}(1)),\theta^*_M\big)+P(X_{iM}=0)\cdot m_{iM}\big((0,y_{iM}(0)),\theta^*_M\big)$, and we do not observe both $y_{iM}(0)$ and $y_{iM}(1)$ at the same time. This observation motivates a simple method to estimate a bound on these extra terms such that the corrected variance estimators are still conservative, but are smaller than the one-way CRVE or CGM2.
The variance estimator depends on the sampling and assignment schemes. If the cluster variance matrix is purely induced by sampling, then $\Delta_{EC,M}=\Delta_{cluster,M}$, which can be consistently estimated, as $\hat{\Delta}_{cluster,N}$ consistently estimates $\rho_{uM} \Delta_{cluster,M}$. Then, it remains to identify a lower bound of $\Delta_{E,M}$. On the other hand, if there is cluster assignment, we have to identify a lower bound of $\Delta_{E,M}+\Delta_{EC,M}$ jointly. Similarly, if there is cluster sampling on the $G$ dimension but cluster assignment on the $H$ dimension, we only need to adjust for the $H$ dimension. If there are multiway cluster assignments, we need to adjust for both dimensions.
We first discuss estimating the $\Delta_{E,M}$ component. We can remove part of $\Delta_{E,M}$ using the regression-based approach below, by eliminating the variation that is linearly predictable from the covariates $z_{iM}$. Consider the estimator,
where $\hat{K}_N=\Big(\sum\limits^M_{i=1}R_{iM}z_{iM}'z_{iM}\Big)^{-1}\bigg[\sum\limits^M_{i=1}R_{iM}z_{iM}'m_{iM}(W_{iM},\hat{\theta}_N)'\bigg]$. With clustered data, we can include cluster dummies as regressors in the linear projection of $m_{iM}(W_{iM},\hat{\theta}_N)$ onto the fixed attributes.
Under one-way clustered sampling but independent assignment, the estimator for an upper bound of the finite population CRAV is
where $G_N$ is the number of clusters in the sample. The composite sampling probability $\rho_{uM}\rho_{gM}$ can be estimated by $N/M$, where the population size $M$ is assumed to be known. If the entire population is observed, $\rho_{uM}\rho_{gM}$ is simply one. We ignore the Hessian matrix as it does not affect the discussion here. This estimator is asymptotically conservative because $\Delta^Z_M \leq \Delta_{E,M}$.
Next, we turn to $\Delta_{E,M}+\Delta_{EC,M}$. We could sum $m_{iM}(W_{iM},\hat{\theta}_N)$ within each cluster, and linearly project $\sum_{i \in \mathcal{N}^G_g}R_{iM} m_{iM}(W_{iM},\hat{\theta}_N)$ onto the fixed attributes. The number of observations in the linear projection is the number of clusters in the sample. To reduce the dimensionality of the regressors, the fixed attributes can also be summed within clusters as one way of aggregation. As a result, $\sum_{i \in \mathcal{N}^G_g}\mathbb{E}\big[m_{iM}(W_{iM},\theta^*_M)\big]$ can be partially estimated by its predicted value from the linear projection. Let
and
Estimate $\Delta_{E,M}+\Delta_{EC,M}$ with
Theorem (ref) proposes an easy way to partially remove $\Delta_{E,M}+\Delta_{EC,M}$ all at once with one-way clustering, and $\hat{\Delta}^Z_{CE,N}$ is positive semidefinite. In this case, one-way cluster sampling is allowed as long as there is no within-cluster sampling. With large samples, even though the limit of the adjusted finite population CRVE is still conservative (as $\Delta^Z_{CE,M}\leq\big(\Delta_{E,M}+\Delta_{EC,M}\big)$), it is still less conservative than the limit of the usual CRVE.
We list all possible cases of sampling and assignment and their corresponding adjusted variance estimators in Table 1 for one-way clustering. Details of the derivation are relegated to Appendix B. Case 4 reduces exactly to the approach in abadie2020sampling for linear regression. While they take the square of the difference between the score and the predicted score as their estimator, their approach is numerically equivalent to taking the difference of the second moments as we propose.\footnote{Our approach is fundamentally different from the proposal in abadie2023should. They split the data into subsamples and directly estimate the finite population variance, while our approach here shrinks the variance using information from covariates. All cluster sizes need to diverge in their approach, which is not required here. We also allow for constant assignment within clusters, which is ruled out by abadie2023should.}
A similar shrinkage procedure can be applied to two-way clustering with CGM2 since the additively separable one-way cluster objects can be shown to converge to their limit even with multiway dependence. The variance estimator after adjustment is still conservative for the finite population two-way CRAV. To be precise, let:
and
When adjusting the variance matrix estimator with two-way clustering, we use the entire population. Theorem (ref) also provides formal guarantees that the cluster estimators for the intersection of $G$ and $H$ either converge to their estimand or are negligible based on the scale of the cluster correlation on different clustering dimensions.
Table 2 summarizes our proposed variance estimators for multiway sampling or multiway assignment (Hessian matrix is ignored here). $H_N$ and $G_N$ are the number of clusters in the sample on the dimensions of $H$ and $G$ respectively. We expect cases 2 and 3 in Table 1 and case 2 in Table 2 to be the leading cases in empirical practice.
In the context of doing inference for functions of M-estimators, we can also apply the same techniques to estimate the two extra terms, $\Delta^f_{E,M}$ and $\Delta^f_{EC,M}$. The only difference is that the dependent variables in the regression-based approach would be the cluster sum of
rather than $m_{iM}(W_{iM},\hat{\theta}_N)$ alone. (See the details of the notation in Appendix A.) The results for multiway clustering can be derived in a similar fashion and hence are omitted here.
In this section, we compare the Monte Carlo standard deviation of the coefficient estimator and the APE estimator of the assignment variable in a binary response model with a set of different standard errors. We are mainly interested in the finite sample performance of the proposed shrinkage variance estimators. As a leading case in empirical practice, we focus on (multiway) clustered assignment with the entire population observed.
In the population generating process, there is a single assignment variable $X_{iM}\in\{0,1\}$ and a single attribute variable $z_{iM}=z_{g(i)M}+z_{h(i)M}$, where $z_{hM}=\pm 1$ with equal probability and $z_{gM}=\pm 1$ with equal probability in the design of two-way clustered assignment, and $z_{gM}=\pm 2$ with equal probability in the design of one-way clustered assignment. The potential outcome of a binary response is generated as
The idiosyncratic unobservable $e_{iM}$ is the residual from regressing random realization of a standard normal distribution on $z_{iM}$. The data of $z_{iM}$ and $e_{iM}$ are generated once and kept fixed in the population $M$.
We partition the population units into 50 clusters each on the two dimensions $G$ and $H$ with one unit for every $(g,h)$ cluster pair. As a result, the population size is 2,500. The results for 100 clusters on both dimensions are similar and hence omitted to save space.
The assignment variable $X_{iM}=A_{g(i)}B_{h(i)}$, where $A_g$ and $B_h$ are binary cluster assignment variables drawn independently with $P(A_g=1)=P(B_h=1)=1/2$, $\forall\ g=1,2,\dots, G$ and $h=1,2,\dots,H$. Therefore, the assignments are clustered at both the $G$ and $H$ dimensions. For the case of one-way cluster assignment, we fix $B_h=1$. There are 10,000 replications for both designs. For each replication, $X_{iM}$ is re-assigned according to the assignment rules above.
Estimates from the pooled probit regression of $Y_{iM}$ on 1, $X_{iM}$, and $z_{iM}$ are displayed in Table (ref) below. Columns (1) and (2) collect results for one-way cluster assignment and columns (3) and (4) show results for two-way cluster assignment.
The first two rows of Table (ref) report the Monte Carlo standard deviation of the point estimates and the coverage rate of the 95% confidence interval based on the oracle standard error, i.e., Monte Carlo standard deviation. The oracle coverage rates are very close to the nominal level of 95%. Thus, normal approximation seems to work well in finite samples. The next two rows report the superpopulation EHW standard errors and the corresponding coverage rate of the 95% confidence interval.\footnote{The $97.5^{th}$ percentile of $t(G-1)$ is used as the critical value in constructing the confidence intervals.} The EHW standard errors are too small and the confidence interval undercovers as expected.
For one-way clustered assignment, we focus on one-way cluster-robust standard errors at the level $G$. Here, we demonstrate how our shrinkage variance estimators work in finite samples. For two-way clustered assignment, we report results on both one-way cluster-robust and two-way cluster-robust standard errors.
The adjusted one-way clustered standard errors at the level $G$ are more than half smaller than the superpopulation one-way clustered standard errors, though still above the Monte Carlo standard deviation under one-way clustered assignment. Switching to two-way clustered assignment, the adjusted one-way clustered standard errors are both too small. The two-way clustered standard errors for both CGM and CGM2 estimators work well in this population generating process. There is slight downward bias of the adjusted CGM standard errors. The adjusted CGM2 standard errors are larger but are guaranteed to be conservative.
We conclude from the simulation results that the usual superpopulation one-way cluster-robust standard errors and the two-way CGM2 cluster-robust standard errors are overly conservative. When there are fixed attributes available, they can be used to estimate an upper bound of the finite population CRAV. Although the adjusted finite population cluster-robust standard error is still conservative, it often improves over the usual cluster-robust standard error. In particular, our adjusted standard errors are 40% smaller than CGM while still maintaining correct coverage.
In this section, we first apply our general theorems to the difference-in-means estimator, which has been studied thoroughly in the literature. Next, we investigate the treatment effect estimand that a cluster fixed effect regression targets. Lastly, we discuss some empirical settings where two-way clustered standard errors have been considered.
To make a direct comparison with abadie2023should, we consider a difference-in-means estimator for a binary assignment variable without covariates with multiway clustering. For treatment variable $X \in \{ 0,1 \}$, we denote the nonstochastic potential outcome as $y_{iM} (x)$. We are interested in the population average treatment effect (ATE):
Multiway assignment is treated in the following way. Data is generated by independently drawing $A_{g}\in[0,1]$, $B_{h}\in [0,1]$ and $e_{i}\sim U[0,1]$, with $X_{iM}=1\left\{ e_{i}<A_{g(i)}B_{h(i)}\right\}$. The random variables $A_{g}$ and $B_{h}$ have means $\mu_{A}, \mu_{B} > 0$ and variances $\sigma_{A}^2, \sigma_{B}^2$ respectively. This process nests several cases. If assignment is one-way clustered, then we can simply set $A_{g}=1$. Another case is where assignment occurs at the intersection level, and we need both dimensions $G$ and $H$ to be assigned treatment for the unit to be treated. Then, $A_{g}, B_{h} \in \{0,1 \}$. If $\mu_{A}$ and $\mu_{B}$ are both non-zero, then even though $X_{iM} = A_{g(i)} B_{h(i)}$ is an interaction model, we do not need to be concerned about non-normality.\footnote{ Let $A_g = \mu_A + \epsilon^A_g$ and $B_h = \mu_B + \epsilon^B_h$ where $\mathbb{E}[\epsilon_g] = \mathbb{E}[\epsilon_h]=0$. Then, $\sum_{i=1}^M X_{iM} = \sum_{i=1}^M (\mu_A + \epsilon^A_{g(i)}) (\mu_B + \epsilon^B_{h(i)}) = \sum_{i=1}^M (\mu_A \mu_B+ \epsilon^A_{g(i)} \mu_B + \mu_A \epsilon^B_{h(i)}+ \epsilon^A_{g(i)}\epsilon^B_{h(i)})$ so the dominant stochastic terms are $\epsilon^A_{g(i)} \mu_B $ and $ \mu_A \epsilon^B_{h(i)}$ instead of $ \epsilon^A_{g(i)}\epsilon^B_{h(i)}$. }
With $\alpha_M := (1/M)\sum_{i=1}^{M} y_{iM}(0)$ and error $U_{iM} := Y_{iM} - \alpha_M - \tau_M X_{iM}$, the potential errors are denoted:
For $R_{iM}=1$, we observe $\{ Y_{iM}, X_{iM} \}$, with $Y_{iM}= X_{iM} y_{iM}(1) + (1-X_{iM})y_{iM}(0)$. Let $b_{1}=\mathbb{E}\left[R_{iM} X_{iM}\right]$, $b_{0}=\mathbb{E}\left[R_{iM} (1-X_{iM})\right]$. $b_{1}$ is the probability that an individual is observed and treated; $b_{0}$ is the probability that an individual is observed and untreated. $N_{1} := \sum_{i=1}^{M} R_{iM} X_{iM}$ and $N_{0} := \sum_{i=1}^{M} R_{iM} (1-X_{iM})$. The least squares estimator is:
To state the main result, we first define a few terms. Let $\mathcal{N}_i$ denote the neighborhood of $i$, which is the set of observations that are plausibly correlated with $i$. The score is
In this context, $\mathbb{E}[\eta_{iM}]=u_{iM}(1) - u_{iM}(0) =: \tau_{iM}-\tau_M$ is not zero in general, but $\sum_i \mathbb{E}[\eta_{iM}]=0$. Hence, we define $\xi_{iM}$ as the demeaned residual for individual $i$ that features in the variance of $\hat{\tau}_M$:
This result is a corollary of the existing normality results and law of large numbers. Comparing this context to our framework, $\eta_{iM}= m_{iM}$, $\xi_{iM} = m_{iM} - \mathbbm{E}[m_{iM}]$, and $v_M = V_{TWM}$. Using our framework, we can answer questions on whether multiway clustering matters and whether multiway clustering is appropriate. Using the one-way CRVE on dimension $G$ (without loss of generality) yields the following estimand:
Then, answering the question on whether two-way clustering matters involves comparing $V_{GM}$ with $V_{TWSM}$ and answering the question on whether two-way clustering is appropriate involves comparing $V_{GM}$ with $V_{TWM}$. The comparisons yield:
Consider the environment with constant treatment effects. Then, $u_{iM}(1) - u_{iM}(0) = 0$ and hence $\mathbbm{E}[m_{iM}] =0$ and all the $\Delta_E$ terms in the equation above are 0, so the two comparisons become identical. If there is multi-way clustered assignment only or if there is clustered assignment on $H$ and clustered sampling on $G$, then ((ref)) cannot be simplified further. However, if there is multiway clustered sampling only or if there is clustered assignment on $G$ and clustered sampling on $H$, then $\Delta_{H,M}(\theta^*_M) =0$ in ((ref)) so $V_{TWSM} = V_{TWM} = V_{GM}$. This result implies that, under constant treatment effects, it suffices to cluster on the assignment dimension and the usual CRVE is no longer conservative.\footnote{This result complements Corollary 1 of abadie2017should: they find that one-way clustering is unnecessary if there is both constant treatment effects and no clustering in the assignment. }
With clustered data, empirical papers often include fixed effects (FE) in a regression. A first-order question is: What are these FE estimators approaching to in the limit within our finite population framework? In particular, when can FE estimands be interpreted as the ATE? We focus on an environment where assignment $X_{iM}$ is binary and two-way clustered with $A_g, B_h \in \{ 0,1 \}$ as described in Section (ref), and we observe the entire population, i.e., $R_{iM}=1$ for all $i$.\footnote{Due to how fixed effects estimators transform variables using within-cluster means, any form of sampling leads to random sample sizes within each cluster which introduces additional technical complications that is beyond the immediate scope of our theory, so it is left for future research.} Nonetheless, since we allow for two-way clustered assignment and study estimands for both one-way and two-way fixed effects, the results in this section are new relative to abadie2023should and athey2022design: abadie2023should study one-way fixed effects with one-way clustering, and athey2022design require a random adoption date in the two-way fixed effects environment with one-way clustered assignment.
We start with one-way fixed effect. Let the ATE be defined in the following way.
In a one-way FE estimator, let
Since $\mu_B \in [0, 1]$, all weights $\omega_i$ are positive, so the estimand is a weighted average of treatment effects. In the special case where we have perfectly balanced clusters $M_{(g(i),h(i))}^{G\cap H}=k$ and hence $M_{g(i)}^{G}=Hk$, then $\tau_{OWFE} = \tau_M$. This result suggests that, unlike the difference-in-means example where clustered dependence only affects the variance, clustered dependence here also affects the interpretation of the estimand. Mechanically, the difference is that the estimator here contains products $X_{iM} X_{jM}$ for $i \ne j$ that is not present in Section (ref), so correlation in $X$ affects the estimand.
Moving on to two-way fixed effect (TWFE), we focus on the case with multiway assignment where a unit is treated if and only if both its clusters are treated. This setup is different from one-way assignment at the intersection level (see (ref)). With a linear model
the TWFE estimator is
where $\tilde{X}_{iM}$ is the residual when regressing $X$ on the fixed effects. Due to baltagi2008econometric Chapter 3,
This result is fairly weak in that it requires balanced clusters, but in this idealized situation, the TWFE estimand is the ATE. Unlike Proposition (ref), we cannot interpret the estimand as a weighted average of treatment effects in general, but this result is not surprising considering the large difference-in-differences (DD) literature on how TWFE cannot be interpreted as a weighted average of treatment effects, even with parallel trends.
One scenario in which two-way clustering is often employed is with panel data, where one clusters at the cross-sectional entity level and time level to adjust for both cross-sectional dependence and serial correlation. Clustering at the entity level can be justified by cluster sampling if cross-sectional entities are independently sampled from a finite population and each entity comes with its serial observations. If the entire cross-sectional population is observed, clustering at the entity level can be justified by serially correlated assignments for each entity across time. A leading example is the standard difference-in-differences, where either before or after the initial treatment period, the treatment assignments within entities are perfectly positively correlated. Meanwhile, treatments across the initial treatment period are negatively correlated. Justification for clustering at the time level is more ambiguous. One leading argument is to account for spatial correlation and hence two-way clustering has been used as a substitute for spatiotemporal-correlation robust inference.
Let us ignore the Hessian matrix for now, as it does not matter for the comparison of different variance matrix estimators. For spatial data, we often work with the entire population and hence sampling may not be a particularly important aspect here. Let us suppose that all sampling probabilities are one. Since our asymptotic theory requires that the number of clusters at both dimensions goes to infinity in the sequence of populations, we work with long panel data where $(|D_M|,T)\to\infty$. Notice that $D_M$ denotes lattices in $\mathbb{R}^d$ and $|D_M|$ stands for the cross-sectional population size. Adapting the argument from xu2022design, the finite population spatiotemporal variance matrix is approximately
where $m_{it}=m_{itM}(W_{itM},\theta^*_M)$ to simplify notation. $\nu(i,j)$ measures the cross-sectional distance and $\nu(t,s)$ measures the time lag. $b_M$ and $b_T$ are the cross-sectional and time dimension bandwidth respectively. $\omega(\cdot)$ is the kernel weighting, which satisfies: $\omega (0)=1$; $\omega\Big(\frac{\nu(\cdot,\cdot)}{b}\Big)=0$ for any $\nu(\cdot,\cdot)>b$; $\Big|\omega\Big(\frac{\nu(\cdot,\cdot)}{b}\Big)\Big|\leq 1$, $\forall$ $\nu(\cdot,\cdot)$.
If the dependence within clusters is weak, the usual CRVE is still consistent for one-way clustering within the superpopulation framework. Although without downweighting, the CRVE will have slower rates of convergence than HAC type variance estimators (see hansen2007asymptotic). On the other hand, if there is spatiotemporal correlation in assignments, two-way clustered variance matrix accounts for the first three terms in $V_{PM}$ but ignores cross-cluster correlation in different time periods for different units. The magnitude of two-way clustering and spatiotemporal robust inference cannot be ordered without additional information on the correlation pattern of assignment variables.
There are instances where assignment on multiple variables could be clustered at differing dimensions. Suppose there are two assignment variables $X_{1iM}$ and $X_{2iM}$. The assignment of $X_{1iM}$ is clustered on the dimension of $G$, whereas the assignment of $X_{2iM}$ is clustered on the dimension of $H$. As an empirical example, hersch1998compensating studies wage and injury risk trade-off for men and women. The two assignment variables are injury rates for individual $i$'s industry and for $i$'s occupation respectively. Hence, one assignment variable is clustered at the industry level and the other assignment variable clustered at the occupation level.
When two assignment variables are independent of each other, the Frisch–Waugh theorem implies that we only need to cluster the standard errors on one dimension for each of the assignment variables in a linear regression. We conduct a simple simulation to compare two-way clustered standard errors with one-way clustered standard errors on each dimension for the coefficient estimator on both assignment variables.
The potential outcome function is given below. \[ y_{iM}(x_{g(i)},x_{h(i)})=\tau_{1i} x_{g(i)}+\tau_{2i} x_{h(i)}+e_{iM}, \] where $e_{iM}$ is the nonstochastic individual unobservable. The cluster assignment variables $X_g$ and $X_h$ are binary with probabilities $P(X_g=1)=P(X_h=1)=1/2$. In the first design, we impose heterogeneous treatment effect on cross dimensions. Namely, $\tau_{1i}=\tau_{h(i)}=\pm 1$ with equal probability; $\tau_{2i}=\tau_{g(i)}=\pm 1$ with equal probability. In the second design, we set $\tau_{1i}=\tau_{g(i)}$ and $\tau_{2i}=\tau_{h(i)}$. We regress $Y_{iM}$ on 1, $X_{g(i)}$, and $X_{h(i)}$ and report different superpopulation standard errors as upper bounds for the finite population ones.
As we can see from Table (ref), indeed one-way clustered standard errors are sufficient for the coefficient estimators on the corresponding assignment variables. Namely, we can report one-way clustered standard errors at the level $G$ for $\hat{\tau}_1$ and one-way clustered standard errors at the level $H$ for $\hat{\tau}_2$. In the second design, using the two-way clustered standard errors is harmless, although both one-way and two-way clustered standard errors are conservative. However, when the heterogeneous treatment effects are more asymmetric in the first design, the two-way clustered standard errors can be a lot larger than the one-way clustered standard errors, making the inference unnecessarily overly conservative.
Recently, triple differences method has got more attention among empirical researchers. There has been ample discussion on the correct inference on difference-in-differences; see, e.g., bertrand2004much. Nevertheless, the dicussion on the inference method for triple differences has been limited. In the simulation of olden2022triple, they report one-way cluster-robust standard errors on the treatment level to directly compare with the simulation results in bertrand2004much. Recently, strezhnev2023decomposing advocates for two-way cluster-robust standard errors for triple difference estimators. Both one-way and two-way clustered standard errors have been reported in empirical research according to our survey.
Borrowing the notation from strezhnev2023decomposing, ((ref)) is a common specification for triple differences.
Unit $i$ in stratum $h$ and group $g$ is treated at time $t$ if $D_{ight}=1$. The terms $\alpha_{gh}$, $\gamma_{ht}$, and $\delta_{gt}$ are the group-stratum, stratum-time, and group-time fixed effects respectively. $\epsilon_{ight}$ is the individual idiosyncratic term. Suppose there are multiple strata $h\in\{1,2,\dots,H\}$ and multiple groups $g\in\{1,2,\dots,G\}$. Otherwise, there is not much we can do in clustering the standard errors.
We can rewrite the treatment variable $D_{ight}$ as an interaction of three variables, $D_{ight}=D_{h(i)} \times D_{g(i)} \times Post$, where $Post$ is a time dummy variable that takes the value of one for periods post the initial treatment period. Consequently, it is tempting to compare one-way clustering with two-way clustering. In our view, other than the time variable, the key boils down to the nature of the other two variables within the triple interaction term. If both grouping indicators are stochastic assignment variables on non-nested dimensions, one would like to report the two-way clustered standard errors. An example is marchingiglio2019employment, where $h$ and $g$ represent state and industry respectively. Their assignment variable is the state level adoption of gender-specific minimum wage laws applying to specific industries that employed a larger share of women. Alternatively, if treatment is assigned based on one grouping indicator, and the other grouping indicator is some nonstochastic attribute, one-way clustered standard errors are the most appropriate. As an example, in bau2021can $h$ indicates province and $g$ represents ethnicity. Since the policy they study is the roll-out of pension plans at the province level, ethnicity is nonstochastic within our design-based framework. In their context, the pension plan policy affects the ethnicity groups differentially by nature.
To showcase the differences between one-way and two-way clustered standard errors, we conduct a simple simulation of a triple differences regression based on ((ref)). There are two time periods. The treatment is essentially randomly assigned in the second time period. Both $D_g$ and $D_h$ are cluster binary variables with probabilities $P(D_g=1)=P(D_h=1)=1/2$. In the first design, $D_g$ and $D_h$ are stochastic, whereas in the second design, $D_g$ is the nonstochastic attribute variable but $D_h$ remains as the stochastic assignment variable. We construct $\tau_i$ in a way that the parallel trends assumption for triple differences holds.\footnote{Specifically, we construct $\tilde{\tau}_i=\tau_{g(i)}+\tau_{h(i)}$, where $\tau_g=\pm 2$ with equal probability and $\tau_h=\pm 1/2$ with equal probability. In design 1, $\tau_i$ is the demeaned $\tilde{\tau}_i$. In design 2, we average $\tilde{\tau}_i$ across units with $D_{g(i)}=1$ and denote this average by $\bar{\tilde{\tau}}$. $\tau_i=\tilde{\tau}_i-\bar{\tilde{\tau}}$.} We report the adjusted finite population standard errors for $\hat{\tau}$ and the coverage rate of the 95% confidence interval based on these standard errors.\footnote{We use $\tau_i$ as the attributes in the estimation of the adjusted finite population standard errors.}
As shown in Table (ref), the EHW standard errors underestimate as expected. When both grouping indicators are stochastic assignment variables, one-way clustered standard errors are generally not sufficient. By fluke the one-way clustered standard errors could be larger than the standard deviation, but that is because they are conservative within the design-based framework. The adjusted CGM2 standard error works pretty well with coverage rate of the confidence interval close to its nominal level. Switching to the case when only the grouping indicator $D_h$ is stochastically assigned, clustering the standard errors at the level of $H$ suffices. Two-way clustered standard error is overly conservative and can be more than three times larger than the one-way clustered standard errors clustered on $H$.
The adjusted finite population CRVE proposed in Theorem (ref) is applied to antecol2018equal, who study the effects of tenure clock stopping policies on tenure rates among assistant professors. The unique dataset collected by them contains all assistant professor hires at the top-50 Economics departments from 1980-2005 as pooled cross sections, resulting in 1,392 observations in total. Furthermore, the tenure clock stopping policies are assigned at the university level while the data are collected at the individual level, implying that we have a setting of observing the entire population with cluster assignment.\footnote{This group of assistant professors is treated as the population.} The standard errors in antecol2018equal are clustered at the policy university level, which is the correct level to cluster the standard errors as implied by Remark (ref). As a result, there are 49 clusters in total with cluster sizes ranging from 11 to 57.
Since the dependent variable is a binary response, we analyze the linear probability model (LPM) given in antecol2018equal along with an additional probit model given in ((ref)) below, which adopts the same notation from their paper.
The dependent variable $Y$ is an indicator of obtaining tenure at the university of initial placement. Binary variables $GN$ and $FO$ are indicators of gender-neutral and female-only tenure clock stopping policies respectively. The dummy variable $F$ is the indicator for females. The variable $E$ is an indicator of starting jobs in years zero through three after policy adoption. The vector $X$ contains individual characteristics and the vector $Z$ includes university level controls.\footnote{We refer to antecol2018equal for the details of the variables included as controls.} The parameter $\rho$ captures gender-specific time trend and $\psi$ represents gender-specific university heterogeneity. The subscripts, $u$, $g$, $i$, $t$, are indicators for university, gender, individual, and the year the job started, respectively.
antecol2018equal include gender-specific university dummies to capture different unobserved university heterogeneity for males and females. Adding group dummies in the linear model is equivalent to performing fixed effects with clustered data. However, adding group dummies in a nonlinear model may cause the incidental parameter problem. Since the cluster sizes are unbalanced, we use pooled probit with correlated random effects as suggested by wooldridge2010econometric to allow correlation between the gender-specific university heterogeneity and the covariates. Using Chamberlain-Mundlak device, the cluster size, the gender-specific university averages of individual and university characteristics and policies, and their interactions with cluster sizes are included as additional controls.
Given the probit model above is a nonlinear “difference-in-differences” model, the common trend assumption is imposed on the latent outcome variable following puhani2012treatment and wooldridge2023simple. The treatment effects are defined as the differences in the probit probabilities when the treatment variables equal one or zero. We report the average of the treatment effect for those actually treated by the specific policy. Since our emphasis in this study is on inference, we adhere to the specification presented in antecol2018equal to facilitate a direct comparison with their reported standard errors. Assume that $\psi$ conditional on the sufficient statistics (the additional controls included) follows a normal distribution, APEs can be obtained via pooled probit.
In Table (ref), panel A presents the total effects for men and women hired in years zero through three after policy adoption, and panel B shows the effects for those employed in years four or later. The left panel (columns (1)-(3)) summarizes the results under the LPM. Columns (1) and (2) report the total effects and the standard errors, as shown in column (1) of Table 2 in antecol2018equal, while column (3) reports the adjusted finite population clustered standard errors. The coefficients (APEs) are interpreted as the policy effect on the tenure attainment of the assistant professors compared with those of the same genders at the same university but without any clock stopping policies.
To estimate the adjusted finite population CRAV, we sum all the estimated score functions and control variables within clusters and apply the variance estimator in Case 2 of Table 1 together with the usual estimator of the Hessian matrix. Since the number of control variables exceeds the number of clusters in the data, we only include university characteristics as the fixed attributes in the linear projection, resulting in a linear regression with 49 observations and eight independent variables. Compared with the usual cluster-robust standard errors, the finite population cluster-robust standard errors shrink by about 4% to 21% across the eight treatment groups. In terms of the statistical significance, the effect of gender-neutral policy for women hired three or more years after the policy adoption is significant at the 5% rather than the 10% level based on the adjusted finite population cluster-robust standard error. The same result holds when the critical values from $t(48)$ distribution are used.
In the right panel (columns (4)-(6)), we can see that the APEs from the probit regression are close in magnitudes to those from the linear model. The adjusted finite population CRAV is estimated applying Theorem (ref) and the delta method. Using the same set of university characteristics as the attribute variables, the reduction from the usual clustered standard errors to the finite population clustered standard errors ranges from 4% to 25%. Based on the critical values from $t(48)$, the effect of gender-neutral policy for women hired in later years is significant at the 5% level rather than the 10% level when the finite population clustered standard error is adopted.
To sum up, control variables can help shrink the standard errors when the population is treated as finite in both linear and nonlinear models. The empirical evidence suggests that gender-neutral tenure clock stopping policy is beneficial to men in obtaining tenured positions but detrimental to women.
This paper develops finite population inference methods for M-estimators with data that is potentially clustered multiway. The takeaway for empirical practice is summarized as follows. One should only adjust standard errors for clustering if there is cluster sampling or cluster assignment. Two-way clustered standard errors are justified if there are two-way cluster sampling or two-way cluster assignments, or cluster sampling and cluster assignment on different dimensions. While the standard one-way CRVE from liang1986longitudinal is conservative for the true variance under one-way clustering, the standard two-way variance estimator from cameron2011robust is no longer conservative. Although a subsequent proposal from davezies2018asymptotic is guaranteed to be conservative for two-way clustering, their variance estimator is often too large, so we provide a refinement. Our proposed refinement uses control variables, such as baseline characteristics, and ensures that the estimators remain valid for inference. Evidence from our simulation and empirical illustration suggests that gains from our variance correction can be substantial.
Through a survey of when clustered standard errors are used in empirical work, we offer insights on the appropriateness of clustering in various contexts from our theory on M-estimation with clustering. The results apply straightforwardly to a difference-in-means estimator. In the context with one-way fixed effects, we find that the estimand is interpretable as an ATE only in special cases, but the estimand is still a weighted average of treatment effects in general. With two-way fixed effects, the requirements for the estimand to be interpretable as an ATE is even more restrictive. With spatiotemporal correlation, the magnitude of two-way clustering and spatiotemporal variance estimators cannot be ordered in general. With two assignment variables clustered on different dimensions, we find that it suffices to apply one-way clustering on the respective dimensions rather than to use two-way clustering in certain cases. In the estimation of triple differences, we find that the choice between one-way and two-way clustering depends on the nature of the variables in the triple interaction term.
The current paper focuses on the asymptotics as the number of clusters tends to infinity in the limit. For a small number of clusters or wildly unbalanced clusters, the wild cluster bootstrap\footnote{See, for example, cameron2008bootstrap and mackinnon2017wild.} has been proposed as a better-performing inference method for linear models in the setting of superpopulations. The finite population inference method for few heterogeneous clusters remains an interesting future research topic.