Fixed-Effect Regressions on Network Data

FIXED-EFFECT REGRESSIONS ON NETWORK DATA

\def\spacingset#1 \spacingset{1}

\if11 \fi

\if01 {

} \fi

\spacingset{1.45}

\setcounter{equation}{0}

Introduction

Data on the interaction between agents are in increasing supply. A workhorse technique to analyze such data is a linear regression model with agent-specific parameters. It has been used to investigate a variety of questions. For example, application of a two-way regression model to matched employer-employee data decomposes (log) wages into worker heterogeneity, firm heterogeneity, and residual variation. Following AbowdKramarzMargolis1999 the correlation between the estimated worker and firm effects is regarded as a measure of assortative matching. A positive correlation indicates that high quality workers are employed in more productive firms. Using the same decomposition, CardHeiningKline2013 study to what extent the evolution of wage inequality is due to changes in the variance of worker and firm heterogeneity. Nimczik2018 reports the whole distribution of the estimated worker and firm effects. In a similar fashion, the literature on student achievement backs out student and teacher effects from test score data. The estimated teacher heterogeneity is interpreted as teacher value-added and their variance as a measure of their importance (see JacksonRockoffStaiger2014 for an overview of this literature). These estimates are used to assess teachers and are important inputs to personnel evaluations and merit pay programs Rothstein2010.\footnote{Fixed-effect regressions of this kind are now part of the standard toolkit of many empiricists in a variety of different areas. FinkelsteinGentzkowWilliams2016 and AmitiWeinstein2018 use them to separate supply and demand factors in healthcare utilization from data on patient migration, and in firm investment behavior from financial data on banks loans, respectively. ChettyHendren2018 evaluate the importance of growing up in a specific neighborhood on labor market outcomes later on in life.}

In spite of their widespread use, there is little to no work on the theoretical properties of such fixed-effect approaches. In fact, the few results that are available point to issues of downward bias in the estimation of the correlation between worker and firm effects, finding a spurious negative correlation in many data sets AndrewsGillSchankUpward2008,AndrewsGillSchankUpward2012, and upward bias in the estimator of the variance of teacher effects Rockoff2004. The presence of bias here is not surprising. Indeed, the individual effects are estimated with noise. Their sampling error then introduces bias in the estimator of nonlinear functionals. A more complicated issue is the assessment of the statistical precision with which the fixed effects are estimated and, more generally, the development of distribution theory. This is important as it allows to establish conditions for consistency and rates of convergence, and yields insight into whether standard test statistics can be expected to be approximately size correct and have non-trivial power. None of these issues has been addressed so far. Providing such theory so is not only relevant for inference on the individual effects and their moments, but may also serve as a stepping stone to address related problems. For example, without theory for the fixed-effect estimator the behavior of the falsification test for value-added models of Rothstein2010 remains unknown, and correct standard errors for regressions of outcomes on estimated fixed effects KettemannMuellerZweimuller2017 cannot be derived.

The data structure arising from interactions between agents is different from that of standard cross-section or panel data. It is typically difficult to see how the data carry information about certain parameters. In this paper we present sufficient conditions for consistency and asymptotic normality of least-squares estimators of fixed effects in linear regression models. We see the data as a network and represent it by a graph where agents are vertices and edges between vertices are present if these agents interact. It is intuitive that the structure of this graph should be a key determinant of the accuracy of statistical inference. We formalize this here. Our setup places no a priori restrictions on the graph structure and our results apply to both dense and sparse settings. A data structure of particular importance is that of a bipartite graph. Here, the data concerns two types of individuals and interactions only occur between the types but not within each type. This is the case in our motivating examples above and we treat this bipartite case in detail. In fact, while we deal with general graphs, our regression setup is designed to capture the main features of the prototypical two-way regression model. We focus on inference on the individual effects but our results also serve as a stepping stone for the analysis of estimators of other parameters, such as the variance and other moments of (the distribution of) the individual effects, and we provide some results on these as well. We do not discuss inference on common slope coefficients. In contemporary work, Verdier2016 provides such results for two-way regression models. The fixed-effect model for test scores, for example, can be used to assess the effect of class-size reductions on student achievement while controlling for student and teacher heterogeneity.

The ability to accurately estimate the individual effect of a given vertex depends on how well this vertex is connected to the rest of the network. Our theory involves both global and local measures of network connectivity. The main global connectivity measure we use is the smallest non-zero eigenvalue of the (normalized) Laplacian matrix of the graph.\footnote{The Laplacian matrix is similar to the adjacency matrix as a devise to represent a graph and can be obtained from it. Both matrices are formally defined below. Eigenvalues and eigenvectors of these and related matrices have also been found of use in determining equilibrium conditions in games played on networks BramoulleKrantonDamours2014 and in (statistical) community detection SchiebingerWainwrightYu2015.} It reflects how easy it is to disconnect a network by removing edges from it. The other measures of connectivity that we use are the degrees of the vertices as well as various harmonic means thereof. All of these measures arise naturally when studying the variance of the fixed-effect estimator. We highlight the interplay between them in deriving conditions for consistent estimation and for inference based on standard first-order asymptotics to be possible. As the network grows the smallest eigenvalue may approach zero, and so the graph may become more sparse, provided the relevant harmonic mean grows sufficiently fast. These findings mimic conditions on the bandwidth in nonparametric estimation problems, although they will typically show up in second-order terms here. This explains why estimation at the parametric rate may be feasible even in sparse networks. Our results also show that inference on averages over the individual effects is more demanding on the network structure and, even after bias reduction, may only be feasible in quite dense networks.

Our analysis shows that it is useful to inspect measures of global and local connectivity when interpreting estimation results from network data. We do so here for two data sets. The first is a large network of teachers in elementary schools in North Carolina, where the object of interest would be teacher value-added. This is arguably one of the most important applications of the two-way regression model. This graph is only very weakly connected and our theory does not support the use of large-sample arguments. When a simple model with homoskedastic errors is applied to these data standard errors based on conventional first-order approximations for teacher value-added are, on average, about $40\%$ smaller than the actual standard deviations. Further, the sample variance of the estimated teacher effects has a substantial upward bias. This bias translates into an overly optimistic view on the ability of teacher value-added to explain variation in test scores. To provide an example of a data set that yields a much stronger connected graph we also construct an occupational network from the British Household Panel Survey (BHPS). This graph would arise in the context of wage regressions with occupational dummies, for example. Here, our connectivity measures are much more supportive of standard inferential approaches and, indeed, again in a simple model, we find that conventional first-order approximationa are quite accurate.

The remainder of this paper is organized as follows. Section 2 details the structure of the data under study and introduces the regression model of interest. Special attention is given to the bipartite graph and the two-way regression model. Section 3 provides distribution theory for the least-squares estimator of the individual effects and also discusses estimation of their moments. Section 4 contains details on our two numerical illustrations. Section 5 concludes. The supplementary material to this paper contains some additional results and illustrations, as well as the proofs of all theorems.

Regression analysis of network data

Data structure

Consider an undirected graph $\mathcal{G}:=\mathcal{G}(V,E)$ where $m := \lvert E \rvert$ edges are placed between $n := \lvert V \rvert$ vertices. We allow for multiple edges between vertices (i.e., $\mathcal{G}$ can be a multigraph) and the edges may be assigned a weight. We do not consider loops (i.e., no edge connects a vertex with itself). Without loss of generality we label the vertices by natural numbers, so that $V$ is $\{1,\ldots, n\}$. The multiset $E$ contains the $m>0$ unordered pairs $(i,j)$ from the product set $V\times V$ that are linked by an edge, possibly with repetition. The same pair $(i,j)$ will appear multiple times in $E$ if they share more than one edge; we let ${E}_{(i,j)}\subset E$ denote the set of edges between them. We have $E_{(i,j)}=E_{(j,i)}$ and may have $E_{(i,j)}=\varnothing$. We will label the edges by natural numbers; so, each edge $e$ edge has assigned to it an integer $\varepsilon_e\in \lbrace 1,\ldots,m \rbrace$. For later use we note that vertices $i$ and $j$ are said to be connected if $\mathcal{G}$ contains a path from $i$ to $j$, and that the graph $\mathcal{G}$ is said to be connected if every pair of vertices in the graph is connected.

For an edge $e\in E$ let $w_e>0$ be its weight. An unweighted graph has $w_e=1$ for all $e\in E$. The graph $\mathcal{G}$ may be represented by its $m\times n$ (oriented) incidence matrix $\mathbold{B}$, with entries

Here, the choice of sign gives each edge $e$ an orientation. As will become apparent, our analysis below is invariant to this choice of orientation. The graph may also be represented through its $n\times n$ adjacency matrix $\mathbold{A}$, which has elements $$ (\mathbold{A})_{ij} := \sum_{e\in E_{(i,j)}} w_e . $$ The incidence matrix and adjacency matrix are related through the $n\times n$ Laplacian matrix $\mathbold{L}$, as $$ \mathbold{L}:= \mathbold{B}^\prime\mathbold{B} = \mathbold{D}-\mathbold{A}, $$ for $\mathbold{D}:=\mathrm{diag}(d_1,\ldots,d_n)$ the diagonal $n\times n$ (weighted) degree matrix, where the degree of vertex $i$ is $$ d_i : = \sum_{j=1}^n (\mathbold{A})_{ij}. $$ When $\mathcal{G}$ is an unweighted graph, for example, $d_i$ equals the number of edges that involve vertex $i$. For a vertex $i$ we will let $ [i] := \lbrace j\in V: E_{(i,j)}\neq \varnothing \rbrace $ denote the set of its direct neighbors. Observe that $d_i$ may be large even if $i$ has few neighbors---i.e, when $\lvert [i] \rvert$ is small---as the edge weights $(\mathbold{A})_{ij}$ for $j\in [i]$ may be large. An example is a multigraph where many edges exist between $i$ and some $j\in [i]$.

