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.
75,497 characters · 24 sections · 64 citation commands
Weak Identification in Peer Effects Estimation
The study of peer effects and social interactions has grown immensely in recent years, driven by the increasing availability of social network data. Researchers seek to understand the extent to which interactions between individuals drive correlations in outcomes. For example, smoking behavior has been shown to spread through social networks, with individuals more likely to smoke if their friends do saari2014smoking. One possible explanation for this is direct transmission of behavior: individuals may mimic their peers' behaviors due to factors such as social pressure. This is known as an endogenous peer effect. However, there are other possible sources of correlation: friends may tend to belong in the same socioeconomic class, and individuals in the same class could behave similarly. This is known as a contextual peer effect.
The seminal work of manski1993identification formalized the challenge of separately identifying endogenous and contextual peer effects. In his framework, endogenous effects are defined as the influence of average peer outcomes, while contextual effects reflect the influence of average peer baseline characteristics. Manski showed that in this setting, the two effects are not separately identifiable, making it statistically impossible to distinguish the source of observed outcome correlations. He termed this the “reflection problem.”
While Manski's work focused on dense group interactions where individuals are presumed to be affected by everyone in their group, it has been recognized that social network data can provide more detailed accounts of which interactions are actually taking place. Borrowing from ideas in the spatial econometric literature, bramoulle2009identification proposed the following network linear-in-means model, which defines endogenous and contextual peer effects as average outcomes and baseline covariates over an individual's social connections, respectively:
where $Y\in \mathbb{R}^{n}$ are outcomes, $X\in \mathbb{R}^{n\times p}$ are baseline $p$-dimensional covariates, $G\in \mathbb{R}^{n\times n}$ is a row-normalized adjacency matrix of social interactions, and $\varepsilon\in \mathbb{R}^n$ is an error term. The parameters of interest are $\theta:=(\alpha, \beta, \delta, \rho)$, where $\delta$ and $\rho$ represent contextual and endogenous peer effects, respectively. Since the term $GY$ is endogenous, identification usually requires a valid instrument. The key insight of bramoulle2009identification is that higher-order lags of $X$, such as $G^2X$, serve as valid instruments; consequently, if the collection $[1, X, GX,G^2X]$ are not collinear, then the peer effect parameters are identified, making Manski's non-identification result a special case where the network is complete.
Most existing applications of this model assume strong identification, where the relevant covariance matrices remain well-conditioned in the limit. This assumption is often justified by an asymptotic regime in which data are drawn from many small, independent networks, such as school classrooms sacerdote2011peer or villages banerjee2013diffusion, cai2015social. Such an assumption mirrors the “increasing-domain” asymptotics used in spatial statistics anselin2022spatial, and standard asymptotic theory for independent and identically distributed (i.i.d.) data holds by considering each network to be an independent sample.
However, many modern social settings do not fit this framework. For example, online social media platforms typically involve a single network of interconnected users who may have a large number of connections aral2012information. In such environments, identification strength may deteriorate as network density increases. This is particularly problematic for the linear-in-means model, where key covariates, such as neighborhood averages, become increasingly similar across individuals in dense networks. In the weak instruments literature, this loss of variation leads to slower convergence rates antoine2009efficient,hahn2002discontinuities or non-standard limiting distributions staiger1994instrumental.
While weak identification has been observed in simulations of peer effects bramoulle2009identification, startz2017improved, the phenomenon remains theoretically underexplored. In this paper, we seek to answer the following question: under what natural conditions does the linear-in-means model suffer from weak identification, and how does this affect standard estimators of peer effects?
Motivated by intuition on neighborhood averages, we analyze peer effects estimation under two simple assumptions: that network density increases at any rate, and that covariates are i.i.d. We show that these two assumptions are sufficient for weak identification to occur. Specifically, we prove that ordinary least squares is generally biased, in contrast to prior results showing consistency under strong identification lee2002consistency. On the other hand, we establish that two-stage least squares estimation is consistent for networks with density $o(\sqrt{n})$, although the estimation rate degrades as network density increases. This non-standard asymptotic behavior occurs as a result of vanishing signal in network averages and asymptotic collinearity of the design matrix.
These results naturally motivate the search for alternative models that are less vulnerable to weak identification. We consider the related linear-in-sums model, which replaces the row-normalized adjacency matrix $G$ with the unnormalized matrix $A$. We demonstrate that this model generally avoids the identification issues seen in the linear-in-means model. In particular, for random graphs, we show that as long as the underlying network model has sufficient structural variation, the model is strongly identified, and the two-stage least squares estimator converges at a standard rate. Linear-in-sums models have been applied more recently in the literature and can be interpreted as local aggregation instead of averaging ballester2006s,liu2014endogenous. We argue that this model can be a useful alternative in discerning endogenous and contextual effects.
This paper is organized as follows: in Section (ref), we survey some related literature on the linear-in-means model. Section (ref) formally introduces the network linear-in-means model along with our asymptotic framework and baseline assumptions. In Section (ref), we state and prove our main results on weak identification in the linear-in-means model. In Section (ref), we show generic strong identification for the linear-in-sums model. We provide Monte Carlo simulations in Section (ref) and conclude in Section (ref). All omitted proofs can be found in the \hyperref[appn]{Appendix}.
Throughout this paper, we use the standard probability notation, such as $O_P(\cdot)$, $\Omega_P(\cdot)$, etc. We write $f(n)\ll g(n)$ if $f = o(g)$, as well as $f(n)\lesssim g(n)$ if $f(n)=O(g(n))$. We also write $\|\cdot\|_F$ for the Frobenius norm and $\|\cdot\|_2$ for the spectral norm on matrices. We further use the notation $a\vee b=\max(a,b)$ and $a\wedge b=\min(a,b)$.
bramoulle2009identification proposed to solve the reflection problem using the network linear-in-means model. They show that in the many-group setting, endogenous and contextual effects are identified if the networks satisfy a linear independence constraint that $I, G, G^2$ are independent in the population, where $G$ is the row-normalized adjacency matrix. This identification strategy has roots in the spatial statistics literature: the network linear-in-means model is therein known as the Spatial Durbin model anselin2022spatial, and the two-stage least-squares approach for estimating this model is due to kelejian1998generalized.
In their simulations, bramoulle2009identification also recognized the potential problem of weak instruments in large, dense network settings. While this issue is not widely acknowledged in the peer effects literature, it has been explored by a few authors. In particular, lee2004asymptotic,lee2007method analyzed the properties of quasi-MLE and GMM for the linear-in-means model without contextual effects. The author recognized a case where asymptotic collinearity occurs, namely the setting of large group interactions where group size increases with the population size. In this setting, the author shows the same estimation rate $\sqrt{n/h_n}$, where $1/h_n$ is the growth rate of the entries of $G$. The rate we show under asymptotic collinearity is the same, though in our context, the collinearity is driven by local averaging instead of group interactions.
startz2017improved also acknowledge the weak instruments problem and propose an alternative estimator based on approximating the series expansion of the endogenous term. While their estimator outperforms standard two-stage least-squares in simulation, they do not show formal guarantees.
lee2002consistency considers estimation of spatial autoregressive models using OLS in a setup with only $X$ and $GY$ in the regression, where $X$ is presumed fixed. While there is endogeneity under this setup, the endogeneity is shown to be vanishing in magnitude if the network is dense enough, and thus OLS can be consistent at a slower rate, with error on the order of $O\left(\sqrt\frac{d}{n}\right)$. Our work contrasts with this by showing that if endogenous and contextual effects are both included in the model, then the OLS estimator is not consistent in general.
The regularity conditions in our work require the network degrees to be near-regular, a property which is satisfied in many statistical models of dense networks frieze2015introduction. On the other hand, lee2023qml consider estimation of spatial autoregressive models when there is a bounded amount of degree heterogeneity. They employ a novel central limit theorem to show that estimates of the model parameters are still typically consistent and asymptotically normal at a $\sqrt{n}$ rate, so long as the fraction of “dominant units” is bounded. The results in our current work could extend to such settings. However, it is possible that estimators would have non-standard asymptotic distributions when the regressors are considered random, since the limiting variance depends on the specific $X$ value of the dominant units.
In concurrent work, hayes2024peer also derive results on the asymptotic behavior of OLS for the full linear-in-means model with endogenous and contextual effects and random regressors. Under an additional subgamma condition on the covariates and error terms, they show that if the degree of the graph is $\Omega(\log n)$ then $GX$ and $GY$ converge uniformly to constants, and they further show that the OLS estimation rates for the linearly dependent terms of the model depend on $\operatorname{Tr}(G^TG)$. Our results for OLS agree with and refine some of their findings, and we are able to provide a more specific description of the reason for inconsistency (i.e., by decomposing the error into a sum of an asymptotically normal term plus bias). Additionally, we are able to provide new results for two-stage least squares. Our proof method is somewhat conceptually simpler than that of Hayes and Levin: while their work produces an approximate singular value decomposition, we directly compute the matrix inverse, which is tractable by invoking the Frisch-Waugh-Lovell theorem. By carefully considering the scaling of the moments as the network structure varies, we can derive estimation rates.
The linear-in-sums model has become more popular in recent years. liu2014endogenous formulate a model which incorporates both endogenous average and aggregate components, which they estimate jointly. Linear-in-sums models have seen applications in treatment effect estimation cai2015social,chandrasekhar2011econometrics and in network causal inference as a parametric model of interference chin2019regression. Our work provides some additional support in favor of these models by arguing that they generically do not suffer from weak identification.
For the peer effects model we describe, the observed data includes nodewise covariates $X_i\in \mathbb{R}^p$ where $p$ is the number of covariates, as well as nodewise outcomes $Y_i\in \mathbb{R}$ for each of $n$ individuals. We collect these individual covariates and outcomes into a matrix $X\in \mathbb{R}^{n\times p}$ and vector $Y\in \mathbb{R}^n$ respectively. Furthermore, we observe a binary symmetric adjacency matrix of social interactions, $A\in [0,1]^{n\times n}$. The asymptotics in this work are taken with respect to the population size, $n$: that is, if $X^{(n)}$, $Y^{(n)}$, and $A^{(n)}$ are the covariates, outcomes, and adjacency matrix for a population of size $n$, then we observe a sequence $(X^{(n)},Y^{(n)},A^{(n)})_{n\in \mathbb{N}}$, and asymptotics are taken in reference to this sequence. We will generally omit the superscript for clarity.
For an adjacency matrix $A$, we denote $D$ to be the diagonal matrix where the $i$-th entry on the diagonal is $d_i:=\sum_{j=1}^n A_{ij}=(A1)_i$, where $1$ is the $n$-dimensional vector of ones. We then define $G=D^{-1}A$ to be the row-normalized version of $A$. To ensure this is well-defined if $d_i=0$, we re-define $d_i^{-1}:=1$ if $d_i=0$,
With this notation, the classical linear-in-means model bramoulle2009identification is defined as follows:
where $\varepsilon\in \mathbb{R}^n$ is an unobserved error term and $(\alpha,\beta,\delta,\rho)$ are parameters of the model. $\alpha\in \mathbb{R}$ and $\beta\in \mathbb{R}^p$ are the parameters of a standard linear model, while $\delta\in \mathbb{R}^p$ and $\rho\in \mathbb{R}$ are network effect parameters modeling contextual and endogenous peer effects, respectively. In this paper, we will take $p=1$ so that $\delta$ is a scalar, and we will primarily be interested in estimation of the network parameters $\delta$ and $\rho$, rather than on estimation of the coefficients $\alpha,\beta$.
To simplify notation, we make the following definitions for instruments $Z$ and endogenous variables $W$, respectively:
If $\theta:= (\alpha,\beta,\delta,\rho)$, then we can rewrite the peer effects model as
The OLS estimator is defined as $$\hat{\theta}_{OLS} := (W^TW)^{-1}(W)^TY,$$ and the two-stage least squares estimator is then $$\hat{\theta}_{2SLS} := ((P_ZW)^TP_ZW)^{-1}(P_ZW)^TY,$$ where $P_Z=Z(Z^TZ)^{-1}Z^T$ is the orthogonal projection onto the column space of $Z$. Since we are working in the just-identified regime, this is equivalent to the IV estimator, $$\hat{\theta}_{2SLS} := (Z^TW)^{-1}Z^TY.$$ We will also introduce the following notation for some important network-related quantities:
In this notation,
We now describe our key assumptions on $X,\varepsilon$, and the network.
Throughout our analysis, we impose the following regularity conditions, which are standard in the literature kelejian1998generalized:
The homoskedasticity assumption on $\varepsilon$ is for convenience but can be relaxed. The assumption that $G$ and $S(\rho)^{-1}$ have bounded column sums is used to ensure the law of large numbers holds and is almost always enforced in the literature, though there is some recent work which allows it to be relaxed lee2023qml. Some assumption on the growth of the column sums of these matrices is necessary to avoid “concentration of influence” phenomena, where moments such as $\frac{1}{n}1^TGX$ would be dominated by the value of a single well-connected dominant unit; see, relatedly, discussions of the implications of this phenomenon for network economies and governance in acemoglu2012network, halpern2021defense.
Assumption (ref) is used to argue that appropriate sample moments converge; however, contrary to much existing literature, we do not make the common assumption that the relevant limiting moment matrices, such as $\Gamma_{ZW}$, are invertible.
Our main assumption for the network sequence $(A^{(n)})_{n\in \mathbb{N}}$ is the following near-degree regularity condition:
This assumption is strictly stronger than Assumption (ref). For example, Assumption (ref) allows networks where nodes of degree one are attached to a clique of growing size, while our Assumption (ref) would require that these low-degree nodes match the degree of the clique asymptotically. Nevertheless, it is well-known to be satisfied by a number of graph models, such as random graphs with degree $\omega(1)$ frieze2015introduction, and is also assumed in some existing analyses of spatial autoregressive models lee2002consistency.
An important consequence of Assumption (ref) are the following bounds, which will govern the rates of convergence for the estimators we analyze:
The last assumption we will use for our analysis is an i.i.d. assumption on the covariates:
This assumption appears benign in general and would hold if $X_i$ were a randomly administered treatment in a network experiment cai2015social,banerjee2013diffusion. However, when these two assumptions hold, $\Gamma_{ZW}$ is provably low-rank:
The proposition holds via simple computation of the limiting moment matrix. The rank-deficient limit is caused by asymptotic collinearity between $1$, $GX$, and $GY$. Intuitively, as the neighborhood of each individual is denser, more of the variation in $X$ which is useful for identification gets averaged out, so that in the limit, $GX$ approaches a constant. hayes2024peer provide a uniform convergence result of this type.
While standard arguments for least-squares consistency cannot be applied directly due to the rank deficiency, it is not immediately clear that standard estimators are inconsistent. Indeed, in similar non-network settings where an invertible covariance matrix has a singular limit, standard estimators are often consistent but at a slower rate knight2008shrinkage,phillips2016inference. In the following section, we analyze both OLS and 2SLS estimators of $\theta$ and show that while OLS is generally inconsistent, 2SLS can be consistent provided the network is not too dense.
In this section, we analyze the behavior of the OLS and 2SLS estimators for the linear-in-means model. We maintain the assumptions defined in the previous section. In addition, our results will require the following mild regularity assumption on the sequence of network moments:
The term $\operatorname{Tr}((GA)^TGB)$ is a signal term which arises from the expected covariance matrix. It is $O\left(\frac{n}{d}\right)$, which can be seen from the following:
using Cauchy-Schwarz, the inequality $\|AB\|_F\leq \|A\|_F\|B\|_2$, Lemma (ref) in the \hyperref[appn]{Appendix}, and Proposition (ref) under our near-degree regularity assumption. In many cases which we describe below, it is also $\Omega(\frac{n}{d})$, so $\frac{1}{\|G\|_F^2}\operatorname{Tr}((GA)^TGB)=\Theta(1)$. In this case, this assumption is simply asserting that this boundedness translates to a deterministic limit.
Since $\|G\|_F^2=\Theta(\frac{n}{d})$, we can alternatively assume that $\frac{d}{n}\|G\|_F^2=\frac{d}{n}\sum_{j=1}^n \frac{1}{d_j}$ converges, and that $\frac{d}{n}\operatorname{Tr}((GA)^TGB)$ converges. In this case, since $\frac{1}{n}\sum_{j=1}^n \frac{1}{d_j}$ is the reciprocal of the harmonic mean of degree sequence, the assumption that $\frac{d}{n}\|G\|_F^2$ converges is equivalent to a natural assumption that the harmonic mean is asymptotically equivalent to a constant factor times $d$. For example, it excludes sequences such as Erd\H{o}s--R\'enyi random graphs with a randomly chosen connection probability for each $n$. We choose to state the assumption using normalization by $\|G\|_F^2$ instead of the network-independent quantity $\frac{n}{d}$ since it will generalize more easily to the analysis of two-stage least squares, which has a rate that depends on the structure of the network beyond the first-order degrees.
In general, since the terms $\operatorname{Tr}((GA)^TGB)$ are upper bounded in magnitude by $O(\frac{n}{d})$, Assumption (ref) is a mild condition which requires that this upper bound translates to deterministic limits. Depending on the parameter values, some of these limits could be zero. For example, if $\beta = 0$, then
In random graphs, $\|G^2\|_F=O_P(1\vee \frac{\sqrt{n}}{d})$ and $\|G\|_F=\Theta_P(\sqrt\frac{n}{d})$, so $\operatorname{Tr}(G^TGH)=o_P(\frac{n}{d})$ and thus $m_{G,GH}=0$. On the other hand, the following lemma records which limits are guaranteed to be nonzero, as well as their sign:
We now analyze the behavior of the OLS estimator of the network effects. Recall that this estimator is defined as
Since we are primarily interested in the network effects $\delta$ and $\rho$, we partial out the intercept and $X$, which allows us to concentrate our analysis on just two parameters. Specifically, denote $W_1=(1,X)$, $W_2=(GX,GY)$, $P_1=W_1(W_1^TW_1)^{-1}W_1^T$, $M=I-P_1$, $\tilde{W}_2=MW_2$, $\theta_1=(\alpha, \beta)$, $\theta_2=(\delta, \rho)$ and $\hat\theta_2=(\hat{\delta}_{OLS},\hat{\rho}_{OLS})$. Then by the Frisch-Waugh-Lovell theorem ding2024linear, we can write \[ \hat{\theta}_2 -\theta_2 = (\tilde{W}_2^T\tilde{W}_2)^{-1}\tilde{W}_2^T\varepsilon. \] To formally state our main result, we introduce the following quantities involving network moments which are well-defined using Assumption (ref):
and $\bar\kappa:=\det(\tilde\Gamma_{WW}).$ With this notation defined, we have the following theorem, which decomposes the regression error into the sum of an asymptotically normal term and a residual bias term:
Since $\|G\|_F=\Theta(\sqrt\frac{n}{d})$ as asserted in Proposition (ref), this theorem decomposes the error into the sum of an an asymptotically normal component with variance $O(\frac{d}{n})$ and a typically non-vanishing bias. The non-convergence result contrasts with the results on OLS derived by lee2002consistency. That work assumes there is no collinearity among the regressors, and shows that OLS converges at the standard rate $\frac{1}{\sqrt{n}}$ if $d=\Omega(n^{1/2})$, or at a slower rate $O_P\left(\sqrt\frac{d}{n}\right)$ if $1\ll d\ll n^{1/2}$. The key driver of consistency is the vanishing of endogeneity when the network density increases. In our setting, we still have vanishing of endogeneity, but the baseline rate of convergence is slower due to the collinearity between endogenous and contextual effects. Thus the increasing degree of the graph simultaneously alleviates endogeneity and induces collinearity, which produces bias. It is not clear if there is a procedure to correct this bias, as it depends on the true value of $\rho$ through $m_{I,GS^{-1}}$.
Our results are qualitatively similar to those of hayes2024peer, who analyze OLS for the linear-in-means model under very similar assumptions. They show that if there is a nonzero endogenous effect, then $\hat\rho$ is inconsistent and $\hat\delta$ has a convergence rate lower bounded by $\Omega\left(\sqrt\frac{d}{n}\right)$. Our results provide more information on the bias and variance of the estimator. In particular, we show that inconsistency is due to the convergence to a bias term of constant order, and that the stochastic component of the error is asymptotically normal with a slower rate $\Theta_P\left(\sqrt{\frac{d}{n}}\right)$. Additionally, our result provides a more refined description of the behavior of $\hat\delta$, namely that if $\rho$ and $\beta$ are nonzero, then $\hat\delta$ is inconsistent as well.
Additionally, we can determine the sign of the bias for $\rho$: since $m_{G,G}$ is non-negative and $m_{I,GS^{-1}}$ is of the same sign as $\rho$ from the proof of Lemma (ref), it follows that the limiting bias for $\hat\rho$, $\frac{1}{\kappa}\sigma^2\sigma_\varepsilon^2m_{I,GS^{-1}}m_{G,G}$, is of the same sign as $\rho$, meaning that the peer effect is overestimated in magnitude.
In Section (ref) of the \hyperref[appn]{Appendix}, we discuss various cases where OLS is consistent, as implied by our results. These include the cases where $\rho=0$, or when $\beta=0$ and the network has a sufficiently small triangle density. Such examples further demonstrate the complexity of linear-in-means estimation under weak identification, as the consistency and estimation rates depend on both the underlying parameters and network structure.
In summary, this theorem provides a precise description of the behavior of OLS when contextual and endogenous effects are assumed to both be present in the model. The OLS estimator is upward-biased if there is an endogenous effect, and this bias is non-vanishing. The estimator generally converges to this biased quantity at a rate $\Theta_P\left(\sqrt\frac{d}{n}\right)$.
For 2SLS estimation, we have instruments $Z=(1,X,GX,G^2X)$. The endogenous variables $W$ are projected onto the instruments via the projection matrix $P_Z:=Z(Z^TZ)^{-1}Z^T$, resulting in the 2SLS estimator,
Since we are working in the just-identified case where the number of instruments equals the number of endogenous variables, this is equivalent to the IV estimator, that is,
A similar strategy by partialing out $(1,X)$ simplifies the analysis. The following extension of the Frisch-Waugh-Lovell theorem for linear instrumental variable estimation was first derived by giles1984instrumental and expanded upon by basu2023yule:
To analyze the two-stage least squares estimator, we will need an extra assumption on the strength of identification. Specifically, for 2SLS, identification is known to depend on the level of intransitivity in the network - for example, the identification result from bramoulle2009identification requires the existence of triplets of nodes which do not form a connected triangle. We introduce a slightly stronger assumption which essentially posits a lower bound on the density of such triplets. Let $C_3$ and $C_4$ denote the sets of $3$-cycles (i.e., triangles) and 4-cycles in the graph, respectively, where 4-cycles can either be simple cycles or three-node cycles of the forms $i\rightarrow j\rightarrow i \rightarrow k\rightarrow i$ or $i\rightarrow j\rightarrow k \rightarrow j\rightarrow i$. We assume the following:
This assumption limits the density of complete triangles in the network; it rules out networks such as unions of fully connected graphs on $d$ nodes, but holds trivially for bipartite graphs. In random graphs with expected degree $\Theta(d)$ where $d\rightarrow \infty$, the assumption is satisfied when the degree is sufficiently small. In particular, for random graphs in this asymptotic regime, standard computations of expectation and variance of $|C_3|$ show that $|C_3|=\Theta_P(E\left[|C_3|\right])=\Theta(d^3)$. Similarly, by considering the types of 4-cycles, it follows that $|C_4|=\Theta_P(nd^2+d^4)$. When $d\ll \sqrt n$, the size of $C_4$ is governed by the $\Theta(nd^2)$ term and is thus the assumption is satisfied. If $d$ is too large, i.e., $d\gtrsim \sqrt{n}$, then the graph has too high a concentration of triangles, and the assumption is not satisfied.
Assumption (ref) is related to the global clustering coefficient, a popular measure of clustering in social network analysis wasserman1994social. Denoting $C_3^{open}$ to be the set of ordered triples $i,j,k$ where $i,j,k$ are distinct and two edges exist between the three nodes, the global clustering coefficient is defined as $$\frac{|C_3|}{|C_3|+|C_3^{open}|}.$$ When the degrees are growing, boundedness of the global clustering coefficient implies Assumption (ref):
Aside from this new identification assumption, we will also need an analogue of Assumption (ref) on the convergence of network quantities involving $G^2$:
As discussed previously, this assumption is a mild regularity condition which enforces that limiting network moments exist. The additional subtracted term corresponds to variation that is lost from the partialling out step of the proof. It is typically constant, as we argue in Section (ref) of the \hyperref[appn]{Appendix}. The term $\eta$ is generally zero, except for the special case $\|G^2\|_F=\Omega(\sqrt\frac{n}{d})$, as is the case for complete bipartite graphs.
Similar to the OLS proof, this assumption allows us to define the limiting quantities which will be relevant:
Here, the quantities in the denominators are guaranteed to be nonzero by the same logic as the proof of Lemma (ref). Additionally, note that $\tilde\Gamma_{ZW}^{-1}$ is rank-one, which is a consequence of the asymptotic collinearity between $G$ and $G^2$.
We can now state a positive result for two-stage least squares:
We can see that the convergence rate can be different for networks with the same degree. In the best case, the network can be a disjoint union of complete bipartite $d$-regular graphs, which achieves a convergence rate of $O\left(\sqrt\frac{d}{n}\right)$. On the other hand, if the network is from a random graph, the rate is $\frac{d}{\sqrt{n}}$ when $d\ll n^{1/2}$.
While our positive result for random graphs is only true for $d\ll \sqrt{n}$, we conjecture that $d\gtrsim \sqrt{n}$ is the threshold for inconsistency, and that $\hat\rho -\rho$ converges to a possibly non-normal distribution. This would be analogous to what occurs in the weak instruments literature; for example, the weak IV framework of staiger1994instrumental provides a signal strength below which the asymptotic distribution of the IV estimator is a mixture of normals.
The 2SLS estimator has a higher variance for random, possibly non-bipartite graphs than the OLS estimator. This reflects the general phenomenon that IV estimators exhibit larger variances than their OLS counterparts, in spite of achieving consistent estimation in the presence of endogeneity.
In the previous section, we showed that estimation of the linear-in-means model is subject to slow convergence or inconsistency when the network degree is increasing. In this section, we introduce the linear-in-sums model as an alternative and argue that it can usually be estimated at a standard rate via two-stage least squares.
In the linear-in-means model, the estimation rate is slow because of the asymptotic collinearity between $1$ and regressors of the form $G^kX$ or $GY$. In this section, we show that the linear-in-sums (LIS) model, defined by replacing the row-normalized adjacency matrix $G$ with the original adjacency matrix $A$, generally does not suffer from the same issues as long as there is sufficient variation in the graph structure.
We first outline the baseline analysis for the LIS model, ignoring weak identification concerns. The model is defined as $$Y=\alpha 1 + X\beta + AX\delta_n + \rho_n AY + \varepsilon.$$ We emphasize that $\delta_n$ and $\rho_n$ are parameter sequences depending on $n$, which is necessary to make $(I-\rho_n A)^{-1}$ well-defined. Henceforth we remove the subscript for ease of notation. We maintain the same assumptions as with the LIM model, including near-degree regularity, i.e., the degrees $d_i$ are uniformly in the interval $(c_1d,c_2d)$ for some constants $c_1,c_2$. To ensure that the model is well-defined and to keep $Y$ on a constant scale, we enforce that $|\rho|<1/\lambda_1(A)$, which implies that the inverse $(I-\rho A)^{-1}$ exists and can be written via the Neumann series, $(I-\rho A)^{-1}=\sum_{k=0}^\infty (\rho A)^k$. We further assume that $d\rho \rightarrow_n \rho_0$ for some constant $\rho_0$; since $\lambda_1(A)=\Theta(d)$ under our near-degree regularity assumption, this implies $\rho$ is on the same scale as $1/\lambda_1(A)$. We impose an analogous restriction on $\delta$: $d\delta \rightarrow_n \delta_0$ for some constant $\delta_0$. We do not need an extra boundedness condition on $\delta$ since the contextual effect $\delta$ does not show up in the inverted term, $(I-\rho A)^{-1}$.
Since the regressors are on different scales, the baseline convergence rate is different for the four parameters, as $A^kX$ typically scales as $O(d^k)$. To analyze this, we thus consider rescaling matrices $F=\mathrm{diag}(1,1,1/d, 1/d^2)$ and $H=\mathrm{diag}(1,1,1/d, 1/d)$ for the instruments $Z$ and endogenous $W$ respectively. For the 2SLS estimator, we have
Along with the same conditions on $\varepsilon$ imposed in the previous section, the standard 2SLS assumptions are now imposed on the rescaled covariates, $ZF$ and $WH$:
Under these assumptions, standard arguments yield the following:
That is, the network components of $\hat\theta$ converge at the rate $O_P\left(\frac{1}{d\sqrt{n}}\right)$ - equivalently, the rate for the estimating the $d$-rescaled coefficients is standard root-n.
We can provide conditions under which this invertibility assumption holds and thus the collinearity issues for linear-in-means are avoided. First, since $X$ is not responsible for collinearity as long as its variance is positive, we only need to consider collinearity among the intercept and network features (this is equivalent to partialling out $X$). We first examine $\frac{1}{n}F(Z^TZ)F$:
To provide a broad class of networks where we can show convergence and invertibility of the matrix, we consider $A$ sampled from a graphon, which is commonly used to model dense networks bickel2011method:
where $U_i$ are i.i.d., $A$ is symmetric, and $f:[0,1]\times [0,1]\mapsto \mathbb{R}_{\geq 0}$ is a symmetric, integrable function known as a graphon. In this setting, $p_n$ governs the sparsity of the graph, and we can take $d=np_n$. The following convergence result is adapted fro avella2018centrality:
Consequently, $\left(\frac{1}{n}F(Z^TZ)F\right)_{ij}\rightarrow m_{i+j-2}$. In the language of graph limits, $m_k$ is the homomorphism density of the length-$(k{-}1)$ path in the graphon $f$. With convergence of moments settled, we can derive conditions for invertibility of the limiting moment matrix. In the following, we provide two conditions, the first of which only uses integrals of the graphon $f$, and the second of which is based on the spectral decomposition of the graphon.
We assume here that the instruments are $1,AX,A^2X$. Define the degree function $g_1(u)=\int_v f(u,v)$, and the codegree function $g_2(u)=\int_{v,w} f(u,v)f(v,w)$. Then the moment matrix has the following structure: \[ \Gamma_{ZZ}=
. \] The decomposition can be verified by re-arranging the integrals: for example, the bottom right entry is
From this representation, the following necessary and sufficient condition is immediate:
The definitions of $g_1$ and $g_2$ are related to the link and codegree functions used by auerbach2022identification to identify a partially linear model with graphon-controlled semiparametric components. Our $g_1$ and $g_2$ are the integrals of the agent link and codegree functions as defined in Auerbach's work. The intuition is that variation in the expected degrees of each node is required for $1$ and $g_1$ to be linearly independent; furthermore, there needs to be excess variation in the number of two-hop neighbors to achieve identification. The condition is also similar to the standard assumption usually employed in linear-in-means models for identification, namely that $I,G,G^2$ are linearly independent. However, in the growing network setting, the corresponding identification assumption for linear-in-means does not imply that the limiting $\Gamma_{ZZ}$ is invertible, even if we impose the graphon assumption. This is because the row-normalization forces the constant vector to be an eigenvector of $G$, so that $GX$ and $G^2X$ limit to constant vectors when $X$ is IID.
By appealing to the spectral decomposition of a graphon, we can provide alternative conditions for identification which are natural for low-rank network models, such as the stochastic block model lee2019review. Following the discussion from avella2018centrality, we define the graphon operator:
Graphon operators admit a spectral decomposition:
With this decomposition in mind, we can write $$m_k=\int_u (T_f^k1)(u)du=\left\langle 1,\sum_{i=1}^\infty \lambda_i^k\langle \phi_i,1\rangle \phi_i\right\rangle=\sum_{i=1}^\infty \langle \phi_i,1\rangle^2\lambda_i^k.$$ Consequently, we can write the moment matrix in the form $\sum_{i=1}^\infty a_i^2b_ib_i^T$, where $a_i=\langle \phi_i,1\rangle$ and $b_i=(1,\lambda_i,\lambda_i^2)$. We can see that this matrix is full rank if and only if there are three linearly independent $b_i$ corresponding to $a_i\neq 0$. Using the Vandermonde structure of $M$ horn2012matrix, it follows that any set of three $b_i$'s are linearly independent if and only if the eigenvalues $\lambda_i$ are distinct. Combining these observations yields the following identification condition:
The theorem generalizes to the case where multiple instruments are used: to distinguish between $k+1$ values $m_0,...,m_k$, we now require $k+1$ distinct eigenvalues.
In the following section, we provide an example of identification in stochastic block models.
The stochastic block model is a common model of communities in social networks lee2019review. The model is parameterized by a matrix $P\in [0,1]^{K\times K}$ and vector $\pi\in [0,1]^K$, where $K$ is the number of communities, $P_{ij}$ is the probability of connection between individuals in communities $i$ and $j$, and $\pi_i$ is the probability an individual is in community $i$. For each individual, a single community membership is drawn independently, and conditioned on the memberships, edges are sampled independently of one another.
We further define $E_{SBM}=PQ$, where $Q=\mathrm{diag}(\pi)$. The spectral properties of the graphon $f_{SBM}$ can be derived from those of the associated SBM, as shown in the following proposition:
Thus for the SBM, we have identifiability if and only if $E_{SBM}$ has at least 3 distinct eigenvalues corresponding to eigenvectors not orthogonal to $1$. We can expect this to be true for most values of $E_{SBM}$ when the number of communities, $K$, is at least $3$. Identification can fail because there are fewer than 3 distinct eigenvalues, but it can also fail because the number of eigenvectors which have nonzero inner product with the constant vector, $1$, is too small. For example, consider a disconnected community matrix, \[ E_{SBM}=\frac{1}{3}
. \] The eigenvalues, $1/3$, $5/6$, and $1/6$, are distinct, and $1$ is not an eigenvector, but the last eigenvector $(0,1,-1)^T$ is orthogonal to 1. It turns out that even by enforcing strict positivity of entries and positive-definiteness of $E_{SBM}$, one can numerically construct values of $E_{SBM}$ that satisfy the distinct eigenvalue condition but not the orthogonality condition, such as the following matrix for which the eigenvalues are distinct, but the bottom eigenvector is orthogonal to $1$: \[
\] However, violation of the identifiability condition is still an edge case, as it is generally restored after a small perturbation of $E_{SBM}$.
Aside from the condition that $\Gamma_{ZZ}$ is full-rank, we also need the condition that $\Gamma_{ZW}$ is full-rank. This is the instrument relevance condition, namely, that the instrument $A^2X$ is correlated with the endogenous $AY$. We argue that the conditions above are sufficient. First, we note that from calculations involving the sampled graphon (i.e., Theorem 2 from avella2018centrality), moments of the form $X^TA^mAY$ behave as
where $h(\lambda)=\frac{1}{1-\rho_0 \lambda}\lambda(\alpha+\mu\beta+\mu\delta_0\lambda)$, and $a_k$ and $\lambda_k$ are as defined in Section (ref). Thus with the rescaling matrix definitions of $F$ and $H$ provided earlier, we have
This limiting matrix can be written as $BC^T$, where $B$ has $b_k:=(1,\lambda_k, h(\lambda_k))$ as columns and $C$ has $(1,\lambda_k,\lambda_k^2)$ as columns. Under the conditions in Theorem (ref), $C$ is full rank. Subsequently, identification depends on the rank of $B$. We now show that the same conditions generally suffice for invertibility of $\Gamma_{ZW}$, and so under the conditions of Theorem (ref), the 2SLS estimator of the LIS model converges at a standard rate.
In this section, we validate our theory through a variety of numerical experiments.
We compare the behaviors of OLS and 2SLS for an Erd\H{o}s--R\'enyi graph. For each $n=100,200,...,2000$, we sample an Erd\H{o}s--R\'enyi random graph with expected degree $d$, where $d=O(n^{1/2}),O(n^{1/4})$. The constant is chosen so that for $n=100$, the expected degree is $10$. We fix $\alpha=1$, $\beta=1.5$, $\delta=0.6$, $\rho=0.3$. The $X_i$ are drawn from a normal distribution with mean $2$ and variance $1$, and the $\varepsilon_i$ are drawn from a normal distribution with standard deviation $0.1$. The resulting values of $\hat{\delta}-\delta$ and $\hat{\rho}-\rho$, along with 95 percent pointwise confidence bands, are plotted in Figures (ref) and (ref).
From these figures, it is evident that there is a nonzero bias in the OLS estimates for $\delta$ and $\rho$, as predicted by Theorem (ref). The sign of the bias is predictable from the expressions we derived, i.e., it is negative for the contextual effect and positive for the endogenous effect when $\rho>0$. We can see in Figure (ref) that the estimator error is converging to the bias. We expect this to be true as well in Figure (ref) for the denser graph $d=n^{1/2}$, but we believe the sample size is too small to see the convergence, since the bias is still increasing significantly at $n=2000$.
In contrast, the 2SLS estimator does not have bias, i.e., the orange line meets the dashed black line at zero. For $d=n^{0.5}$, the estimate for $\rho$ appears to maintain a constant amount of variation. This agrees with our analysis, which suggested that for random graphs, $d\gtrsim \sqrt{n}$ is the threshold where $\hat{\rho}_{2SLS}$ is inconsistent and has a non-degenerate limiting distribution. For moderate values of density $d=n^{1/4}$, we can see that the estimates from 2SLS converge at a rate slower than $\sqrt{n}$.
In the discussion after Theorem (ref), we suggested a surprising result that when $\beta=0$, OLS could be asymptotically unbiased for $\delta$ if the graph structure is not highly clustered. In Figure (ref), we validate that this is indeed the case: for random graphs, as we described, OLS appears to be asymptotically unbiased for $\delta$ and closely matches the behavior of 2SLS. However, the estimate for $\rho$ still admits a non-zero bias.
For this experiment, we generated a disjoint union of $d$-regular complete bipartite graphs, with $n$ nodes total. This union of bipartite cliques has different properties from random graphs: for random graphs, $\operatorname{Tr}(A^3)=\Theta(d^3)$ and $\operatorname{Tr}(A^4)=\Theta(d^4)$, while for this ensemble, $\operatorname{Tr}(A^3)=0$ and $\operatorname{Tr}(A^4)=\Theta(nd^3)$, i.e., it has a minimal number of triangles and maximal number of 4-cycles. Our theory suggests that when $\beta=0$, the OLS estimate of $\delta$ is still biased, while the 2SLS estimate of $\rho$ enjoys an improved convergence rate of $O\left(\sqrt\frac{d}{n}\right)$.
Figure (ref) demonstrates that this is indeed the case. Although OLS estimation of $\delta$ was asymptotically unbiased for random graphs when $\beta=0$, for the union of complete bipartite graphs, there is a non-vanishing bias term. On the other hand, $\hat{\rho}_{2SLS}$ enjoys an improved convergence rate, as compared visually with Figure (ref). This example highlights the complexities of the linear-in-means model: the parameter values and graph structure can affect OLS and 2SLS estimators differently.
Finally, we consider the case where the network is a disjoint union of $d$-cliques. In this case, our theory for OLS estimation still holds, and we expect that the estimator will converge to a biased solution. However, for the 2SLS theory, we prohibited such networks in our identifiability assumption. In Figures (ref) and (ref), we plot the results of estimating the linear-in-means model in this regime. We can see that the OLS estimates, while biased, are stable. However, the 2SLS estimates are wildly unstable. This is to be expected: it is known that Manski's reflection problem can be solved from an identification perspective by excluding oneself from their own neighborhood for group averaging bramoulle2020peer, but standard identification does not imply strong identification in our asymptotic framework. The weak instrument problems thus lead to highly unreliable estimates.
In this work, we investigated the phenomenon of weak identification in the popular linear-in-means model of peer effects. By concentrating our analysis on the network contextual and endogenous peer effects, we showed that the behavior of OLS and 2SLS estimators can be complex when the degrees of the network are allowed to grow: the model parameters and network structure can contribute jointly to produce biases or slower rates of convergence. We also showed that this weak identification is generally not a problem in the linear-in-sums model: for a large class of random graphs, the linear-in-sums model is strongly identified. Our results indicate that the linear-in-sums model can be efficiently estimated in settings where linear-in-means models suffer from slow rates.
One broader theme that is reflected in this work is the difference between OLS and instrumental variable estimation methods. The weak identification literature has provided diagnostics to determine whether or not an instrument can reliably used in IV estimation staiger1994instrumental. For the linear-in-means model, while instrumental variables generally outperforms ordinary least-squares, when there truly is weak identification, as in the case of the union of cliques, instrumental variables estimation is unreliable. In such a setting, OLS is biased, but can still potentially be used to say something informative about the parameters of interest if the bias is small or known to be in a particular direction.
There are a few interesting possible extensions to this work. The asymptotic framework we studied corresponds to the infill asymptotic setting, where a single network is observed. In contrast, some applications fit in the increasing-domain asymptotic framework, where many isolated communities, such as villages, are observed. In such cases, it is common to incorporate village-level fixed effects into the model. This creates an incidental parameter problem that is often resolved by differencing. Since this procedure removes useful information from the data, a stronger condition based not only on triangles, but also on fourth-order motifs, is required to ensure identification bramoulle2009identification. Once this differencing is performed, identification has been shown to be impossible in both complete graphs and complete bipartite graphs. Similarly, we would expect that the identification conditions for 2SLS in our setup would need to change in order to limit the number of triangles and 4-cycles; it is likely that the estimation rate would then depend on quantities such as $\|G^3\|_F$ instead. We leave this extension for future work.
Additionally, it would be interesting to see how the results extend to networks with degree heterogeneity, as discussed in the related work. Degree heterogeneity is known to be widespread but problematic in other areas of statistical network analysis, such as community detection, as it slows down concentration of useful random matrix quantities ke2024optimal. However, this could potentially be beneficial for the weak instrument problem, since rapid averaging, or concentration, is the cause of weak instruments in our setup. Nevertheless, degree heterogeneity poses other problems for inference: it may be the case that the limiting covariance matrix depends on the values of covariates for specific individuals due to concentration of influence acemoglu2012network,halpern2021defense. This would make inference non-standard when interpreting regressors as random, as is commonly assumed in econometrics and statistics abadie2020sampling,buja2019models. Studying the impact of degree heterogeneity seems to be a complex but worthwhile endeavor that can help us understand how the theory extends to more practically relevant network models.
Another question is whether or not the estimation rates derived are fundamental lower bounds over the class of all estimators. In particular, we can ask if there are estimators which dominate the two-stage least squares estimator. It is known that in the strongly-identified case, there is a three-stage procedure which is efficient among estimators using $GX,G^2X,...$ as instruments lee2003best,zhang2010more. How do these estimators perform in weakly identified setups? Additionally, it may be possible for regularized estimators to outperform their unregularized counterparts in this setup, as suggested in some related work on weakly identified models knight2008shrinkage,tchuente2019weak.
Finally, it would be worthwhile to explore further connections with the weak instruments literature. A large body of work has been devoted to developing finite-sample approximations of estimator distributions staiger1994instrumental,bound1995problems,andrews2019weak or devising robust testing procedures which have the correct size regardless of the strength of identification mikusheva2010robust. In settings where peer effects are weakly identified, can these procedures be used to improve estimates? Developing improved techniques for inference could help increase the credibility of peer effects research.
This work was supported by ARO MURI W911 NF-19-1-0217 and a Vannevar Bush Fellowship from the Office of the Under Secretary of Defense.