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.
108,373 characters · 17 sections · 75 citation commands
Nonparametric Regression under Cluster Sampling
\address{Department of Economics, University of Wisconsin, Madison. 1180 Observatory Drive, Madison, WI 53706-1393, USA.} \email{[email removed]}
Nonparametric regression is widely used in economics for its flexibility. Typically, data are assumed to be independently and identically distributed; however, in reality, observations may exhibit dependence within a group structure called a cluster. Examples of clusters are classrooms, schools, families, hospitals, firms, industries, villages, regions, and so on. The cluster sampling framework assumes independence between observations from different clusters but allows dependence within each cluster.
The previous literature on nonparametric regression under cluster sampling assumes a bounded and homogeneous number of observations per cluster. This assumption may not hold in real data due to heterogeneous cluster sizes. To fill this gap, this paper studies nonparametric kernel regressions that accommodate heterogeneous cluster sizes, including those that grow to infinity asymptotically. Our approach is general, allowing for both bounded and growing clusters simultaneously, and includes cluster-level regressors.
We develop a comprehensive asymptotic theory for nonparametric density estimation, Nadaraya-Watson kernel regression, and local linear estimation. Our results on asymptotic conditional bias and variance, uniform consistency, and asymptotic normality enable us to propose valid methods for bandwidth selection and inference.
For clusters of growing sizes, the asymptotic variance contains a novel term for within-cluster dependence, which does not appear under the assumption of bounded cluster sizes. This term becomes significant due to the potential for a cluster to contain a growing number of observations within a local neighborhood, making cluster dependence non-negligible asymptotically. We propose consistent estimators of the asymptotic variance that account for cluster dependence and validate its importance through simulation. Our cluster-robust confidence interval achieves improved coverage ratios, while conventional confidence intervals could suffer from under-coverage in our simulated datasets.
Nonparametric regression, while significant on its own, also serves as an intermediate tool for other estimators, such as regression discontinuity design, nonparametric auction estimation, and semiparametric models under cluster sampling. Our results could extend to these areas as well.
There is a substantial body of literature on cluster sampling in econometrics. C. hansen2007asymptotic provides an asymptotic theory for parametric regression with homogeneous cluster sizes. djogbenou2019asymptotic and B. hansen2019asymptotic extend this theory to heterogeneous cluster sizes. bugni2022inference considers heterogeneous and random cluster sizes for cluster-level randomized experiments. For further literature on parametric models under cluster sampling, the reader can refer to cameron2015practitioner and mackinnon2022cluster.
Conversely, the theory on nonparametric regression under cluster dependence, even with homogeneous cluster sizes, is limited. lin2000nonparametric and wang2003marginal examine local polynomial and local linear regressions, assuming fixed and homogeneous cluster sizes and focusing primarily on asymptotic efficiency. bhattacharya2005asymptotic offers an asymptotic theory for local constant estimators under multi-stage samples, analogous to cluster sampling. When the number of first-stage strata is set to one, his setup becomes a standard cluster sampling with fixed and homogeneous cluster sizes. He puts a similar structure on error terms as this paper, but the fixed cluster sizes render the term reflecting within-cluster dependence asymptotically negligible. For the regression discontinuity literature, bartalotti2017regression has derived asymptotic theories for local polynomial regression under bounded and homogeneous cluster sizes.
menzel2024transfer proposes a method for estimating nonparametric regressions in the presence of cluster dependence, aiming to extrapolate treatment effects across clusters. He considers independent but not identical observations between clusters, with a fixed number of clusters exhibiting uniformly growing size. Our approach differs by incorporating general dependence within a cluster and allowing for both bounded and growing cluster sizes simultaneously, leading to distinct asymptotic results and theories.
To the best of our knowledge, there is no literature on nonparametric models with growing and heterogeneous size clusters except for hu2024some, which became available online after the working paper version of our paper was posted. While both our paper and theirs accommodate flexible cluster sizes, their work concentrates on series regression and does not include any theory on inference, model selection, and kernel regression. Their Assumption 2(i) excludes cluster-level covariates. Our paper adopts the same cluster size framework as djogbenou2019asymptotic and hansen2019asymptotic. The presence of clusters with growing sizes complicates the proofs for asymptotic theories, as cluster dependence becomes non-negligible. Consequently, this paper introduces new technical results for nonparametric regressions under cluster sampling, notably developing Bernstein's inequality for cluster sampling to demonstrate uniform consistency. These novel contributions are believed to offer valuable theoretical tools for future research.
This research also sheds new light on the literature regarding nonparametric regressions with dependence. Following the foundational work on i.i.d. datasets (e.g., stone1982optimal, fan1992design, ruppert1994multivariate), the results have been extended to time series (robinson1983nonparametric, hansen2008uniform, kristensen2009uniform, vogt2012nonparametric, vogt2020multiscale) and spatial datasets (robinson2011asymptotic, lee2016series), as well as to the cluster dependence framework discussed above.
The remainder of this paper is organized as follows: Section $\text{\ref{sec:Cluster-sampling}}$ introduces the cluster sampling framework under consideration. Sections $\text{\ref{sec:Nonparametric-density-estimation}}$-$\text{\ref{sec:Local-linear-estimator}}$ discuss asymptotic theories for nonparametric density estimators, Nadaraya-Watson estimators, and local linear estimators, respectively. Section $\text{\ref{sec:Uniform-convergence}}$ demonstrates uniform convergence of these estimators. Section $\text{\ref{sec:Bandwidth-selection}}$ provides guidelines for selecting bandwidth in nonparametric regressions. Section $\text{\ref{sec:Cluster-robust-variance-estimati}}$ addresses cluster-robust inference. Section $\text{\ref{sec:Monte-Carlo-simulation}}$ presents Monte Carlo simulations for bandwidth selections and inference. Section $\text{\ref{sec:Empirical}}$ illustrates our methods with an application in development economics using a dataset by alatas2012targeting. The paper concludes with Section $\text{\ref{sec:Conclusion}}$. All proofs, technical lemmas, technical discussions, and additional simulation results are included in the Appendix.
The researcher observes $\left(Y_{i},X_{i}\right)\in\mathbb{R}\times\mathbb{R}^{d}$ for $i=1,\ldots,n$, with cluster sizes given by $n_{g}\in\{1,2,\cdots\}$ for $g=1,\ldots,G$. Here, $Y_{i}$ represents a dependent variable, and regressors $X_{i}$ are continuous random variables with the Lebesgue density $f(x)$. Assume that each observation can be grouped into one cluster.\footnote{Formally, we assume that for any $i$, we know a function $g(i)\in\{1,\cdots,G\}$.} Thus, the total number of observations is $n=\sum_{g=1}^{G}n_{g}$. To explicitly represent the cluster structure, we also use the notation $\left(Y_{gj},X_{gj}\right)$ for $g=1,\ldots,G$ and $j=1,\ldots,n_{g}$. We treat cluster size $n_{g}$ as nonrandom and possibly heterogeneous across clusters. We assume that observations belonging to different clusters are mutually independent but permit general dependence within the same cluster. We decompose $X_{gj}$ into $X_{gj}=\left(X_{gj}^{(\mathrm{ind})\top},X_{g}^{(\mathrm{cls})\top}\right)^{\top}\in\mathbb{R}^{d}$ where $X_{gj}^{(\mathrm{ind})}\in\mathbb{R}^{d_{\mathrm{ind}}}$ represents individual-level regressors and $X_{g}^{(\mathrm{cls})}\in\mathbb{R}^{d_{\mathrm{cls}}}$ represents cluster-level regressors. We assume that the regressors contain at least one individual-level regressors, $d_{\mathrm{ind}}\geq1$. By construction, $d=d_{\mathrm{ind}}+d_{\mathrm{cls}}$ holds.
We denote $\mathbf{X}_{g}=\left(X_{g1},\dots,X_{gn_{g}}\right)$ and aim to estimate the nonparametric regression model:
We also assume
The model specified through $\text{\eqref{eq:estimand}}$-$\text{\eqref{eq:cond_cov_model}}$ exhibits greater flexibility than initially apparent. The constraint imposed by $\text{\eqref{eq:cond_var_model}}$ is that the conditional variance of the error term for an individual is dependent only on the individual's own regressors, both at the individual and cluster levels. Additionally, $\text{\eqref{eq:cond_cov_model}}$ states that the conditional covariance of the error terms between any two individuals within the same cluster is a function only of their individual-level regressors and shared cluster-level regressors. This framework accommodates the inclusion of cluster random effects in $e_{gj}$ and allows for the dependence of regressors within clusters. In economic applications, cluster dependencies often arise from strategic interactions within the cluster or from cluster-level unobserved shocks, including measurement errors.
In this section, we show the consistency of nonparametric density estimators. In this paper, we will use kernel functions satisfying the following definitions.
The kernel density estimator for $f(x)$ is:
where $h>0$ is a bandwidth.
For the sake of simplicity, our discussion will focus on scenarios where a single bandwidth is used for all components of $X$. However, our theory can be generalized to accommodate multivariate bandwidths by substituting $h$ with a bandwidth matrix, as discussed by ruppert1994multivariate.
In this section, we derive an asymptotic theory for the Nadaraya-Watson estimator (a.k.a. local constant estimator) for estimating the conditional expectation $\mathbb{E}\left[Y_{gj}\mid X_{gj}=x\right]$. The estimator is:
We use the following assumption to derive the asymptotic variance.
In the special case of $\lambda=0$, the asymptotic conditional variance is equivalent to the i.i.d. case. A sufficient condition for $\lambda=0$ is $\left(\max_{g\leq G}n_{g}\right)h^{d_{\mathrm{ind}}}=o(1)$. In a finite sample, it is more precise to consider $\lambda>0$. The sign of the second term of $\text{\eqref{eq:asy_cond_var}}$ depends on the sign of $\sigma\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$. In economic applications, it usually takes a positive value, indicating positive conditional covariance of error terms within clusters. Neglecting this term will lead to under-coverage in empirical applications.
The asymptotic distribution has the same bias and the same convergence rate as in the i.i.d. case. The asymptotic variance is a scaled value of the primal terms of asymptotic conditional variance that include the conditional covariance term due to the cluster dependence. Our simulation in Section $\text{\ref{sec:Monte-Carlo-simulation}}$ shows the importance of considering this term in inference.
The asymptotic variance in the previous literature with bounded cluster sizes (e.g., bhattacharya2005asymptotic) has only the first term of $\text{\eqref{eq:nw_asydist}}$. Under bounded cluster sizes, cluster dependence is asymptotically negligible since an observation in the $g$-th cluster has a negligible number of observations belonging to the same cluster around the local neighborhood. On the other hand, under growing cluster sizes $n_{g}\rightarrow\infty$, the observation could have a non-negligible number of neighboring observations belonging to the same cluster. Thus, the conditional covariance of error terms matters in our general setup.
Theorem $\text{\ref{thm:nw_asy_dist}}$ is an asymptotic result that holds pointwise-in-the-underlying distribution. Consequently, the asymptotic bias pertains to a specific data-generating process, which influences how it should be used for bandwidth selection (see Remark $\text{\ref{rem:uniformity}}$ for derails). Also, Theorem $\text{\ref{thm:nw_asy_dist}}$ does not provide guidance on how the asymptotic bias should be addressed in inference. Remark $\text{\ref{rem:bias_handling}}$ discusses a practical approach for handling this issue in inference.
In this section, we consider the local linear estimator
where
and $K_{h}\left(\cdot\right)=\frac{1}{h^{d}}K\left(\frac{\cdot}{h}\right)$. We will assume an additional condition for the simplicity of proofs.
We can establish similar asymptotic theories for local linear estimators as we derived for Nadaraya-Watson estimators. As in the i.i.d. case, the asymptotic bias of a local linear estimator does not include the term of first-order derivatives.
If we impose further assumptions, our pointwise consistency result can be strengthened to uniform consistency. Before proving uniform consistency for nonparametric estimators, we will show uniform consistency for the generic function
to its expectation, where $X_{gj}\in\mathbb{R}^{d}$ and $W_{gj}\in\mathbb{R}$.
We assume the cluster samples $\left\{ W_{gj},X_{gj}\right\} $ satisfy the following assumptions.
We also require a further assumption on the kernel function.
The proof for Theorem $\text{\ref{thm:psi_unifconv}}$ relies on the following cluster sampling version of Bernstein's inequality, which could be of independent interest.
Based on Theorem $\text{\ref{thm:psi_unifconv}},$ we will show the uniform consistency of the nonparametric density estimator and nonparametric regressions. It requires the following conditions, including uniform smoothness.
The range $\left\{ x:\left\Vert x\right\Vert \leq c_{n}\right\} $ expands slowly to $\mathbb{R}^{d}$ since our condition $\text{\eqref{eq:cn}}$ can cover a sequence $\left\{ c_{n}\right\} $ such that $c_{n}\rightarrow\infty$ slowly as $n\rightarrow\infty$. This expansion is useful to establish asymptotic theories for semiparametric estimation with a nonparametric kernel estimator in the first-stage.
Suppose that $c_{n}=c$ (constant) and $\delta_{n}$ is far away zero. Then, the uniform convergence rate for kernel regressions is $a_{n}+h^{2}=\left(\log n/(nh^{d})\right)^{1/2}+h^{2}$. By choosing $h=\left(\log n/n\right)^{1/(d+4)}$, the optimal rate $\ensuremath{\left(\log n/n\right)^{2/(d+4)}}$ is attained. This convergence rate is equivalent to stone1982optimal's optimal rate in the i.i.d. case.
In this section, we provide guidelines for selecting bandwidth in nonparametric regressions. We suggest three types of methods: the asymptotic integrated mean squared error (AIMSE) optimal bandwidth selection, the cluster-robust rule-of-thumb, and the cluster-robust cross-validation.
Let $B_{*}(x)=B_{\mathrm{nw}}(x)$ or $B_{\mathrm{LL}}(x)$. The asymptotic integrated mean squared error of the estimator $\widehat{m}_{*}\left(x\right)$ is
where $w(x)$ is some integrable weight function which ensures that $\overline{B}\equiv\int_{\mathbb{R}^{d}}B_{*}(x)^{2}f(x)w(x)\mathrm{d}x$, $\overline{\sigma}^{2}\equiv\int_{\mathbb{R}^{d}}\sigma^{2}(x)w(x)\mathrm{d}x$, and \[ \overline{\sigma}_{\mathrm{cls}}\equiv\int_{\mathbb{R}^{d}}\frac{f_{2}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)\sigma\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)}{f(x)}w(x)\mathrm{d}x \] are finite. We define
as an objective function for bandwidth selection since the third term in $\text{\eqref{eq:aimse_derivation}}$ does not depend on $h$ and the fourth term in $\text{\eqref{eq:aimse_derivation}}$ is asymptotically negligible.
Our asymptotic theorems rely on the assumption $\left(\max_{g\leq G}n_{g}\right)h^{d_{\mathrm{ind}}}=O(1)$.\\ When $\left(\max_{g\leq G}n_{g}\right)n^{-d_{\mathrm{ind}}/(d+4)}\rightarrow\infty$, the AIMSE-optimal $h_{0}$ does not satisfy this order. In this case, the AIMSE-optimal bandwidth does not make sense since the AIMSE criterion itself relies on the assumption $\left(\max_{g\leq G}n_{g}\right)h^{d_{\mathrm{ind}}}=O(1)$. Thus, when the largest cluster size is large compared to the sample size $n$, we recommend using the cross-validation criterion (see Section $\text{\ref{subsec:Cross-validation}}$).
In practice, it is not easy to compute the AIMSE-optimal bandwidth since $\text{\eqref{eq:h0}}$ contains unknown parameters. As suggested by fan1996local (fan1996local, Section 4.2) for the i.i.d. case, we provide a cluster-robust Rule-of-Thumb (CR-ROT) bandwidth choice for a one-dimensional individual-level regressor $x\in\mathbb{R}$. This bandwidth could be a crude estimator of the AIMSE-optimal bandwidth, but the primary purpose of it is to give a guess of the bandwidth requiring little computational effort. Let
be a fitted 4th-order global polynomial regression leaving out the $g$-th cluster. Given this parametric model and a user-specified integrable weight function $w(x)$, the CR-ROT bandwidth is calculated by
where
and $\check{e}_{gj}=Y_{gj}-\check{m}_{-g}\left(X_{gj}\right)$. In words, $\check{B}$ and $\check{\sigma}^{2}$ are computed by the parametric model $\text{\eqref{eq:global-g}}$ and the homoskedastic standard error assumption for local linear estimators. For Nadaraya-Watson estimators, we also assume that $X$ has a uniform distribution for simplicity. Then, we have $f^{\prime}(x)=0$ and can compute $\bar{B}$ as for local linear estimators by $B_{\text{nw}}(x)=B_{\text{LL}}(x)$. A common choice of $w(x)$ is an indicator function of some interval.
Equation $\text{\eqref{eq:CR-ROT}}$ is different from the standard Rule-of-Thumb (ROT) bandwidth choice by fan1996local since it uses $\check{m}_{-g}(x)$ instead of $\check{m}(x)$, which is estimated by the full sample. We use $\check{m}_{-g}(x)$ to eliminate dependence between the estimator $\check{m}_{-g}(\cdot)$ and $\left(Y_{gj},X_{gj}\right)$. This modification should provide a better estimation of out-of-sample prediction error.
A heuristic cross-validation function for clustered sampling is
where $\widetilde{e}_{gj}=Y_{gj}-\widetilde{m}_{-g}\left(X_{gj},h\right)$, and $\widetilde{m}_{-g}\left(X_{gj},h\right)$ is the leave-one-cluster-out nonparametric estimator computed with bandwidth $h$ and without cluster $g$. For example, hansen2022econometrics (hansen2022econometrics, Section 19.20) suggests this form of cross-validation, but he does not provide any theoretical guarantees. For Nadaraya-Watson estimators, the leave-one-cluster-out nonparametric estimator is defined by
Similarly, for local linear estimators, the leave-one-cluster-out nonparametric estimator is defined by
where
$\mathbf{X}_{x,-g}$ and $\mathbf{W}_{x,-g}$ are defined by the same way as $\mathbf{X}_{x}$ and $\mathbf{W}_{x}$, but without using the variables in the $g$-th cluster. We will show that this cross-validation criterion works appropriately.
Since $\overline{\sigma}_{w}^{2}$ does not depend on $h$, minimizing $\mathbb{E}\left[\mathrm{CV}(h)\right]$ on $h$ is equivalent to minimizing $\operatorname{IMSE}_{G-1}(h)$, which is a sum of the expected mean squared errors weighted by cluster sizes. Thus, this theorem justifies the use of the leave-one-cluster-out cross-validation. We can choose the bandwidth by minimizing a cluster-robust cross-validation function $\mathrm{CV}(h)$ over some finite grid points $H=[h_{1},\cdots,h_{J}]$,
Note that the decomposition theorem holds for finite samples and does not rely on assumptions such as $\left(\max_{g\leq G}n_{g}\right)h^{d_{\mathrm{ind}}}=O(1)$.
Since the asymptotic variance of $\text{\eqref{eq:nw_asydist}}$ contains the joint density $f\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$, the conditional variance $\sigma^{2}(x)$, and the conditional covariance $\sigma\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$, we need to estimate each of them for inference. Alternatively, calonico2019nprobust and hansen2022econometrics (hansen2022econometrics, Section 19.20) propose to use a finite sample conditional variance of $\widehat{m}\left(X_{gj}\right)$ with estimated error terms as an estimator of the asymptotic variance. To the best of our knowledge, there is no theoretical guarantee of their methods, and this paper is the first research providing asymptotic theories of inference for nonparametric regressions under general cluster sizes.
For the joint density estimation, we propose to use
where $b$ is a bandwidth and $N=\sum_{g:n_{g}\geq2}n_{g}(n_{g}-1)/2$.
The expression $\text{\eqref{eq:joint_density_estimator}}$ can be interpreted as a standard nonparametric density estimator. We estimate the density using $\left(2d_{\mathrm{ind}}+d_{\mathrm{cls}}\right)$-dimensional regressors $\left(X_{gj}^{(\mathrm{ind})\top},X_{g\ell}^{(\mathrm{ind})\top},X_{g}^{\mathrm{(cls)}\top}\right)^{\top}$, thus we have $b^{2d_{\mathrm{ind}}+d_{\mathrm{cls}}}$ in the denominator in $\text{\eqref{eq:joint_density_estimator}}$. For clusters larger than $2$ (i.e., $n_{g}\geq2$), there are $\sum_{1\leq j<\ell\leq n_{g}}1=n_{g}(n_{g}-1)/2$ possible combinations of $X_{gj}^{(\mathrm{ind})}$ and $X_{g\ell}^{(\mathrm{ind})}$. Each cluster has a $n_{g}(n_{g}-1)/2$ effective size observations, and we have the $N=\sum_{g:n_{g}\geq2}n_{g}(n_{g}-1)/2$ effective size sample in total. In these senses, $\text{\eqref{eq:joint_density_estimator}}$ is a standard nonparametric density estimator for $\left(2d_{\mathrm{ind}}+d_{\mathrm{cls}}\right)$-dimensional regressors and $n_{g}(n_{g}-1)/2$ size clusters.
We make the following assumptions to estimate the joint density consistently.
Next, we will consider conditional variance and covariance estimation. We only provide Nadaraya-Watson type estimators, but they can be easily extended to local linear type ones. Since the goal here is to estimate $\sigma^{2}(x)$, we can estimate it as we did for $m(x)$. The infeasible Nadaraya-Watson estimator is \[ \widehat{\sigma}_{\mathrm{nw}}^{2*}\left(x\right)=\frac{\sum_{g=1}^{G}\sum_{j=1}^{n_{g}}K\left(\frac{X_{gj}-x}{h}\right)e_{gj}^{2}}{\sum_{g=1}^{G}\sum_{j=1}^{n_{g}}K\left(\frac{X_{gj}-x}{h}\right)}, \] This estimator is infeasible because $e_{gj}$ is unknown. We can replace it by $\widehat{e}_{gj}=Y_{gj}-\widehat{m}_{*}\left(X_{gj}\right)$ with $\widehat{m}_{*}(x)=\widehat{m}_{\mathrm{nw}}(x)$ or $\widehat{m}_{\mathrm{LL}}(x)$. The feasible variance estimator of the conditional variance is
The following theorem shows $\widehat{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ is a consistent estimator.
Similar to the joint density estimator, we can construct a Nadaraya-Watson type estimator for the conditional covariance using $\left(2d_{\mathrm{ind}}+d_{\mathrm{cls}}\right)$-dimensional regressors $\left(X_{gj}^{(\mathrm{ind})\top},X_{g\ell}^{(\mathrm{ind})\top},X_{g}^{\mathrm{(cls)}\top}\right)^{\top}$:
Because $e_{gj}$ is unknown, it is infeasible as $\widehat{\sigma}_{\mathrm{nw}}^{2*}\left(x\right)$. The feasible version of $\widehat{\sigma}_{\mathrm{nw}}^{2*}$ is estimated by replacing $e_{gj}$ with $\widehat{e}_{gj}=Y_{gj}-\widehat{m}_{*}\left(X_{gj}\right)$,
Corollary $\text{\ref{cor:nw_asydist_condvar}}$ suggests to use
as a standard error. The estimator $\widehat{\lambda}R_{k}^{d_{\mathrm{cls}}}\widehat{f}_{2}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)\widehat{\sigma}_{\mathrm{nw}}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)/\left(\widehat{f}(x)\right)^{2}$ could be too difficult to estimate in practice for the following two main reasons. First, it contains $\widehat{f}_{2}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$ and $\widehat{\sigma}_{\mathrm{nw}}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$, which put most kernel weights for observations that $X_{gj}^{(\mathrm{ind})}$ and $X_{g\ell}^{(\mathrm{ind})}$ are both in the neighborhood of $x^{\mathrm{(ind)}}$. In a finite sample, such observations could be rarely observed, and these estimators could be imprecise. Second, it contains a density ratio $\widehat{f}_{2}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)/\widehat{f}(x)$, which is difficult to estimate, especially nonparametrically.
To overcome these difficulties, we provide a parametric compromise under additional assumptions. We assume that $f_{2}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$ follows a multivariate normal distribution, $x^{\mathrm{(ind)}}$ and $x^{\mathrm{(cls)}}$ are independent or there are no cluster-level regressors (see also Remark $\text{\ref{rem:simple_cov}}$), and the conditional covariance is homoskedastic. Then, we can simplify
where $\check{e}_{gj}=Y_{gj}-\check{m}_{-g}\left(X_{gj}\right)$, $\check{m}_{-g}\left(x\right)$ is estimated by the global polynomial regression as $\text{\eqref{eq:global-g}}$, and $p\left(x_{1}\mid x_{2},\mu,\Sigma\right)$ is a conditional density function of $x_{1}$ given $x_{2}$ with the joint distribution $\left(x_{1}^{\top},x_{2}^{\top}\right)^{\top}\sim\mathrm{N}\left(\mu,\Sigma\right)$. We can estimate $\widehat{\mu}=\left(\widehat{\mu}_{1}^{\top},\widehat{\mu}_{1}^{\top}\right)^{\top}$ and $\widehat{\Sigma}=\left(
\right)^{\top}$ easily by using sample moments. Note that the expectation $\widehat{\mu}_{1}$ and the variance matrix $\widehat{\Sigma}_{11}$ are the same for $x_{1}$ and $x_{2}$ since we initially assumed identical marginal densities in Assumption $(ref)$.
In practice, we can estimate $\sigma^{2}\left(x\right)$ and $\sigma\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$ by using clustered-level jackknife estimators. hansen2022jackknife shows that for parametric linear regressions, clustered-level jackknife variance estimators are better than conventional variance estimators with respect to the worst-case bias. Clustered-level jackknife variance estimators $\widetilde{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ and $\widetilde{\sigma}_{\mathrm{nw}}\left(x^{\mathrm{(ind)}},x^{\mathrm{(ind)}};x^{\mathrm{(cls)}}\right)$ are estimated by replacing $e_{gj}$ with $\widetilde{e}_{gj}=Y_{gj}-\widetilde{m}_{-g}\left(X_{gj}\right)$, where $\widetilde{m}_{-g}\left(\cdot\right)$ is a nonparametric estimator estimated leaving out the $g$-th cluster observations. In the simulation section, we will compare coverage ratios of confidence intervals constructed by the conventional standard error $\widehat{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ and the cluster-robust standard error $\widetilde{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$.
In this section, we will check the validity of bandwidth selections and confidence intervals in simulated datasets under cluster sampling. For both simulation studies, we consider the following setup. We fix the number of clusters $G=100$ and cluster sizes $n_{g}=20$ for $g=1,\dots G-1$. To evaluate the effect of the largest cluster size, we try two cluster sizes $n_{G}\in\left\{ 20,100\right\} $ for cluster $g=G$. Thus, we try two scenarios with $\left(\max_{g\leq G}n_{g}\right)/n\approx\left\{ 0.02,0.09\right\} $, also corresponding to homogeneous or heterogeneous size clusters. We generated 2000 datasets for replication. For the data-generating process, the following two models are considered.
Setup 1 (homoskedastic errors):
\[ Y_{gj}=\sin\left(2X_{gj}\right)+2\exp\left(-16X_{gj}^{2}\right)+0.5e_{gj}, \] where $X_{gj}=\sqrt{\rho_{X}}\left(X_{1}\right)_{g}+\sqrt{1-\rho_{X}}\left(X_{2}\right)_{gj}$, $e_{gj}=\sqrt{\rho_{e}}c_{g}+\sqrt{1-\rho_{e}}u_{gj}$, and we generate $\left(X_{1}\right)_{g}\sim\mathcal{N}\left(0,1\right)$, $\left(X_{2}\right)_{gj}\sim\mathcal{N}\left(0,1\right)$, $c_{g}\sim\mathcal{N}\left(0,1\right)$, $u_{gj}\sim\mathcal{N}\left(0,1\right)$ independently. We set $\rho_{X},\rho_{e}\in\{0.2,0.5\}$. Note that larger $\rho_{X}$ and $\rho_{e}$ imply stronger cluster dependence on the regressor and the error term, respectively.
Setup 2 (heteroskedastic errors):
and $X_{gj}$ and $e_{gj}$ are generated in the same way as Setup 1.
A key feature is that Setup 1 has homoskedastic errors, and Setup 2 has heteroskedastic errors. We adopted the functional form $m(\cdot)$ for Setup 1 from fan1992variable and Setup 2 from kai2010local. The data-generating process for $X_{gj}$ and $e_{gj}$ are standard in the cluster dependence literature (cameron2008bootstrap; bartalotti2017regression). We set the weight function $w(x)$ for cross-validation and IAMSE equals to $w(x)=\mathbb{I}\left\{ \xi_{\mathrm{L}}\leq x\leq\xi_{\mathrm{U}}\right\} $, where we set $\xi_{\mathrm{L}}=-1.5$ and $\xi_{\mathrm{U}}=1.5$ for Setup 1, and $\xi_{\mathrm{L}}=0$ and $\xi_{\mathrm{U}}=1$ for Setup 2, respectively. For nonparametric regression, we use the Epachenikov kernel and local linear estimators. Results when using Nadaraya-Watson estimators are presented in Appendix $\text{\ref{sec:Add_sim}}$ because their values are similar to the ones by local linear estimators.
We will compare four methods of bandwidth choice: (i) rule-of-thumb (ROT), (ii) cluster-robust rule-of-thumb (CR-ROT), (iii) cross-validation (CV), and cluster-robust cross-validation (CR-CV). $h_{\text{CR-ROT}}$ (Equation $\text{\ref{eq:CR-ROT}}$) and $h_{\text{CR-CV}}$ (Equation $\text{\ref{eq:h_cv}}$) are what we suggested. The ROT bandwidth choice $h_{\text{ROT}}$ is proposed by fan1996local for i.i.d. observations. Instead of leave-one-cluster-out global fit as $\text{\eqref{eq:global-g}}$ for $h_{\text{CR-ROT}}$, it uses the global fit using the entire sample. $h_{\text{CV}}$ minimizes the cross-validation function. The difference between $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ is that $h_{\text{CV}}$ minimizes a criterion based on leave-one-out prediction errors, while $h_{\text{CR-CV}}$ minimizes a criterion based on leave-one-cluster-out prediction errors.
In simulation, we first compute $h_{\text{ROT}}$ and $h_{\text{CR-ROT}}$. Then, $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ are found by the grid search for $50$ points over $[h_{\text{CR-ROT}}/3,3h_{\text{CR-ROT}}]$. The performance of the methods of bandwidth selection is evaluated by the average squared error (ASE): \[ \operatorname{ASE}(h)=\frac{1}{n_{\text{grid }}}\sum_{k=1}^{n_{\text{grid }}}\left\{ \widehat{m}_{\mathrm{LL}}\left(u_{k},h\right)-m\left(u_{k}\right)\right\} ^{2}, \] where $\widehat{m}_{\mathrm{LL}}\left(u_{k},h\right)$ is the local linear estimator with the bandwidth $h$, and $\text{\ensuremath{\left\{ u_{1},\dots,u_{n_{\text{grid }}}\right\} } }$ are the grid points to evaluate the performance. We set the number of the grid $n_{\text{grid }}=50$ and $\text{\ensuremath{\left\{ u_{1},\dots,u_{n_{\text{grid }}}\right\} } }$ are evenly distributed over $\left[\xi_{\mathrm{L}},\xi_{\mathrm{U}}\right]$.
Tables $\text{\ref{tab:bw_LL1}}$ and $\text{\ref{tab:bw_LL2}}$ show means of the ASE for the local linear estimator and means of selected bandwidths (in curly brackets) across each simulation draw for Setup 1 and 2, respectively. Each table contains four methods of bandwidth choice in several scenarios. We consider combinations of homogeneous or heterogeneous size clusters, high or low cluster dependence on regressors, and high or low cluster dependence on error terms. In Setup 1 (Table $\text{\ref{tab:bw_LL1}}$, homoskedastic errors), $h_{\text{ROT}}$ and $h_{\text{CR-ROT}}$ have similar values of the ASE and the selected bandwidth, and $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ have the similar values of them, but $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ work better than $h_{\text{ROT}}$ and $h_{\text{CR-ROT}}$ in terms of the ASE. Within the same method of bandwidth choice, heterogeneous size clusters $n_{G}=100$ and high cluster dependence on regressors $\rho_{X}=0.5$ give a slightly larger ASE. Compared to them, high cluster dependence on error terms $\rho_{e}=0.5$ gives a much larger ASE.
In Setup 2 (Table $\text{\ref{tab:bw_LL2}}$, heteroskedastic errors), $h_{\text{ROT}}$ and $h_{\text{CR-ROT}}$ work poorly because they assume homoskedasticity. Different from Setup 1, $h_{\text{ROT}}$ has a larger ASE than $h_{\text{CR-ROT}}$. As Setup 1, $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ work well and have similar values of the ASE and the selected bandwidth. The good performance of $h_{\text{CV}}$ can not be explained by our theoretical results. We probably need an asymptotic analysis of $h_{\text{CV}}$ under the cluster dependence, which is outside of the scope of this paper.
To investigate how close the selected bandwidths are to the bandwidth that minimizes the ASE, we plot two figures for a scenario with $n_{G}=100$ and $\rho_{X}=\rho_{e}=0.5$. Figures $\text{\ref{fig:bw_LL1}}$ and $\text{\ref{fig:bw_LL2}}$ have values of bandwidth $h$ in the $x$-axis and means of the function $\mathrm{ASE}(h)$ in the $y$-axis, which are calculated from simulation draws for Setup 1 and 2, respectively. These figures also contain means of selected bandwidths by four selection methods and $h_{\mathrm{argmin}}$ minimizing $\mathrm{ASE}(h)$. We find that $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ are very close to $h_{\mathrm{argmin}}$ in both setups.
We recommend $h_{\text{CR-CV}}$ because it has a theoretical guarantee (Theorem $\text{\ref{thm:cr-cv}}$) and because it performs the best in our simulation, although the difference of ASEs between $h_{\text{CV}}$ and $h_{\text{CR-CV}}$ is subtle. In terms of the computational cost, $h_{\text{CR-CV}}$ is also better than $h_{\text{CV}}$ since leave-one-cluster-out estimators use smaller sample sizes than leave-one-out estimators do. $h_{\text{CR-ROT}}$ is useful for a rough estimation and for choosing the range of the grid search in cross-validation. We recommend $h_{\text{CR-ROT}}$ over $h_{\text{ROT}}$ for these purposes because it has a smaller ASE.
We will compare three methods to calculate 95% confidence intervals: (i) using the conventional standard error as for i.i.d. datasets ($CI$), (ii) using the cluster-robust standard error without the term related to the conditional covariance ($CI_{\mathrm{CR}}$), and (iii) using the cluster-robust standard error with the term related to the conditional covariance ($CI_{\lambda}$). More precisely, we calculate $CI$ with the standard error $\sqrt{R_{k}^{d}\widehat{\sigma}_{\mathrm{nw}}^{2}\left(x\right)/\left(nh^{d}\widehat{f}(x)\right)}$, $CI_{\mathrm{CR}}$ with the standard error $\sqrt{R_{k}^{d}\widetilde{\sigma}_{\mathrm{nw}}^{2}\left(x\right)/\left(nh^{d}\widehat{f}(x)\right)}$, and $CI_{\mathrm{\lambda}}$ with the standard error \[ \sqrt{\frac{1}{nh^{d}}}\sqrt{\frac{R_{k}^{d}\widetilde{\sigma}_{\mathrm{nw}}^{2}\left(x\right)}{\widehat{f}(x)}+\frac{\widehat{\lambda}\widehat{f}_{2}\left(x,x\right)\widehat{\sigma}_{\mathrm{nw}}\left(x,x\right)}{\left(\widehat{f}(x)\right)^{2}}}, \] where $\widehat{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ and $\widetilde{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ are nonparametrically estimated with $\widehat{e}_{gj}=Y_{gj}-\widehat{m}_{\mathrm{LL}}\left(X_{gj}\right)$ and $\widetilde{e}_{gj}=Y_{gj}-\widetilde{m}_{\mathrm{LL},-g}\left(X_{gj}\right)$, and $\widehat{\lambda}\widehat{f}\left(x,x\right)\widehat{\sigma}_{\mathrm{nw}}\left(x,x\right)$ is calculated parametrically as $\text{\eqref{eq:cov_simplified}}$. Note that in our data-generating processes, we have no cluster-level regressor $x^{\mathrm{(cls)}}$. In nonparametric regressions, bandwidths are selected as follows. The bandwidth $h_{f}$ for $\widehat{f}\left(x\right)$ is calculated by the reference bandwidth of the Epanechnikov kernel $h_{f}\approx1.049\cdot S_{X}\cdot n^{-1/5}$ where $S_{X}$ is a standard deviation of $X$ (e.g., see li2007nonparametric, Section 1.2). The bandwidth $h_{\sigma^{2}}$ for $\widehat{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ and $\widetilde{\sigma}_{\mathrm{nw}}^{2}\left(x\right)$ is set to $h_{f}$. Choosing $h_{\sigma^{2}}=h_{f}$ is a conventional choice, for example, used by imbens2012optimal. In this simulation, $\ensuremath{\widehat{\lambda}=20\cdot h_{m}}$ for $n_{G}=20$ and $\widehat{\lambda}\approx23.846\cdot h_{m}$ for $n_{G}=100$. To ignore the asymptotic bias, we set the bandwidth for $\hat{m}_{\mathrm{LL}}(x)$ using an undersmoothing approach defined as $h_{m}=h_{\text{CR-CV}}\times n^{1/5}\times n^{-2/7}$, where $h_{\text{CR-CV}}$ is computed using the CR-CV method as outlined in Section $\text{\ref{subsec:sim_bw}}$.
Appendix $\text{\ref{sec:Add_sim}}$ contains results without undersmoothing, where we set $h_{m}=h_{\text{CR-CV}}$. In the appendix, we explore two additional strategies on the bias: an infeasible analytical bias correction using knowledge of the true data-generating process, and simply ignoring the bias term.
The CIs are constructed at $x=0.75$ for Setup 1 and at $x=0.8$ and $0.4$ for Setup 2. Performances of confidence intervals are measured by the coverage ratio across each simulation draw.
Tables $\text{\ref{tab:CI_LL1}}$-$\text{\ref{tab:CI_LL2_0.4}}$ show the coverage ratio for local linear estimators and means of the length of confidence intervals (in curly brackets) across each simulation draw for Setup 1, Setup 2 with $x=0.8$, and Setup 2 with $x=0.4$, respectively. Each table contains results for three types of confidence intervals in several scenarios. As in Section $\text{\ref{subsec:sim_bw}}$, we consider 8 different scenarios with all possible combinations of $n_{G}\in\left\{ 20,100\right\} $, $\rho_{X}\in\{0.2,0.5\}$ and $\rho_{e}\in\{0.2,0.5\}$.
In Setup 1 (Table $\text{\ref{tab:CI_LL1}}$, homoskedastic errors), $CI_{\mathrm{CR}}$ has slightly better coverages than $CI$ does although both confidence intervals have severe under-coverage values when $\rho_{e}=0.5$. These confidence intervals work more poorly for the case $\max_{g\leq G}n_{g}=100$. On the other hand, $CI_{\lambda}$ performs the best among the three methods. It has accurate coverage (94.5%-96%) for every data-generating process.
For Setup 2 (heteroskedastic errors), we consider two different points (Table $\text{\ref{tab:CI_LL2_0.8}}$ for $x=0.8$ and Table $\text{\ref{tab:CI_LL2_0.4}}$ for $x=0.4$). Table $\text{\ref{tab:CI_LL2_0.8}}$ shows that $CI_{\lambda}$ improves the accuracy greatly, and it attains coverage ratios close to 95%. $CI_{\mathrm{CR}}$ and $CI$ fail to reach even 90% coverage ratios for almost all cases. However, Table $\text{\ref{tab:CI_LL2_0.4}}$ shows that all three methods have 95% coverage ratios, and $CI_{\lambda}$ has over-coverage values at $x=0.4$. Differences between Table $\text{\ref{tab:CI_LL2_0.8}}$ and Table $\text{\ref{tab:CI_LL2_0.4}}$ come from the functional form of the error term $\sigma\left(X_{gj}\right)e_{gj}$. Since $\sigma\left(x\right)=\left(2+\cos\left(2\pi x\right)\right)/5$ takes a large value at $x=0.8$ and a small value at $x=0.4$, the conditional variance and covariance of error terms also do so. Overall, $CI_{\lambda}$ is the most conservative choice among the three methods. Our proposed confidence interval $CI_{\lambda}$ performs well even without undersmoothing (see Appendix $\text{\ref{sec:Add_sim}}$).
We recommend $CI_{\lambda}$ because it works the best for homoskedastic errors, and it provides a conservative interval for heteroskedastic errors in our simulation.
In this section, we will apply our methods to a dataset from alatas2012targeting,\footnote{Their replication package, including datasets, is available on the AEA website.} which ran an experiment in 640 Indonesian subvillages with heterogeneous cluster sizes from 17 to 72. The purpose is to investigate a good way to target people with low incomes. In their subvillage-level randomized assignments, they compare three different ways of targeting: using demographic characteristics as proxies of income, using the community knowledge on the ranking of wealth (community targeting), and using a hybrid of them. The wealth ranking for community targeting was measured as follows. In each subvillage, people were asked to rank everyone in the community from the richest to the poorest. A facilitator used randomly ordered index cards, each representing a household. Starting the first two cards, the facilitator asked the community which household was better off in terms of wealth. Based on the community\textquoteright s response, the cards were placed with the wealth order. By sequentially adding one more index card to the comparison, the facilitator continued the process until all the households had been ranked.
One concern for this ranking process is that human errors could happen since it took 1.68 hours on average. alatas2012targeting investigated this concern by running a nonparametric regression of the mistarget rate ($Y_{gj}$) on the card order in the ranking process ($X_{gj}$). The mistarget rate is calculated based on the household\textquoteright s per capita consumption. The card orders in the ranking process are scaled from $0$ to $1$. In nonparametric regression, error terms may exhibit dependence within the same cluster due to unobserved subvillage-level shocks during the ranking process (e.g., instances where some individuals leave the room, distracting others) or strategic interactions (e.g., groups colluding to appear poorer than they actually are to secure future aid). These dependencies are captured by the subvillage-level cluster random effects.\footnote{alatas2012targeting state that “Since the targeting methods were assigned at the subvillage level, the standard errors are clustered to allow for arbitrary correlation within a subvillage.” in a parametric regression analysis.} We revisit alatas2012targeting with theoretically justified methods for cluster sampling. We will use the local linear regression with the Epachenikov kernel while alatas2012targeting used the local linear regression with the quartic kernel (what they call nonparametric Fan regression). Since the regressor is an observed variable rather than a treatment assignment, the model-based approach is taken.\footnote{Although cluster sizes were randomly selected in their study, our theoretical results remain applicable when conditioning on the realized cluster sizes.}
By the random card order, it is reasonable to assume that the regressor $X_{gj}$ is independent within the cluster (subvillage). Since the distribution of $X_{gj}$ does not follow from $U[0,1]$ due to the lack of observations on the mistarget rate $Y_{gj}$, we also estimate it nonparametrically. Thanks to the independence of the regressor, we can estimate the joint density by the product of marginal densities. Other detailed calculations for the bandwidth selection and standard errors are done in the same way as in Section $\text{\ref{sec:Monte-Carlo-simulation}}$.
The sub-dataset for the above regression contains $n=3784$ observations, $G=431$ subvillages, and each subvillage has from 4 to 9 observations. Thus, $\max_{g\leq G}n_{g}$ is $9$. The bandwidth selected by CR-CV was $h_{\text{CR-CV}}=0.1301$, and the undersmoothing version, obtained by multiplying by the factor $n^{1/5}\times n^{-2/7}$ was $h_{\text{undersmoothing}}=0.0642$. In contrast, alatas2012targeting heuristically choose a bandwidth of $(\max(X_{gj})-\min(X_{gj}))/5=0.1979$. We plot the cluster-robust cross-validation function in Figure $\text{\ref{fig:cv_Alatas}}$.
We calculated three 95% confidence intervals: $CI$, $CI_{\mathrm{CR}}$, and $CI_{\lambda}$. We calculate $\widehat{\lambda}\approx1.148$. Since $CI$ and $CI_{\mathrm{CR}}$ are almost identical, we only draw $CI_{\mathrm{CR}}$ on the plot. Figures $\text{\ref{fig:Alatas}}$ and $\text{\ref{fig:Alatas-undersmooth}}$ show the estimated nonparametric regression values and estimated pointwise confidence intervals with $h_{m}=h_{\text{CR-CV}}$ and $h_{m}=h_{\text{undersmoothing}}$, respectively. The asymptotic bias exists under $h_{m}=h_{\text{CR-CV}}$, whereas it is asymptotically dominated under $h_{m}=h_{\text{undersmoothing}}$. The estimated nonparametric function appears more oscillatory and has wider CIs in Figure $\text{\ref{fig:Alatas-undersmooth}}$ due to the larger variance. We found that $CI_{\lambda}$ is slightly wider than $CI_{\mathrm{CR}}$ in both figures, however the difference is smaller in Figure $\text{\ref{fig:Alatas-undersmooth}}$, as the smaller bandwidth reduces within-cluster correlation in smaller neighborhoods. We recommend trying both bandwidth choices as a robustness check in practice. We still have significant pointwise differences between the first few households and the household in the middle of the ranking process (mistargeting rate rises 5-10%) even under wider confidence intervals $CI_{\lambda}$. The conclusions are similar to alatas2012targeting.
This article has developed a comprehensive theoretical framework for nonparametric regression analysis under cluster sampling. Our contributions are threefold, addressing critical aspects of cluster-dependent data analysis that have significant implications for econometric methodologies and applied research. First, we allow both growing and bounded size clusters. This extension is crucial, as growing cluster sizes introduce a non-negligible within-cluster dependence, necessitating the inclusion of an additional term in the asymptotic variance to capture this phenomenon accurately. Second, we cover the case where regressors contain common variables within the same clusters. These cluster-level regressors are the extreme case of cluster-dependent regressors, and they require the careful estimation of the joint density function. Third, our proposed inference is valid with heterogeneous and growing cluster sizes. The simulation studies illustrate the critical role of accounting for within-cluster dependence, affirming the practical relevance of our theoretical insights.
While this article establishes a foundation for nonparametric regression analysis under cluster sampling, several avenues for future research emerge. Theoretical work on other nonparametric estimators, such as local polynomial regressions and series regressions, would be an interesting extension. Investigating boundary analysis is crucial due to its impact on estimator bias. Additionally, developing cluster bootstrap inference methods for nonparametric regressions is important since it would provide more practical statistical inference for clustered data. Lastly, deriving honest and adaptive uniform confidence bands, as done by chernozhukov2014anti and chen2024adaptive in the i.i.d. case, represents an important extension. This direction would depend on future advancements in empirical process theory under cluster sampling.