Regression model and least-squares estimator

Now, given a graph $\mathcal{G}$, for each edge $e\in E$ we observe an outcome $y_{\varepsilon_e}$ and a $p$-vector of covariates $\mathbold{x}_{\varepsilon_e}$. Allowing $\mathcal{G}$ to be a multigraph covers the (unbalanced) panel data case, where multiple outcomes are available for some vertex pairs. Collect all outcomes in the $m$-vector $\mathbold{y}$ and all covariates in the $m\times p$ matrix $\mathbold{X}$. Let $\mathbold{\alpha}:=(\alpha_1,\ldots,\alpha_n)^\prime$ be an $n$-vector of vertex-specific parameters and let $\mathbold{\beta}:=(\beta_1,\ldots,\beta_p)^\prime$ be a $p$-vector of regression slopes. Our interest lies in estimating the model

where $\mathbold{u}$ is an $m$-vector of regression errors.\footnote{A change of edge orientation corresponds to a sign flip in the corresponding outcome and regressor matrices. This does not affect the least squares estimator of (ref).} We will treat $\mathbold{B}$ and $\mathbold{X}$ as fixed throughout. This implies that we consider the network as non-random and exogenous.\footnote{Exogeneity of the network is the standard assumption in the literature building on AbowdKramarzMargolis1999. Accounting for endogenous network formation requires more complicated models and has started to receive some attention; see BonhommeLamadonManresa2015 and LentzPiyapromdeeRobin2018. } Our focus is on the vector $\mathbold{\alpha}$. In the two-way regression model of our motivating examples these are the worker and firm effects or the student and teacher effects, respectively. By the definition of $\mathbold{B}$ one of these effects will enter (ref) with a minus sign. While this may appear to be an unusual convention from an applied perspective, it is convenient for our theoretical analysis. As will be explained below, this sign convention is without loss of generality in the two-way model, where the underlying graph is bipartite.

In (ref) the outcomes for a given vertex pair $(i,j)$ depend on the individual effects through their difference $\alpha_i - \alpha_j$. This implies that our model is overparameterized. Indeed, we have $\mathbold{B} \mathbold{\iota}_n =0$, where $\mathbold{\iota}_n:=(1,\ldots,1)^\prime$ is the $n$-vector of ones, as each row of $\mathbold{B}$ sums up to zero. It follows that the mean of the vertex-specific parameters cannot be learned from the data and a normalization is required. We impose that

which will prove a convenient choice for our purposes. Denoting the degree vector by $\mathbold{d}:=(d_1,\ldots, d_n)^\prime$ we may write the constraint (ref) compactly as $\mathbold{d}^\prime \mathbold{\alpha}=0$. A normalization can be dispensed with if interest lies in parameter differences, i.e, $\alpha_i-\alpha_j$, as in FinkelsteinGentzkowWilliams2016, for example. Results for such differences that parallel those developed below are given in the supplementary material.

The standard estimator of $\mathbold{\alpha}$ is the constrained least-squares estimator

where $\lVert\cdot \rVert$ denotes the Euclidean norm, $ \mathbold{M}_{\mathbold{X}}:= \mathbold{I}_m - \mathbold{X}(\mathbold{X}^\prime \mathbold{X})^{-1} \mathbold{X}^\prime, $ and $\mathbold{I}_m$ is the identity matrix of dimension $m \times m$ . The following theorem gives conditions under which this estimator exists and is unique. For any matrix $\mathbold{C}$ we denote its Moore-Penrose pseudoinverse by $\mathbold{C}^\dagger$. When $\mathbold{C}$ is $n\times n$ we let $\mathbold{C}^\star := \mathbold{D}^{-1/2} \left( \mathbold{D}^{-1/2} \mathbold{C} \mathbold{D}^{-1/2}\right)^\dagger \mathbold{D}^{-1/2}$. It is easily shown that $\mathbold{C} \mathbold{C}^\star \mathbold{C} = \mathbold{C}$ and $ \mathbold{C}^\star \mathbold{C} \mathbold{C}^\star = \mathbold{C}^\star $. Therefore, $\mathbold{C}^\star$ is a pseudoinverse of $\mathbold{C}$.

The need for a pseudoinverse arises because $\mathbold{B}^\prime \mathbold{M}_{\mathbold{X}} \mathbold{B}$ is singular, which follows from the fact that $\mathbold{B} \mathbold{\iota}_n =0$. The use of the particular pseudoinverse $(\mathbold{B}^\prime\mathbold{M}_{\mathbold{X}}\mathbold{B})^{\star}$ is a consequence of our normalization $\mathbold{d}^\prime \mathbold{\alpha} = 0$. A change of normalization would imply a different pseudoinverse in the statement of Theorem (ref). The result of the theorem is intuitive and generalizes results in the literature on matched employer-employee data (AbowdCreecyKramarz2002). When the graph $\mathcal{G}$ is disconnected a separate normalization of the form in (ref) is needed for each connected component of $\mathcal{G}$. Our results then apply to each of these subgraphs. In practice, the analysis is typically confined to the largest connected component of $\mathcal{G}$ (see, for example, CardHeiningKline2013).

While $\check{\mathbold{\alpha}}$ is routinely used its statistical properties are not well understood. Our aim here is to shed light on how the structure of the network $\mathcal{G}$ affects its sampling behavior and, with it, the reliability of standard inferential procedures based on $\check{\mathbold{\alpha}}$. For our analysis edge-specific covariates mostly complicate notation and presentation. It will on occasion be convenient to first analyze (ref) when $\mathbold{\beta}$ is treated as known and the outcome vector is redefined as $\mathbold{y}-\mathbold{X}\mathbold{\beta}$. Then $$ {\mathbold{\hat{\alpha}}} := (\mathbold{B}^\prime\mathbold{B})^{\star} \mathbold{B}^\prime (\mathbold{y}-\mathbold{X}\mathbold{\beta} ) $$ is the least-squares estimator of $\mathbold{\alpha}$ subject to (ref). To appreciate how the structure of $\mathcal{G}$ relates to our problem of estimating the parameter $\mathbold{\alpha}$ suppose first that $\mathbold{u}\sim(\boldsymbol{0},\sigma^2 \mathbold{I}_m)$. Then

So, up to a scale factor, the variance of $\hat{\mathbold{\alpha}}$ is completely determined by the Laplacian of $\mathcal{G}$. If, in addition, we were to assume that $\mathbold{u}\sim N(\boldsymbol{0},\sigma^2\mathbold{I}_m)$ we would be in the classical regression setting and, given unbiasedness of $\hat{\mathbold{\alpha}}$, size-correct inference could be performed for any sample size. It is not clear, however, how one should proceed with non-classical regression errors.

The validity of standard large-sample arguments is not immediate here. From (ref) we have

where $$ \mathbold{S}:=\mathbold{D}^{-1/2} \mathbold{L} \mathbold{D}^{-1/2} = \mathbold{I}_n - \mathbold{D}^{-1/2} \mathbold{A} \mathbold{D}^{-1/2} $$ is the normalized Laplacian. Equation (ref) follows from the fact that $\mathbold{L}^*=\mathbold{D}^{-1/2}\mathbold{S}^\dagger\mathbold{D}^{-1/2}$, such that $(\mathbold{L}^\star)_{ii} = (\mathbold{S}^\dagger)_{ii}/d_i$.\footnote{Our choice of normalization (ref) guarantees the appearance of the Moore-Penrose pseudoinverse of $\mathbold{S}$ in $\mathbold{L}^*$, which is the main reason for that choice. We are grateful to Nadine Geiger for pointing out an inconsistency in our normalization in an earlier version of this paper. } While (ref) shows the importance of the sample size in the variance of $\hat{\alpha}_i$ through the presence of the degree $d_i$, it does not imply that ${\rm var}( \hat{\alpha}_i )$ shrinks as $d_i\rightarrow \infty$, nor would it give a convergence rate if it did. This is because the normalized Laplacian $\mathbold{S}$ also changes when $\mathcal{G}$ grows.

Let $\lambda_1\leq \lambda_2\leq \cdots \leq \lambda_n$ be the eigenvalues of $\mathbold{S}$. The spectrum of Laplacian matrices is well studied (see, e.g., Chung1997). We have $0 \leq \lambda_i \leq 2$ for all $i$. We always have that $\lambda_1=0$, with $\mathbold{\iota}_n$ as eigenvector. The number of zero eigenvalues $\lambda_i$ equals the number of connected components in $\mathcal{G}$. Hence, if $\mathcal{G}$ is connected, then $\lambda_2>0$ is the smallest non-zero eigenvalue of the normalized Laplacian.

