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.
126,180 characters · 15 sections · 3 citation commands
Leverage, Influence, and the Jackknife in Clustered Regression Models: Reliable Inference Using summclust
\vskip -12pt
\onehalfspacing
It is now standard in many fields of economics and other disciplines to employ cluster-robust inference for the parameters of linear regression models. In the most common case, each of the $N$ observations is assigned to one of $G$ disjoint clusters, which might correspond to, for example, families, schools, villages, hospitals, firms, industries, years, cities, counties, or states. The assignment of observations to clusters is assumed to be known, and observations in different clusters are assumed to be independent, but any pattern of heteroskedasticity and/or dependence is allowed within each cluster. Under these assumptions, a cluster-robust variance matrix, or CRVE, yields asymptotically valid $t$-tests, Wald tests, and confidence intervals. However, even when $N$ is very large, the resulting inferences may be unreliable when $G$ is not large or the clusters are not sufficiently homogeneous.
The literature on cluster-robust inference has grown rapidly in recent years. \citet*{CM_2015} is a classic survey article. \citet*{CGH_2018} surveys a broader class of methods for dependent data. \citet*{MNW-guide} is a comprehensive guide to empirical practice. As it discusses, there are two situations in which cluster-robust $t$-tests and Wald tests are at risk of over-rejecting to an extreme extent, even when $G$ is not small. The first is when one or a few clusters are much larger than the rest, and the second is when the only “treated” observations belong to just a few clusters; \citet*{DMN_2019} discusses the first case, and \citet*{MW-JAE,MW-TPM,MW-EJ} discuss the second. In both of these cases, one cluster (or a few of them) has high leverage, in the sense that omitting this cluster has the potential to change the OLS estimates substantially. When that actually happens, a cluster is said to be influential.
The concepts of leverage and influence are normally applied at the observation level \citep*{BKW_1980}, but they are equally applicable at the cluster level. Just as high-leverage observations can make heteroskedasticity-robust inference unreliable Chesher_1989, so too can high-leverage clusters make cluster-robust inference unreliable. Just as highly influential observations may lead us to suspect that there is something wrong with the model or the data, so too may highly influential clusters. Any situation in which a few clusters have high leverage or high influence should be worrying.
There are at least two different concepts of leverage. The usual one focuses on fitted values or, equivalently, residuals. A cluster is said to have high leverage if removing it has the potential to change the fitted values for that cluster by a lot. The second concept is partial leverage \citep*{CW_1980}. A cluster is said to have high partial leverage for the \th{j} coefficient if removing that cluster has the potential to change the estimate of the \th{j} coefficient by a lot. We discuss both concepts in (ref).
Whether a cluster has high leverage, high partial leverage, or is influential can depend on the sample in rather complicated ways. We provide a new Stata package called summclust that implements computationally-efficient ways to identify high-leverage and influential clusters and provides a number of statistics that collectively summarize the cluster structure of the dataset. These can be useful for detecting cases in which cluster-robust inference may not be reliable. Our leverage and influence calculations also allow us to compute two cluster jackknife variance matrix estimators, which we refer to as CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$, at little additional cost. These estimators are already available in Stata by using either the vce(jackknife) option or the jackknife prefix. Recent work \citep*{Hansen-jack,MNW-bootknife} suggests that CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ generally perform better in finite samples than more widely-used CRVEs; see also (ref).
The remainder of the paper is organized as follows. The next section begins with a brief review of cluster-robust inference for linear regression models. Then (ref) introduces our new measures of leverage, partial leverage, and influence at the cluster level. (ref) shows how our results can be used to compute the CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ jackknife variance matrix estimators. (ref) discusses what quantities are reported by summclust and should, at least in some cases, be reported by the investigator.
(ref) provides a detailed description of the summclust package which computes these variance estimators and diagnostic measures. The command uses the syntax:
The package has quite a few options and can even be used by itself to estimate a linear regression model with clustered disturbances. The last few sections of the paper illustrate the use of summclust and provide evidence on the value of the measures that it calculates. (ref) presents an empirical illustration in which measures of leverage, partial leverage, and influence are highly informative. (ref) discusses several special cases in which some or all of these measures can be determined analytically. (ref) briefly discusses two\kern 0.04167em-way clustering, where summclust can be valuable even though it is not explicitly designed to handle this case. (ref) describes some simulation experiments which suggest that it may be desirable to report many of the quantities calculated by summclust, and (ref) concludes.
We focus on the linear regression model
where the data have been divided into $G$ disjoint clusters. The \th{g} cluster has $N_g$ observations, so that the sample size is $N = \sum_{g=1}^G N_g$. In \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, ${\bm{X}}_g$ is an $N_g\times k$ matrix of regressors, ${\bm\beta}$ is a $k$-vector of coefficients, ${\bm{y}}_g$ is an $N_g$-vector of observations on the regressand, and ${\bm{u}}_g$ is an $N_g$-vector of disturbances (or error terms). The ${\bm{X}}_g$ may of course be stacked into an $N\times k$ matrix ${\bm{X}}$\kern -0.08333em, and likewise the ${\bm{y}}_g$ and ${\bm{u}}_g$ may be stacked into $N$-vectors ${\bm{y}}$ and ${\bm{u}}$, so that \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} can be rewritten as ${\bm{y}} = {\bm{X}}\kern -0.08333em{\bm\beta} + {\bm{u}}$.
Dividing the sample into clusters only becomes meaningful if we make assumptions about the disturbance vectors ${\bm{u}}_g$ and, consequently, the score vectors ${\bm{s}}_g = {\bm{X}}_g^\top {\bm{u}}_g$. For a correctly specified model, ${\rm E}({\bm{s}}_g)={\bm{0}}$ for all $g$. We further assume that
where ${\bm{\Sigma}}_g$ is the symmetric, positive semidefinite variance matrix of the scores for the \th{g} cluster. The second assumption in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}} is crucial. It says that the scores for every cluster are uncorrelated with the scores for every other cluster. We take the number of clusters $G$ and the allocation of observations to clusters as given. The important issue of how to choose the clustering structure, perhaps by testing for the correct level of clustering, is discussed in detail in \citet*{MNW-testing}.
The OLS estimator of ${\bm\beta}$ is
where the second equality depends on the assumption that the data are actually generated by \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} with true value ${\bm\beta}_0$. It follows that
From the rightmost expression in \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}}, we see that the distribution of $\hat{\bm\beta}$ depends on the disturbance subvectors ${\bm{u}}_g$ only through the distribution of the score vectors ${\bm{s}}_g$. Asymptotic inference commonly uses the empirical score vectors $\hat{\bm{s}}_g = {\bm{X}}_g^\top \hat{\bm{u}}_g$, in which the ${\bm{u}}_g$ are replaced by the residual subvectors $\hat{\bm{u}}_g$, to estimate the variance matrix of the ${\bm{s}}_g$. This should work well if the sum of the ${\bm{s}}_g$, suitably normalized, is well approximated by a multivariate normal distribution with mean zero, and if the ${\bm{s}}_g$ are well approximated by the $\hat{\bm{s}}_g$. However, asymptotic inference can be misleading when either of these approximations is poor.
It follows immediately from \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}} that an estimator of the variance of $\hat{\bm\beta}$ may be based on the usual sandwich formula,
The natural way to estimate \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}} is to replace the ${\bm{\Sigma}}_g$ matrices by their empirical counterparts, that is, the $\hat{\bm{s}}_g\hat{\bm{s}}_g^\top$. If, in addition, we multiply by a correction for degrees of freedom, we obtain the cluster-robust variance estimator, or CRVE,
This is by far the most widely used CRVE in practice, and it is the default one implemented in Stata; alternatives to this estimator will be discussed in (ref). When $G=N$\kern -0.08333em, the CV$_{\kern -0.08333em1}$ estimator reduces to the familiar HC$_1$ estimator \citep*{MW_1985} that is robust only to heteroskedasticity of unknown form.
The fundamental unit of inference for clustered observations is not the observation but the cluster; this is evident from \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}}, \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}}, and \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}. The asymptotic theory for cluster-robust inference has been analyzed by \citet*{DMN_2019} and HansenLee_2019 under the assumption that $G\to\infty$. The quality of the asymptotic approximation depends on the number of clusters $G$ and the heterogeneity of the score vectors \citep*{MNW-guide}. The more the distributions of the scores vary across clusters, the worse the asymptotic approximation will likely be. Heterogeneity can arise from variation in cluster sizes and/or from variation in the distributions of the disturbances, the regressors, or both. As we discuss in (ref), leverage, partial leverage, and summary statistics based on them provide useful measures of heterogeneity across clusters.
Inference about ${\bm\beta}$ is typically based on cluster-robust $t$-statistics and Wald statistics. If $\beta_j$ denotes the \th{j} element of ${\bm\beta}$ and $\beta_{0j}$ is its value under the null hypothesis, then the appropriate $t$-statistic is
where $\hat\beta_j$ is the OLS estimate, and $\textrm{s.e.}(\hat\beta_j)$ is the square root of the \th{j} diagonal element of \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}. Under extremely strong assumptions \citep*{BCH_2011}, it can be shown that $t_j$ asymptotically follows the $t(G-1)$ distribution. Conventional inference in Stata and other programs is based on this distribution.
As the articles cited in the second paragraph of (ref) discuss, inference based on $t_j$ and the $t(G-1)$ distribution can be unreliable when $G$ is small and/or the clusters are severely heterogeneous. This is true to an even greater extent for Wald tests of two or more restrictions \citep*{PT_2018}. The measures of leverage and partial leverage at the cluster level that we introduce in the next section may help to identify the sort of heterogeneity that is likely to make inference unreliable.
Instead of using the $t(G-1)$ distribution, we can obtain both $P$ values for $t_j$ and confidence intervals for $\beta_j$ by employing the wild cluster restricted (or WCR) bootstrap \citep*{CGM_2008}. It can sometimes provide much more reliable inferences than the conventional approach; see (ref). \citet*{RMNW} describes a computationally efficient implementation of this method in the Stata package boottest. \citet*{MNW-bootknife} proposes new versions of the wild cluster bootstrap that involve transforming the empirical scores. When $G$ is reasonably large and the clusters are not very heterogeneous, inferences based on the WCR bootstrap and inferences based on CV$_{\kern -0.08333em1}$ $t$-statistics combined with the $t(G-1)$ distribution will often be very similar. When they differ noticeably, neither should be relied upon without further investigation.
(ref) discusses two CRVEs, which we refer to as CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$, that are both based on the cluster jackknife. In practice, these estimators are often extremely similar. CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ tend to yield more reliable inferences in finite samples than does CV$_{\kern -0.08333em1}$, especially when the clusters are quite heterogeneous; see (ref) and \citet*{MNW-bootknife}. Based on this simulation evidence, we recommend computing either CV$_{\kern -0.08333em3}$ or CV$_{\kern -0.08333em3{\rm J}}$ essentially all the time. This is easy to do using summclust.
At the observation level, there are three classic measures of heterogeneity, namely, leverage, partial leverage, and influence \citep*{BKW_1980,CH_1986}. In this section, we propose analogous measures at the cluster level.
Measures of leverage at the observation level are based on how much the residual for observation $i$ changes when that observation is omitted from the regression. If $h_i$ denotes the \th{i} diagonal element of the “hat matrix” ${\bm{H}} = {\bm{P}}_{\bm{X}} = {\bm{X}}({\bm{X}}^\top\!{\bm{X}})^{-1}{\bm{X}}^\top$\kern -0.08333em, then omitting the \th{i} observation changes the \th{i} residual from $\hat u_i$ to $\hat u_i/(1-h_i)$. Because $0 < h_i < 1$, this delete\kern 0.04167em-one residual is always larger in absolute value than $\hat u_i$. The factor by which the delete\kern 0.04167em-one residual exceeds $\hat u_i$ increases with $h_i$. Since the average of the $h_i$ is $k/N$\kern -0.08333em, observations with values of $h_i$ substantially larger than $k/N$ may reasonably be said to have high leverage.
Dropping the \th{g} cluster when we estimate ${\bm\beta}$ yields the delete\kern 0.04167em-one\kern 0.04167em-cluster estimate $\hat{\bm\beta}^{(g)}$\kern -0.08333em. Using $\hat{\bm\beta}^{(g)}$ in place of $\hat{\bm\beta}$ changes the residual vector for the \th{g} cluster from $\hat{\bm{u}}_g$ to $\hat{\bm{u}}_g^{(g)}$. These delete\kern 0.04167em-one\kern 0.04167em-cluster residual vectors can be written in two ways:
In the rightmost expression above,
is the $N_g\times N_g$ diagonal block of ${\bm{H}}$ that corresponds to cluster $g$. The matrix ${\bm{H}}_g$ is the cluster analog of the scalar $h_i$. Of course, it is not feasible to report the ${\bm{H}}_g$. In fact, when any of the clusters is sufficiently large, even computing and storing these matrices may be challenging. As a measure of leverage, we therefore suggest using a matrix norm of the ${\bm{H}}_g$. Specifically, we suggest the scalar
When the \th{g} cluster contains just one observation, say the \th{i}, then $L_g = h_i$. Thus, in this special case, the leverage measure that we are proposing reduces to the usual measure of leverage at the observation level.
The trace in \hyperref[{eq:traceXX}]{\tagform@{\ref*{eq:traceXX}}} is the nuclear norm of the matrix ${\bm{H}}_g$. In general, the nuclear norm of a matrix ${\bm{A}}$ is the sum of the singular values of ${\bm{A}}$. When ${\bm{A}}$ is symmetric and positive semidefinite, the singular values are equal to the eigenvalues, which are non-negative. Since the trace of any square matrix is equal to the sum of the eigenvalues, the trace of a symmetric and positive semidefinite matrix is also its nuclear norm. In principle, we could report any norm of the ${\bm{H}}_g$ matrices, but the nuclear norm is particularly easy to compute. Also, because it is linear, we can sum over $g$ and take the sum inside the norm just as if the ${\bm{H}}_g$ were scalars. Since $\sum_{g=1}^G {\bm{X}}_g^\top\!{\bm{X}}_g = {\bm{X}}^\top\!{\bm{X}}$, this result means that $G^{-1} \sum_{g=1}^G \operatorname{Tr}({\bm{H}}_g) = k/G$, which is analogous to the result that the average of the $h_i$ over all observations is $k/N$\kern -0.08333em.
High-leverage clusters can be identified by comparing the $L_g$ to $k/G$, their average. If, for some cluster $h$, $L_h$ is substantially larger than $k/G$, then cluster $h$ may be said to have high leverage. Just how much larger $L_h$ has to be is a matter of judgement. A cluster with $L_h=2k/G$ probably does not qualify, but a cluster with $L_h=5k/G$ probably does. Cluster $h$ can have high leverage either because $N_h$ is considerably larger than $G/N$ or because the matrix ${\bm{X}}_h$ is somehow extreme relative to the other ${\bm{X}}_g$ matrices, or both. We can compare the leverage of any two clusters by forming ratios. For example, if $L_1 = 3$ and $L_2 = 1$, then we can say that the first cluster has three times the leverage of the second cluster.
The leverage measure we suggest in \hyperref[{eq:traceXX}]{\tagform@{\ref*{eq:traceXX}}} shows the potential impact of a specified cluster on residuals and fitted values, but not on any particular regression coefficient. When interest focuses on just one such coefficient, say the \th{j}\kern -0.08333em, it may be more interesting to calculate the partial leverage of each cluster. The concept of partial leverage was introduced, for individual observations, in \citet*{CW_1980}. Let
where ${\bm{x}}_j$ is the vector of observations on the \th{j} regressor, and ${\bm{X}}_{[j]}$ is the matrix of observations on all the other regressors. Thus $\acute{\bm{x}}_j$ denotes ${\bm{x}}_j$ after all the other regressors have been partialed out. The partial leverage of observation $i$ is simply the \th{i} diagonal element of the matrix $\acute{\bm{x}}_j(\acute{\bm{x}}_j^\top\acute{\bm{x}}_j)^{-1}\acute{\bm{x}}_j^\top$, which is just $\acute x_{ji}^2/(\acute{\bm{x}}_j^\top\acute{\bm{x}}_j)$, where $\acute x_{ji}^2$ is the \th{i} element of $\acute{\bm{x}}_j$.
The analogous measure of partial leverage for cluster $g$ is
where $\acute{\bm{x}}_{gj}$ is the subvector of $\acute{\bm{x}}_j$ corresponding to the \th{g} cluster. This is what \hyperref[{eq:traceXX}]{\tagform@{\ref*{eq:traceXX}}} reduces to if we replace ${\bm{X}}$ and ${\bm{X}}_g$ by $\acute{\bm{x}}_j$ and $\acute{\bm{x}}_{gj}$, respectively. It is easy to calculate the partial leverage for every cluster for any coefficient of interest. The average of the $L_{gj}$ is evidently $1/G$, so that if cluster $h$ has $L_{hj} >\!> 1/G$, it has high partial leverage for the \th{j} coefficient. Moreover, as we will see in (ref), the empirical distribution of the $L_{gj}$ across clusters seems to provide useful diagnostic information.
AY-IV derives a measure of cluster-level leverage for the first-stage regression used to obtain a linear instrumental variables estimator. The paper calls $L_{gj}$ the group $g$ “share of coefficient leverage” for instrument $j$ and then uses the maximum of the $L_{gj}$ over all the instruments excluded from the structural equation as a measure of the leverage of cluster $g$. Using simulations based on 1309 IV regressions from 30 published papers, the paper finds that inference is much less reliable for models where one or two clusters have high leverage in the first-stage regression than for models where no clusters do so.
One possible consequence of heterogeneity is that the estimates may change a lot when certain clusters are deleted. It can therefore be illuminating to delete one cluster at a time, so as to see how influential each cluster is. To do this in a computationally efficient manner, summclust first computes the cluster-level matrices and vectors
These are then used to construct ${\bm{X}}^\top\!{\bm{X}}$ and ${\bm{X}}^\top{\bm{y}}$, and the vector of least squares estimates when cluster $g$ is deleted is computed as
Unless $k$ is extremely large, it should generally not be expensive to compute $\hat{\bm\beta}^{(g)}$ for every cluster using \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}}. summclust simply has to invert $G$ matrices, each of them $k\times k$, and then multiply each of those matrices by a $k$-vector.
Especially when they vary a lot, the $\hat{\bm\beta}^{(g)}$ can reveal a great deal about the sample. In addition, as we shall see in (ref), they may be used to calculate jackknife variance matrices. When there is a parameter of particular interest, say $\beta_j$, it may be a good idea to report the $\hat\beta_j^{(g)}$ for $g=1,\ldots,G$ in either a histogram or a table. By default, summclust creates several figures with these and other cluster-level statistics. If $\hat\beta_j^{(h)}$ differs greatly from $\hat\beta_j$ for some cluster $h$\kern -0.08333em, then cluster $h$ is evidently influential.
In a few extreme cases, there may be a cluster $h$ for which it is impossible to compute $\hat\beta_j^{(h)}$\kern -0.08333em. This will happen, for example, when the regressor corresponding to $\beta_j$ is a treatment dummy and cluster $h$ is the only treated one. This is an extreme example of the problem of few treated clusters, and inferences based on either the $t(G-1)$ distribution or the WCR bootstrap are likely to be completely unreliable in this case \citep*{MW-TPM,MW-EJ,MW-RI}.
Identifying influential clusters by comparing the $\hat{\bm\beta}^{(g)}$ with $\hat{\bm\beta}$ is very similar to identifying influential observations using the classic methods discussed in \citet*{BKW_1980} and \citet*{CH_1986}; for an interesting recent extension, see \citet*{BGM_2021}. Unlike the leverage measures, the $\hat\beta^{(g)}_j$ may be either positive or negative, must depend on the ${\bm{y}}_g$, and necessarily vary across clusters. They may sometimes reveal features of the model or dataset that require further investigation. Perhaps the model does not seem to apply to some clusters, or perhaps there are measurement errors or observations that have been miscoded.
Regression models often include cluster fixed effects. When one of the regressors is a fixed-effect dummy for cluster $g$, the matrices ${\bm{X}}_g^\top\!{\bm{X}}_g$ and ${\bm{X}}^\top\!{\bm{X}} - {\bm{X}}_g^\top\!{\bm{X}}_g$ are singular. This will cause the calculation in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} to fail unless a generalized inverse routine, such as the invsym routine in Mata, is used. Although summclust uses this routine, it also provides options to avoid the problem, and save some computer time, by partialing out the fixed-effect dummies prior to computing the cluster-level matrices and vectors in \hyperref[{eq:subthings}]{\tagform@{\ref*{eq:subthings}}}; see (ref).
Partialing out cluster fixed effects means replacing ${\bm{X}}$ and ${\bm{y}}$ by $\tilde{\bm{X}}$ and $\tilde{\bm{y}}$, the deviations from their cluster means. For example, the element of $\tilde{\bm{y}}$ corresponding to the \th{j} observation in the \th{g} cluster is $y_{g,j} - N_g^{-1}\sum_{i=1}^{N_g} y_{g,i}$. The \th{g} subvector of $\tilde{\bm{y}}$ is $\tilde{\bm{y}}_g$, and the \th{g} submatrix of $\tilde{\bm{X}}$ is $\tilde{\bm{X}}_g$. Since there is just one fixed effect per cluster, $\tilde{\bm{y}}_g$ depends solely on ${\bm{y}}_g$, and $\tilde{\bm{X}}_g$ depends solely on ${\bm{X}}_g$. The calculations in \hyperref[{eq:traceXX}]{\tagform@{\ref*{eq:traceXX}}} and \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} are now based on $\tilde{\bm{X}}^\top\!\tilde{\bm{X}}$\kern -0.08333em, $\tilde{\bm{X}}^\top\tilde{\bm{y}}$, the $\tilde{\bm{X}}_g^\top\!\tilde{\bm{X}}_g$, and the $\tilde{\bm{X}}_g^\top\tilde{\bm{y}}_g$. Importantly, the sum of the $L_g$ is now equal to the number of columns in $\tilde{\bm{X}}$ instead of the number of columns in ${\bm{X}}$\kern -0.08333em.
Although the CV$_{\kern -0.08333em1}$ variance estimator defined in \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}} is very widely used, it often does not have good finite\kern 0.04167em-sample properties. Two alternative CRVEs, which are usually known as CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$, were proposed in \citet*{BM_2002}. They are the cluster analogs of the hetero\-skedasticity-consistent estimators HC$_2$ and HC$_3$, which are appropriate when the $u_i$ are independent. These names were coined in \citet*{MW_1985}, which proposed HC$_3$ as a jackknife variance estimator. In the remainder of this section, we focus on CV$_{\kern -0.08333em3}$, because CV$_{\kern -0.08333em2}$ is not a jackknife estimator and is not amenable to the computational methods that we propose; for more on it, see \citet*{Imbens_2016}, \citet*{PT_2018}, and \citet*{NAAMW_2020}. Stata 18 added the ability to rapidly calculate CV$_{\kern -0.08333em2}$ standard errors, using the option vce(hc2 clustvar). Simulations in \citet*{MNW-bootknife} suggest that CV$_{\kern -0.08333em2}$ is preferred to CV$_{\kern -0.08333em1}$, but that CV$_{\kern -0.08333em3}$ is almost always preferred to CV$_{\kern -0.08333em2}$.
CV$_{\kern -0.08333em3}$ can be written in several ways. One of them is
where the modified score vectors $\ddot{\bm{s}}_g$ are defined as
Here ${\bm{M}}_{gg} = {\bf I}_{N_g} - {\bm{X}}_g({\bm{X}}^\top{\bm{X}})^{-1}\kern -0.08333em{\bm{X}}_g^\top$ is the diagonal block corresponding to the \th{g} cluster of the projection matrix ${\bm{M}}_{\bm{X}}$, which satisfies $\hat{\bm{u}} = {\bm{M}}_{\bm{X}}{\bm{u}}$. Although computing CV$_{\kern -0.08333em3}$ using \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} works well when all the $N_g$ are very small, it becomes expensive, or perhaps computationally infeasible, when one or more of the $N_g$ is large. The problem is that an $N_g\times N_g$ matrix needs to be stored and inverted for every cluster. \citet*{NAAMW_2020} proposes a method that is much faster for large clusters, versions of which apply to both CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$. However, recognizing that CV$_{\kern -0.08333em3}$ is a jackknife estimator makes a method that is even simpler and usually faster available.
There are actually two cluster jackknife estimators of $\operatorname{Var}(\hat{\bm\beta})$. The simplest is probably
where $\bar{\bm\beta}$ is the sample mean of the $\hat{\bm\beta}^{(g)}$, which were defined in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}}. The expression in \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} is the cluster analog of the usual jackknife variance matrix estimator given in \citet*[eqn.\ (11)]{MW_1985}. Each of the $\hat{\bm\beta}^{(g)}$ is obtained by deleting a cluster instead of an observation, and the summation is over clusters instead of observations. If $\bar{\bm\beta}$ in \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} is replaced by $\hat{\bm\beta}$, we obtain instead
This version of CV$_{\kern -0.08333em3}$ is numerically identical to the one in \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} \citep*[Section 3]{MNW-bootknife}. Unless all the clusters are very small, computing CV$_{\kern -0.08333em3}$ using \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}} is much faster than using \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}; timings are reported in \citet*{MNW-bootknife}.
Many discussions of jackknife variance estimation follow \citet*{Efron_79} and use $\bar{\bm\beta}$ as in \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}}, but others, including \citet*{BM_2002}, use $\hat{\bm\beta}$ as in \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}}. Although these two jackknife variance estimators are asymptotically the same, they are rarely equal, since CV$_{\kern -0.08333em3}$ minus CV$_{\kern -0.08333em3{\rm J}}$ is a positive semi-definite matrix. In practice, however, they tend to be very similar \citep*{MNW-bootknife}, and there seems to be no good reason to expect either CV$_{\kern -0.08333em3}$ or CV$_{\kern -0.08333em3{\rm{J}}}$ to perform better in general. Interestingly, the original HC$_3$ estimator proposed in \citet*{MW_1985} is actually the analog of CV$_{\kern -0.08333em3{\rm J}}$. The modern version of HC$_3$, which is the analog of CV$_{\kern -0.08333em3}$, seems to be due to \citet*[Chapter 16]{DM_1993}. This version of HC$_3$ is normally computed by dividing each residual by the corresponding diagonal element of ${\bm{M}}_{\bm{X}}$, and the factor of $(N-1)/N$ is usually (but incorrectly) omitted.
The factor of $(G-1)/G$ in both \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} and \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}} is designed to compensate for the tendency of the $\hat{\bm\beta}^{(g)}$ to be too spread out. This factor is the analog of the usual factor of $(N-1)/N$ for a jackknife variance matrix at the individual level. It implicitly assumes that all clusters are the same size and perfectly balanced, with disturbances that are independent and homoskedastic. In this special case, the estimators CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm{J}}}$ would be identical and unbiased \citep*{BM_2002}. These estimators are already available in Stata. When used with the cluster option, the vce(jackknife) option computes CV$_{\kern -0.08333em3{\rm{J}}}$ standard errors, and the vce(jackknife,mse) option computes CV$_{\kern -0.08333em3}$ standard errors. Because it is specialized for linear regression models, the implementation in \texttt{summclust} is quite a bit faster.
Both jackknife estimators may readily be used to compute cluster-robust $t$-statistics. Because there are $G$ terms in the summation, it seems natural to compare these with the $t(G-1)$ distribution, as usual. These procedures should almost always be more conservative than $t$-tests based on the widely-used CV$_{\kern -0.08333em1}$ estimator. In an important recent paper, \citet*{Hansen-jack} shows that CV$_{\kern -0.08333em3}$ has much better worst-case theoretical properties than CV$_{\kern -0.08333em1}$. This strongly suggests that $t$-statistics based on CV$_{\kern -0.08333em3}$ are likely to yield lower rejection frequencies than ones based on CV$_{\kern -0.08333em1}$. The simulation results in (ref) and in \citet*{MNW-bootknife} are consistent with this conjecture.
When a model includes fixed effects, some care needs to be taken when computing CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$. As noted in (ref), it is computationally attractive to partial out fixed effects prior to calculating $\hat{\bm\beta}$. However, if we were to partial out any arbitrary regressors prior to computing the delete\kern 0.04167em-one\kern 0.04167em-cluster estimates in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}}, then the computed $\hat{\bm\beta}^{(g)}$ would depend on the values of the partialed-out regressors for the full sample, including those in the \th{g} cluster. Consequently, the values of CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ will be incorrect if we partial out any regressor that affects more than one cluster (such as industry-level fixed effects with firm-level clustering). The regressors that are partialed out must be cluster fixed effects or fixed effects at a finer level (such as firm-level fixed effects with industry-level clustering), because each of them affects only one cluster. See the discussion of the absorb and fevar options in (ref).
It is possible that the vector ${\bm\beta}$ is identified for the full sample but not when one cluster is deleted. For example, consider the coefficient on a dummy variable that takes on non-zero values only for observations in the \th{g} cluster. This coefficient cannot be identified when cluster $g$ is omitted. In such a case, the matrix ${\bm{X}}^\top\!{\bm{X}} - {\bm{X}}_g^\top\!{\bm{X}}_g$ in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} is singular, and CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ cannot be computed using an ordinary matrix inverse. However, because summclust uses the invsym function in Stata, which implements a generalized inverse, the offending element of $\hat{\bm\beta}^{(g)}$ is simply replaced by $0$. The package therefore checks whether any of the $\hat\beta^{(g)}$ coefficients of interest are equal to $0$ and issues a warning when they are; see (ref).
There may be more than one set of fixed effect that are invariant at the cluster level. For example, imagine an analysis of students' test scores where the researcher wants to control for both school and neighborhood fixed effects and cluster the standard errors at the state level. In this case, neither of Stata's built-in regress and areg commands can produce an estimate of CV$_{\kern -0.08333em3}$, because the fixed effects for schools and neighborhoods in state $g$ cannot be identified when state $g$ is omitted. However, summclust can produce such an estimate.
We believe that investigators should routinely compute the $L_g$. They should also compute the $L_{gj}$ for any coefficient(s) of particular interest. In some cases, the $L_g$ and the $L_{gj}$ will be roughly proportional to the $N_g$ (the cluster sizes). That in itself would be informative. It may be even more interesting, however, to find that the relative size of $L_g$ and/or $L_{gj}$ for some cluster(s) $g$ is much larger, or much smaller, than the relative size of $N_g$.
When the number of clusters is small, it is easy enough to look at all the $N_g$, $\hat\beta^{(g)}_j$, $L_g$, and $L_{gj}$ to see whether any clusters are unusually large, unusually influential, or have unusually high leverage or partial leverage. Once $G$ exceeds 10 or 12, however, it may be more informative to report summary statistics or to plot these quantities. The summclust package always reports the minimum, first quartile, median, mean, third quartile, and maximum of the $N_g$ and the $L_g$. It also reports these quantities for the $L_{gj}$ and the $\hat\beta^{(g)}_j$ for the specified regressor $j$, and by default it provides a figure containing four scatterplots of the $L_g$ and the $L_{gj}$ against the $N_g$ and the $\hat\beta^{(j)}$; see (ref).
Another possibility is to report a few summary statistics, as summclust also does. Consider a generic (positive) quantity $a_g$, which might denote any of $N_g$, $L_g$, or $L_{gj}$ for $g=1,\ldots,G$. It seems plausible that inference may be unreliable when any of the $a_g$ vary substantially across clusters, and we provide some evidence to support this conjecture in (ref).
There are many measures of how much the distribution of the $a_g$ differs from what it would be in the perfectly balanced case. One of these is the scaled variance
where the argument $a_{\bullet}$ is to be interpreted as the entire set of $a_g$ for $g=1,\ldots,G$, and $\bar a$ denotes the arithmetic mean, which is $N/G$ for the $N_g$, $k/G$ for the $L_g$, and $1/G$ for the $L_{gj}$. These are all positive numbers, so it is reasonable to scale by their squares. Larger values of $V_s$ imply that the $a_g$ are more variable across clusters, relative to their mean. We could report either $V_s$ or its square root, which is often called the coefficient of variation. In the perfectly balanced case, $V_s=0$. By default, summclust reports the coefficient of variation for the cluster sizes, the leverages, the partial leverages, and the $\hat\beta_j^{(g)}$\kern -0.08333em.
Another possibility, which is only valid for positive quantities, is to report one or more alternative sample means. The more these differ from the arithmetic mean, the more heterogeneous must be the clusters. Three common alternatives to the arithmetic mean are the harmonic, geometric, and quadratic means:
Unless all the $a_g$ are the same, the harmonic and geometric means will be less than the arithmetic mean $\bar a$, and the quadratic mean (which has the same form as the root mean squared error of an estimator) will be greater than $\bar a$. summclust optionally reports all three of these alternative means, along with the ratio of each of them to $\bar a$. The three ratios provide scale\kern 0.04167em-free measures of cluster heterogeneity; the closer they are to one, the more homogeneous are the clusters.
Another way to quantify the heterogeneity of the cluster sizes and the regressors is to calculate $G^*$\kern -0.08333em, the “effective number of clusters,” as proposed in \citet*{CSS_2017}. The value of $G^*$ depends on the coefficient $j$ for which it is being computed and on a parameter $\rho$ to be discussed below, so we denote it $G^*_j(\rho)$. It is defined as
where $0\le\rho\le1$, and the $\gamma_{gj}(\rho)$ are given by
Here ${\bm{e}}_j$ is a $k$-vector with 1 in the \th{j} position and 0 everywhere else, so that ${\bm{e}}_j^\top({\bm{X}}^\top{\bm{X}})^{-1}$ is the \th{j} row of $({\bm{X}}^\top{\bm{X}})^{-1}$\kern -0.08333em, and ${\bm{\Omega}}_g(\rho)$ is an $N_g\times N_g$ matrix with 1 on the principal diagonal and $\rho$ everywhere else. It is easy to see that
where ${\bm\iota}$ is an $N_g$-vector of 1s, and ${\bf I}$ is an $N_g\times N_g$ identity matrix. Notice that $\Gamma_j(\rho)$ is just the scaled variance of the $\gamma_{gj}(\rho)$; compare \hyperref[{eq:coefvar}]{\tagform@{\ref*{eq:coefvar}}}.
The parameter $\rho$ may be interpreted as the intra-cluster correlation coefficient for a model with cluster-level random effects. Since $\rho$ is unknown, \citet*{CSS_2017} suggests calculating $G^*_j(1)$ as a sort of worst case. However, when there are cluster-level fixed effects, or fixed effects at a finer level nested within clusters, they will absorb all of the intra-cluster correlation. Thus it does not make sense to specify $\rho>0$ in either of these cases. It does seem natural to use $G^*_j(0)$, however, because the amount of intra-cluster correlation that remains in models with cluster fixed effects is often quite small.
From \hyperref[{eq:gammag}]{\tagform@{\ref*{eq:gammag}}} and \hyperref[{eq:omega}]{\tagform@{\ref*{eq:omega}}}, we see that
This result makes it inexpensive to compute the $\gamma_{gj}(\rho)$ for any value of $\rho$ by first computing them for $\rho=0$ and $\rho=1$. The needed equations are
where ${\bm{w}}_j$ is the \th{j} column of $({\bm{X}}^\top{\bm{X}})^{-1}$. After obtaining the $\gamma_{gj}(\rho)$ from \hyperref[{eq:gammarho}]{\tagform@{\ref*{eq:gammarho}}}, it is trivial to compute $G^*_j(\rho)$ using \hyperref[{eq:gstar}]{\tagform@{\ref*{eq:gstar}}}. Evidently, $G^*_j(\rho)$ is always less than $G$. When it is much smaller than $G$, it can provide a useful warning.
Suppose that we have partialed out cluster fixed effects prior to computing $G^*_j(\rho)$. Then the first term on the right-hand side of \hyperref[{eq:gammiddle}]{\tagform@{\ref*{eq:gammiddle}}} should in theory be a zero matrix, because every column of ${\bm{X}}_g$ should add to zero. In practice, however, the limitations of floating-point arithmetic mean that this matrix will actually contain extremely small positive numbers. This will cause the computation of $G^*_j(\rho)$ to be numerically unstable. When the fixed effects are not partialed out, similar but more complicated numerical issues arise.
The Stata package clusteff discussed in \citet*{LS_2018} is designed to calculate $G^*_j(\rho)$, with $\rho=0.9999$ rather than $\rho=1$ by default to avoid numerical instabilities. However, the only version of this package that we have used does so in a computationally inefficient way that does not use \hyperref[{eq:gammarho}]{\tagform@{\ref*{eq:gammarho}}}. When any of the $N_g$ is large, it can take a very long time, or even fail because Stata runs out of memory. For example, it failed with some of the samples in \citet*{MNW-guide}.
Like $V_s(a_{\bullet})$ and the alternative sample means for measures of leverage and partial leverage discussed above, $G^*_j(\rho)$ is sensitive not only to variation in cluster sizes but also to other features of the ${\bm{X}}_g$ matrices. But it is not sensitive to heteroskedasticity or to any other features of the disturbances. summclust computes $G^*_j(0)$, $G^*_j(1)$, and (optionally) $G^*_j(\rho)$ for a specified covariate. However, when there are cluster fixed effects, or fixed effects nested within clusters, it only computes $G^*_j(0)$. For example, it will not compute $G^*_j(\rho)$ for $\rho\ne0$ whenever there are state\kern 0.04167em-level fixed effects and clustering at the region level.
The quantity $G^*_j(0)$ is very closely related to $V_s(L_{{\bullet}j})$, where $L_{{\bullet}j}$ denotes the entire set of $L_{gj}$, for $g=1,\ldots,G$. It is not hard to see that the $\gamma_g(0)$ defined in \hyperref[{eq:gammag}]{\tagform@{\ref*{eq:gammag}}} and \hyperref[{eq:gammarho}]{\tagform@{\ref*{eq:gammarho}}} are equal to the $L_{gj}$ defined in \hyperref[{eq:partial}]{\tagform@{\ref*{eq:partial}}} divided by $\acute{\bm{x}}_j^\top\acute{\bm{x}}_j$. Since this makes the $\gamma_g(0)$ proportional to the $L_{gj}$, $V_s(L_{{\bullet}j})$ is numerically identical to $\Gamma(0)$; compare \hyperref[{eq:coefvar}]{\tagform@{\ref*{eq:coefvar}}} and the middle equation in \hyperref[{eq:gstar}]{\tagform@{\ref*{eq:gstar}}}. Thus we see from the first equation in \hyperref[{eq:gstar}]{\tagform@{\ref*{eq:gstar}}} that $G^*_j(0)$ is simply a monotonically decreasing function of the scaled variance of our measures of partial leverage at the cluster level. When $V_s(L_{{\bullet}j})$ is large, $G^*_j(0)$ is necessarily much smaller than $G$.
The summclust package may be obtained from SSC or \url{https://github.com/mattdwebb/summclust}. It implements the summclust command, which calculates a large number of statistics to help assess cluster heterogeneity and also provides CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ standard errors. The package does not rely on any other Stata packages, but it does require a version of Stata that provides Mata's \texttt{panelsum()} function (Version 13 or later).
We first present an overview of the summclust command, followed by a simple illustration using the webuse dataset nlswork.
Syntax
\noindentvarlist: the dependent variable, the independent variable of interest, and other (binary or continuous) independent variables. At least one additional regressor must be specified. Time\kern 0.04167em-series operators and factor variables are not permitted.
\noindentcluster: the clustering variable, for which the number of unique values equals $G$.
\vskip -12pt \hrule
\noindentDescription
summclust is a stand-alone command for summarizing cluster variability in several ways. It always calculates measures of cluster-level influence and leverage, and it optionally calculates the effective number of clusters. It also always reports CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ standard errors for a single coefficient, and it optionally reports a CV$_{\kern -0.08333em3{\rm J}}$ standard error as well. If requested, it can calculate additional measures of cluster-level heterogeneity. Unless it is told not to, it produces a figure which can help identify the source of cluster level heterogeneity. Finally, it can optionally produce a full table of regression results with CV$_{\kern -0.08333em3}$ standard errors.
By default, summclust calculates the CV$_{\kern -0.08333em3}$ standard error based on \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}. With well-behaved samples, this should match the standard error calculated using either Stata's native jackknife: \textbf{\texttt{reg y x, cluster(group)}} or \textbf{\texttt{reg y x, cluster(group) vce(jackknife)}} commands. However, many samples are not well-behaved, in the sense that the regressor matrices for some of the omit-one-cluster subsamples may not have full rank. We will refer to such subsamples, rather informally, as “singular subsamples.”
Whenever there are singular subsamples, summclust calculates two standard errors. One of these drops the singular subsamples, as the native Stata commands do. The other uses a generalized inverse. summclust provides guidance as to which standard error is likely to be more reliable. When regtable is specified, and singular subsamples are present, two versions of the regression table are displayed. Similarly, if \textbf{\texttt{jackknife}} is specified and there are singular subsamples, four different standard errors are shown, either CV$_{\kern -0.08333em3}$ or CV$_{\kern -0.08333em3{\rm J}}$, combined with either the generalized inverse or one that drops the singular subsamples.
nograph suppresses creation of the figure, which is otherwise shown by default. The figure shows four scatter plots: leverage against observations per cluster, partial leverage against observations per cluster, leverage against omit-one-cluster coefficients, and partial leverage against omit-one-cluster coefficients. This figure can be quite informative, but it is computationally costly to produce. We recommend invoking this option after the figure has been inspected.
When jackknife is specified, regtable uses the CV$_{\kern -0.08333em3{\rm J}}$ estimates to produce the regression table. Otherwise, CV$_{\kern -0.08333em3}$ estimates are used.
To illustrate summclust's functionality and syntax, we consider a simple example using the online dataset nlswork, which contains a sample of women who were 14--26 years of age in 1968 from the National Longitudinal Survey of Young Working Women. For the purposes of this exercise, we restrict the sample to individuals who are 20 to 40 years old.
We estimate a simple Mincer regression using the nlswork dataset to see whether there is a marriage premium for wages. The variable msp is equal to 1 if the person is married and cohabits with their spouse, and equal to 0 otherwise. For the purposes of this example, we cluster by industry. The following code opens the dataset and estimates the regression using Stata's regress command:
The Stata output from the command above provides CV$_{\kern -0.08333em1}$ standard errors. Alternatively, we can estimate CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ standard errors using this code:
When either of these commands is run, Stata displays the warning “Note: One or more parameters could not be estimated in 2 jackknife replicates; standard-error estimates include only complete replications.”
The coefficient on msp and two or three standard errors can also be obtained using summclust. The basic command is:
This code results in the default output from summclust, which is mostly contained in two tables. The first one includes the coefficient on the second variable in the varlist (in this case msp), the CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ standard errors for this coefficient, and the associated $t$-statistics, $P$ values, and confidence intervals. In this case, summclust also displays a warning about singular subsamples and thus produces two “Regression Output” tables. The standard errors in the table which drops singular subsamples match those produced natively in Stata.
In the first table for this example, the CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ standard errors are noticeably different, with the latter being considerably larger. However, in the second table, where the two singular subsamples are dropped, the CV$_{\kern -0.08333em3}$ standard error becomes much smaller.
The “Cluster Variability” table from summclust (below) provides insight into what is happening. It reports summary statistics for $N_g$, $L_g$, $L_{gj}$, and $\hat\beta_{j}^{(g)}$. Whenever singular subsamples are dropped, two sets of statistics are shown for $\hat\beta_{j}^{(g)}$. The first (second-last column) uses all the jackknife subsamples with a generalized inverse standard error. The second (final column) uses only the non-singular subsamples. We can see that the largest value of $\hat\beta_{j}^{(g)}$ is considerably smaller (and therefore more different from the other values) when none of the subsamples is dropped. This explains why the CV$_{\kern -0.08333em3}$ standard error is larger in the first table above than in the second one.
\eject
It is evident from this table that the clusters are extremely heterogeneous. The largest cluster contains almost one\kern 0.04167em-third of the sample and is 167 times the size of the smallest. There are also extreme differences in both leverage and partial leverage across clusters. The ratio of the largest to the smallest value is 198 for leverage and 504.5 for partial leverage. The sum of the leverages is $12 \times 4.583333 = 55$, which is the number of estimated coefficients. Although both sets of $\hat\beta_{j}^{(g)}$ vary quite a bit, dropping one cluster never changes the sign of the coefficient.
The option fevar is used when there are factor variables, which would be specified as i.varname in conventional Stata syntax. In the above example, the arguments to fevar are grade, age, and \texttt{birth_yr}. For each argument, a set of temporary dummy variables is created. These dummy variables are included in the regression, and there is no constant term if they are present.
The sample code above does not illustrate several additional options. The most important of these is the absorb option, which operates like fevar. It treats its argument, a single variable, as an additional factor variable to include in the set of regressors. absorb(varname) can be used when including i.varname in a regression would result in many fixed effects. Speed can often be increased, perhaps substantially, by partialing out the absorbed fixed effects from the dependent and all the independent variables. It is advisable to use absorb rather than fevar whenever their argument corresponds to a set of cluster fixed effects, since the elements of $\hat{\bm\beta}^{(g)}$ that correspond to the fixed effects cannot be identified in that case; see (ref).
The absorb option should be used with care. Partialing out fixed effects is valid for the measures of leverage and influence and for the jackknife variance matrices only when the absorbed variable yields fixed effects that can be partialed out on a cluster-by-cluster basis. That is, absorb should only be used for straight cluster fixed effects or for fixed effects at a finer level, such as state $\times$ year fixed effects for a panel with clustering at the state level. It is not valid to partial out fixed effects that are not limited to a single cluster. In that case, the $\hat{\bm\beta}^{(g)}$ and quantities based on them would be different for the original data and the data after partialing out, because the partialed-out observations for the \th{g} cluster would depend on other clusters as well. Accordingly, summclust checks to ensure that the clustering variable is invariant within each value of the absorbed variable. When it is not invariant, a warning is displayed, and the values of $L_{g}$, $L_{gj}$, $\hat\beta_{j}^{(g)}$\kern -0.08333em, CV$_{\kern -0.08333em3}$, and CV$_{\kern -0.08333em3{\rm J}}$ are not available.
To see the difference between fevar and absorb, we can estimate an expanded regression that includes industry fixed effects. Consider the following two commands:
For the command which uses fevar for all the categorical variables, some of the output is
Because every one of the jackknife subsamples is singular, only the results based on the generalized inverse are reported. In contrast, when absorb is used for the industry fixed effects, the corresponding output is instead
These two tables highlight a key reason for using absorb. Because only two of the jackknife subsamples are singular, summclust is able to report both standard errors. Observe that, when all 12 jackknife samples are used, the standard errors are the same regardless of whether industry fixed effects are specified using fever or absorb.
Whether we use fevar or absorb leads to somewhat different output for the measures of cluster variability.
The $\hat\beta_j^{(g)}$ when all clusters are retained are identical for both options. But since there are two singular subclusters, there are two versions of the $\hat\beta_j^{(g)}$ for the fevar results.
The leverage estimates are also smaller when we use the absorb option. Recall that, for the original model with no industry fixed effects, the leverages summed to 55. In the first case just above, where the industry fixed effects are included as regressors in fevar, the regression has 66 coefficients, and the leverages therefore sum to $12 \times 5.5 = 66$. In the second case, where the industry fixed effects are partialed out using absorb, the regression has 54 coefficients, and the leverages therefore sum to $12 \times 4.5 = 54$. Thus for the first case, each of the leverages is larger than the corresponding one for the second case by precisely 1.
\noindentExamples
In the examples that follow, we include the nograph option to reduce computational time.
This example illustrates the jackknife and table options:
In addition to the two standard tables, it displays the following table:
This table makes it easy to see whether the high leverage clusters are also the largest clusters. That is clearly the case here. After running the program, this table is stored as the Mata matrix scall.
To obtain summary statistics on the four (or five) measures of cluster variability, we can use the addmeans option:
This command produces the following table:
Once again, we see that there is extreme variability across the clusters. This is particularly noticeable for the ratio of the harmonic mean to the arithmetic mean, which is between 0.125 and 0.143 for the cluster size, leverage, and partial leverage measures. Recall that these ratios would be close to one if the clusters were relatively homogeneous. This table is stored in Mata's memory as bonus.
To obtain estimates of the effective number of clusters, we can use either the gstar option or the rho() option. The former displays $G_j^*(0)$ and $G_j^*(1)$. The latter requires a specified value of $\rho$ and displays $G_j^*(0)$ and $G_j^*(1)$ along with $G_j^*(\rho)$. When there are fixed effects at the cluster or subcluster level, only $G_j^*(0)$ is reported.
For the nlswork example, the first option may be called as:
This yields:
The second option, using $\rho = 0.5$ as an illustration, may be called as:
This yields:
In this example, it is clear that the effective number of clusters is substantially less than the actual number of clusters. This provides more evidence that inference using the CV$_{\kern -0.08333em1}$ standard error together with the $t(G-1)$ distribution is likely to be unreliable. These three quantities can be accessed in Mata's memory as gstarzero, gstarrho, and gstarone, respectively.
By using the regtable option, one can display a modified version of the regression table, which is similar to the default output from Stata's regress command. The command is:
When there are singular subsamples, two versions of this table will be displayed. In this example, the table is quite long, so we do not reproduce it here.
All the results that are displayed as output can also be found in Mata's memory. To access one of these after running summclust, simply add the following line:
The \verb+object_name+ can take one of the following values:
\hangindent=\parindent\hangafter=1 {\tt cvstuff}: This matrix stores the table with the title “Regression Output”. It is $2\times6$ when the jackknife option is not used (the default), and $3\times6$ when jackknife is used.
\hangindent=\parindent\hangafter=1 {\tt scall}: This matrix stores the $G\times4$ table created by the table option with the title “Cluster by Cluster Statistics”.
\hangindent=\parindent\hangafter=1 {\tt bonus}: This $6\times4$ matrix contains the alternative sample means and their ratios to the arithmetic mean created by the addmeans option.
\hangindent=\parindent\hangafter=1 {\tt gstarzero}: This scalar contains $G^*(0)$ created by the gstar or rho options.
\hangindent=\parindent\hangafter=1 {\tt gstarone}: This scalar contains $G^*(1)$ created by the gstar or rho options.
\hangindent=\parindent\hangafter=1 {\tt gstarrho}: This scalar contains $G^*(\rho)$ created by the rho option.
\hangindent=\parindent\hangafter=1 {\tt regresstab}: This matrix contains the table shown when the regtable option is specified.
Scalars within matrices can be referenced on a cell-by-cell basis. For example, the CV$_{\kern -0.08333em3}$ standard error is stored in the second row and second column of cvstuff, and to display it one can enter the following command:
Additionally, several results are available as scalars or matrices in return memory using r(). The available scalars are:
\hangindent=\parindent\hangafter=1 {\tt beta}: The estimate $\hat\beta$ for the coefficient of interest.
\hangindent=\parindent\hangafter=1 {\tt cv1se}: The CV$_{\kern -0.08333em1}$ standard error for the coefficient of interest.
\hangindent=\parindent\hangafter=1 {\tt cv1t}: The CV$_{\kern -0.08333em1}$ $t$-statistic for the coefficient of interest.
\hangindent=\parindent\hangafter=1 {\tt cv1p}: The $P$ value for the null hypothesis that $\beta = 0$ for the coefficient of interest using the CV$_{\kern -0.08333em1}$ standard error.
\hangindent=\parindent\hangafter=1 {\tt cv1lci}: The lower bound of the 95% confidence interval for $\beta$ using the CV$_{\kern -0.08333em1}$ standard error.
\hangindent=\parindent\hangafter=1 {\tt cv1uci}: The upper bound of the 95% confidence interval for $\beta$ using the CV$_{\kern -0.08333em1}$ standard error.
\hangindent=\parindent\hangafter=1 {\tt gstarzero}: The effective number of clusters for the coefficient of interest using $\rho=0$.
\hangindent=\parindent\hangafter=1 {\tt gstarone}: The effective number of clusters for the coefficient of interest using $\rho=1$.
\hangindent=\parindent\hangafter=1 {\tt gstarrho}: The effective number of clusters for the coefficient of interest using the value of $\rho$ specified in rho($\rho$).
The standard error, $t$-statistic, $P$ value, and confidence interval bounds are also available for the CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ standard errors. To access these, replace “1” in the above with either “3” or “3J”; for example, the $P$ value using CV$_{\kern -0.08333em3{\rm J}}$ is available in cv3Jp. In the event of singular subsamples, there are two versions of the CV$_{\kern -0.08333em3}$ or CV$_{\kern -0.08333em3{\rm J}}$ results. The ones where singular subsamples have been dropped have a suffix of `drop'. For instance, cv3sedrop is used instead of cv3se.
The available matrices are:
\hangindent=\parindent\hangafter=1 {\tt ng}: This $G\times1$ matrix contains the number of observations, $N_g$, for each cluster.
\hangindent=\parindent\hangafter=1 {\tt leverage}: This $G\times1$ matrix contains the leverage, $L_g$, for each cluster.
\hangindent=\parindent\hangafter=1 {\tt partlev}: This $G\times1$ matrix contains the partial leverage, $L_{gj}$, for each cluster.
\hangindent=\parindent\hangafter=1 {\tt betanog}: This $G\times1$ matrix contains the $\hat\beta_j^{(g)}$ for each cluster.
We consider an empirical example from \citet*{BG_2019}, which studies an experiment where retail firms were randomly assigned to enter one of 72 different geographic markets (in Spanish, mercados), within the Dominican Republic. After randomization, 21 markets had no entrants and so were in the control group, 18 had one entrant, another 18 had two, and the remaining 15 had three. The primary analysis only distinguishes between the 51 treated markets and the 21 control markets. The number of observations (stores) per market varies from 20 to 55.
This example is interesting because conventional wisdom \citep*[e.g.,][]{MNW-guide} suggests that, with 72 clusters that do not vary much in size, and with neither few treated nor few control clusters, inference based on CV$_{\kern -0.08333em1}$ standard errors and the $t(71)$ distribution should work well. However, our leverage measures suggest otherwise, and alternative inference methods yield noticeably different results.
The model we estimate is
Here $s$ indexes stores, and $d$ indexes markets. The treatment variable $Z_d$ equals 1 if market $d$ is treated (there was entry) and 0 if it was a control (there was no entry). The coefficient of interest is $\gamma$, which measures the causal effect of increased competition on an outcome $Y$\kern -0.08333em. We focus on just one of several outcomes, namely, the log of demeaned prices after treatment. The results from this regression are found in Table 5, Panel B, column 4, row 1 of \citet*{BG_2019}. The table states that there are 72 clusters and $2,\kern -0.08333em025$ observations; however, the replication dataset that we use contains just $1,\kern -0.08333em926$ observations.
Regression \hyperref[{mainreg}]{\tagform@{\ref*{mainreg}}} includes 17 control variables in the row vector ${\bm{X}}_{sd}$. These are the first lag of the outcome variable, the number of retailers in each district pre\kern 0.04167em-treatment, a lagged quality index, eight province fixed effects, total district beneficiaries of a conditional cash transfer program, percent beneficiaries of that program, average income in the market, two market education measures, and a binary indicator for the urban status of the market. Thus the total number of regressors is 19.
\setcounter{table}{0}
The OLS estimate of $\gamma$, its CV$_{\kern -0.08333em1}$ standard error, the $P$ value for a test that $\gamma=0$, and a .95 confidence interval are shown in the first row of (ref). Allowing for different numbers of reported digits, these estimates accord with the ones in \citet*{BG_2019}. The estimate of $-0.01469$ has the expected sign (average prices declined). However, the $P$ value is just slightly less than 0.05, and the confidence interval barely excludes 0.
We next use the summclust package to calculate the cluster-level characteristics of the model and dataset. Some key ones are reported in (ref). It is evident that cluster sizes are well balanced, varying from 20 to 55, with the first and third quartiles equal to 24 and 27. However, both the leverages $L_g$ and the partial leverages $L_{g1}$ vary considerably. The former range from 0.1308 to 0.7378, and the latter from 0.0001 to 0.0642. The coefficients of variation are 0.3887 and 1.0598, respectively. The latter is moderately large, although not enormous. The two values of $G^*$ are slightly smaller than $G/2$, which also suggests that the sample is not well balanced.
Most of the $\hat\gamma^{(g)}$ do not vary much, and thus their coefficient of variation is small. However, the most extreme values are notable. The estimate of $\gamma$, which is $-0.01469$, could be as small as $-0.01755$ or as large as $-0.01072$ if just one out of 72 clusters were dropped.
These results suggest that CV$_{\kern -0.08333em1}$, the default CRVE, may not be particularly reliable in this case. We therefore consider five alternative procedures. The second, third, and fourth rows of (ref) report the CV$_{\kern -0.08333em2}$, CV$_{\kern -0.08333em3}$, and CV$_{\kern -0.08333em3{\rm J}}$ standard errors, along with the $P$ values and confidence intervals associated with them. The CV$_{\kern -0.08333em2}$ $P$ value is noticeable larger than the CV$_{\kern -0.08333em1}$ one and suggests that the estimate is not significant at the .05 level. The CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ rows are almost identical. At $0.1105$, the CV$_{\kern -0.08333em3}$ $P$ value does not even allow us to reject the null at the .10 level. The fifth and six rows of (ref) report two WCR bootstrap $P$ values and the associated .95 confidence intervals. At 0.0891 and 0.0913, these are a bit smaller than the jackknife ones, but they clearly do not allow us to reject the null hypothesis at the .05 level.
In view of the reasonably large number of clusters and the fact that cluster sizes do not vary much, the large discrepancy between the results for CV$_{\kern -0.08333em1}$ and the other procedures may seem surprising. However, it is not all that surprising when we note how much the leverages and, especially, the partial leverages vary.
By default, summclust produces a figure like Figure (ref), with its title created by the program using the name of the clustering variable, in this case “mercado”. This figure plots both leverage and partial leverage against the number of observations per cluster and also against the omit-one-cluster coefficients. These four subfigures may help to reveal the source of cluster-level heterogeneity. For this example, neither the large leverages nor the large partial leverages come exclusively from clusters with large numbers of observations or extreme omit-one-cluster coefficients.
To explore what is driving the differences in partial leverage, we create an additional scatter plot. (ref) plots partial leverage against the number of observations per cluster, with different colors and symbols depending on whether or not a given market (cluster) was treated. The figure has two interesting features. The first is that the three rather large clusters have fairly small partial leverage. The second is that the 12 clusters with the highest partial leverage are all control markets. The first result is quite surprising, since large clusters often tend to have high leverage. But (ref) makes it clear that there is, in general, no simple relationship between cluster sizes and partial leverage. The second result is not so surprising, because only 21 out of the 72 clusters are controls. Many of the control clusters presumably have high partial leverage because control clusters are relatively rare. See \hyperref[{eq:treatlev2}]{\tagform@{\ref*{eq:treatlev2}}} in (ref) in the next section for an explanation.
In this section, we discuss a number of simple examples in which it is possible to calculate our measures of leverage and influence analytically. These examples are quite revealing.
Up to this point, we have focused on one\kern 0.04167em-way clustering. However, it is also important to compute measures of leverage, partial leverage, and influence when there is clustering in two or more dimensions \citep*{CGM_2011}. In the simplest and most commonly-encountered case, where there is two\kern 0.04167em-way clustering, we recommend computing the usual one\kern 0.04167em-way measures of leverage, partial leverage, and influence for each of the two clustering dimensions. This requires calling summclust twice.
When the number of clusters in either dimension is small, or when the data are seriously unbalanced in either dimension, conventional inference based on a two\kern 0.04167em-way version of CV$_{\kern -0.08333em1}$, together with the $t(\min(G-1,H-1))$ distribution, can be seriously unreliable. \citet*{MNW_2021} therefore suggests using the usual two\kern 0.04167em-way CV$_{\kern -0.08333em1}$ estimator and applying the original WCR bootstrap to the dimension with the fewest clusters or the most unbalanced clusters. Simulation evidence suggests that this often provides more reliable inferences than the $t$ distribution, but these inferences may still be problematic.
It may also be interesting to calculate measures of leverage, partial leverage, and influence for the intersection of the two clustering dimensions, especially when the number of non-empty intersections is not large. This means calling summclust a third time. Suppose there are two clustering dimensions, with $G$ clusters in the first dimension and $H$ clusters in the second. Then the number of intersection clusters is at most $GH$\kern -0.08333em, but it can be smaller if some of the intersection clusters are empty. In order to use summclust for the intersections, it is necessary to create a new variable that uniquely identifies each of the non-empty intersection clusters. Running summclust for this case may be expensive when the number of non-empty intersections is large, especially if $k$ is also large.
It is important to remember that, when summclust is invoked three times for each of two clustering dimensions and their intersection, the CV$_{\kern -0.08333em3}$ standard error that it reports for each of the three cases is based on a different pattern of one\kern 0.04167em-way clustering. When two\kern 0.04167em-way clustering is appropriate, none of these standard errors is valid. However, what summclust reports can be used to compute an asymptotically valid variance as
Here $\hat\beta_j$ is the OLS estimate of a coefficient of interest, and the three estimated variances on the right-hand side of \hyperref[{eq:twoway}]{\tagform@{\ref*{eq:twoway}}} are the squares of the CV$_{\kern -0.08333em3}$ or CV$_{\kern -0.08333em3\rm{J}}$ standard errors reported by summclust for clustering dimension $G$, clustering dimension $H$\kern -0.08333em, and the intersection of the two clustering dimensions, respectively.
Asymptotically, the two\kern 0.04167em-way variance $\operatorname{Var}_{\rm 2W}(\hat\beta_j)$ should not be less than either of the one\kern 0.04167em-way variances. Therefore, if $\widehat\operatorname{Var}_{\rm 2W}(\hat\beta_j)$ is less than either $\widehat\operatorname{Var}_G(\hat\beta_j)$ or $\widehat\operatorname{Var}_H(\hat\beta_j)$, it makes sense to replace it by the larger of those two variance estimates. Doing this also eliminates the risk of having to take the square root of a negative number. The appropriate $t$ distribution has $\min(G-1,H-1)$ degrees of freedom if $\widehat\operatorname{Var}_{\rm 2W}(\hat\beta_j)$ is used and $G-1$ or $H-1$ degrees of freedom if it is replaced by either $\widehat\operatorname{Var}_G(\hat\beta_j)$ or $\widehat\operatorname{Var}_H(\hat\beta_j)$, respectively. We conjecture that, especially when this is done, the two\kern 0.04167em-way standard error based on either jackknife estimator will yield more conservative, and generally more reliable, inferences than the usual two\kern 0.04167em-way standard error based on CV$_{\kern -0.08333em1}$.
As we discuss in (ref), it is often invalid to partial out fixed effects when computing a jackknife CRVE. This can be particularly tricky in the case of two\kern 0.04167em-way clustering. For example, suppose there are $G$ states and $H$ years. Then it may be desirable to partial out the state fixed effects when computing $\widehat\operatorname{Var}_G(\hat\beta_j)$ but invalid to partial out the year fixed effects. Similarly, it may be desirable to partial out the year fixed effects when computing $\widehat\operatorname{Var}_H(\hat\beta_j)$ but invalid to partial out the state fixed effects. Finally, it is invalid to partial out either set of fixed effects when computing $\widehat\operatorname{Var}_{GH}(\hat\beta_j)$. The absorb option of summclust normally detects cases where partialing out is invalid and refuses to display jackknife standard errors and several other quantities.
One of the reasons for calculating leverages and partial leverages is to identify cases in which inference may be problematical. The objective of the simulation experiments in this section is to see whether the rejection frequencies for cluster-robust $t$-tests can be predicted from the features of the ${\bm{X}}$ matrix reported by summclust. There are 3000 cases, each corresponding to a particular ${\bm{X}}$ matrix. For each case, we generate 10,000 values of ${\bm{y}}$ and use them to estimate rejection frequencies for $t$-tests or bootstrap tests at the .05 level.
In the experiments, there are either 20 clusters and 2000 observations or 30 clusters and 3000 observations. The cluster sizes $N_g$ are determined by a parameter $\gamma\ge0$, as follows:
where $[\cdot]$ denotes the integer part of its argument, and $N_G = N - \sum_{j=1}^{G-1} N_g$. As $\gamma$ increases, the cluster sizes become increasingly unbalanced. The value of $\gamma$ is chosen randomly from the ${\rm U}[2,4]$ distribution, so that the cluster sizes tend to vary quite a lot. When $G=20$, the smallest cluster has between 8 and 32 observations, and the largest has between 229 and 378. When $G=30$, the smallest cluster has between 7 and 32 observations, and the largest has between 237 and 396.
There are five regressors, one of which is the test regressor, plus a constant term. The regressors equal either 0 or 1. With probability $1-p_c$, all the observations in a cluster are 0. With probability $p_c$, they randomly equal either 0 or 1, both with probability 0.5. Thus, when $p_c=1$, all variation is at the individual level, and leverage tends to be proportional to cluster sizes. As $p_c$ declines, the samples become more unbalanced. In the experiments, the values of $p_c$ are chosen to be 0.25, 0.30, 0.35, 0.40, 0.50, and 0.60, each for one\kern 0.04167em-sixth of the cases. Smaller values of $p_c$ tend to be associated with larger discrepancies between actual rejection frequencies and .05, the nominal level of the tests.
For each experiment, we obtain $12,\kern -0.08333em000$ estimated rejection frequencies. One\kern 0.04167em-quarter of these are based on CV$_{\kern -0.08333em1}$ and the $t(G-1)$ distribution, one\kern 0.04167em-quarter on CV$_{\kern -0.08333em3}$ and the $t(G-1)$ distribution, and one\kern 0.04167em-quarter on each of the WCR-C and WCR-S bootstraps. To predict these rejection frequencies, we use a generalized additive model based on smoothing splines; see \citet*[Section 7.7]{JWHT}. The base model can be written as
where $r_i$ is the rejection frequency for case $i$. Here $V_{si}$ denotes $V_s(L_{{\bullet}j})$, the scaled variance of the partial leverages $L_{gj}$ for the test regressor for case $i$, $G^*_{i0}$ denotes $G^*_j(0)$ for the test regressor for case $i$ (recall from (ref) that it is a monotonically decreasing function of the $L_{gj}$), and $f_1(\cdot)$ and $f_2(\cdot)$ are smoothing splines with five degrees of freedom. Since everything on the right-hand side of \hyperref[{eq:ssmodel}]{\tagform@{\ref*{eq:ssmodel}}} is a function of $V_{si}$, this model is simply using the $V_{si}$ to predict the $r_i$ in a potentially nonlinear way.
(ref) shows the fitted values from \hyperref[{eq:ssmodel}]{\tagform@{\ref*{eq:ssmodel}}}, which are predicted rejection frequencies, plotted against the scaled variance of the partial leverages $L_{gj}$ for four methods of inference and two sample sizes. Panel (a) shows them for $t$-tests based on both CV$_{\kern -0.08333em1}$ (solid lines) and CV$_{\kern -0.08333em3}$ (dashed lines) for $G=20$ and $G=30$, and Panel (b) shows them for WCR-C and WCR-S bootstrap tests for the same two cases. The model seems to fit quite well, at least for the asymptotic tests, as can be seen from the values of $R^2$ reported for each of the curves. It also fits well for the bootstrap tests, and in fact it has smaller residuals for them than for the asymptotic tests. The lower $R^2$ values for the bootstrap tests simply reflect the fact that there is much less variation to explain.
We can see from (ref) that $t$-tests based on CV$_{\kern -0.08333em1}$ often over-reject to an extreme degree. For the very smallest values of $V_s(L_{{\bullet}j})$, the tests tend to over-reject modestly, with predicted rejection frequencies of 0.058 for $G=20$ and 0.055 for $G=30$. However, these then rise quite rapidly and almost linearly. For $G=30$, there are four cases (out of 3000) for which $V_s(L_{{\bullet}j})>15$. These are not shown in the figure, but the approximately linear relationship continues to hold, and the fit for these extreme cases is reasonably good.
In contrast, the $t$-tests based on CV$_{\kern -0.08333em3}$ tend to under-reject for small values of $V_s(L_{{\bullet}j})$. For the very smallest values, the predicted rejection frequencies are 0.033 for $G=20$ and 0.039 for $G=30$. Although it is not obvious from the figure, the CV$_{\kern -0.08333em3}$ tests are predicted to under-reject somewhat more than half the time, because, in our experiments, most values of $V_s(L_{{\bullet}j})$ are quite small. As $V_s(L_{{\bullet}j})$ increases, rejection frequencies increase, although for $G=20$ they start to decline again once $V_s(L_{{\bullet}j})$ exceeds about 9.6. The predicted rejection frequencies never exceed 0.105 for $G=20$ and 0.118 for $G=30$. In a few cases (74 for $G=20$ and 5 for $G=30$), the matrix that is inverted in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} was singular for at least one omit-one-cluster subsample. This happened whenever one of the regressors took the same value for all observations in $G-1$ of the clusters. These cases were dropped.
Panel (b) of (ref) shows the fitted values from \hyperref[{eq:ssmodel}]{\tagform@{\ref*{eq:ssmodel}}} for WCR-C and WCR-S bootstrap $t$-tests plotted against the scaled variance of the $L_{gj}$. Notice that the scale of the vertical axis differs greatly from the one in Panel (a). All tests, especially the WCR-S ones, perform quite well for smaller values of $V_s(L_{{\bullet}j})$. Except for WCR-S with $G=30$, however, the rejection-frequency curves are not even close to being linear. This is also the only case for which the fitted values do not deviate greatly from 0.05 for large values of $V_s(L_{{\bullet}j})$. In every other case, a large value of $V_s(L_{{\bullet}j})$ tends to be associated with substantial levels of over-rejection or under-rejection.
It is natural to ask whether we can improve the fit of \hyperref[{eq:ssmodel}]{\tagform@{\ref*{eq:ssmodel}}} by adding additional explanatory variables that are not simply functions of the $V_s(L_{{\bullet}j})$. The answer is that we can. In particular, the variables $\bar a_{\rm geo}(L_{{\bullet}j})$ and $G^*_j(1)$ are often significant when they are added. However, the spline $f_1(V_{si})$ always remains highly significant, even when many other regressors are included. Thus, at least in these experiments, the scaled variance of the partial leverages, which is the square of their coefficient of variation, seems to be particularly revealing.
Based on these results, which are of course extremely dependent on the way in which the regressors are generated, it seems sensible for investigators to look at a number of different summary measures for both leverage and partial leverage. That is why summclust reports several of them. In this case, the most informative summary measure appears to be the scaled variance, defined in \hyperref[{eq:coefvar}]{\tagform@{\ref*{eq:coefvar}}}, of the partial leverage measures $L_{gj}$, defined in \hyperref[{eq:partial}]{\tagform@{\ref*{eq:partial}}}, for the regressor of interest. summclust reports the square root of this in the “Coefvar” line of the “Cluster Variability” table. In general, cluster-robust inference seems to be most reliable when the partial leverages do not vary greatly across clusters.
We have discussed a new Stata package called summclust that is designed to summarize the cluster structure of the dataset for a linear regression model with clustered disturbances. Since the key unit of observation is the cluster, it makes sense to examine measures of influence, leverage, and partial leverage at the cluster level. These are easy to compute and are conceptually very similar to the corresponding classic measures at the observation level \citep*{BKW_1980,CH_1986}. The summclust package calculates all of them and also reports a number of summary statistics.
Our measure of influence at the cluster level can provide valuable information about how empirical results depend on the data in the various clusters. Investigators should be wary if dropping one or two clusters changes the results dramatically. However, apart from such cases, the most interesting quantities that summclust calculates generally seem to be the partial leverages and measures that summarize their distribution.
It has long been known that cluster-robust inference can be unreliable when the number of clusters is small. More recent work, including \citet*{MW-JAE,MW-EJ} and \citet*{DMN_2019}, has shown that it can also be severely unreliable when cluster sizes vary a lot or when few clusters are treated in the context of difference-in-differences and other treatment models. In both of these cases, leverage and partial leverage tend to vary greatly across clusters. It therefore seems natural to use our measures of leverage and partial leverage as diagnostic tools to identify datasets and regression designs in which cluster-robust inference is likely to be challenging. Simulation results in (ref) suggest that the extent to which partial leverage varies across clusters can be particularly informative. We believe that investigators should always look at the summary statistics reported by summclust and exercise caution whenever they indicate substantial variation across clusters.
As we discuss in (ref), the computations needed for leverage and influence are very similar to the ones needed to compute cluster jackknife variance matrix estimators. The summclust package therefore computes two very similar jackknife estimators, which we refer to as CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$, almost as a byproduct of other computations. These are the same estimators that Stata can produce using the vce(jackknife,mse) and vce(jackknife) options. However, because summclust is designed explicitly for linear regression models estimated by OLS, it is faster than using these vce options. Moreover, when \texttt{summclust} is already being used to obtain cluster-level measures of influence and leverage for diagnostic purposes, the additional cost of computing the jackknife variance estimators is minimal.
When the number of clusters is reasonably large and the variation of leverage and partial leverage across clusters is small, we would expect conventional inference based on CV$_{\kern -0.08333em1}$ standard errors to perform well. If so, the CV$_{\kern -0.08333em3}$ standard errors reported by summclust should be very similar to the CV$_{\kern -0.08333em1}$ standard errors reported by one of Stata's regression commands. When this is the case, there is probably no need for investigators to worry further about the reliability of their inferences. In many cases, however, the CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em1}$ standard errors will differ noticeably. This happens for the empirical example in (ref), where there are 72 clusters but partial leverage varies a lot. In such cases, the CV$_{\kern -0.08333em3}$ standard errors are almost certain to be more conservative, and very likely to be more reliable, than the CV$_{\kern -0.08333em1}$ ones.
$P$ values and confidence intervals that are even more reliable can often be obtained by using the restricted wild cluster bootstrap, which is implemented natively with wildbootstrap in Stata 18 and in the package boottest \citep*{RMNW}. Recent versions of that package implement the WCR-S bootstrap \citep*{MNW-bootknife} in addition to the classic WCR-C bootstrap. We strongly recommend that both variants be calculated whenever the CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em1}$ standard errors disagree. When the two bootstrap $P$ values agree, as they do for the empirical example in (ref), then it is probably safe to rely on either of them. When they disagree, then neither of them may be entirely reliable, but we would be inclined to use the one given by the WCR-S bootstrap.
Up to this point, everything in this section has been based on the assumption that there is one\kern 0.04167em-way clustering with a known clustering structure. When more than one level of clustering is plausible, investigators need to choose among them, and this can be challenging; see the discussions in \citet*{MNW-guide,MNW-testing}. The measures of leverage and influence produced by summclust may be helpful in deciding at what level to cluster.
The current version of summclust is not explicitly designed to handle two\kern 0.04167em-way clustering. However, as we discuss in (ref), it can be called for each clustering dimension so as to produce two sets of diagnostic statistics. If it is called three times, once for each dimension and once for their intersection, then it can also be used to compute two\kern 0.04167em-way cluster jackknife variance matrix estimators. At present, however, little is known about the properties of these estimators.
To install the software files as they exist at the time of publication of this article, type
The command summclust can be installed from the Statistical Software Components archive by typing
or from GitHub by typing
{1pt} \addcontentsline{toc}{section}{\refname}