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.
164,029 characters · 22 sections · 58 citation commands
-2cm Optimal Estimation of Two-Way Effects under Limited Mobility
\and Sheng Chao Ho \and Frank Schorfheide.} }
JEL CLASSIFICATION: C11, C13, C23, C55, I21
KEY\ WORDS: Bipartite Graphs, Compound Loss, Empirical Bayes Methods, Limited Mobility, Matched Data, Shrinkage Estimation, Teacher Value-Added, Two-Way Fixed Effects, Unbiased Risk Estimation, Weak Identification.
\thispagestyle{empty} \setcounter{page}{0}
With the availability of large-scale linked datasets, researchers increasingly estimate high-dimensional two-way effects in interaction-based models following Abowd1999, see e.g., Card2013 for worker and firm effects, Kramarz2008 for student and school effects, and Finkelstein2016 for patient and area effects, among many others. Identification of unit-specific parameters in a linear two-way model crucially depends on forming a connected set through movements AbowdEtal2002, e.g., workers who move from one firm to another. Without any workers that connect two firms directly or indirectly, the difference in their firm effects cannot be identified. Similarly, teachers who move from one school to another establish the links to common students that are needed to identify their teacher value-added. In empirical data sets mobility, however, is often limited, which causes econometric challenges. Estimation becomes fragile if removal of a small number of movers could lead to disconnected sets and therefore the loss of identification. Based on graph theory, Jochmans2019 show how weak connectivity of the bipartite graph induced by the matched data adversely affects the commonly used least squares (LS) estimate. They also document weak connectivity among teachers across schools, suggesting the difficulty with teacher value-added estimation in the two-way model.
This paper takes on the challenge of limited mobility and proposes an optimal estimation of high-dimensional effects in a decision-theoretic framework. Specifically, we develop an empirical Bayes (EB) estimator based on a Gaussian prior distribution for the two sets of unit-specific parameters, which are indexed by hyperparameters that control mean, variance, and covariance. We treat the resulting posterior mean vector as a class of shrinkage estimators. Shrinkage estimators introduce bias in exchange for a reduction in variance, which is particularly appealing when the variance is large due to limited mobility. The shrinkage estimators are evaluated under a compound loss function that averages squared estimation errors across the large number of units. An unbiased estimate of the expected compound loss (risk) is minimized with respect to the hyperparameters to obtain an empirical Bayes estimator of the unit-specific parameters.
Our first contribution is to provide a novel prior on the joint distribution of the two sets of heterogeneous parameters by incorporating the observed matching patterns in the data. Assortative matching is an important feature in this interaction-based model, i.e., productive firms tend to employ high-skilled workers; see KLINE2024 for a review. Empirical Bayes methods leverage the common prior to improve estimation of unit-specific parameters. As such, it is important to have a prior that allows for stronger correlation of two unit-specific parameters associated with pairs that match more often in the dataset.\footnote{See Bonhomme2020 for a discussion of the conditionally independent random-effect model and the difficulty to specify conditional random-effects distributions consistent with link formation. Bonhomme2024 survey the random effect approach to network regression and its applications. } Our prior achieves this in a parsimonious fashion by introducing a hyperparameter that controls the scale of correlation and coupling it with the normalized matched frequency in the data. The flexibility to choose this correlation hyperparameter significantly enlarges the class of shrinkage estimator we consider. It yields remarkably improved estimates by pulling information across two sides of the interaction, even when the researcher is only interested in one of the two unit-specific parameters, e.g., firm effects or teacher-value added.
Our second important contribution is to establish frequentist asymptotic optimality in an asymptotic framework that allows for weak connectivity. The benchmark is an oracle estimator that determines the hyperparameters based on knowledge of the true parameters. We show that, conditional on the graph induced by the matched data and the true parameters, the compound loss differential between the proposed empirical Bayes estimator based on the unbiased risk estimate (URE) and the oracle estimator goes to zero asymptotically. Although the optimality is defined only for the unit specific parameters and their linear functions, the estimator provides a stepping stone for their nonlinear functions and subsequent regression analysis where these unit specific parameters are used as dependent variables. The frequentist optimality is robust to misspecification of the prior described above, i.e., when the prior does not align with the distribution of parameters across units conditional on the observed matches. The prior simply determines the class of shrinkage estimators considered. The optimality also does not depend on the Gaussian distribution of the regression errors.
Crucially, the asymptotic analysis is conducted in a framework that allows for weak connectivity, a scenario where the connectivity measure is small in finite-sample and is modeled to converge to zero asymptotically akin to weak instruments in instrumental variable regressions. We consider a sequence of graphs where the number of matches is bounded for each unit as the sample sizes of both sides go to infinity proportionally. As in Jochmans2019, we consider non-random graphs and use the small non-zero eigenvalues of the Laplacian matrix to measure its (global) connectivity. For asymptotic optimality of the URE based EB estimator, we derive conditions under which the URE criterion uniformly converges to the loss function over a large parameter space of the hyperparameters. In a setup where we only allow a finite-number of near-zero eigenvalues of the Laplacian matrix, this uniform convergence holds as long as each eigenvalue converges to zero more slowly than the inverse of the square root of the sample size.\footnote{This rate is associated with the parameter space of the hyperparameter and therefore the scope of the shrinkage estimator class where optimality is defined. If we confine the variance components of the hyperparameter to some bounded set, we can improve this rate from square root of the sample size to the sample size. } The asymptotic optimality is robust when all eigenvalues are bounded away from zero for well connected graphs.
Third, we show that our class of shrinkage estimators includes many commonly used estimators as special cases, and document in a Monte Carlo study that the proposed EB estimator dominates these competitors. By setting hyperparameters to specific values, we obtain not only the LS estimate of two-way effects, but also some one-way EB estimators widely used in the study of value-added in labor economics, see e.g., Kane2008, Chetty2014, and the review in WALTERS2024. With assortative matching and unobserved heterogeneity on both sides, the proposed estimator shows substantial improvement over the one-way shrinkage estimators. The class also includes the conventional EB estimator that determines the hyperparameters based on marginal likelihood of the data. Compared to the conventional EB estimator, the proposed URE-based EB estimator is less sensitive to the specification of the prior and demonstrate sizable improvement when the prior is misspecified to a large degree.
Fourth, we apply the proposed estimator in an empirical application based on a matched student-teacher data set from the North Carolina Education Research Data Center (NCERDC). Our estimators suggest that there is significant student-teacher assortative matching, and further that the rankings of teachers are sensitive to the estimator employed. Besides this specification, our method also apply to studies with assortative matching between students and schools and between teachers and schools, e.g., Kramarz2008 and Mansfield2015.
Shrinkage estimation has a long tradition in the statistics literature, dating back to Stein1956 and James1961 in the context of a vector of unknown means. The empirical Bayes approach can be traced back to Robbins1956. The recent statistics literature, e.g., Brown2009 and Jiang2009, and the recent panel data literature, e.g., Gu2017a,Gu2017b, Liu2020 have emphasized a nonparametric treatment of the prior distribution. Setting up a nonparametric prior for the joint distribution of two types of unit-specific parameters in the presence of matching pattern is an open question. The parametric prior we propose is a step in this direction.
Hyperparameter choice based on URE dates back to Stein1981 and has been used more recently, for instance, in Xie2012 and Kwon2021 for hyperparameter determination in models with fixed and time-varying one-dimensional heterogeneity, respectively, and in Brown2018 for a model with two-dimensional heterogeneity. The setup in Brown2018 is similar to ours but they are interested in optimal estimation of the cell means, which is the sum of two unit-specific parameters associated with the two sides, rather than the unit-specific parameters themselves. Our results do not follow from theirs directly. Furthermore, although they allow for missing cells, the amount of missing cells allowed is not suitable for the sparsely matched data we consider.\footnote{To measure the imbalance of the design matrix due to missing cells, Brown2018 use a diverging eigenvalue comparable to the inverse of the smallest non-zero eigenvalues of our Laplacian matrix. They require this divergence is slower than $n^{1/8}/(log(n))^2$, where $n$ denotes the sample size. This rate is suitable for occasional missing data rather than sparsely matched data.}
Our paper builds on some recent developments on the estimation of two-way effects with sparsely matched data. We adopt the bipartite graph representation and the measure of graph connectivity proposed by Jochmans2019. Our results focus on the compound decision that involves high-dimensional parameters rather than the estimation of an individual parameter, and we extend the study from the LS estimator to a large class of shrinkage estimator. Verdier2020 views the two-way effects as nuisance parameters and provides inference for the homogeneous slope coefficients in the context of sparsely matched data. We suggest using his estimator to concentrate out observed covariates before our shrinkage estimator is constructed.
In a contemporaneous paper, HeRobin2025 independently propose to use a ridge estimator for high-dimension two-way fixed effects estimation. Ridge estimators belong to the class of shrinkage estimator we consider.\footnote{We can obtain the ridge estimator by setting the correlation parameter and the location parameter in the hyperparameters to zero.} Although both papers find shrinkage estimation is highly desirable, our analysis focuses on different objects and studies them in a distinct asymptotic setups. HeRobin2025 focus on a regularized estimator of the Laplacian matrix, its inverse, and the bias and variance of the fixed effects estimator. In contrast, we focus on estimation of the fixed effects themselves (or its subvector). Their asymptotic analysis is based on a degree-corrected stochastic block model and closes the gap between the aforementioned estimators and their counterparts in a deterministic equivalent graph. Our asymptotic analysis conditions on the graph, hence it is deterministic, and closes the gap between a feasible decision rule and the optimal but infeasible decision rule to choose the regularization parameter. Through these complementary analyses, these two papers offer different perspectives to justify the adoption of shrinkage methods in two-way effect models.
Our empirical investigation complements the estimation of teacher value-added in the education economics literature see, for instance, Kane2008, Chetty2014, and Kwon2021 with various parametric EB methods and Gilraine2020 with a nonparametric EB methods. Our two-way EB estimator is robust to unobserved student heterogeneity not fully controlled by past test scores and other observable and assortative matching between students and teachers.\footnote{See Graham2008 for discussions of matching and sorting in the study of student achievement in different classrooms.} After obtaining the EB estimator, we also use it to compute various nonlinear functions of the high-dimensional effects. Some alternative methods are available in the literature, e.g., see Kline2020 for unbiased estimation of the variance and covariance, Armstrong2022 for EB confidence intervals, and Gu2023 for ranking and selection with EB methods
The remainder of the paper is organized as follows. The model and the LS estimator are presented in Section (ref). Section (ref) introduces our novel prior distribution, derives the posterior mean estimator of the unit-level effects, and discusses the URE based hyperparameter determination. The optimality theory for the proposed EB estimator is developed in Section (ref). Results from our Monte Carlo experiments are summarized in Section (ref) and the empirical analysis is presented in Section (ref). Finally, Section (ref) concludes. Proofs of the main results appear in the Appendix to this paper. Additional derivations, computational details, simulation results, and further information about the empirical analysis are relegated to an Supplementary Online Appendix.
We consider a two-way fixed effects model for matched data
where $i=1,\dots r$, $t = 1,\dots,T$, $j(i,t) \in \left\{1,\dots,c\right\}$ is the unit matched to $i$ in period $t$, $\alpha_i$ and $\beta_{j(i,t)} $ are fixed effects, $ \boldsymbol{x} _{it}\in\mathbb{R}^k$ is a set of observed regressors that could vary with $i$, $t$, and $u_{it}$ is an exogenous error with mean zero and known variance $\sigma^{2}$, and is i.i.d. across $i$ and $t$. When $\sigma^{2}$ is unknown, we can plug in an unbiased and consistent estimator. We discuss an extension to the heteroskedastic case when introducing the unbiased risk estimate. For the asymptotic optimality of the URE-based estimator, we focus on the homoskedastic case only.
In an application with student-teacher matched data, $y_{it}$ could be a student test score, which is a function of student ability $\alpha_i$ and teacher value-added $\beta_{j(i,t)}$, as well as other student and teacher characteristics. Here the function $j(i,t)$ characterizes the teacher $j$ assigned to student $i$ in period $t$, with the understanding that each student has only one teacher in any given time period. We will explain more carefully below that our econometric analysis conditions on the observed $j(i,t)$. Alternatively, in an application with employer-employee matched data $y_{it}$ could be the wage of employee $i$, $\alpha_i$ is the worker's skill or effort, $\beta_j$ is the firm-specific productivity, and $j(i,t)$ is the firm for which employee $i$ works in period $t$. Importantly, $\alpha_i$ and $\beta_{j(i,t)}$ are allowed to be correlated to account for assortative matching. We consider settings where the sample size of the two sides $r$ and $c$, symbolizing rows and columns of a matrix, both grow asymptotically and the time period $T$ is finite.
To facilitate the subsequent analysis, we stack the time series for each student in the order from $i=1$ to $r$ and write the model in ((ref)) in matrix notation as
where $ \boldsymbol{Y} ^*\in \mathbb{R}^{rT}$, $X\in \mathbb{R}^{rT\times k}$, $ \boldsymbol{U} \in \mathbb{R}^{rT}$, $ \boldsymbol{\theta} := ( \boldsymbol{\alpha} ^\prime, \boldsymbol{\beta} ^\prime)^\prime$ with $ \boldsymbol{\alpha} :=(\alpha_1,\dots,\alpha_r)^\prime$ and $ \boldsymbol{\beta} :=(\beta_1,\dots,\beta_c)^\prime$, and $B=[B_1,B_2]$ selects the unit $\alpha_i$ and $\beta_j$ out of the corresponding vectors for the matched outcome.\footnote{Specifically, in the student-teacher example, the entries of the matrix $B_1$ are indicators denoting the student associated with the corresponding outcome $y_{it}^*$ in $ \boldsymbol{Y} ^*$, and entries of the matrix $B_2$ are indicators denoting the teacher matched to student $i$ at time $t$.} The econometrician observes $ \boldsymbol{Y} ^*$, $B$, and $X$. Because $\alpha$ and $\beta$ appear additively in the model, a normalization is required for their separate identification. To this end, we impose the additional restriction $ \boldsymbol{1} _c^\prime \boldsymbol{\beta} =0$, which is conducive to the subsequent analysis on the subvector $ \boldsymbol{\beta} $.
Our goal is the optimal estimation of the fixed effects $ \boldsymbol{\theta} $ or its sub-vector instead of the regression coefficients $ \boldsymbol{\gamma} $.\footnote{General results in the Appendix cover linear functions of $ \boldsymbol{\theta} $.} To focus on the optimal shrinkage estimator and how it is related to the structure in $B$, we first study the case where $ \boldsymbol{\gamma} $ is known and consider the simplified model:
We subsequently show in Corollary (ref) that the optimality results can be extended to the case in which $ \boldsymbol{\gamma} $ is consistently estimated.
The identification of $ \boldsymbol{\theta} $ is determined by the link information captured in the matrix $B$, or equivalently the function $j(\cdot,\cdot)$. To study this relation in graph theoretic terms, the links are viewed as edges in a bipartite graph between two sets of nodes, $\mathcal{S}:=\left\{1,\dots,r\right\}$ (students or workers) and $\mathcal{T}:=\left\{1,\dots,c\right\}$ (teachers or firms). The set of edges between $i$ and $j$ is the set $\mathcal{E}_{ij}:=\left\{t\leq T: \ j(i,t)=j\right\}$. The graph implied by the matrix $B$ is then fully defined as $\mathcal{G}:=\left(\mathcal{S},\mathcal{T},\left\{\mathcal{E}_{ij}\right\}\right)$. We assume that the graph $\mathcal{G}$ is connected to ensure the identification of $\theta$.\footnote{A pair of nodes $i$ and $j$ is said to be connected if there exists a path between them -- that is, a sequence of unique edges that begins with one of $(i,j)$ and ends with the other. A graph is said to be connected if there exists a path for every pair of distinct vertices.}. Subsequently, we will be concerned with weak identification of $ \boldsymbol{\theta} $ in finite samples when the connectivity is weak.
We consider the compound loss of a high-dimensional parameter estimator $\hat{ \boldsymbol{\theta} }$. Specifically, the estimates of the unit-specific parameters are evaluated under the conventional quadratic loss function
where $W$ is a positive semi-definite matrix that serve both the role of selection and normalization. In particular, we consider two leading cases for $W$:
In the first case where $W=W_{a+b}$, we are interested in the full vector $ \boldsymbol{\theta} $. In the second case where $W=W_{b}$, we are interested in the subvector $ \boldsymbol{\beta} $ only. This weighting matrix appears subsequently in the loss-based evaluation of estimators.
The conventional two-way fixed effects estimator estimates $ \boldsymbol{\theta} $ by least squares (LS) which requires an inverse of the Laplacian matrix
The matrix $L$ is singular and the number of zero eigenvalues corresponds to the number of connected components in the graph. Therefore, a generalized inverse of $L$ is employed to obtain the LS estimator. We show in Lemma (ref) in the Supplementary Online Appendix that, under the normalization $ \boldsymbol{1} _c^\prime \boldsymbol{\beta} =0$ and the assumption of a single connected component, the unique (restricted) LS estimator is given by
where $L^{-}$ is a generalized inverse of $L$ that takes the form
and $ L^\dagger$ is the Moore-Penrose inverse of $L=B'B$. Using the partitions of ${\cal R}$, we can express the LS estimators of $ \boldsymbol{\alpha} $ and $ \boldsymbol{\beta} $ as $ \hat{ \boldsymbol{\alpha} }^{\text{ls}}=\mathcal{R}_a \hat{ \boldsymbol{\theta} }^{\text{\text{ls}}}$ and $\hat{ \boldsymbol{\beta} }^{\text{ls}}=\mathcal{R}_b \hat{ \boldsymbol{\theta} }^{\text{\text{ls}}}. $ The rotation matrix $\mathcal{R}$ ensures that the LS estimator satisfies the normalization condition: by construction, $ \boldsymbol{1} _c^\prime\hat{ \boldsymbol{\beta} }^{\text{ls}}=0$.
While $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ is an unbiased estimator of $ \boldsymbol{\theta} $, it is generally not optimal under the quadratic compound loss function in ((ref)). In particular, it is associated with a large compound loss when the graph $\mathcal{G}$ is weakly connected. Because $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ relies on the presence of movers, such as teachers/workers moving between different sets of students/firms, for estimation, it is imprecise when there are only a few movers between some isolated clusters, implying that the removal of a small number of edges is sufficient to separate $\mathcal{G}$ into disconnected components. In this case, some non-zero eigenvalues of the Laplacian matrix $L$ are close to zero, which implies that the generalized inverse $L^{\dagger}$ that appears in $\hat{ \boldsymbol{\theta} }$ has large eigenvalues. Below we study an asymptotically optimal estimator that minimize the compound loss in an asymptotic framework where a subset of the eigenvalues of $L$ converge to $0$ asymptotically, modeling estimation with a weakly connected graph.
To obtain an optimal estimator that minimizes the quadratic compound loss function, we pursue an empirical Bayes strategy. Our starting point is a hierarchical Bayes model with a family of prior distributions indexed by a vector of hyperparameters. Based on the class of priors we derive a posterior mean estimator and then subsequently determine the hyperparameters in a data-driven manner, which leads to the empirical Bayes estimator.
A generative economic model would start from the vector of unit-specific coefficients $ \boldsymbol{\theta} $ characterizing agents “technologies and preferences.” A matching mechanism that groups students into classes and assigns teachers to classes, or a search model that connects workers and firms, determines the links encoded in the matrix $B$ and generates a conditional distribution $p(B| \boldsymbol{\theta} )$. Finally, outcomes are determined based on $p( \boldsymbol{Y} |B, \boldsymbol{\theta} )$ which is characterized by ((ref)). In a Bayesian setting, one specifies a prior distribution for $ \boldsymbol{\theta} $, which we denote by $p( \boldsymbol{\theta} )$. This leads to a joint distribution of data and parameters:
According to Bayes Theorem the posterior distribution of $ \boldsymbol{\theta} $ given the observables $( \boldsymbol{Y} ,B)$ is
where $\propto$ denotes proportionality. To obtain the second line, we replaced $p(B| \boldsymbol{\theta} ) p( \boldsymbol{\theta} )$ by $p( \boldsymbol{\theta} |B) p(B)$ and absorbed $p(B)$ into the constant of proportionality because it does not depend on $\theta$. Rather than deriving $p( \boldsymbol{\theta} |B)$ from a structural model of link formation $p(B| \boldsymbol{\theta} )$, we directly specify a prior $p( \boldsymbol{\theta} |B, \boldsymbol{\lambda} )$, where $ \boldsymbol{\lambda} :=(\mu,\lambda_a,\lambda_b,\phi)$ is a vector of hyperparameters that incorporates the assortative matching patterns reasonably expected in practice. We emphasize that the prior only serves to define the estimator class, and is not assumed to be correctly specified for our subsequent optimality results to hold.
The prior distribution $p( \boldsymbol{\theta} |B, \boldsymbol{\lambda} )$ is described in Section (ref) and the mean of the posterior distribution $p( \boldsymbol{\theta} |Y,B, \boldsymbol{\lambda} )$ is derived in Section (ref). Hyperparameter determination is discussed in Section (ref) and Section (ref) shows how to covariates can be used to center the prior for $ \boldsymbol{\theta} $, which may be useful in applications.
We proceed by introducing our proposed prior $p( \boldsymbol{\theta} |B, \boldsymbol{\lambda} )$. To build the observed link information in the prior for $ \boldsymbol{\theta} $, we first consider the Laplacian matrix $L=B'B$ and its interpretation. Because $B_1$ and $B_2$ are selection matrices by construction, we know that $B_1'B_1$ is a $r \times r$ diagonal matrix whose diagonal elements counts the number of times student $i$ appears in the data, $B_2'B_2$ is a $c \times c$ diagonal matrix whose diagonal elements count the number of times teacher $j$ appears in the data, $B_1'B_2$ is a $r \times c$ matrix whose $(i,j)$ element equals the number of times student $i$ and teacher $j$ are matched over the $T$ time periods. Define
which are the degree, adjacency, and normalized adjacency matrices of $\mathcal{G}$, respectively. Because $B_1'B_1$ and $B_2'B_2$ are both diagonal, we have
where the $(i,j)$ element of the upper right submatrix $\mathcal{A}_{12}=(B_1'B_1)^{-1/2}B_1'B_2(B_2'B_2)^{-1/2}$ measures the normalized matched frequency between $i$ and $j$. In the proposed prior, we model assortative matching based on the observed matched frequency in $\mathcal{A}$.
Conditional on $B$, we introduce a normal prior for $ \boldsymbol{\theta} $ that depends on a vector of hyperparameters $ \boldsymbol{\lambda} :=(\mu,\lambda_a,\lambda_b,\phi)$, where $\mu$ models the mean value of $\alpha_i$, $\lambda_a$ and $\lambda_b$ model the precision (inverse of variance) of $\alpha_i$ and $\beta_j$, respectively, and $\phi$ models the correlation between $\alpha_i$ and $\beta_j$ given $\mathcal{A}$. The mean of $\beta_j$ is $0$ in the prior to be consistent with the normalization. Specifically, the the proposed prior takes the form
In the special case where $\phi=0$, the prior in (ref) reduces to
Furthermore, the prior under $\phi=0$ makes the simplifying assumption that $\alpha_i$ and $\beta_j$ are all independent and the prior does not depend on $B$:
On the other hand, when $\phi\neq0$, Lemma (ref) shows that the prior incorporates the assortative matching patterns present in the bipartite graph in a parsimonious way through $\phi$.
Lemma (ref) shows that the correlation between $\alpha_i$ and $\beta_j$ is entirely determined by the parameter $\phi$ and their observed matched frequency in the normalized adjacency matrix $\mathcal{A}$. A positive value of $\phi$ implies positive assortative matching. For example, a higher value of $\beta_j$, i.e., a teacher with higher value-added or a firm with higher productivity, increases the conditional mean of $\alpha_i$ connected to it, implying on average better students taught by teacher $j$ or more productive workers in firm $j$.
For $\phi \not=0$ our prior (ref) avoids the joint independence assumption of (ref) and adopts only a conditional independence assumption:
with a correspondence for $p( \boldsymbol{\beta} |B, \boldsymbol{\lambda} , \boldsymbol{\alpha} )$. This bears similarity to the one-sided correlated random effects specification employed by Bonhomme2019.
Under the assumption that $u_{it} | (B, \boldsymbol{\theta} ) \sim_{iid} \mathcal{N}(0,\sigma^2)$, the posterior mean under the prior in ((ref)) is given by
where
Due to the normalization $ \boldsymbol{1} _c^\prime \boldsymbol{\beta} =0$, this posterior mean formula involves the rotation matrix $\mathcal{R}$, previously defined in ((ref)). A derivation is provided in Lemma (ref) in the Supplementary Online Appendix.
The hyperparameter $ \boldsymbol{\lambda} $ is chosen over the set $ \boldsymbol{\lambda} =(\mu, \lambda_a, \lambda_b, \phi)\in \bar{\mathcal{J}} := [-\bar{\mu},\bar{\mu}] \times [0,\infty] \times [0,\infty] \times [\-\bar{\phi},\bar{\phi}] $ for some finite $\bar{\mu}$ and $0<\bar{\phi}<1$.\footnote{In many applications bounds for $\mu$ can be obtained from bounds on the outcome variable, e.g., test scores or wages. Shrinking to a location outside these bounds does not further reduce the loss. In contrast, we do not impose upper and lower bounds for $\lambda_a$ and $\lambda_b$ in order to consider optimal estimation in a wider class. As we show below, setting $\lambda_a$ or $\lambda_b$ to $0$ or $\infty$ leads to some commonly used one-way fixed effect estimators.} The hyperparameter $ \boldsymbol{\lambda} $ affects the posterior mean through $S$ and $ \boldsymbol{v} $, which determine the shrinkage weight and the shrinkage location, respectively. Specifically, $S$ depends on $(\lambda_a,\lambda_b,\phi)$ and $ \boldsymbol{v} $ depends on $\mu$. When $\lambda_a=0$ and $\lambda_b=0$, the posterior mean is the LS estimator. Moreover, when $\lambda_a \rightarrow \infty$ and $\lambda_b \rightarrow \infty$, the posterior mean becomes the prior mean. We will regard $\{\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ) \, | \, \lambda \in {\bar{\mathcal{J}} } \}$ as a class of shrinkage estimators and study the problem of optimally choosing an estimator within this class.
Given a data-driven choice of the hyperparameter $\hat{ \boldsymbol{\lambda} }$, the resulting estimator $\hat{ \boldsymbol{\theta} }(\hat{ \boldsymbol{\lambda} })$ is an empirical Bayes (EB) estimator; see Robbins1956 or, for instance, Robert1994 for a textbook treatment. Although this class of shrinkage estimator is motivated by the Gaussian prior and the Gaussian regression error, the the EB estimator proposed below is asymptotically optimal within this large class of estimators even when these Gaussian distributions are misspecified.
{\bf Unbiased Risk Estimation.} We consider an EB based estimator that chooses the hyperparameter $ \boldsymbol{\lambda} $ by minimizing an unbiased estimate of the risk. To define the estimation risk, we consider a frequentist approach that conditions on the unit-specific effects $ \boldsymbol{\theta} $. Because we deliberately abstracted from a theory of how matches are formed, we also condition on the observed matrix $B$. In applications in which the matching mechanism is deterministic and $p(B| \boldsymbol{\theta} )$ is a point mass, there would be no other choice than to condition on $B$. For a given hyperparameter $ \boldsymbol{\lambda} $, we define the risk as
The expectation is taken with respect to $ \boldsymbol{Y} $ and the subscript indicates that we condition on $( \boldsymbol{\theta} ,B)$.
We now construct an unbiased estimate of $R_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$, denoted by $\text{URE}( \boldsymbol{\lambda} )$, that only depends on the observations $( \boldsymbol{Y} ,B)$ and not on the unknown parameter vector $ \boldsymbol{\theta} $. Using the formula for $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} )$ in (ref), the risk function can be written as
Under the quadratic loss, the second term on the right-hand side of (ref) comes from the variance of $S_1\hat{ \boldsymbol{\theta} }^{\text{ls}}$, which is known, and the first term is a quadratic function of the bias induced by the prior. A natural unbiased estimator of the first term is
where $\sigma^2\operatorname*{tr}[ S'WS L^- ]$ is used for bias correction because
Combining (ref) and (ref) , we obtain
The unbiasedness of $\text{URE}( \boldsymbol{\lambda} )$ is summarized in the following Lemma.
We denote the URE hyperparameter estimate by
where $ \boldsymbol{\lambda} =(\mu, \lambda_a, \lambda_b, \phi)\in \mathcal{J}:= [-\bar{\mu},\bar{\mu}] \times (0,\infty) \times (0,\infty) \times [\-\bar{\phi},\bar{\phi}] $.\footnote{Here we rule out $0$ and $\infty$ from the parameter space of $\lambda_a$ and $\lambda_b$ to facilitate showing uniform convergence of some functions over $\mathcal{J}$, which is needed for the asymptotic analysis below. Importantly, we do allow $\lambda_a$ and $\lambda_b$ to be arbitrarily close to the boundaries.} The proposed estimator is $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$. Justification for $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ goes beyond the unbiased property in Lemma (ref). In the next section, we show that $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ is the asymptotic optimal estimator in a decision theoretic framework that allows for weakly connected graphs.
{\bf Treatment of $\sigma^2$.} To implement the hyperparameter determination one can replace $\sigma^2$ with its LS squares estimate
Lemma (ref) in the Supplementary Online Appendix shows that this estimator is unbiased, and it is also consistent under weak conditions. A practitioner might be concerned about heteroskedasticity. In this case $\sigma^2\operatorname*{tr}[ S'WS L^- ]$ and $\sigma^2\operatorname*{tr}[ S_1'WS_1 L^- ]$ in (ref) can be, respectively, replaced with the leave-out estimates of $\operatorname*{tr}\big[ S'WS\mathbb{V}_{ \boldsymbol{\theta} , B}[\hat{ \boldsymbol{\theta} }^{ \text{ls}}] \big]$ and $\operatorname*{tr} \big[ S_1'WS_1 \mathbb{V}_{ \boldsymbol{\theta} , B}[\hat{ \boldsymbol{\theta} }^{\text{ls}}] \big]$ provided by Kline2020.\footnote{We thank Patrick Kline for this insightful suggestion.} This leads to a heteroskedastic-robust $\text{URE}( \boldsymbol{\lambda} )$. In addition, one could consider shrinkage estimation for $\sigma^2_i$ as well, as it is commonly done in full Bayesian implementations of one-way effect models, e.g., Liu2023 and Liu2023a. In the remainder of this paper we will establish asymptotic optimality within the homoskedastic framework. Extensions to the heteroskedastic case are left for future research.
{\bf Marginal Likelihood.} Alternatively, one could follow the tradition in the EB literature and estimate $ \boldsymbol{\lambda} $ by maximizing the marginal likelihood function
which leads to
and $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{mle}})$. We expect $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{mle}})$ to perform well whenever the prior (ref) is correctly specified. If $ \boldsymbol{\theta} $ is generated from a different distribution than (ref), then the marginal likelihood will be misspecified and $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{mle}})$ is likely to be sub-optimal. Comparisons in the Monte Carlo studies in Section (ref) confirm that $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{mle}})$ is more sensitive to the prior distribution than $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$.
In practice it might be useful to shrink the fixed effects to an index constructed from observable covariates that are informative about the latent heterogeneity. Define $Z_a\in\mathbb{R}^{r\times k_a}$, $Z_b\in\mathbb{R}^{c\times k_b}$, $ \boldsymbol{\delta} _a\in\mathbb{R}^{k_a}$ and $ \boldsymbol{\delta} _b\in\mathbb{R}^{k_b}$, where $k_a$ and $k_b$ are the dimensions of the regressors used to construct covariate-based prior means for the $\alpha_i$s and $\beta_i$s. For student-teacher data sets, examples of potential $Z_a$ and $Z_b$ include student demographic information, such as gender, ethnicity, and teacher information such as years of teaching experience, respectively. Define
The location of shrinkage for $ \boldsymbol{\theta} $ is now $Z \boldsymbol{\delta} $, where $ \boldsymbol{\delta} $ is determined jointly with the other hyperparameters in minimizing the URE. This is operationalized by replacing $ \boldsymbol{v} $ of (ref), (ref) and (ref) by $Z \boldsymbol{\delta} $ and then performing the minimization as before -- now minimizing the URE over $ \boldsymbol{\delta} $ as well. Section (ref) of the Supplementary Online Appendix shows how to concentrate out $ \boldsymbol{\delta} $ when performing the URE-minimization.
We next provide some asymptotic justifications for the empirical Bayes estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$. We study a sequence of growing graphs $\left\{\mathcal{G}_{r,c}\right\}$ as $r,c \rightarrow \infty $, where we use $r$ and $c$ in the subscript to indicate that the sizes of the nodes as well as other characteristics, such as the connectivity measures, could change with the sample sizes. The number of periods $T$ is fixed in our asymptotic analysis. For notational simplicity, we omit the dependence on $r, c$ in other graph characteristics such as $B$ and the Laplacian matrix $L$. Section (ref) discusses our assumption and the convergence of the URE objective function. The optimality result is stated in Section (ref). A comparison of the proposed two-way shrinkage estimator to other popular estimators is presented in Section (ref) and the extension to the case of additional estimated regression coefficients is provided in Section (ref).
We first state some regularity assumptions on the regression model. Throughout all the assumptions and proofs, we use $\epsilon$ and $M$ to represent generic constants for lower and upper bounds. They do not have to take the same value when they appear in different contexts.
Condition (i) restricts the heterogeneity of fixed effects, and is a condition widely used in the empirical Bayes literature, see e.g., Xie2012 and Brown2018. Condition (ii) imposes only a finite fourth moment on the regression error rather than requiring its normal distribution as often seen in these papers. The finite fourth moment is used to bound the variance of some quadratic forms of the regression errors. Next, we consider some conditions on the sequence of graphs $\{\mathcal{G}_{r,c}\}$.
Condition (i) requires the sizes of the two sets of nodes $\mathcal{S}$ and $\mathcal{T}$ to be proportional and grow asymptotically at the same order. Condition (ii) restricts each node to having a finite degree, characterizing its finite number of matches. In a student-teacher application both conditions could be satisfied by assuming that each student only takes one class per time period, the number of time periods $T$ is fixed, and the class-size stays asymptotically bounded from below and above. In a worker-firm application the conditions allow for a distribution of firm sizes (in terms of number of workers), but rule out asymptotics along which, for instance, one firm employs 50% of the workers. Under Assumption (ref) it is not possible to estimate $ \boldsymbol{\theta} $ consistently.
To model the limited mobility issue in finite sample, we consider a sequence of graphs $\left\{\mathcal{G}_{r,c}\right\}$ whose connectivity gets weaker asymptotically. This connectivity measures are the eigenvalues of the Laplacian matrix $L=B'B$ when we are interested in the full vector $ \boldsymbol{\theta} $. When we are interested in the subvector $ \boldsymbol{\beta} $, the connectivity measures become the eigenvalues of the Laplacian matrix of the one-mode projected graph\footnote{In the one-mode projected graph, the vertices are $\mathcal{T}$ and the nodes from $\mathcal{S}$ moving between them act as the edges. Two nodes from $\mathcal{T}$ are thus connected whenever they are connected to a common node from $\mathcal{S}$ in the original graph $\mathcal{G}$, with the strength of their connection determined by the number of common matches, see Jochmans2019 for detailed discussions of this one-mode projection graph.}
For clarity of the graph-theoretic interpretation, we focus on these two leading case in the main paper. Proofs of all theoretical results allow for a general subvector selection and weighting under the high-level conditions in Assumptions (ref) and (ref) in the Appendix. The Assumptions below are their sufficient conditions in these two leading cases.
Let $\rho_{\ell}(A)$ denote the $\ell^{th}$ smallest eigenvalue of a symmetric matrix $A$. For a connected graph, it is known that $\rho_1(L)=0$. Jochmans2019 show $\rho_2(L)$ is crucial in the analysis of the LS estimate. Here, we allow $\rho_2(L)$ as well as some other small eigenvalues to converge to $0$ as $r,c \rightarrow \infty$, in order to study optimal estimation in an asymptotic framework for weakly connected graphs. When the parameter of interest is the subvector $ \boldsymbol{\beta} $, we allow $\rho_2(L_{2,\perp})$ and some other small eigenvalues of $L_{2,\perp}$ to converge to $0$ asymptotically.
Below we present some assumptions on these connectivity measures in order to obtain desirable properties of $\text{URE}( \boldsymbol{\lambda} )$ as a selection criterion for $ \boldsymbol{\lambda} $ and the proposed EB estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$. Beyond the unbiasedness in Lemma (ref), we show the distance between the compound loss function $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and the selection criterion $\text{URE}( \boldsymbol{\lambda} )$ is small under some stronger measures, which eventually leads to asymptotic optimality of $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$. To this end, the conditions typically involve how weak the connectivity can be in order for $\text{URE}( \boldsymbol{\lambda} )$ to be informative. As such, we impose bounds on how fast $\rho_2(L)$ or $\rho_2(L_{2,\perp})$ can go to zero and how many of these eigenvalues can go to zero. It is important to note that although the empirically relevant asymptotics are obtained when $\rho_2(L)$ or $\rho_2(L_{2,\perp})$ converge to zero, none of our assumptions require that the second smallest eigenvalue vanishes. These desirable properties of $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ also hold for well connected graphs where all non-zero eigenvalues are bounded away from zero as the sample size increases.
Assumptions (ref) and (ref) present two versions of this connectivity assumption, one for $W=W_{a+b}$ for the full vector $ \boldsymbol{\theta} $ and one for $W=W_{b}$ for the subvector $ \boldsymbol{\beta} $. This set of Assumptions on graph connectivity ensures that not only the compound loss function $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and the selection criterion $\text{URE}( \boldsymbol{\lambda} )$ share the same expectation by construction, but also the variance of their difference converge to zero asymptotically.
Define
and note that Assumptions (ref) and (ref) imply that $\delta_{r,c} \rightarrow \infty$.
Theorem (ref) shows pointwise convergence between $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and $\text{URE}( \boldsymbol{\lambda} ) $ in the $L_2$ norm. It also shows that the convergence depends both on the sample size and the smallest non-zero eigenvalue of the associated graph. As long as the smallest non-zero eigenvalue converges to $0$ slower than the square root of the sample size, we have the convergence in the $L_2$ norm.
Below we impose stronger conditions on the connectivity measure to ensure uniform convergence between $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and $\text{URE}( \boldsymbol{\lambda} ) $ over $ \boldsymbol{\lambda} \in \mathcal{J}$.
This assumption allows for at most $k$ eigenvalues converge to zero, requiring the full graph to have at most $k$ weakly connected components. This condition is particularly useful to obtain uniform convergence, without imposing stronger conditions on the rate in Assumptions (ref) and (ref). In this large-scale estimation problem, $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and $\text{URE}( \boldsymbol{\lambda} ) $ both involve trace of matrices whose dimensions grow with the dimension of $ \boldsymbol{\theta} $. In particular, they involve the inverse of the non-zero eigenvalues of $L$, which are explosive for weakly connected components. Having only a finite number of weakly connected components ensures that we only work with a finite number of explosive ones. Our proof in the Appendix can accommodate a slowly growing number of weakly connect components, if we restrict the weakness of the connectivity. We show this trade off through a general trace condition in Assumption (ref) in the Appendix. In the empirical application with student-teacher linked data, this finite component assumption requires that students and teachers are generally well connected within each school, and that there are only a finite number of schools across which movements of students and teachers is limited.
Assumption (ref) is intended for the estimation of the full vector $ \boldsymbol{\theta} $, i.e., $W=W_{a+b}$. When the interest is on the subvector $ \boldsymbol{\beta} $, i.e., $W=W_{b}$, we have a similar finite component assumption on the one-mode projected graph.
In addition to this requirement, uniform convergence between $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and $\text{URE}( \boldsymbol{\lambda} ) $ for $W=W_{b}$ also require some minimum assumption on the full graph because we conduct all the estimation simultaneously, despite our interest is on the subvector only.
This is a very weak condition on the full graph by considering $\epsilon$ close to $0$. It is only relevant when we are interested in the subvector. When the full vector is of interested, this condition is already implied by Assumption (ref).
Building on the uniform convergence between $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ and $\text{URE}( \boldsymbol{\lambda} )$, now we show that $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ is asymptotically optimal in the sense that its compound loss is comparable to that of an oracle estimator that choose $ \boldsymbol{\lambda} $ based on $l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} )$ directly. To this end, define
and the oracle estimator oracle $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ol}})$. The superscript “ol” stands for oracle loss because it is with knowledge of the parameters $ \boldsymbol{\theta} $. Note that this oracle estimator targets the actual compound loss rather than the risk and it achieves the lowest possible compound loss within the class of shrinkage estimator considered. Under parameter heterogeneity, this oracle loss is strictly positive. We say that an estimator is asymptotically optimal if it achieves the oracle loss.
The main theoretical result of the paper is that $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ is asymptotically optimal.
Theorem (ref) follows immediately from the uniform convergence in Lemma (ref) and $\text{URE}( \boldsymbol{\lambda} )$ is minimized by $ \boldsymbol{\lambda} ^{ure}$. It shows that the proposed estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ is asymptotically as good as the oracle estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ol}})$ in the sense that they achieve the same compound loss in probability. Crucially, this holds even under weak connectivity permitted by Assumptions (ref) to (ref). The asymptotic optimality established in Theorem (ref) is a statement that conditions on the fixed effects $ \boldsymbol{\theta} $. It therefore holds even if the priors of (ref) are misspecified, i.e., when $ \boldsymbol{\theta} $ is generated from a distribution other than (ref). This makes the estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ robust to misspecification, which is not guaranteed for other empirical Bayes estimators that rely on distributional assumptions in selecting $ \boldsymbol{\lambda} $, such as $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{mle}})$.
The class of two-way estimators defined by (ref) trivially covers the LS estimator by setting $\lambda_a=\lambda_b=0$. In addition, it covers two widely-used estimators that only shrink in the $\beta$ dimension. We denote the first estimator by $\hat{ \boldsymbol{\beta} }^{\text{1-way}}(\tilde \lambda_b)$. It is constructed by assuming that the $\alpha_i$s are homogeneous. We initially impose the normalization on $\alpha_i=\alpha=0$ and derive the posterior mean based on the prior $\beta_j \sim_{iid} {\cal N}(\underline{\beta},\sigma^2/\tilde \lambda_b)$. {\em Ex post} we demean the estimator to make it comparable to the one that we used previously. This leads to
Notice that the {\em ex post} demeaning eliminates the effect of the prior mean $\underline{\beta}$. The estimator $\hat{ \boldsymbol{\beta} }^{\text{1-way}}(\tilde \lambda_b)$ can be obtained as the limit of our two-way shrinkage estimator as the precision of the prior for $\alpha_i$ goes to infinity.
If we let $\hat{ \boldsymbol{\beta} }^{\text{1-way}}(\lambda_b^{\text{mom}})$, where $\lambda_b^{\text{mom}} := \tfrac{\sigma^2}{\text{var}(\beta_j)}$ and $\text{var}(\beta_j)$ is an empirical estimate of the variance within $ \boldsymbol{\beta} $, we obtain the one-way EB estimator of Kane2008. This estimator is widely used for the estimation of teacher-value added and is generalized to time-varying $ \boldsymbol{\beta} $ in Chetty2014 and Kwon2021. Compared to $\hat{ \boldsymbol{\beta} }^{\text{1-way}}(\cdot)$, our two-way EB estimator uses a data-dependent choice of $\lambda_a$ rather than setting it to infinity. When $\alpha_i$ is indeed heterogeneous, due to unobserved heterogeneity not fully controlled by observed covariates, the omitted variable bias could be large in the one-way estimator under assortative matching.
The second estimator actually maintains the heterogeneity of the $\alpha_{i}$s but conducts shrinkage in the $ \boldsymbol{\beta} $ dimension only. Specifically, this estimator first projects out $ \boldsymbol{\alpha} $ using LS and subsequently estimates $ \boldsymbol{\beta} $ with a one-way EB estimator that respects the normalization restriction $ \boldsymbol{1} _c' \boldsymbol{\beta} =0$. We use $\hat{ \boldsymbol{\beta} }^{\text{2-way}}(\tilde \lambda_b)$ to denote the posterior of $ \boldsymbol{\beta} $ under a Gaussian prior, after $ \boldsymbol{\alpha} $ is projected out. Its formula is given by
Here $\hat{ \boldsymbol{\beta} }^{\text{ls}}$ is the $\beta$ component of $\hat{\bs \theta}^{\text{ls}}$ defined in ((ref)), $L_{2,\perp}$ is the Laplacian matrix of the one-mode projected graph in (ref), and $M_c$ replaces ${\cal R}$ in the definitions of $S_1$ and $S$ in the posterior mean equation (ref). When $\tilde{\lambda}_b= 0$, define $\hat{ \boldsymbol{\beta} }^{\text{2-way}}= M_c\cdot \hat{ \boldsymbol{\beta} }^{\text{ls}}$. This estimator $\hat{ \boldsymbol{\beta} }^{\text{2-way}}(\tilde \lambda_b)$ is our two-way shrinkage estimator with the precision of the prior for $\alpha_i$ set to zero.
Chetty2018 use a one-way shrinkage estimator akin to $\hat{ \boldsymbol{\beta} }^{\text{2-way}}(\lambda_b)$ to estimate the causal effects of growing up in different neighborhoods on the future outcomes of resident children. In this application, the mover-based identification\footnote{More precisely, the identification scheme uses a more specific variation, namely differences in movers' timings from one neighborhood to another to identify the relative per-year causal effects of the two neighborhoods.} and LS estimation eliminates unobserved assortative matching between individuals and neighborhoods. Then, a one-way shrinkage is conducted to estimate the causal effect for different neighborhoods.\footnote{While the shrinkage estimators in these applications typically use only the variances of $\hat{ \boldsymbol{\beta} }^{\text{ls}}$ instead of $L_{2,\perp}$ in the shrinkage weight, we expect the latter to be more efficient as it takes into account the entire correlational structure.} Compared to $\hat{ \boldsymbol{\beta} }^{\text{2-way}}(\lambda_b)$, the two-way shrinkage estimator benefit from pooling information in the $ \boldsymbol{\alpha} $ dimension, and importantly utilizing information in the observed assortative matching pattern with a data-dependent choice of the correlation parameter $\phi$.
In Theorem (ref) we established that choosing $ \boldsymbol{\lambda} = \boldsymbol{\lambda} ^{ure}$ is optimal. It is tempting to deduce that as soon as a researcher chooses $ \boldsymbol{\lambda} \not= \boldsymbol{\lambda} ^{ure}$, as she most likely would when using $\hat{ \boldsymbol{\beta} }^{\text{ls}}$, $\hat{ \boldsymbol{\beta} }^{\text{1-way}}(\tilde{\lambda}_b)$, or $\hat{ \boldsymbol{\beta} }^{\text{2-way}}(\tilde \lambda_b)$, the resulting estimator is no longer optimal in the sense of Definition (ref). Unfortunately, this form of suboptimality does not follow directly from Theorem (ref). We proceed with a formal definition of dominance that is akin to inadmissibility.
{\bf Example:} To find parameter values $(B, \boldsymbol{\theta} )$ under which $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ is dominated by $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$, it is sufficient to find $(B, \boldsymbol{\theta} )$ under which $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ has strictly larger risk than a trivial estimator $ \boldsymbol{0} _{r+c}$, which can be generated as a limit of $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ by setting $\mu=0$, $\phi=0$, and letting $\lambda_a,\lambda_b \rightarrow \infty$. The risk of the estimator $ \boldsymbol{0} _{r+c}$ is $\tfrac{1}{r+c} \boldsymbol{\theta} ' \boldsymbol{\theta} $, and the risk of $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ is\footnote{The inequality in (ref) follows from von Neumann's trace inequality, the fact that ${\cal R'R}$ has a single eigenvalue of $r/c$ and eigenvalues of 1 with multiplicity of $r+c-1$, and an assumption that $r>c$ that corresponds to the setting of our empirical application.}
Thus, for DGPs satisfying
$\hat{ \boldsymbol{\theta} }^{\text{ls}}$ has larger risk than $ \boldsymbol{0} _{r+c}$. Under these DGPs, $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ also has a larger risk than $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ and is dominated in the sense of Definition (ref). Condition ((ref)) is satisfied in settings in which the graph is weakly connected (some eigenvalues of $L=B'B$ are close to zero) and the deviations of the two-way effects from zero are small. $\nobreak\hfill$ $\qed$
Rather than conducting a formal analysis for the comparison of $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$ with other estimators, we provide numerical illustrations in the Monte Carlo experiments in Section (ref).
Next, we study conditions for the asymptotic optimality in Theorem (ref) to hold when the regressor coefficients $ \boldsymbol{\gamma} $ are estimated. Let $\tilde{ \boldsymbol{\gamma} }$ be an estimator of $ \boldsymbol{\gamma} $. Define all the estimators exactly the same as before with $ \boldsymbol{Y} := \boldsymbol{Y} ^*-X \boldsymbol{\gamma} $ replaced by $ \boldsymbol{Y} := \boldsymbol{Y} ^*-X\tilde{ \boldsymbol{\gamma} }$. In particular, let $\tilde{ \boldsymbol{\theta} }( \boldsymbol{\lambda} )$, $\widetilde{\text{URE}}( \boldsymbol{\lambda} )$, and $ \boldsymbol{\lambda} ^{\widetilde{\text{ure}}}$ be the counterparts of $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} )$, $\text{URE}( \boldsymbol{\lambda} )$, and $ \boldsymbol{\lambda} ^{\text{ure}}$, respectively. We provide conditions under which the impact of $\tilde{ \boldsymbol{\gamma} }$ is negligible uniformly on both the compound loss and the URE, e.g., $ \sup_{ \boldsymbol{\lambda} \in\mathcal{J}} | l_w(\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} ) - l_w(\tilde{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ), \boldsymbol{\theta} ) | \rightarrow_p 0,$ and $ \sup_{ \boldsymbol{\lambda} \in\mathcal{J}} | \text{URE}( \boldsymbol{\lambda} ) - \widetilde{\text{URE}}( \boldsymbol{\lambda} ) | \rightarrow_p 0 $ in order to extend the asymptotic optimality to the EB estimator $\tilde{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\widetilde{\text{ure}}})$. To this end, we impose the following assumption on the regressors and the estimator $\tilde{ \boldsymbol{\gamma} }$. In this case, $\text{URE}( \boldsymbol{\lambda} )$ can be interpreted as an asymptotic unbiased risk estimate.
We suggest using the $r^{1/2}$-consistent estimators of $\gamma$ studied by Verdier2020 in sparsely connected two-way model. With strictly exogenous regressors, this holds for the standard OLS estimator. When regressors are only sequentially exogenous, e.g., the lagged test score, Verdier2020 propose a $r^{1/2}$-consistent estimator of $ \boldsymbol{\gamma} $ by extending the recursive orthogonal transformation from the one-way fixed effect model to the two-way fixed effects model. The following Corollary is the generalization of Theorem (ref) to estimated regressors.
The estimators are introduced in Section (ref), the simulation designs are described in Section (ref) and the results are summarized in Section (ref).
We study the properties of the proposed estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ure}})$, henceforth EB-URE, vis-\`{a}-vis the infeasible oracle estimator $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{ol}})$, henceforth OL, the least squares estimator $\hat{ \boldsymbol{\theta} }^{\text{ls}}$, henceforth LS, and $\hat{ \boldsymbol{\theta} }( \boldsymbol{\lambda} ^{\text{mle}})$, henceforth EB-MLE, that selects $ \boldsymbol{\lambda} $ to maximize the marginal likelihood of $ \boldsymbol{Y} $, i.e., the likelihood that integrates out $ \boldsymbol{\theta} $ according to (ref). The computational details for these estimators are detailed in Section (ref) of the Supplementary Online Appendix. For simulations that compare the estimation of only $ \boldsymbol{\beta} $, we also include the one-way fixed effects estimator $\hat{ \boldsymbol{\beta} }^{\text{1-way}}(\lambda_b^{\text{mom}})$ defined in ((ref)), henceforth EB-1way.
The Monte Carlo designs are based on a prototypical application of estimating teacher value-added based on a linked student-teacher data set. We assume that there are $r=40,000$ students, $c=4,000$ teachers, and students and teachers are allocated to $s=200$ schools, numbers that are of the same order as those in the empirical application. The time dimension of the panel is $T=2$. The designs begin with the generation of the unit-specific parameters $ \boldsymbol{\theta} $, then a bipartite graph ${\cal G}$ is generated conditional on $ \boldsymbol{\theta} $, which in turn determines the matrix $B$.
We start from perfect assortative matching and then make random re-assignments controlled by the parameter $\pi_{\text{match}}$ to break the matching: $\pi_{\text{match}}=1$ ($\pi_{\text{match}}=0$) corresponds to perfect positive matching (no matching). If we skip the period $t=2$ step of teachers switching schools, then the resulting graph ${\cal G}$ will be disconnected. While there would be a potentially connected subgraph ${\cal G}_l$ for each school, these subgraphs are disconnected from each other. Thus, we can use $\pi_{\text{mob}}$ to control the connectivity of ${\cal G}$. The outcomes for periods $t=1,2$ are determined according to
The simulation designs are summarized in Table (ref). Design 1 is the reference design, where we calibrate the DGP parameters such that the empirical moments of $\hat{ \boldsymbol{\theta} }^{\text{ls}}$ and connectivity of projected teacher graph match those estimated from the empirical application in Section (ref). The subsequent designs then each perturb Design 1 along a single dimension: Design 2 assumes a relatively high level of student-teacher assortative matching and Design 3 assumes a relatively high level of teacher mobility across schools. Finally, Design 4 inverts the heterogeneity in student abilities and teacher value-added, so that students are relatively homogeneous. Under all designs, the estimators are evaluated according to their RMSEs defined as
where $w_b$ represents the usage of $W_b$ in the loss as defined in (ref).\footnote{Section (ref) of the Supplementary Online Appendix repeats the simulations for Design 1 using the RMSE $ \sqrt{ l_{w_{a+b}}(\hat{ \boldsymbol{\theta} }, \boldsymbol{\theta} ) }$ and also considers alternative specifications of the distributions governing the unit-specific parameters $ \boldsymbol{\alpha} $ and $ \boldsymbol{\beta} $ beyond Gaussianity, such as with skewness or fat tails.}
Figure (ref) shows the distribution of RMSEs of estimators across Monte Carlo repetitions under different designs. EB-URE tracks the RMSE of the infeasible benchmark OL very closely across all designs, regardless of connectivity strength, matching intensity, or heterogeneity magnitudes. In particular, this is true for Design 1 which is calibrated towards the empirical application and provides justification for implementing the EB-URE in the application. Table (ref) lists the selected hyperparameters and shows that EB-URE and OL are based on similar hyperparameter values. This is a reflection of the URE objective function providing a good approximation to the infeasible loss function. The EB-MLE, on the other hand, which relies on a correct specification of the prior (ref), is suboptimal in most designs. The suboptimality is particularly pronounced under strong matching (Design 2). This suggests a higher degree of assortative matching renders (ref), which motivates our class of estimators $ \boldsymbol{\beta} ( \boldsymbol{\lambda} )$, less representative of the true distribution of $ \boldsymbol{\theta} | B$. The proposed EB-URE instead targets a purely frequentist criterion, and thus continues to perform well relative to the benchmark OL in the presence of a misspecified prior.
We also see in Design 1 that the LS estimator performs poorly relative to OL. This is unsurprising, because the weak connectivity of the projected teacher graph results in a large variance of the LS estimators. It is therefore desirable to induce some level of shrinkage to optimize the bias-variance trade-off. However, Figure (ref) also suggests that the method used in optimizing this trade-off is crucial. The feasible EB-URE attains at least 50% RMSE reduction relative to LS across most designs. On the other hand, the performance of EB-1way is highly dependent on the underlying DGP. For instance, with even a weak positive student-teacher matching as in Design 1, the EB-1way estimator only performs marginally better than the LS. If we have a strong level of matching as in Design 2, then its performance may be even worse than the LS. In fact, its RMSE there is uniformly greater than $0.5$. By assuming $\lambda_a=\infty$, EB-1way is forced by the student-teacher matching to wrongly attribute variability in $ \boldsymbol{\alpha} $ as variability in $ \boldsymbol{\beta} $. It thus selects a smaller $\lambda_b$ than is optimal; see Table (ref). Only if the student heterogeneity is small relative to the error variance $\sigma^2$ as in Design 4, the EB-1way performs nearly optimally.
We now examine the effects of weak connectivity on the three types of estimates at a more granular level. Figure (ref) shows scatter plots of the teacher value-added $\beta_j$ ($x$-axis) and the average matched student ability for each teacher ($y$-axis), defined as
where $d_{b,j}$ is the total number of matches for teacher $j$. We provide scatter plots of the true pairs $(\beta_j,\mu_j)$, $j=1,\ldots,c$ and their estimates (LS, EB-URE, and EB-MLE) for Designs 1 and 2. The true values depicted in the top-left panels for both designs have positive slopes, reflecting the positive assortative matching patterns underlying the DGPs. The top-right panels show the same scatter plot for the LS estimates. The weak connectivity results in a bias large enough that the implied correlation between $(\beta_j,\mu_j)$ is actually highly negative for Design 1. The bottom-left panels show the EB-URE estimates of $(\beta_j,\mu_j)$. Even though the hyperparameter selection does not target the estimation loss of $ \boldsymbol{\alpha} $, the estimates indicate a strong correlation between $\beta_j$ and $\mu_j$, meaning that the EB-URE estimator captures a significant portion of matching.
Table (ref) provides further evidence through empirical moments computed from the true effects and their estimates. We report the variances of $\alpha_i$ and $\beta_j$, respectively, and the correlations between $\mu_j$ and $\beta_j$. For instance, under Design 1 the empirical correlation between the true effects is 0.15, whereas the correlation between the EB-URE estimates is 0.19. For Design 2 these values change to 0.46 and 0.50, respectively. The empirical variances of $\alpha_i$ and $\beta_j$ computed from the true and their estimates are also very close. The empirical moments computed from the EB-MLE estimates exhibit slightly larger discrepancies. For instance, for Designs 2 and 3 the correlations among the true effects are 0.46 and 0.14, whereas the correlation among the EB-MLE estimates are 0.56 and 0.28, respectively. We conclude from the simulations that the EB-URE estimates not just generate low RMSEs, but they also are able to reproduce key empirical moments of the underlying true effects.
We utilize a matched student-teacher dataset from the North Carolina Education Research Data Center with observations on students from grades three to five for the years 2017 and 2018 period to estimate teacher value-added. For identification purposes, we restrict the dataset to the largest connected component of the student-teacher graph. We remove students that appear only once in the data set, since they are not relevant for the identification of the teacher value-added parameters, students that repeat grades, and students with special accommodations. This leaves us with $r=41,243$ students, $c=5,332$ teachers, and $s=258$ schools. Because we have test scores for two consecutive years, this leads to 82,486 observations in total.
Starting point of the empirical analysis is model (ref). We take the outcome variable $y_{it}$ to be a math test score, standardized within each (year, grade) cell in accordance with the teacher value-added literature to have a mean of zero and variance of one. We only include a single regressor $ \boldsymbol{x} _{it}$, namely the student's lagged test score $y_{it-1}$.\footnote{We decided to keep the model simple. As a robustness exercise we repeat the analysis by also including cubic polynomials in age and class size; see Section (ref) of the Supplementary Online Appendix. The resulting teacher value-added estimates have a correlation of $0.989$ with the ones reported here.} We apply the transformation in ((ref)) by subtracting $X\tilde{ \boldsymbol{\gamma} }$ from the outcome variable. Because $y_{it-1}$ is only sequentially and not strictly exogenous, we use Verdier2020's estimator to obtain $\tilde{ \boldsymbol{\gamma} }=0.067$. The estimator intuitively extends the recursive orthogonal transformation in one-way to two-way effects models for consistent estimation; see Section (ref) of the Supplementary Online Appendix for more details. The subsequent analysis also conditions on the error variance estimate $\hat{\sigma}^2 = 0.12$. In Section (ref) we report measures of connectivity for our data set and discuss some features of the LS estimates. The EB-URE estimates are presented in Section (ref) and in Section (ref) we compare teacher rankings based on value-added estimates from different estimators.
We previously emphasized that the performance of estimators for the two-way effects model depends on the connectivity of the student-teacher graph. Table (ref) provides information about the distribution of eigenvalues of the Laplacian $L_{2,\perp}$ of the projected teacher graph. The smallest eigenvalue is 0.013, the first percentile is 0.069, and the fifth percentile is 0.23. This distribution appears to be broadly in line with our Assumption (ref) Section (ref) that only a finite number of the eigenvalues can be small, but not too small. When scaled by $c^{1/2}$ as in Assumption (ref) the smallest non-zero eigenvalue is 0.96 and the first percentile is 5.1. We interpret these numbers as evidence of weak connectivity: many subsets of teachers share few common students, in particular teachers from different schools.
Weak connectivity generates instability and imprecision of LS estimates. To illustrate the instability, we compare the LS estimator $\hat{ \boldsymbol{\beta} }^{\text{ls}}$ to an estimator $\tilde{ \boldsymbol{\beta} }^{\text{ls}}$ that is obtained by replacing $\alpha_i$ by an index function $ \boldsymbol{\alpha} _x' \boldsymbol{x} _i$, where the regressors $ \boldsymbol{x} _i$ are time-invariant demographic variables: gender, ethnicity, economically disadvantaged, English learner; see Section (ref) of the Supplementary Online Appendix for more details. The left panel of Figure (ref) shows a scatter plot of the two estimators and documents that $\hat{\beta}_j^{\text{ls}}$ is much more variable than $\tilde{\beta}_j^{\text{ls}}$. The empirical variances of $\hat{\beta}_j^{\text{ls}}$ and of $\tilde{\beta}_j^{\text{ls}}$ are $0.15$ and $0.063$, respectively. The right panel shows a scatter plot of $\hat{\beta}_j^{\text{ls}}$ and the average ability of students taught by teacher $j$, previously denoted by $\hat{\mu}_j^{\text{ls}}$ and defined in (ref). The slope of a linear regression line is -0.61. While this pattern seems to suggest negative assortative matching between students and teachers, it is likely an artifact of the graph's weak connectivity. In fact, in the Monte Carlo simulation we observed a similar pattern in the second and fourth top-row panels of Figure (ref), where it was a limited mobility bias.
To construct our proposed EB estimator, we adopt the covariate-based prior of Section (ref) where the student effects $ \boldsymbol{\alpha} $ are shrunk to the location $Z_a \boldsymbol{\delta} _a$ and $Z_a$ contains the same student demographics as $ \boldsymbol{x} _i$ in the definition of $\tilde{ \boldsymbol{\beta} }^{\text{ls}}$ above. We shrink the teacher effects to $0$, i.e., $ \boldsymbol{\delta} _b=0$ under the normalization. The parameters $(\lambda_\alpha,\lambda_\beta,\phi, \boldsymbol{\delta} _a)$ are jointly determined via URE-minimization using $W=W_b$. This produces our proposed EB-URE estimates that targets the $ \boldsymbol{\beta} $ estimation loss. The estimated scale hyperparameters are $(\lambda_a^{\text{ure}},\lambda_b^{\text{ure}},\phi^{\text{ure}})=(0.075,1.67,0.7)$. The estimate of $ \boldsymbol{\delta} _a$ is reported in Section (ref) of the Supplementary Online Appendix.
The left panel of Figure (ref) plots $\hat{\beta}_j^{\text{ure}}$ against $\hat{\mu}^{\text{ure}}_j$. The regression coefficient is now $0.533$, which implies positive student-teacher matching as opposed to the negative matching implied by the LS estimates. The magnitude of matching patterns is obscured by the difference in the scale of variances in the $\alpha$ and $\beta$ dimension. Thus, we also provide the empirical correlation of $\hat{\alpha}_i^{\text{ure}}$ and $\hat{\beta}_{j(i,t)}^{\text{ure}}$ in Table (ref), which is computed to be $0.12$. In combination with the right panel of Figure (ref) the empirical results are consistent with the Monte Carlo results in Section (ref), demonstrating that the EB-URE is able to accurately pick up on matching patterns otherwise undetected by the LS estimates due to weak connectivity. We also plot the school averages of $\hat{\beta}_j^{\text{ure}}$ against $\hat{\alpha}_j^{\text{ure}}$ in the right panel of Figure (ref) with a regression coefficient of $0.698$ and empirical correlation of $0.188$, suggesting that part of the observed positive student-teacher matching pattern arises from assortative matching between schools.
Figure (ref) plots the EB-URE estimates against LS estimates and against $\hat{ \boldsymbol{\beta} }^{\text{1-way}}$, the standard shrinkage estimator from a one-way model that assumes $\alpha_i = \boldsymbol{\alpha} _x' \boldsymbol{x} _i$. The latter is commonly used in the literature and a simplified version was previously defined in ((ref)). Section (ref) of the Supplementary Online Appendix provides a detailed explanation of its construction for this application. We refer to it as EB-1way consistent with the naming in the simulations of Section (ref). While the EB-URE and LS estimators are positively correlated (left panel of the figure), by construction the EB-URE estimator is more stable and has a five-times smaller variance. Numerical values for the sample variances are reported in Table (ref).
The right panel of Figure (ref) shows that the EB-1way estimates are less volatile than the LS estimates, but still more volatile than the two-way EB-URE estimates. To explore the extent to which the demographic regressors $ \boldsymbol{x} _i$ explain the student heterogeneity, we project $\hat{\alpha}_i^{\text{ls}}$ onto $ \boldsymbol{x} _i$, where we retrieve an $R^2$ of $0.21$. This is relatively small and indicates that most of the variation in student heterogeneity is not captured by the $ \boldsymbol{x} _i$ regressors. In conjunction with the usual concern regarding assortative matching, this provides an additional justification for adopting the two-way effects model.
A frequent use of value-added estimates is to rank teachers for remuneration purposes. We conclude the empirical analysis with a comparison of the rankings based on $\hat{ \boldsymbol{\beta} }^{\text{ls}}$, $\hat{ \boldsymbol{\beta} }^{\text{ure}}$ and $\hat{ \boldsymbol{\beta} }^{\text{1-way}}$. For each set of teacher value-added estimates, teachers are split into quintiles. We then take two sets of estimates and compute the number of teachers for each of the 25 quintile pairs. The results are reported in Table (ref). For instance, the number of teachers who are in the first quntile of the $\hat{ \boldsymbol{\beta} }^{\text{ls}}$ distribution and also in the second quintile of $\hat{ \boldsymbol{\beta} }^{\text{ure}}$ distribution is 253. Note that the number of teachers per quintile is approximately 1,066. The diagonals thus represent the number of teachers for which two estimators agree in terms of quintile rank, the principal off-diagonals represent the numbers of disagreements by one rank, and so on.
There is substantial variation in how $\hat{ \boldsymbol{\beta} }^{\text{ls}}$ and $\hat{ \boldsymbol{\beta} }^{\text{ure}}$ rank teachers. In particular, there are 39 teachers (the bottom left and top right corners of the left panel) who are simultaneously ranked at opposite quintiles by $\hat{ \boldsymbol{\beta} }^{\text{ure}}$ and $\hat{ \boldsymbol{\beta} }^{\text{ls}}$. While the rankings within the first and fifth quintiles are relatively more similar across $\hat{ \boldsymbol{\beta} }^{\text{1-way}}$ and $\hat{ \boldsymbol{\beta} }^{\text{ure}}$ at the extreme quintiles, there is still substantial variation in the middle quintiles. In view of the theoretical optimality of $\hat{ \boldsymbol{\beta} }^{\text{ure}}$ and its adaptivity that were demonstrated in the simulations calibrated to fit this application (Design 1), we believe that EB-URE is an attractive alternative to these estimators.
We develop an empirical Bayes estimator for two-way effects models using a novel prior that incorporates assortative matching patterns based on the graph normalized Laplacian matrix and hyperparameter selection based on the minimization of an unbiased risk estimate. We prove its asymptotic optimality, allowing for weakly connected graphs due to limited mobility. We show in a Monte Carlo study that the estimator dominates, in terms of RMSE, a number of competitors and is able to capture assortative matching. An application of our estimator to the NCERDC data set suggests that there is substantial student-teacher assortative matching which is ignored by both the common one-way shrinkage estimator and the two-way LS estimator. We also find that teacher rankings based on value-added estimates are quite sensitive to the methodology employed.
\setstretch{1.1} \setstretch{1.3}
\setstretch{1.3}
\setcounter{page}{1}
Consider the LS estimate of the error variance in (ref). The following Lemma shows that it is unbiased and consistent.
The following lemma shows that the impact of the estimated coefficient is negligible uniformly over the hyperparameter. It is the key result to deliver asymptotic optimality with the estimated coefficient in Corollary (ref).