Our theory involves conditions on $\lambda_2$ and on the degree structure of the network through various harmonic means thereof. $\lambda_2$ can be seen as a measure of global connectivity of $\mathcal{G}$. To see this we note that it can be linked to the Cheeger constant, $$ C := \min_{U\in\left\{ U\subset V: \, 0<\sum_{i\in U} d_i \leq \sum_{i\notin U} d_i \right\}} \frac{\sum_{i\in U} \sum_{j\notin U } (\mathbold{A})_{ij}}{\sum_{i\in U} d_i} . $$ The constant $C\in [0,1]$ measures how difficult it is to separate $\mathcal{G}$ into two disconnected components by removing edges from it. The numerator in the definition of $C$ is the total weight of the removed edges, the denominator is the total degree in the smallest of the two components. A larger value of $C$ implies a more strongly-connected graph, and it is linked to $\lambda_2$ through the inequalities

which are due to FriedlandNabben2002. Thus, like the Cheeger constant, $\lambda_2$ is a measure of global connectivity of the graph $\mathcal{G}$. Our results below allow for $\lambda_2\rightarrow 0$ as $\mathcal{G}$ grows, and so cover situations where the graph becomes increasingly more sparse. We will give explicit rates on $\lambda_2$ for consistent estimation to be possible.

Two-way regression model on bipartite graph

To relate our model to our main motivating examples consider the case of a bipartite graph $\mathcal{G}$, i.e., $V=V_1\cup V_2$ and $V_1\cap V_2 = \varnothing$, and edges are formed only between the subsets $V_1$ and $V_2$ but not within. So, for an edge $(i,j)$ we necessarily have that $i\in V_1$ and $j\in V_2$. A bipartite graph describes the interaction between two types of units, such as workers and firms or students and teachers. The outcome of interest here would typically be (log) wages or earnings and test scores, respectively. If we have panel data, so $\mathcal{G}$ is a multigraph, we may observe workers match with different firms over time and observe students in different classrooms or across multiple subjects. In fact, in these applications, such longitudinal data are necessary for $\mathcal{G}$ to be connected. A two-way regression model for such data takes the form

where $\mathbold{\mu}:=(\mu_1,\ldots,\mu_{n_1})^\prime$ and $\mathbold{\eta}:=(\eta_1,\ldots,\eta_{n_2})^\prime$ are the $n_1:=\lvert V_1 \rvert$ and $n_2:=\lvert V_2 \rvert$ parameter vectors for the two types of units, and the $m\times n_1$ and $m\times n_2$ matrices $\mathbold{B}_1$ and $\mathbold{B}_2$ have entries

This is a workhorse specification to capture heterogeneity across units in linked data sets. It can be cast into (ref) by setting $$ \alpha_i = \left\lbrace

\right. $$ sorting the units in $V$ by type so that we can write $\mathbold{\alpha}=(\mathbold{\mu}^\prime,-\mathbold{\eta}^\prime)^\prime$, and constructing the $m\times n$ matrix $\mathbold{B}=(\mathbold{B}_1,-\mathbold{B}_2)$ by concatenation. Choosing the sign in front of $\eta_i$ is without loss of generality because links are only formed between, but never within, the subsets $V_1$ and $V_2$. The need for a normalization built-in in our general specification arises here from the fact that \eqref{eq:twoway} is invariant to reparametrizations of the form $(\mu_i,\eta_j)\mapsto (\mu_i+c, \eta_j-c)$ for any $c$.

The two-way regression model provides an interesting example where a weighted graph arises naturally. In many applications the researcher is primarily interested in learning the parameters of one type, say those $i\in V_2$. This is so in teacher value-added models, for example. There, interest lies in estimating the $n_2$ teacher effects while controlling for unobserved student-specific heterogeneity through the inclusion of $n_1$ student effects (see, e.g., JacksonRockoffStaiger2014). Partialling-out the vector $\mathbold{\mu}$ from the two-way model in (ref) gives

From standard partitioned-regression theory, the least-squares estimator of $\mathbold{\eta}$ from this equation is numerically identical to the one obtained from joint estimation of $\mathbold{\mu}$ and $\mathbold{\eta}$ in (ref). However, the formulation in (ref) is helpful in understanding the behavior of the estimator of $\mathbold{\eta}$. The properties of the matrix $\mathbold{B}_2^\prime\mathbold{M}_{\mathbold{B}_1}\mathbold{B}_2$ drive the sampling behavior of $\check{\mathbold{\eta}}$. This matrix is the Laplacian of a weighted one-mode projection Newman2010 of the bipartite graph $\mathcal{G}$ on the $n_2$ vertices in $V_2$.

It is instructive to discuss this one-mode projection in more detail and to formalize how it fits the general setup in (ref). Projecting the bipartite graph $\mathcal{G}=\mathcal{G}(V_1\cup V_2,E)$ on $V_2$ is done by suppressing the vertices in $V_1$. This gives a new (unipartite) graph, say $\mathcal{G}^\prime=\mathcal{G}(V_2,E^\prime)$. Each edge pair $(e_1,e_2)$ with $e_1\in E_{(i,j)}$ and $e_2\in E_{(i,j^\prime)}$ in $\mathcal{G}$ for some $i\in V_1$ and $j, j^\prime\in V_2$ gives rise to a single edge $e=(e_1,e_2)\in E^\prime_{(j,j^\prime)}$ in $\mathcal{G}^\prime$. In the student-teacher example, two teachers $j$ and $j^\prime$ are connected by an edge in $\mathcal{G}^\prime$ if and only if there exists at least one student $i$ that they have both taught. Alternatively, the edge $e=(e_1,e_2)\in E^\prime$ exists because $e_1\in E$ and $e_2\in E$ both connect to the same vertex $i\in V_1$. Given the edges $e_1,e_2$ this connecting vertex $i$ is unique; for later use we denote it by $c((e_1,e_2))$. In $\mathcal{G}^\prime$ we have $$ m^\prime := \lvert E^\prime \rvert = \sum_{i\in V_1} \frac{d_i(d_i-1)}{2} $$ edges; $m^\prime$ need not equal $m$ and, indeed, may be much larger. We again label $e^\prime\in E^\prime$ by natural numbers $\varepsilon^\prime_{e^\prime}$. The process of concatenating edges in $E$ to form the new edge set $E^\prime$ can be described by the $m^\prime\times m$ matrix $\mathbold{Q}$ with entries $$ (\mathbold{Q})_{\varepsilon^\prime_{e^\prime} \, \varepsilon_{e_1}} := \left\lbrace

\right. $$ Choosing the orientation of the rows of $\mathbold{Q}$ is without loss of generality. The matrix $\mathbold{Q}$ has a first-differencing interpretation. Indeed, when applied to the two-way regression model \eqref{eq:twoway} we get $ \mathbold{Q}\mathbold{y} = \mathbold{Q}\mathbold{B}_2\mathbold{\eta} + \mathbold{Q}\mathbold{X}\mathbold{\beta} + \mathbold{Q}\mathbold{u} , $ because $\mathbold{Q}\mathbold{B}_1 = 0$. Thus, $\mathbold{Q}$ sweeps out the $n_1$ nuisance parameters $\mathbold{\mu}$ and transforms original outcomes into first differences. The matrix $\mathbold{Q}\mathbold{B}_2$ is the (oriented) incidence matrix of an unweighted graph, and this first-differenced regression equation fits \eqref{eq:regression1}.\footnote{ $\mathbold{Q}\mathbold{B}_2$ will contain rows with only zero entries if there are differenced outcomes that do not depend on $\mathbold{\eta}$. This is at odds with the definition of an incidence matrix. Dropping these differences from $\mathbold{Q}\mathbold{y}$, however, restores the incidence matrix interpretation of $\mathbold{Q}\mathbold{B}_2$. This operation does not affect estimation of $\mathbold{\eta}$ and so is irrelevant for our purposes. However, the differenced outcomes may provide information on $\mathbold{\beta}$, which is why we prefer to work with $\mathbold{Q}\mathbold{B}_2$ as defined here. } Applying least-squares directly to the first differences is inefficient and is not equivalent to estimation of the two-way regression model. Ordinary least-squares estimation of \eqref{eq:profiled_eq} is numerically equivalent to weighted least-squares estimation of the first-differenced equation. The relevant $m^\prime\times m^\prime$ diagonal weight matrix $\mathbold{W}$ has entries $ (\mathbold{W})_{\varepsilon^\prime_e\, \varepsilon^\prime_e} := {1}/\sqrt{d_{c(e)}}. $ Ordinary least-squares applied to (ref) and to

yields the same result. Here, ${\mathbold{W}}\mathbold{Q}\mathbold{B}_2$ is the incidence matrix of a weighted one-mode projection of $\mathcal{G}$. This $\mathcal{G}^\prime$ determines the properties of the least-squares estimator. Its Laplacian is $$ \mathbold{L}^\prime := \mathbold{B}_2^\prime (\mathbold{Q}^\prime \mathbold{W}^2 \mathbold{Q}) \mathbold{B}_2 = \mathbold{B}_2^\prime\mathbold{M}_{\mathbold{B}_1}\mathbold{B}_2, $$ where we use the fact that $ \mathbold{Q}^\prime \mathbold{W}^2 \mathbold{Q}=\mathbold{M}_{\mathbold{B}_1}$.

