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.
42,528 characters · 11 sections · 0 citation commands
Unbiased estimation of the OLS covariance matrix when the errors are clustered
Within-cluster dependence presents a considerable challenge for reliable inference. Even with large data sets, a small number of clusters induces substantial finite sample bias in the estimated variance of the regression coefficients. Several options are available to mitigate this bias. Stata uses a scalar correction to the Liang and Zeger (1986) cluster-robust variance estimator, while Bell and McCaffrey (2002) develop cluster extensions of the MacKinnon and White (1985)\nocite{McKW1985} heteroskedasticity-robust variance estimators. See Cameron and Miller (2015) and MacKinnon, Nielsen, and Webb (2022)\nocite{McKOrrWe22} for recent surveys on the topic. However, with the exception of some special cases, none of these variance adjustments completely eliminates the bias.
In this paper, we develop variance estimators that are unbiased under progressively more complicated dependence structures. Our aim is to investigate whether removing the bias in the variance estimators leads to improved inference, in particular by delivering hypothesis tests with more accurate size control. The key idea underlying the unbiased variance estimator is a cluster extension of the variance estimator by Hartley, Rao, and Kiefer (1969), which is unbiased under heteroskedasicity. In its original form, this variance estimator has the drawback that it requires inverting a matrix that grows quadratically with the sample size. We show how the underlying structure of this matrix can be exploited to make the computation feasible even with large microeconometric data sets.
With a large number of clusters, test statistics based on cluster-robust variance estimators have a standard normal distribution, see for instance Hansen and Lee (2019)\nocite{HaLe19}. With a small number of clusters, the use of the normal distribution to obtain confidence intervals and critical values can lead to substantial size distortions as discussed in Cameron and Miller (2015), Section VI.D, unless the within-cluster dependence is restricted as in Ibragimov and M\"uller (2016)\nocite{IbMu10}. The use of a $t$-distribution reduces the size distortion, but this requires selecting the appropriate degrees of freedom (d.f.). For our proposed variance estimators, we derive a data-driven estimator for the d.f.\ following the approach based on an independence assumption on the errors as in Bell and McCaffrey (2002) as well as the generalization to a random-effects (RE) structure studied in Imbens and Koles{\'a}r (2016).
We focus on three dependence structures of increasing generality. First, we assume that in each cluster the errors follow the same RE structure. In this case, the covariance structure depends on two (unknown) parameters. Second, we extend this setting by allowing the RE parameters to be cluster dependent, increasing the number of parameters to two times the number of clusters. Finally, we consider a fully unrestricted setting where each cluster has an arbitrary covariance matrix. This captures for example a setting with conditional heteroskedasticity where the covariance matrix depends via an unknown functional form on a set of continuous regressors. Our approach can be readily adapted to a panel data setting.
For each of the three dependence structures, we numerically evaluate the size properties of hypothesis tests based on the unbiased variance estimators. We compare their performance with the default Stata option as well as the HC2 variance estimator by Bell and McCaffrey (2002) with d.f.\ as in Imbens and Koles{\'a}r (2016). The model we consider includes a treatment dummy and a continuous variable. For each covariance structure, we vary the number of treated clusters and consider both a balanced design, where each cluster has the same number of observations, as well as an unbalanced design.
Under the specification where the RE covariance structure is the same across clusters, we find that the corresponding unbiased variance estimator performs remarkably well. Even with only a single treated cluster, hypothesis tests provide accurate size control on both the treatment dummy and the continuous variable. When the number of observations differs between clusters, we find that the d.f.\ calculated under the more general RE assumption improve substantially over those calculated under independence assumptions. In a more general setting where the RE structure is cluster dependent, we find that using the corresponding variance estimator improves over the benchmarks particularly when the design is unbalanced. Finally, we consider a setting where there is conditional heteroskedasticity that depends on the continuous variable. The most general unbiased variance estimator continues to control size in this set-up.
After these simulations with fully artificial data we compare methods using real-life data with an artificial element added. That is, we estimate a wage equation on the basis of U.S. data, clustered by state. To the real-life data we added an artificial state-wide policy dummy variable. We study the size of an hypothesis test on the effect of this dummy variable by sampling subsets of states either at random or based on their number of observations.
The paper is organized as follows. In Section (ref) we start by deriving the general form of unbiased estimators for error covariance matrices with a linear structure. We then specify this for clusters in Section (ref). We first consider in Section (ref) a simple structure with just two parameters, one for the overall error and one for the within-cluster error. In Section (ref) we generalize this and make these parameters specific per cluster. In Section (ref) we generalize this one more step and allow all covariances within clusters to vary freely. We proceed to compare the performance of the various unbiased variance estimators, first by simulation and then through an application to real-life data. Our performance measure is the size of the $t$-test. The d.f.\ of the $t$-tests play an important role, and in Section (ref) we discuss how we set them. Section (ref) describes the set-up of the simulations, while the results are presented in Section (ref). The results for the real-life data are given in Section (ref). Section (ref) concludes.
We consider linear regression $\by=\bX\bbeta+\bepsi$, with $\bX$ exogenous of order $n\times k$. We follow the usual notation $\bM\equiv\bI_n-\bX(\bX'\bX)^{-1}\bX'$ and $\bP\equiv\bX(\bX'\bX)^{-1}\bX'$. The errors are distributed according to $\bepsi\sim(\bzero,\bSigma)$ and we consider the case where $\bSigma$ is linear in parameters, \[ \mbox{vec}\bSigma=\bD\bpi, \] with $\bpi$ of order $r\times 1$ and the design matrix $\bD$ of order $n^2\times r$. We are interested in unbiased estimation of the covariance matrix $\bV$ of the OLS estimator $\hat\bbeta$ of $\bbeta$, \[ \bV=(\bX'\bX)^{-1}\bX'\bSigma\bX(\bX'\bX)^{-1}. \] With \[ \bR'\equiv\left((\bX'\bX)^{-1}\bX'\otimes(\bX'\bX)^{-1}\bX'\right)\bD, \] we have in stacked form, which is more convenient for our analysis,
We base our estimator on a function of the residuals $\hat{\bepsi}\equiv\bM\bepsi$ that is aligned with the structure of $\bSigma$. We hence project the squared residuals on the space spanned by $\bD$, so we use $\bD(\bD'\bD)^{-1}\bD'(\hat{\bepsi}\otimes\hat{\bepsi})$, leading to the estimator
However, this estimator is biased; with $\E\left(\bD'(\hat{\bepsi}\otimes\hat{\bepsi})\right)=\bD'(\bM\otimes\bM)\bD\bpi$ there holds \[ \E(\tilde\bv)=\bR'(\bD'\bD)^{-1}[\bD'(\bM\otimes\bM)\bD]\bpi\ne\bR'\bpi=\bv. \] The bias is easily removed by replacing the term $(\bD'\bD)^{-1}$ by $[\bD'(\bM\otimes\bM)\bD]^{-1}$. For the special case of heteroskedasticity, this idea is due to Hartley, Rao, and Kiefer (1969)\nocite{HaRK69}. The adapted, unbiased estimator of $\bv$ then is
For computational purposes ((ref)) is unattractive as the matrix $\bM\otimes\bM$ is huge with large data sets. However, we show below how the simple structure of $\bM$, being the sum of the unit matrix and a matrix of low rank, can be exploited to avoid computational difficulties. A relatively common issue with unbiased estimation of variance components, see for instance Kline, Saggio, and S{\o}lvsten (2020)\nocite{KlSo21}, is that the estimator is not guaranteed to be positive definite. However, corrections that make the estimator positively biased are readily available and avoid overrejection.
Below we will consider three cases, with different design matrices $\bD$. In the third case the number of columns of $\bD$ can be very large. Then we can use an adapted version of ((ref)). Let
Then
Since \[ (\bW+\bF'\bA^{-1}\bF)\bW^{-1}\bF'=\bF'\bA^{-1}(\bA+\bF\bW^{-1}\bF') \] there holds \[ \bW^{-1}\bF'(\bA+\bF\bW^{-1}\bF')^{-1}=(\bW+\bF'\bA^{-1}\bF)^{-1}\bF'\bA^{-1}. \] Substitution in ((ref)) yields
This expression still contains the inverse of the matrix $\bA$, which has the same number of columns as $\bD$. It will appear, though, that $\bA^{-1}$ occurs only in the form $\bF'\bA^{-1}$, which appears to have a simple expression in this case.
We now turn to the cluster structure. We denote the number of clusters by $C$ and index them by $c=1,\ldots,C$. Cluster $c$ has $n_c$ observations, so $\sum_cn_c=n$. We let
Let $\biota_c$ an $n_c$-vector of ones. With a slight abuse of notation we will write $\bI_c$ for $\bI_{n_c}$ and let
The regressors for cluster $c$ are collected in $\bX_c\equiv\bG_c'\bX$ and their sum over the cluster in the row vector $\tilde{\bx}_c'\equiv\bb_c'\bX$. The $\tilde{\bx}_c'$s are collected in the $C\times k$ matrix $\tilde{\bX}\equiv\bB'\bX$. Likewise, $\hat{\bepsi}_c\equiv\bG_c'\hat{\bepsi}$ and $\tilde{\hat{\bepsi}}_c\equiv\bb_c'\hat{\bepsi}$ so $\tilde{\hat{\bepsi}}=\bB'\hat{\bepsi}$.
Below we will frequently perform matrix operations using
for conformable generic $\bA,\bB,\bC$ and $\bD$. A piece of notation that is useful in the third case that we will study is the Kronecker product with a dot on top. With $\be_c$ be the $c$th unit vector, we write \[ \textstyle\sum_c\be_c'\;\dot\otimes\;\bA_c=\left(\bA_1,\ldots,\bA_C\right) \] for matrices $\bA_1,\ldots,\bA_C$ with the same number of rows but possibly different number of columns. The use of $\dot\otimes$ is as straightforward as the use of $\otimes$.
In this section we consider three, increasingly general structures for $\bSigma$ and present the variance estimator ((ref)) for each case. The derivations are given in Appendix A.
We first consider the case where the errors are equicorrelated within clusters, so \[ \bSigma=\sigma^2\bI_n+\tau^2\bB\bB', \] with $\bB$ as given in ((ref)). Let \[ \bPsi=\left(
\right), \] with
Then
is an unbiased estimator of $\bv$.
Two remarks are in order here. The first one concerns symmetry. The $k\times k$ covariance matrix $\hat\bV$, obtained by rearranging $\hat\bv$ into a matrix, should be symmetric. The derivation of ((ref)) did not take this requirement into consideration. However, it is easy to show that $\hat\bV$ is symmetric, by employing the commutation matrix $\bK_k$, with properties $\bK_k(\bA\otimes\bB)=(\bB\otimes\bA)\bK_k$ for any $k\times k$ matrices $\bA$ and $\bB$ and $\bK_k\mbox{vec}\bC=\mbox{vec}\bC$ for any symmetric $k\times k$ matrix $\bC$. Symmetry of $\hat\bV$ is equivalent to $\bK_k\mbox{vec}\hat\bv=\mbox{vec}\hat\bv$. By using $\bK_k=\bK_k^{-1}$ this readily follows. The same holds for the other two variance estimators derived below.
The second remark concerns the role played by the regressors. When they would have been neglected in the derivation, that is, estimating $\bv$ by ((ref)) instead of by ((ref)), we would have obtained
We can then write \[ \hat{\bv}=(\bX'\bX\otimes\bX'\bX)^{-1}\left(\mbox{vec}\;\bX'\bX,\mbox{vec}\;\tilde{\bX}'\tilde{\bX}\right)(\hat{\sigma}^2,\hat{\tau}^2)' \] or
with $\hat\bSigma=\hat{\sigma}^2\bI_n+\hat{\tau}^2\bB\bB'$, where
In this form, $\hat\bSigma$ is the estimator for $\bSigma$ used by Imbens and Koles{\'a}r (2016) in their d.f.\ derivation, to be discussed below in Section (ref).
We next let $\sigma^2$ and $\tau^2$ vary over clusters, so now \[ \bSigma=\textstyle\sum_c(\sigma_c^2\bG_c\bG_c'+\tau_c^2\bb_c\bb_c'). \] Let \[ \bPhi=\left(
\right), \] with
while $\bA, \bL$ and $\bQ$ are matrices of order $C\times C$ with typical elements
Then
is the unbiased estimator for the variance $\bv$.
Also here it is interesting to consider the result when the regressors are neglected. Then \[ \bPhi=\left(
\right). \] By permuting rows and columns we can rearrange $\bPhi$ into a block-diagonal matrix with $c$th block equal to \[ \bPhi_c=\left(
\right). \] from ((ref)) it is clear that this leads to a generalization of ((ref)) to the case of cluster-specific parameters, with obvious adaptations of ((ref)) and ((ref)).
The third case we consider has errors correlate freely within clusters, in a way that differs over clusters. Thus,
where the $\bLambda_c$ are $n_c\times n_c$ matrices of parameters. With \[ \bS_c\equiv\bI_{k^2}-\bI_k\otimes\bX_c'\bX_c(\bX'\bX)^{-1}-\bX_c'\bX_c(\bX'\bX)^{-1}\otimes\bI_k \] we now obtain
as the unbiased estimator of $\bv$ for this case.
Again it is interesting to consider the version of $\hat\bv$ that neglects the regressors. Rearranged into matrix format, it appears to be
This estimator directly generalizes the White (1980)\nocite{Whit80} for cross-sections to clusters and was introduced in the context of panel data analysis by Liang and Zeger (1986), where it underlies the widely used panel-robust standard errors allowing for both heteroskedasticity and correlation over time, see e.g. Cameron and Trivedi (2005)\nocite{CaTr05}.
When the interest shifts from clustered data to panel data one might like to consider the counterpart of ((ref)) that is homogeneous over the observational units. The setting then is the panel data model with $N$ units and $T$ waves, so $n=NT$, and the covariance structure is $\bSigma=\bI_N\otimes \bLambda$, with $\bLambda$ of order $T\times T$. We discuss unbiased estimation for this case in Appendix B.
The various expressions for $\hat\bV$ or $\hat\bv$ may be of interest by themselves but their main use will be in inference on one particular regression coefficient, $\beta_\ell$, say. For large $C$, the critical values from a standard normal distribution can be used. However, in practice $C$ is often small, and using a $t$-distribution is to be preferred. For instance, Stata uses a $t(C-1)$-distribution after the command regress y x, vce(cluster clustvar).
Bell and McCaffrey (2002) proposed a refinement by making the d.f.\ in the $t$-distribution data-dependent. The idea is as follows. Let $v^2_\ell$ be the variance of the OLS estimator $\hat{\beta}_\ell$ and $\hat{v}^2_\ell$ an estimator of $v^2_\ell$. Let \[ T=\frac{\hat{\beta}_\ell}{v_\ell}/{\frac{\hat{v}_\ell}{v_l}}. \] Under normality of the regression errors, the numerator is $N(0,1)$ when $\beta_\ell=0$. Letting $\hat{v}^2_\ell$ be the usual OLS-based estimator of $v^2_\ell$, the denominator is distributed according to
leading to the $t(n-k)$-distribution for $T$. This classical result gets lost when we employ another estimator $\hat{v}^2_\ell$ than the usual one, like one of the cluster-robust estimators discussed in Section (ref). The proposal of Bell and McCaffrey (2002) is to stay close to ((ref)), by seting the d.f.\ $d_\ell$ such that \[ d_\ell\;\frac{\hat{v}^2_\ell}{v^2_\ell}\stackrel{\mbox{\tiny{app}}}{\sim}\chi^2_{d_\ell}, \] where “app” stands for “approximately” in the sense that the first two moments of $d_\ell\hat{v}^2_\ell/v^2_\ell$ match those of a $\chi^2$-distribution with $d_\ell$ d.f.\ Using unbiased estimators of the variance as derived in the preceding section proves its usefulness here since then the first moments left and right match. Letting the second moments match means $\mbox{var}(d_\ell\hat{v}^2_\ell/v^2_\ell)=2d_\ell$ or
Obviously, $d_\ell$ is not known and needs to be estimated. There are two issues with this. One is that $d_\ell$ may depend on parameters, which have to be estimated. A second issue is that evaluating $v^2_\ell$ and $\mbox{var}(\hat{v}^2_\ell)$ requires the distribution of $\bepsi$. As a practical solution to obtain a reasonable value of $\hat{d}_\ell$, Bell and McCaffrey (2002) propose to take $\bepsi\sim N(\bzero,\sigma^2\bI_n)$ as the “reference distribution.” Imbens and Koles{\'a}r (2016) suggested to take the RE model as the reference distribution, $\bepsi\sim N(\bzero,\sigma^2\bI_n+\tau^2\bB\bB)$, with $\bB$ as defined in ((ref)). We will now derive expressions for $d_\ell$ for both cases. Given our focus on unbiased estimation, we extend previous results by using an unbiased estimator for $\mbox{var}(\hat{v}^2_\ell)$ and by using an unbiased estimator of any parameter that we meet in $d_\ell$.
So, first following Bell and McCaffrey (2002), we let $\bepsi\sim N(\bzero,\sigma^2\bI_n)$. As $\hat{v}^2_\ell$ is quadratic in $\hat\bepsi$, we can write $\hat v^2_\ell=\hat\bepsi'\bA_\ell\hat\bepsi$ for some symmetric $n\times n$ matrix $\bA_\ell$ whose particular form follows from ((ref)), ((ref)) or ((ref)), depending on the case under consideration. For notational simplicity we will omit the subscript $\ell$ to $\bA$ from now on and denote $\ba\equiv\mbox{vec}\bA$, so
hence
From ((ref)), ((ref)) and ((ref)), $\bA$ readily appears to be block-diagonal, with $c$th block $\bA_c$ given by
respectively, with $\blf_\ell\equiv\be_\ell\otimes\be_\ell$ and $\br_1\equiv(r_{11},\ldots,r_{1C})'$ and likewise for $\br_2$. Since $v_\ell^2=\sigma^2\be_\ell'(\bX'\bX)^{-1}\be_\ell$, we obtain
with
Computational gains can be had by exploiting the structure of $\bA_c$. Notice that the expression for $d_\ell$ does not depend on unknown parameters since the factor $\sigma^4$ in the numerator and the denominator cancel out.
Next, following Imbens and Koles{\'a}r (2016), we let $\bSigma = \sigma^2\bI_{n} +\tau^2\bB\bB'$, with $\bB$ as defined in ((ref)). Instead of ((ref)) we now have $\mbox{var}(\hat{v}^2_\ell)=2\sigma^4\mbox{tr}\bA\bM\bSigma\bM\bA\bM\bSigma\bM$ , and ((ref)) generalizes to
Here, both numerator and denominator depend on the parameters $\sigma^4, \tau^4$ and $\sigma^2\tau^2$, which do not cancel out and hence have to be replaced by estimators. The lengthy expression in the denominator posses another complication. Both complications are addressed in Appendix C.
We take the simulation design of MacKinnon and Webb (2018)\nocite{McKWe18} as our point of departure. The data generating process includes a treatment dummy and a continuous variable. For $c=1,\ldots,C$ it is
with $\biota_c$ the intercept, $\bd_{c}$ the treatment dummy equal to 1 in clusters $1,\ldots,C_{1}$, which we will vary from $1$ to $C-1$, and $\bx_c$ the continuous regressor, whose elements are independent $N(0,1)$. The regression errors $\bepsi_c$ within cluster $c$ are normally distributed with their covariance matrix $\bSigma_{c}$ specified below. The errors are independent across clusters. We set the parameters $\alpha=\beta=\gamma=0$, the number of clusters $C=14$, and the total number of observations $n=2800$. The results below are based on 200,000 draws of (ref). We draw the continuous variable $\bx_c$ only once.
\paragraph{Error covariance matrix} To generate the data, we consider three increasingly complicated designs for the covariance matrix of the $\bepsi_{c}$.
\paragraph{Balance} An important design choice is the number of observations per cluster. We first consider a balanced design, where the number of observations per cluster is equal to $n/C=200$, and next an unbalanced design, where the number of observations depends on the cluster index according to
We set $\gamma = 2$, which implies cluster sizes ranging from 67 to 438 observations.
\paragraph{Variance estimators and reference distributions} We consider the following methods to obtain $t$-values for the OLS estimate for $\beta$ in (ref).
Notice that LZ2 does not exist when the number of (un)treated clusters is smaller than two, and that UV2($\cdot$), UV3($\cdot$) do not exist when the number of (un)treated clusters is smaller than three. We then set the size to zero.
The main results of the simulations are presented in Figure (ref), (ref) and (ref), based on data simulated with error covariance matrix as in ((ref)), ((ref)) and ((ref)), respectively. They show the size of the $t$-test for $H_0:\beta=0$, with $\beta$ the coefficient of the dummy variable in ((ref)). The number of treated clusters is on the horizontal axis. The upper panel of each figure is for the balanced case and the lower panel for the unbalanced case as described in ((ref)). Each figure shows eight curves, for STATA, LZIK, UV1(RV0), UV1(RV1), UV2(RV0), UV2(RV2), UV3(RV0), and UV3(RV3). Notice that three variances are involved: the reference variance to obtain $d_\ell$; the variance whose unbiased estimator was used; and the variance used in the simulation. For clarity, Table (ref) summarizes.
The most relevant curves in all three figures are the ones labeled UV1(RV1) in Figure (ref), upper and lower panels. The homogeneous RE design can be considered the more or less generic case in the clustered-error literature and, as is apparent from Table (ref), this particular curve is maximally based on this design as it underlies the data generation SV1, the variance estimator UV1, and $d_\ell$ based on RV1.
\paragraph{SV1} Inspecting Figure (ref) we see, for the balanced design in the upper panel, excellent size control for UV1($\cdot$). This holds even when there is only a single treated cluster. It does not appear to matter whether the d.f.\ are calculated under the more restrictive i.i.d.\ assumption, UV1(RV0), or the RE structure, UV1(RV1). By contrast, UV2($\cdot$), UV3($\cdot$) and LZIK are slightly conservative when we have a small or large number of treated clusters. The STATA variance estimator performs quite poorly, especially when the number of treated clusters is small or large.
Moving to the unbalanced set-up in the lower panel of Figure (ref), we see that UV1(RV0) no longer provides accurate size control. However, UV1(RV1), the most relevant case as argued above, still exhibits excellent performance. The additional computational complexity of this approach appears to pay off. We also see that, unlike in the balanced case, the results for UV2(RV0) and UV3(RV0) differ from the benchmark variance estimator LZIK. The unbiased variance estimators are more conservative for a small number of treated clusters, while becoming slightly oversized for 9-11 treated clusters. UV2(RV1) and UV3(RV1) are again very close to LZIK. The STATA variance estimator again is found not to accurately control size.
\paragraph{SV2} In Figure (ref) we show the size for $t$-tests based on the various variances estimators under the restricted heterogeneous design where each cluster has its own variance and covariance parameter. This set-up is more general than the homogeneous design in which each cluster has the same variance and covariance parameter. As expected, the performance of UV1(RV0) and UV1(RV1) somewhat deteriorates in this set-up, with size slightly below 0.10 for the case of a single treated cluster and balanced design. The same is observed for an unbalanced design, with the size obtained under UV1(RV1) being just over 0.10.
For UV2 and UV3, under both d.f., and LZIK, we see the test slightly overrejects for a small number of treated clusters. When the number of treated clusters increases, the tests become progressively more conservative. Again, a difference emerges between UV2, UV3 and LZIK in the unbalanced case presented in the lower panel of Figure (ref). Here, size control is more accurate for UV2 and UV3 compared to LZIK. Especially UV1(RV0) and UV2(RV0) perform well in this set-up, providing accurate size control up to roughly 8 treated clusters. With more treated clusters, they tend to be conservative, although not as much as LZIK.
\paragraph{SV3} The results for the unrestricted heterogeneous design are nearly identical to those in the homogenous design for the STATA variance, LZIK and UV2 and UV3 under both d.f.\ corrections. For UV1, we find reasonable performance when clusters are balanced. When the clusters are unbalanced, UV1(RV0) becomes oversized for a small number of treated clusters, and undersized when the number of treated clusters is large. The more general d.f.\ correction in UV1(RV1) partly corrects these size distortions.
So far for the test on $\beta$, the coefficient of the cluster-specific dummy variable. We can be much more concise as to $\gamma$, the coefficient of the continuous variable. For SV1 and SV2 the size control is almost perfect. This no longer holds for SV3, where the size is still almost perfect for STATA, LZIK, UV3($\cdot$) but appears to be double the nominal size for $UV1(\cdot$) and UV2($\cdot$); the latter methods are apparently sensitive when the data are generated according to more general scheme SV3.
\paragraph{Degrees of freedom in SV1} Given the notable differences in performance when using degrees of freedom based on RV0 or RV1, we analyze the degrees of freedom under SV1 in Figure (ref). For a balanced design, we see that the degrees of freedom for UV1 are equal to $C-2$. Donald and Lang (2007)\nocite{DoLa07} show that if the design is balanced and if all regressors are invariant within clusters, the $t$-statistic is $t(C-k)$ distributed, where $k$ is the number of regressors in the model. We can expect the same result to apply here since the continuous variable is uncorrelated with the treatment dummy.
Under a balanced design, the degrees of freedom for the other methods are nearly identical. They are low when the number of treated clusters is low and increase to their maximum value when half of the clusters is treated. This maximum appears to coincide numerically with $C-k$ as well.
When the design is unbalanced, we see a strong deviation from the degrees of freedom under RV0 to those under RV1. This is especially true for UV1 and a small number of treated clusters. For the remaining variance estimators, we see that under RV0 the degrees of freedom are asymmetric in the number of treated clusters, while those under RV1 are symmetric.
To analyze the performance of the unbiased variance estimators in an empirical setting, we consider an application similar to that in Cameron and Miller (2015)\nocite{CaMi15}. We use the Current Population Survey (CPS) 2012 data set that can be obtained from {\tt https://cps.ipums.org/cps/}. The data consist of 51 clusters: the fifty American states and the District of Columbia. The number of observations in each cluster varies from 519 (Montana) to 5866 (California).
For observation $h$ in cluster $i=1,\dots,C$, we define the model
Here, $\mbox{\tt policy}$ is a fake policy variable that is randomly assigned to $C_{1}=1,\ldots,C-1$ sampled clusters and constant within each cluster. Since the policy variable is fake, we expect 5% rejections across the replications when we test the hypothesis ${\mbox{H}}_0:\beta_4=0$ at the 5% level.
In line with the simulations in the previous section, we sample a subset of $C=14$ clusters from the 51 available clusters. We consider two different ways of sampling this subset. In the first, we randomly sample clusters with replacement. To test the methods in an unbalanced set-up, we also consider using the $3$ states with the most observations and the $11$ states with the fewest observations. To preserve the relative share of observations in each cluster, we randomly sample with replacement 20% of the observations within each sampled cluster.
Figures (ref)--(ref) show the empirical size (upper panel) and the degrees of freedom (lower panel) averaged over 10,000 replications for the four different designs. The $x$-axis again depicts the number of treated clusters.
In line with the Monte Carlo results from the previous section, we see that the Stata variance estimator with $C-1$ degrees of freedom is severely oversized. In contrast, we find remarkably good size control for UV1(RV1) across the designs. The degrees of freedom drop considerably when moving from RV0 to RV1. This shows that the use of RV1 is of empirical relevance, especially in the settings with higher imbalance and a small number of (un)treated clusters. The LZIK variance estimator also performs well, although it is oversized in the highly unbalanced “3-11" setting. There the unbiased variance matrix estimators control size more accurately.
The point of departure in this paper has been to drive unbiased estimators of the covariance matrix of the OLS estimator when the data are clustered. We considered three cases, the leading one being the RE model. This led to our main research question, which is to assess the performance of these estimators in the $t$-test for a particular regression coefficient, both among each other and vis-\`a-vis two oft-used alternatives.
We addressed this question by simulation, in a regression model with a two regressors, one being continuous and distributed equally in all clusters, while the other regressor represented a cluster-specific treatment dummy. The main finding of the simulation study was the excellent behavior of the $t$-test based on the unbiased estimator for the RE model, for the case that the data actually have been generated according to this model and the degrees of freedom have been based on it. So the three variances that play a role are aligned. This result holds for the coefficient of the cluster-specific dummy variable; there is hardly a noticeable difference in performance for between the other variance estimators underlying the $t$-test.
A next step is to see if this excellent behavior also shows up in the case where the three variances are still aligned but now pertain to the more flexible RE model where the two error-components parameters differ over clusters. While by itself this is eminently doable, the question arises to test this heterogeneous RE structure against the homogeneous one. An obvious starting point is the score test context proposed by Breusch and Pagan (1980)\nocite{BrPa80}. Deriving the relevant expression is straightforward but deriving the (limiting) distribution of the test statistic is not since the number of parameters grows with the number of clusters.
The results in the paper on the quality of unbiased estimators in the $t$-test is based on simulation only. We are not aware of any theory that might help giving these results a theoretical basis. There is certainly a research challenge here.