The adjacency matrix of $\mathcal{G}^\prime$ is the $n_2\times n_2$ matrix ${\mathbold{A}}^\prime$ with entries

Here, $[j] \cap [j^\prime]$ is the set of all vertices in $V_1$ that are connected to both $j\in V_2$ and $j^\prime \in V_2$ in the original bipartite graph $\mathcal{G}$. In the student-teacher example two teachers are connected by an edge if there is at least one student who was taught by both teachers. The weight $({\mathbold{A}}^\prime)_{jk}$ of the edge is larger the more students there are connecting teachers $j$ and $k$, and the more courses they have taken from these teachers. $\mathcal{G}^\prime$ determines the accuracy with which teacher value-added can be estimated.

The matrix $\mathbold{A}^\prime$ is also the adjacency matrix of the simple graph obtained from $\mathcal{G}^\prime$ by replacing all edges $e\in E_{(j,j^\prime)}^\prime$ by one weighted edge, with weight $(\mathbold{A}^\prime)_{jj^\prime}$. Figure (ref) provides an illustration of a simple bipartite graph for students (circular vertices) and teachers (square vertices), given in the left plot, and its induced weighted graph featuring only teachers, given in the right plot. The thickness of the edge between $(j,j^\prime)$ in the latter plot reflects the magnitude of the weight $({\mathbold{A}}^\prime)_{jj^\prime}$.

The device of a one-mode projection highlights the importance of having movers in panel data. In matched worker-firm data sets workers do not frequently switch employer over the course of the sampling period. This lack of mobility is one cause of the substantial bias that is observed in the correlation coefficient between (estimated) worker and firm effects (AbowdKramarzLengermannPerezDuarte2004, AndrewsGillSchankUpward2008,AndrewsGillSchankUpward2012). While this is now well recognized, limited mobility has consequences more broadly. Indeed, it implies that few workers connect firms in the one-mode firm projection. Therefore, the induced graph may be only weakly connected (and $\lambda_2$ will be close to zero) and the variance of the estimator of the firm effects may be large. This is not only detrimental for identifying sorting between workers and firms but, indeed, complicates estimation and inference of the firm effects as well as all their moments, such as their variance. Restricting attention to large firms need not resolve this problem. An analogous argument holds for teacher effects and their estimated variance, and so for our ability to infer the contribution of teacher value-added to observed variation in test scores. We illustrate this in our data below.

Variance bound and asymptotic analysis

Finite-sample bound

To work towards general distribution theory it is instructive to start with a finite-sample bound on the variance of the fixed-effect estimator when the errors are homoskedastic and uncorrelated. Let

This is a (weighted) harmonic mean of the (weighted) degrees $d_j / (\mathbold{A})_{ij}$ of all $j\in [i]$. Note that, for a given vertex $i$, $h_i$ is increasing in the degree of its direct neighbors.

Theorem (ref) states that, for a given degree $d_i$ and global connectivity measure $\lambda_2$, the upper bound on the variance of $\hat{\alpha}_i$ is smaller if the direct neighbors of vertex $i$ are themselves more strongly connected to other vertices in the network. The theorem provides insight into how the local connectivity structure of the network, around vertex $i$, affects statistical precision.

Additional calculations for analytically-tractable cases where $\lambda_2 \rightarrow 0$ as the network grows are provided in the supplementary material.

Theorem (ref) highlights the importance of $\lambda_2 h_i \rightarrow \infty$ as a sufficient condition for the parametric rate ${d_i}^{-1/2}$ to be attainable for estimation of $ \hat{\alpha}_i$. This result carries over to the model with covariates. Let $$ \rho : = \left\lVert (\mathbold{X}^\prime \mathbold{X})^{-1/2} (\mathbold{X}^\prime \mathbold{M}_{\mathbold{B}} \mathbold{X} ) (\mathbold{X}^\prime \mathbold{X})^{-1/2} \right\rVert_2, $$ where $\lVert \cdot \rVert_2$ denotes the spectral norm, and $\mathbold{M}_{\mathbold{B}}: = \mathbold{I}_m - \mathbold{B}(\mathbold{B}^\prime \mathbold{B})^\star \mathbold{B}^\prime$. Note that $\rho \in [0,1]$ is a measure of non-collinearity between the columns of $\mathbold{X}$ and $\mathbold{B}$, with $\rho$ close to zero indicating near-collinearity. Indeed, while $\mathbold{X}^\prime \mathbold{X}$ measures the total variation in $\mathbold{X}$, $\mathbold{X}^\prime \mathbold{M}_{\mathbold{B}} \mathbold{X}$ captures the residual variation in $\mathbold{X}$ after its linear dependence on $\mathbold{B}$ has been partialled out. For $i \in V$ let $\mathbold{b}_i$ be $i^{\text{th}}$ column of $\mathbold{B}$, and let $ \overline{\mathbold{x}}_i : = {\mathbold X}' \mathbold{b}_i / d_i $ and $\mathbold{\Omega}:= \mathbold{X}^\prime \mathbold{X}/m$ in the following theorem.

This result shows that, if $\rho$ is bounded away from zero, introducing covariates only has a higher-order effect on the statistical precision of the fixed-effect estimator. In particular we have

provided that $\lambda_2 h_i \rightarrow \infty$ as $\mathcal{G}$ grows. Furthermore, the parametric rate is achievable even if $\lambda_2$ is not treated as fixed, and $\mathcal{G}$ becomes less dense as more vertices are added to the network.

Large-sample analysis

We now discuss asymptotic results under more general conditions on the regression errors. The following theorem provides a first-order representation of $ \check{\alpha}_i$. Let $\mathbold{\Sigma}:=\mathbbmsl{E}(\mathbold{u}\mathbold{u}^\prime)$.

From the definition of the incidence matrix in (ref), the $m$-vector $\mathbold{b}_i$ has as many non-zero entries as there are edges involving vertex $i$. Further, $\mathbold{b}_i^\prime \mathbold{b}_i = d_i$. Hence, the term ${\mathbold{b}^\prime_i \mathbold{u}} /{d_i}$ is a (weighted) sample mean of the regression errors associated with the edges that involve vertex $i$.

We next consider sequences of growing networks such that

These are relatively weak conditions that ensure that the fact that $\mathbold{\beta}$ is estimated can be ignored in large samples. Moreover, they imply that $$ \epsilon_i = O_p\left(\frac{1}{\sqrt{d_i (\lambda_2\, h_i)}}\right) $$ and that $\tilde \epsilon_i = O_p(1/\sqrt{m}) = o_p(1/\sqrt{d_i})$. The main implication of the theorem is that, then, under the now familiar condition $\lambda_2 h_i\rightarrow\infty$, as $d_i\rightarrow\infty$,

This result allows the errors to be heteroskedastic and correlated.

With Theorem (ref) in hand the limit distribution of $\check{\alpha}_i$ can be deduced under conventional conditions. As an example we do so next for independent but heterogeneously distributed (i.n.i.d.) regression errors.

When the errors $\mathbold{u}$ are independent and homoskedastic we have $\mathbold{b}_i^\prime{\mathbold{\Sigma}}\mathbold{b}_i = \sigma^2\, \mathbold{b}_i^\prime\mathbold{b}_i = \sigma^2 d_i$ and the variance in the theorem reduces to $\sigma^2/d_i$, which agrees with (ref).

A plug-in estimator of $\mathbold{b}_i^\prime{\mathbold{\Sigma}}\mathbold{b}_i$ is $\mathbold{b}_i^\prime\check{\mathbold{\Sigma}}\mathbold{b}_i$, where $\check{\mathbold{\Sigma}}:= \mathrm{diag}(\check{\mathbold{u}}\check{\mathbold{u}}^\prime)/{m}$ and $\check{\mathbold{u}}$ are the residuals from the least-squares regression. This involves estimation of $\alpha_j$ for all $j\in[i]$. We have that $$ (\mathbold{b}_i^\prime\check{\mathbold{\Sigma}}\mathbold{b}_i - \mathbold{b}_i^\prime{\mathbold{\Sigma}}\mathbold{b}_i)/d_i \overset{p}{\rightarrow} 0 $$ as $d_i\rightarrow \infty$, provided that, in addition to the conditions of Theorem (ref) holding, we have that $\lambda_2 H_i$ is bounded away from zero, where $$ H_i : = \left( \sum_{j\in[i]} \frac{(h_i/d_i)/d_j}{h_j} \right)^{-1} $$ is a weighted harmonic mean. At the heart of this result lies (a local version of) a global convergence rate on ${\lVert \check{\mathbold{\alpha}}-\mathbold{\alpha} \rVert}$, which is interesting in its own right. More precisely, letting $$ h:= \left(\frac{1}{n} \sum_{i=1}^n \frac{1}{d_i} \right)^{-1} \qquad \text{ and } \qquad H : = \left( \sum_{i=1}^n \frac{(h/n)/d_i}{h_i} \right)^{-1}, $$ it is easy to see that $$ {\lVert \check{\mathbold{\alpha}}-\mathbold{\alpha} \rVert} = O_p({\textstyle{\sqrt{n/h}}}), $$ provided that $\lambda_2 H$ is bounded away from zero.

Estimation of moments

Suppose that the $\alpha_i$ are sampled from some distribution. One might be interested to learn the variance of this distribution---as in, say, Rockoff2004 or CardHeiningKline2013---or some other moment. The typical estimator is the corresponding sample moment of the estimated effects. Sampling noise in the estimated individual effects will introduce bias in the moment estimator, however. To see this, consider estimation of the variance in a simple model without regressors. The sample variance of the estimated effects in this case is $$ {\hat{\mathbold{\alpha}}^\prime \mathbold{M}_{\mathbold{\iota}_n} \hat{\mathbold{\alpha}}}/{(n-1)}, $$ where $ \mathbold{M}_{\mathbold{\iota}_n} := \mathbold{I}_n -{\mathbold{\iota}_n\mathbold{\iota}_n^\prime}/{n}$ is the usual demeaning matrix. When $\mathbold{u}\sim(\boldsymbol{0},\sigma^2 \mathbold{I}_m)$ its bias is $$ \sigma^2 \, \frac{\mathrm{tr}(\mathbold{L}^\star)}{(n-1)} = (n-1)^{-1} \sum_{i=1}^n \mathrm{var}(\hat{\alpha}_i), $$ which clearly shows how imprecise estimation of $\alpha_i$ contributes to the bias in the variance estimator.

It is difficult to derive an exact expression for the bias for more general functionals. Theorem (ref) is instrumental here. Suppose that $\tau:=\mathbbmsl{E}(\varphi(\alpha_i))$ is of interest. Its plug-in estimator is $$ \check{\tau}:= n^{-1}\sum_{i=1}^n \varphi(\check{\alpha}_i). $$ Under the conditions of the theorem we can calculate the leading bias in this estimator as $$ b: = n^{-1} \sum_{i=1}^n {\mathbbmsl{E}\left(\frac{\varphi^{\prime\prime}(\alpha_i)}{2} \, \frac{\mathbold{b}_i^\prime\mathbold{\Sigma}\mathbold{b}_i/d_i}{d_i}\right)} = O(h^{-1}), $$ where $\varphi^{\prime\prime}$ denotes the second derivative, provided $\lambda_2 \, H \rightarrow \infty$. Simple regularity conditions on $\varphi$ for this bias result to hold are that it is differentiable with $\mathbbmsl{E}(\varphi^{\prime\prime}(\alpha_i)^2) < \infty$ and bounded third derivative. So, quite generally, the bias will shrink like $h^{-1}$. Therefore, for the bias to vanish and $\check{\tau}$ to be consistent, we need that the degrees of the individual vertices grow with $n$ for an increasing fraction of the vertices.

If the functional of interest is the variance, an exact bias correction can be performed (see AndrewsGillSchankUpward2008 and KlineSaggioSolvsten2018). For functionals like $\tau$, a plug-in estimator of the leading-order bias $b$ is easily formed and so an adjusted estimator is readily constructed. Its effectiveness as a bias-correction device will again depend on the connectivity structure of the graph. We postpone a detailed analysis to future work. In a recent contribution KlineSaggioSolvsten2018 present limit theory for quadratic forms in $\mathbold{\alpha}$.

Empirical illustrations

Teacher value-added

We construct a graph connecting teachers as the (weighted) one-mode projection from matched student-teacher data from the North Carolina Education Research Center. The projection of interest is the one discussed in Section 2.3. The full data set includes scores for a standardized test in reading in elementary schools in North Carolina and was used by Verdier2016 to estimate the effect of class-size reduction on student performance. The analysis conducted here is useful to assess the precision with which teacher value-added can be estimated. The data concern pupils in Grades 4 and 5 of elementary school over the period 2008--2012. The full teacher graph (with a single weighted edge between neighboring teachers, as in Figure (ref)) has 12,057 vertices and 53,741 edges and is disconnected. The largest connected component involves 41,612 edges between 11,945 teachers and we work with this subgraph. With $\lambda_2=.0039$ the projected teacher graph is weakly connected. Its local connectivity is summarized in Table (ref). The table contains the mean, standard deviation, and deciles of the relevant degree distributions. Inspection reveals that the degrees are small for all teachers.

The weak connectivity suggests that inference on teacher value-added will be difficult. To get a sense of the precision of a first-order asymptotic approach we can look at the ratio $$ \frac{(\mathbold{L}^\star)_{ii}}{1/d_i} = (\mathbold{S}^\dagger)_{ii}. $$ This is the exact variance of $\hat{\alpha}_i$ to its large-sample approximation in a regression model with homoskedastic and uncorrelated errors. This ratio is free of $\sigma^2$ and can be computed directly from the graph. The left plot of Figure (ref) shows the deciles of the distribution of $(\mathbold{S}^\dagger)_{ii}$. The asymptotic approximation is revealed to be widely inaccurate. On average, the actual variance is about 2.5 times larger than its approximation. Even the first decile equals $1.29$. This implies that confidence intervals based on the large-sample arguments in Theorem (ref) are overoptimistic. To illustrate this the right plot in Figure (ref) gives the distribution of the width of 95% confidence intervals for the $\alpha_i$ using both the exact variance (solid line) and its large-sample approximation (dashed line) for the case $\sigma^2=1$. The former stochastically dominates the latter.

The large variability in the estimators of teacher value-added implies a large bias in their estimated variance. We calculate $$ \mathrm{tr}(\mathbold{L}^\star)/(n-1) = .3545, $$ so the bias in the plug-in estimator of the variance is about one-third of the error variance when $\mathbold{u}\sim(\boldsymbol{0},\sigma^2\mathbold{I}_m)$. The large-sample approximation to the bias here is proportional to $h^{-1}$. With $h=5.4554$ this yields a bias of about 18%, roughly half the size of the exact bias.

One reason for the global connectivity of the teacher graph to be low is limited mobility of teachers between schools. If interest lies in comparing teacher effectiveness within a given school it suffices to restrict attention to that subgraph. Of course, the effective sample size from which teacher value-added is estimated will remain small unless additional years of data are collected. Making accurate comparisons between schools is more complicated as it requires many teachers teachers to switch schools during the sampling period. Collecting additional years of data will not automatically lead to more precise estimates. Mansfield2015 discusses the feasibility of ranking teachers within and between schools (see also MihalyMcCaffreySassLockwood2013 for related discussions). A previous draft of this paper contains versions of our main theoretical results specialized to within and between decompositions of graphs.

Occupational network

Wage regressions on worker and occupational dummies (as in KambourovManovskii2009, for example) provide an interesting example of a situation where more accurate results can be obtained. We use all 18 available waves from the BHPS (for a total of 132,097 observations) to construct the induced (weighted) occupational network. The Standard Occupational Classification (SOC90) in the BHPS distinguishes (at the three-digit level) between 374 occupations. We again focus on the largest connected component, which contains 365 occupations with 14,825 weighted edges between them. As a measure of global connectivity here we find $\lambda_2=.3289$. Compared to traditional matched employer-employee data our occupational network does not suffer as much from limited mobility. One reason is that the number of occupations is relatively small compared to the number of workers. Another is that workers may switch occupation also if they remain employed by the same firm, for example due to internal promotions. Finally, as we are dealing with self-reported occupations there is also the possibility of spurious mobility due to misreporting. A look at the distributions summarized in Table (ref) reveals that the degrees and harmonic means tend to be larger here than in the teacher graph.

The distribution of $(\mathbold{S}^\dagger)_{ii}$ now places most of its mass in close vicinity of unity. Its mean and standard deviation are $1.034$ and $.0521$. The median is $1.0202$ while the first and ninth decile are $1.0046$ and $1.0759$, respectively. This suggests that, here, the large-sample approximation to the variance is a much more accurate reflection of actual estimation uncertainty. Similarly, we may again calculate $\mathrm{tr}(\mathbold{L}^\star)/(n-1) = .0577$, which is about 7 times smaller than in the previous example. Further, as $h^{-1} = .0566$, here, the bias approximation is quite accurate.

Conclusion

We have presented inference results on individual effects in a linear fixed-effect regression model when the underlying data structure constitutes a (weighted) graph. An important example is a two-way regression model on a bipartite graph. The main contribution of this paper is to quantify the dependence of statistical precision of the estimator on the connectivity structure of the graph. A key measure of global connectivity is the smallest non-zero eigenvalue of the (normalized) Laplacian matrix of the graph. It reflects the intuitive notion of mobility in the network. A small eigenvalue captures the presence of bottlenecks, which is detrimental to statistical precision. Several measures of local connectivity, such as the degree structure and various harmonic means thereof, also arise naturally in our analysis.

Our theoretical work highlights the importance of and the interplay between global and local measures of network connectivity for conventional inferential approaches to be reliable. The analysis points to a set of simple statistics that can be inspected to evaluate whether the network is sufficiently well connected in a given application. In an application to teacher value-added we find that this is not the case. We further find that conventional standard errors on teacher value-added estimates are much too small, resulting in a false sense of (statistical) precision on these parameter estimates. In an occupational network, on the other hand, we find much higher measures of connectivity and support for our large-sample approximations.

\newtheorem{theoremS}{Theorem} \newtheorem{lemmaS}[theoremS]{Lemma} \newtheorem{corollaryS}[theoremS]{Corollary} \newtheorem{exampleS}{Example}

\setcounter{section}{0} \setcounter{equation}{0} \setcounter{theoremS}{0} \setcounter{exampleS}{0} \setcounter{figure}{0}

\setcounter{section}{0} \setcounter{footnote}{0} \setcounter{page}{1} \pagenumbering{roman}

Additional illustrations

Recall that our measure of global connectivity of the graph $\mathcal{G}$ is $\lambda_2$, the second smallest eigenvalue of the normalized Laplacian matrix. In the following we provide some concrete examples of graphs for which $\lambda_2$ can be explicitly calculated, and we discuss the implications of our variance bound in Theorem (ref)

Our first example illustrates that, even if $\lambda_2\rightarrow 0$ with the sample size, we may still have that $\mathrm{var}(\hat{\alpha}_i) \asymp d_i^{-1}$.

Theorem (ref) allows to establish the convergence rate for the hypercube, but the conditions are too stringent to obtain (ref). The reason is that $h_i$ does not increase fast enough to ensure that $\lambda_2\,h_i\rightarrow \infty$. The following example deals with an extended hypercube and illustrates that, despite $\lambda_2 \rightarrow 0$, we still have $\lambda_2\,h_i\rightarrow \infty$ in this case.

The next example shows that our bound can still be informative if $h_i$ is finite.

The previous example also illustrates that $\lambda_2$ can be large despite having many vertices with small degrees. It is largely due to this property that we prefer to measure global connectivity by $\lambda_2$ and not by the “algebraic connectivity” (the second smallest eigenvalue of ${\mathbold L}$; see, e.g., Chung1997), which has been studied more extensively.

Our last example shows the effect on the upper bound in Theorem (ref) when neighbors themselves are more strongly connected.

The last two examples also illustrate that adding edges to the graph (in this case, to obtain the Wheel graph from the Star graph) can result in a decrease of our measure of global connectivity $\lambda_2$. This is not a problem, however, for our results as we only require that $\lambda_2$ be sufficiently different from zero. The Wheel graph with $\lambda_2 \geq 1/3$, for example, clearly describes a very well globally connected graph by that measure.

Variance bounds for differences

\setcounter{equation}{0}

Our focus in the main text has been inference on the $\alpha_i$, under the constraint in (ref), $\sum_i d_i \alpha_i = 0$. An alternative to normalizing the parameters that may be useful in certain applications is to focus directly on the differences $\alpha_i-\alpha_j$ for all $i\neq j$. An example where this is the case is FinkelsteinGentzkowWilliams2016. We give a corresponding version of Theorem (ref) here.

Let $d_{ij}:= \sum_{k\in V} (\mathbold{A})_{ik} \, (\mathbold{A})_{jk}$. for an unweighted graph $d_{ij} = \left| [i] \cap [j] \right|$, the number of vertices that are neighbors of both $i$ and $j$. Write

for the corresponding harmonic mean of the degrees of the vertices $ k\in [i] \cap [j]$. We have the following theorem.

For a simple graph $\mathcal{G}$, when $[i]=[j]$ but $i\notin [j]$ and $i\notin [j]$, that is, when vertices $i$ and $j$ share exactly the same neighbors and are not connected themselves, the theorem implies

as, in that case, both $(\mathbold{A})_{ij}$ and the second term in the upper bound in Theorem (ref) are zero.

Alternative normalization

If we change the normalization constraint in the least-squares minimization problem (ref) to $$ \sum_{i=1}^n \alpha_i = 0, $$ we obtain the estimator $\hat{\mathbold{\alpha}}^{\, \diamond} = \mathbold{M}_{\mathbold{\iota}} \hat{\mathbold{\alpha}}$, where $ \mathbold{M}_{\mathbold{\iota}} = \mathbold{I}_n - n^{-1} \mathbold{\iota}_n \mathbold{\iota}_n'$ is the projector orthogonal to $\mathbold{\iota}_n$. We then have $ {\rm var}\left( \hat{\mathbold{\alpha}}^{\, \diamond} \right) = \sigma^2 \mathbold{L}^\dagger$, because this variance needs to satisfy $ {\rm var}\left( \hat{\mathbold{\alpha}}^{\, \diamond} \right) \mathbold{\iota}_n = 0$, and the Moore-Penrose pseudoinverse guarantees that the nullspace of $\mathbold{L}$ equals the nullspace of $\mathbold{L}^\dagger$. Thus, changing the normalization corresponds to changing the particular pseudoinverse of $\mathbold{L}$ that features in the expression for the variance. From $\hat{\mathbold{\alpha}}^{\, \diamond} = \mathbold{M}_{\mathbold{\iota}} \hat{\mathbold{\alpha}}$ we find

which thus also shows that $ \mathbold{L}^\dagger = \mathbold{M}_{\mathbold{\iota}} \mathbold{L}^\star \mathbold{M}_{\mathbold{\iota}}$. We have $\mathbold{L}^\star \leq \lambda_2^{-1} \mathbold{D}^{-1}$, and therefore $ \mathbold{L}^\dagger \leq \lambda_2^{-1} \mathbold{M}_{\mathbold{\iota}} \mathbold{D}^{-1} \mathbold{M}_{\mathbold{\iota}} $. We thus find ${\rm var}( \hat{\alpha}^{\, \diamond}_i ) = \sigma^2 \mathbold{e}_i' \mathbold{L}^\dagger \mathbold{e}_i \leq \lambda_2^{-1} \sigma^2 \mathbold{e}_i' \mathbold{M}_{\mathbold{\iota}} \mathbold{D}^{-1} \mathbold{M}_{\mathbold{\iota}} \mathbold{e}_i$, and evaluating the last expression gives the following theorem.

Notice that $d_i / (n \, h) \leq 1/h \leq 1$, and therefore ${\rm var}( \hat{\alpha}^{\, \diamond}_i ) \leq \frac{2}{d_i} \frac{\sigma^2}{\lambda_2} $. For the estimator $\hat \alpha_i$ obtained under the normalization in the main text we immediately find from (ref) and $\left( \mathbold{S}^\dagger \right)_{ii} \leq \lambda_2^{-1}$ that ${\rm var}( \hat{\alpha}_i ) \leq \frac{1}{d_i} \frac{\sigma^2}{\lambda_2} $. Thus, for sequences of growing networks we find the pointwise consistency results $(\hat{\alpha}^{\, \diamond}_i-\alpha_i) \overset{p}{\rightarrow} 0$ and $(\hat{\alpha}_i-\alpha_i) \overset{p}{\rightarrow} 0$ for both estimators, under the sufficient condition $\lambda_2 \, d_i \rightarrow \infty$.

Analogously one can extend Theorem (ref) from $\hat \alpha_i $ to $\hat{\alpha}^{\, \diamond}_i $ as follows.

Analogous to (ref) in the main text we thus find

provided that $\lambda_2 h_i \rightarrow \infty$ and $n h / d_i \rightarrow \infty $ and $n h_i^{(2)} / d_i \rightarrow \infty $ and $\lambda_2 \, h \, H / d_i \rightarrow \infty $ as $n\rightarrow\infty$. Therefore, under plausible assumptions on the sequence of growing networks we find the same asymptotic properties for $ \hat \alpha_i^{\, \diamond} $ as for $\hat \alpha_i $. The particular choice of normalization in the main text is not necessary for our main results, but it makes all derivations as well as the presentation of the results more convenient.

Proofs

\paragraph{PROOF OF THEOREM (ref) (EXISTENCE)} \addcontentsline{toc}{subsection}{ Proof of Theorem (ref)} \newline The estimator is defined by the constraint minimization problem in (ref). For convenience we express the constraint in quadratic form, $ \left( \mathbold{a}^\prime \mathbold{d} \right)^2 = 0$. By introducing the Lagrange multiplier $\lambda > 0$ we can write

Solving the corresponding first-order condition we obtain

where $\mathbold{S}_{\mathbold{X}} := \mathbold{D}^{-1/2} \mathbold{B}' \mathbold{M}_{\mathbold{X}} \mathbold{B} \mathbold{D}^{-1/2} $ and $\mathbold{\psi} := \mathbold{D}^{1/2} \mathbold{\iota}_n = \mathbold{D}^{-1/2} \mathbold{d}$. Since we assume that the graph is connected we have $d_i>0$ for all $i$, that is, $ \mathbold{D}$ is invertible. Our assumption $\mathrm{rank}((\mathbold{X},\mathbold{B}))=p+n-1$ implies that $\mathrm{rank}(\mathbold{B}' \mathbold{M}_{\mathbold{X}} \mathbold{B})=n-1$, that is, the zero eigenvalue of $\mathbold{B}' \mathbold{M}_{\mathbold{X}} \mathbold{B}$ has multiplicity one. By construction of $\mathbold{B}$ we have $\mathbold{B} \mathbold{\iota}_n=0$, that is, the zero eigenvector of $\mathbold{B}' \mathbold{M}_{\mathbold{X}} \mathbold{B}$ is given by $ \mathbold{\iota}_n$. It follows that the zero eigenvalue $\mathbold{S}_{\mathbold{X}}$ has multiplicity one and eigenvector $ \mathbold{\psi} $. This explains why the matrix $\mathbold{S}_{\mathbold{X}} + \lambda \, \mathbold{\psi} \, \mathbold{\psi}' $ is invertible, which we already used in (ref). Furthermore, the matrices $ \mathbold{S}_{\mathbold{X}} $ and $ \mathbold{\psi} \, \mathbold{\psi}' $ commute, and by properties of the Moore-Penrose inverse we thus have

We furthermore have

where $m = \mathbold{\psi}' \mathbold{\psi}$ is the total number of observations. Because $ \mathbold{B} \mathbold{\iota}_n = 0$, the contribution from $ \left( \mathbold{\psi} \, \mathbold{\psi}' \right)^\dagger$ drops out of (ref), and we obtain

according to the definition of the pseudoinverse $^\star $ in the main text. Notice that $\check{\mathbold{\alpha}}$ given in the last display does not depend on $\lambda$, and automatically satisfies the constraint $\mathbold{d}'\check{\mathbold{\alpha}}=0$, that is, any value of $\lambda$ can be chosen in the above derivation. \qed

\paragraph{PROOF OF THEOREMS (ref) AND (ref) (VARIANCE BOUNDS)} \addcontentsline{toc}{subsection}{ Proof of Theorems (ref) and (ref)} \newline We first show that, if $\mathcal{G}$ is connected, then

Theorems (ref) and (ref) will then follow readily. Analogous to (ref) we also hav $\left(\mathbold{S} + \lambda \, \mathbold{\psi} \, \mathbold{\psi}' \right)^{-1} = \mathbold{S}^\dagger + \lambda^{-1 } \left(\mathbold{\psi} \, \mathbold{\psi}' \right)^\dagger$. Using this and (ref) we find

and since $\mathbold{S} \mathbold{\psi} = 0$ and $ \mathbold{\psi}' \mathbold{\psi} = m$ we thus find that $ \mathbold{S}^\dagger \mathbold{S} = \mathbold{I}_n - m^{-1} \mathbold{\psi} \mathbold{\psi} ' , $ which is simply the idempotent matrix that projects orthogonally to $\mathbold{\psi}$. We thus find $ \mathbold{L}^\star \mathbold{L} = \mathbold{D}^{-1/2 } \mathbold{S}^\dagger \mathbold{S} \mathbold{D}^{1/2 } = \mathbold{I}_n - m^{-1} \mathbold{\iota}_n \mathbold{d}' . $ Plugging in $\mathbold{L} = \mathbold{D} - \mathbold{A}$, and then solving for $\mathbold{L}^\star$ gives

The Laplacian is symmetric, and so transposition gives

Replacing $\mathbold{L}^\star$ on the right-hand side of (ref) by the expression for $ \mathbold{L}^\star$ given by (ref), and also using that $ \mathbold{D}^{-1} \mathbold{A} \mathbold{\iota}_n = \mathbold{\iota}_n$, yields

Re-arranging this equation allows us to write

From $ \mathbold{L}^* = \mathbold{D}^{-1/2} \mathbold{S}^\dagger \mathbold{D}^{-1/2} $ and $\mathbf{0} \leq \mathbold{S}^\dagger \leq \lambda_2^{-1} \mathbold{I}_n$ we obtain $\mathbf{0} \leq \mathbold{L}^\star \leq \lambda_2^{-1} \mathbold{D}^{-1}$, and therefore $$ \mathbf{0} \leq \mathbold{D}^{-1} \mathbold{A} \mathbold{L}^\star \mathbold{A} \mathbold{D}^{-1} \leq \lambda_2^{-1} \mathbold{D}^{-1} \mathbold{A} \mathbold{D}^{-1} \mathbold{A} \mathbold{D}^{-1}. $$ Put together this yields

and multiplication with $\sigma^2 $ gives the bounds stated in (ref).

To show Theorems (ref) and (ref) we calculate, for $i\neq j$,

where $\mathbold{e}_i$ is the vector that has one as its $i^{\mathrm{th}}$ entry and zeros elsewhere. Combining these results with (ref) gives the bounds on, respectively, ${\rm var}( \hat \alpha_i ) = {\mathbold e}_i' {\rm var}( \hat {\mathbold \alpha} ) {\mathbold e}_i$ and ${\rm var}( \hat \alpha_i - \hat \alpha_j ) = ( {\mathbold e}_i - {\mathbold e}_j )' {\rm var}( \hat {\mathbold \alpha} ) ( {\mathbold e}_i - {\mathbold e}_j )$ stated in the theorems. \qed

\paragraph{PROOF OF THEOREMS (ref) AND (ref)} \addcontentsline{toc}{subsection}{ Proof of Theorems (ref) and (ref)} \newline Using that $ \mathbold{L}^* \leq \lambda_2^{-1}\mathbold{D}^{-1}$ we find that

and we calculate

Combing those results gives the statement of Theorem (ref)

Next, multiplying $\mathbold{M}_{\mathbold{\iota}}$ from the left and right to the matrix bounds (ref) and using ${\rm var}\left( \hat{\mathbold{\alpha}}^{\, \diamond} \right) = \mathbold{M}_{\mathbold{\iota}} {\rm var}\left( \hat{\mathbold{\alpha}} \right) \mathbold{M}_{\mathbold{\iota}}$ gives

and therefore

We already calculated $ \mathbold{e}_i' \mathbold{M}_{\mathbold{\iota}} \mathbold{D}^{-1} \mathbold{M}_{\mathbold{\iota}} \mathbold{e}_i$ in (ref) above. We furthermore have

and by applying the Cauchy-Schwarz inequality we find $\sum_{j,k} \frac {( \mathbold{A} )_{jk} } {d_j d_k} \leq \sum_{j,k} \frac {( \mathbold{A} )_{jk} } {d_j^2} = \sum_j \frac 1 {d_j}$, and therefore

Similarly, $ \mathbold{e}_i' \mathbold{M}_{\mathbold{\iota}} \mathbold{D}^{-1} \mathbold{A} \mathbold{D}^{-1} \mathbold{A} \mathbold{D}^{-1} \mathbold{M}_{\mathbold{\iota}} \mathbold{e}_i \geq 0$ contains three terms, for which we have

where in the last line we again applied the Cauchy-Schwarz inequality, and the definitions of the harmonic means $h$ and $H$ in the main text. Combining the above gives the statement of Theorem (ref).

\paragraph{PROOF OF THEOREM (ref) (COVARIATES)} \addcontentsline{toc}{subsection}{ Proof of Theorem (ref)} \newline Define the $n \times n$ matrix

Let $\lambda_i(\mathbold{C})$ denote the $i$th eigenvalue of $\mathbold{C}$, arranged in ascending order. $\mathbold{C}$ is similar to the positive semi-definite matrix $$ \left( {\mathbold X}' {\mathbold X} \right)^{-1/2} {\mathbold X}' \mathbold{B} \left(\mathbold{B}' \mathbold{B} \right)^\star \mathbold{B}' {\mathbold X} \left( {\mathbold X}' {\mathbold X} \right)^{-1/2}, $$ and since similar matrices share the same eigenvalues we have $\lambda_1(\mathbold{C}) \geq 0$. $\mathbold{C}$ is also similar to the matrix $$ \mathbold{B} \left(\mathbold{B}' \mathbold{B} \right)^\star \mathbold{B}' {\mathbold X} \left( {\mathbold X}' {\mathbold X} \right)^{-1} {\mathbold X}', $$ which is the product of two projection matrices, whose spectral norm is thus bounded by one. Hence, $\lambda_n(\mathbold{C})\leq 1$. In addition, we must have $\lambda_i( \mathbold{C} ) \neq 1$ for any $1<i<n$ because, otherwise, ${\rm rank}\left( {\mathbold I}_n - \mathbold{C} \right) < n$, which implies that ${\rm rank}( \mathbold{B}' {\mathbold M}_{\mathbold X} \mathbold{B}) < n-1$, contradicting our non-collinearity assumption (since the graph is connected we have ${\rm rank}( \mathbold{B}'\mathbold{B} ) = n-1 $, which together with the non-collinearity assumption $\mathrm{rank}((\mathbold{X},\mathbold{B}))=p+n-1$ implies that ${\rm rank}( \mathbold{B}' {\mathbold M}_{\mathbold X} \mathbold{B}) = n-1$). We therefore have $\lVert \mathbold{C}\rVert_2 < 1$, implying that $ {\mathbold I}_m - \mathbold{C}$ is invertible.

Using (ref) and (ref) with $\lambda=m^{-1}$ we find that $ \left( \mathbold{B}' {\mathbold M}_{\mathbold X} \mathbold{B} + m^{-1} \mathbold{D} {\mathbold \iota}_n {\mathbold \iota}_n' \mathbold{D} \right)^{-1} = (\mathbold{B}' {\mathbold M}_{\mathbold X} \mathbold{B})^\star + m^{-1} \, {\mathbold \iota}_n {\mathbold \iota}_n' , $ or equivalently

and analogously we have

Subtracting the expressions in the last two displays gives

and by multiplying with $ \left[ (\mathbold{B}' \mathbold{B})^\star + m^{-1} \, {\mathbold \iota}_n {\mathbold \iota}_n' \right]$ from the left and $ \left[ (\mathbold{B}' {\mathbold M}_{\mathbold X} \mathbold{B})^\star + m^{-1} \, {\mathbold \iota}_n {\mathbold \iota}_n' \right]$ from the right, and using $ \mathbold{B} {\mathbold \iota}_n = 0 $, we obtain

which can equivalently be expressed as $ \left( {\mathbold I}_m - \mathbold{C} \right) (\mathbold{B}' {\mathbold M}_{\mathbold X} \mathbold{B})^\star = \left(\mathbold{B}' \mathbold{B} \right)^\star$. We have already argued that $\left( {\mathbold I}_m - \mathbold{C} \right) $ is invertible, and therefore

Since $\left\| \mathbold{C} \right\|_2 < 1$ we can expand $ \left( {\mathbold I}_m - \mathbold{C} \right)^{-1}$ in powers of $ \mathbold{C} $, as

Defining the $p \times p$ matrix

we can rewrite (ref) as

The parameter $\rho$ defined in the main text satisfies

that is, we have $\lVert \widetilde {\mathbold C} \rVert_2 = 1 - \rho$, and since $ \widetilde {\mathbold C}$ is symmetric and semi-definite this can equivalently be written as $ \widetilde {\mathbold C} \leq (1 - \rho) {\mathbold I}_p$. Therefore,

We thus have

and, therefore,

Using the expression (ref) and (ref) for $ \left(\mathbold{B}' \mathbold{B} \right)^\star = \mathbold{L}^\star $ we obtain

where

and we used the Cauchy-Schwarz inequality to bound the mixed term. Again, because similar matrices have the same eigenvalues we have

and therefore,

where in the last step we used $ {\mathbold e}_i' \mathbold{D}^{-1} \mathbold{A} \mathbold{D}^{-1} \mathbold{A} \mathbold{D}^{-1} {\mathbold e}_i = (d_i h_i)^{-1}$. Using our definitions $ \overline {\mathbold x}_i = {\mathbold X}' \mathbold{B} \mathbold{D}^{-1}{\mathbold e}_i $ and $\mathbold{\Omega} = \mathbold{X}^\prime \mathbold{X}/m$ we obtain

Combining the above results we find

For any $a,b \geq 0$ we have $a + b + 2 \sqrt{a b} \leq 2 (a+b)$. Thus, a slightly cruder but simpler bound is given by

where we also used that $ {\rm var}\left( \check\alpha_i \right) \geq {\rm var}\left( \hat \alpha_i \right) $, because adding regressors can only increase the variance of the least squares estimator under homoskedasticity. \qed

\paragraph{PROOF OF THEOREM (ref) (FIRST ORDER REPRESENTATION)} \addcontentsline{toc}{subsection}{ Proof of Theorem (ref)} \newline Remember that we treat $\mathbold{B}$ and $\mathbold{X}$ as fixed (i.e.\ non-random) throughout. Let $\check{\mathbold{\beta}} := \left( \mathbold{X}' \mathbold{M}_{\mathbold{B}} \mathbold{X} \right)^{-1} \mathbold{X}' \mathbold{M}_{\mathbold{B}} \mathbold{y} $. Using the model for $\mathbold{y} $ we find $\check{\mathbold{\beta}} - {\mathbold{\beta}} = \left( \mathbold{X}' \mathbold{M}_{\mathbold{B}} \mathbold{X} \right)^{-1} \mathbold{X}' \mathbold{M}_{\mathbold{B}} \mathbold{u} $. Using our assumptions $ \mathbbmsl{E} (\mathbold{u}) = 0$ and $ \mathbold{\Sigma} \leq \mathbold{I}_m \overline{\sigma}^2$ we find $ \mathbbmsl{E}( \check{\mathbold{\beta}} - {\mathbold{\beta}} ) = 0$ and

The result in (ref) can be rewritten as

The constrained least-squares estimator in (ref) can be expressed as

and analogous to Theorem (ref) we then find $ \check{\mathbold{\alpha}} = \mathbold{L}^\star \mathbold{B}' \left( \mathbold{y} - \mathbold{X} \check{\mathbold{\beta}} \right) = \left( \mathbold{L} + m^{-1} \mathbold{d} \mathbold{d}^\prime \right)^{-1} \mathbold{B}' \allowbreak \left( \mathbold{y} - \mathbold{X} \check{\mathbold{\beta}} \right) $. Multiplying by $\left( \mathbold{L} + m^{-1} \mathbold{d} \mathbold{d}^\prime \right)$ from the left and using our normalization $ \mathbold{d}^\prime \check{\mathbold{\alpha}} = 0$ gives $$ \mathbold{L} \, \check{\mathbold{\alpha}} = \mathbold{B}' \left( \mathbold{y} - \mathbold{X} \check {\mathbold{\beta}} \right). $$ Plugging $\mathbold{L}=\mathbold{D}-\mathbold{A}$ and $\mathbold{y} = \mathbold{B} \mathbold{\alpha} + \mathbold{X} \mathbold{\beta} + \mathbold{u}$ into the last display, multiplying from the left with $\mathbold{D}^{-1}$, and rearranging terms, we obtain

where

We have $ \mathbbmsl{E} (\check{\mathbold{\beta}} - {\mathbold{\beta}} ) = 0$ and $ \mathbbmsl{E} ( \check{\mathbold{\alpha}} - \mathbold{\alpha} ) = 0$, and, therefore, also $\mathbbmsl{E} (\mathbold{\epsilon}) = \mathbf{0}$ and $\mathbbmsl{E}( \tilde{\mathbold{\epsilon}} ) = \mathbf{0}$. The definition $ \rho = \left\lVert (\mathbold{X}^\prime \mathbold{X})^{-1} \mathbold{X}^\prime \mathbold{M}_{\mathbold{B}} \mathbold{X} \right\rVert_2$ can equivalently be written as $ \rho \mathbold{X}^\prime \mathbold{X} \geq \mathbold{X}^\prime \mathbold{M}_{\mathbold{B}} \mathbold{X} $, and therefore $ \rho^{-1} \left( \mathbold{X}^\prime \mathbold{X} \right)^{-1} \leq \left( \mathbold{X}^\prime \mathbold{M}_{\mathbold{B}} \mathbold{X} \right)^{-1} $. Using this and (ref) we obtain

Using $ \check{\mathbold{\alpha}} - \mathbold{\alpha} = (\mathbold{B}^\prime\mathbold{M}_{\mathbold{X}}\mathbold{B})^{\star} \mathbold{B}^\prime \mathbold{M}_{\mathbold{X}}\mathbold{u} $ and the assumption $ \mathbold{\Sigma} \leq \overline{\sigma}^2 \mathbold{I}_n$ we calculate

where in the last step we used (ref). Since furthermore ${\mathbold X} \left( {\mathbold X}' {\mathbold X} \right)^{-1} {\mathbold X}' \leq {\mathbold I}_m$ and $ \left(\mathbold{B}' \mathbold{B} \right)^\star = \mathbold{L}^\star \leq \lambda_2^{-1} \mathbold{D}^{-1}$ we obtain

Denote the elements of $\mathbold{\epsilon}$ and $\tilde{\mathbold{\epsilon}}$ by $ \epsilon_i $ and $\tilde \epsilon_i $. Equation (ref) can then be written as

and we have

and

where we used our definitions $\overline{\mathbold{x}}_i = {\mathbold X}' \mathbold{b}_i / d_i = \mathbold{X}' \mathbold{B} \mathbold{D}^{-1} \mathbold{e}_i$ and $\mathbold{\Omega}:= \mathbold{X}^\prime \mathbold{X}/m$. \qed

\paragraph{PROOF OF THEOREM (ref) (ASYMPTOTIC DISTRIBUTION)} \addcontentsline{toc}{subsection}{ Proof of Theorem (ref)} \newline We have $\rho \leq 1$ by definition. Together with the assumptions $\overline{\sigma}^2 = O(1)$, $\lambda_2 h_i\rightarrow\infty$, and the conditions in (ref) this implies that $\mathbbmsl{E}(\epsilon_i^2)\leq \overline{\sigma}^2 (1+\rho) /(\rho \, d_i \, \lambda_2 \, h_i) = o(d_i^{-1})$, and $\mathbbmsl{E}(\tilde \epsilon_i^{\,2})\leq \overline{\sigma}^2 \, \overline{\mathbold{x}}_i' \, \mathbold{\Omega}^{-1} \, \overline{\mathbold{x}}_i / (\rho \, m) = o(d_i^{-1})$. By Markov's inequality we thus have $\epsilon_i = o_p(d_i^{-1/2})$ and $\tilde \epsilon_i = o_p(d_i^{-1/2})$, and applying Theorem (ref) gives, as $d_i\rightarrow\infty$,

The number of terms $\nu_{\varepsilon_{e} i} $ summed over in the last display grows to infinity asymptotically, because we assume that $d_i = \sum_{j \in [i]} \sum_{e \in E_{(i,j)}} w_e \rightarrow \infty$, while the weights $w_e = \left( \mathbold{B} \right)_{\varepsilon_{e} i}^2$ are bounded. Our assumptions furthermore guarantee that the $\nu_{\varepsilon_{e} i}$ are independent and satisfy $\mathbbmsl{E} (\nu_{\varepsilon_{e} i}) = 0$, $\mathbbmsl{E} (\nu_{\varepsilon_{e} i}^2 )\geq c_1 > 0$, and $\mathbbmsl{E} (\left| \nu_{\varepsilon_{e} i} \right|^3) \leq c_2 < \infty$ for constants $c_1,c_2$. Thus, the Lyapunov condition is satisfied, and the statement of the theorem then follows from a standard application of Lyapunov's central limit theorem. \qed