The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
92,407 characters
Regression with Observational Multilayered Network Data
\begin{bibunit}[jpe]
\def\spacingset#1{\renewcommand{1.3}
{#1}\small\normalsize} \spacingset{1}
\if11
{
\title{Regression with Observational Multilayered Network Data}
\author{Juan Estrada\thanks{Analysis Group Economic Consulting, Washington, DC, USA. \faEnvelopeO: [email removed].} \and Kim P. Huynh\thanks{Department of Economics, Indiana University, 100 S Woodlawn, Bloomington, IN 47405, USA. The Laboratoire d’\'Economie d’Orl\'eans, Universit\'e d'Orl\'eans, Orl\'eans, France. \faEnvelopeO: [email removed].} \and David T. Jacho-Ch\'{a}vez\thanks{Corresponding Author: Department of Economics, Emory University, Rich Building 306, 1602 Fishburne Dr., Atlanta, GA 30322-2240, USA. \faEnvelopeO: [email removed].} \and Leonardo S\'{a}nchez-Arag\'{o}n\thanks{Facultad de Ciencias Sociales y Human\'{i}sticas, Escuela Superior Polit\'{e}cnica del Litoral, ESPOL, Campus Gustavo Galindo Km. 30.5 V\'{i}a Perimetral, P.O. Box 09-01-5863, Guayaquil, Ecuador. \faEnvelopeO: [email removed].}}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Regression with Observational Multilayered Network Data}
\end{center}
\medskip
} \fi
\bigskip
\begin{abstract}
\noindent A novel method to estimate social effect coefficients in the popular so-called linear-in-means regression model in the Social Sciences is presented here that utilizes non-experimental multidimensional network data. The procedure can accommodate social interactions that correlate with the error in the model by making use of a different set of network links among the same observations that are exogenous in the traditional sense. In particular, the full observability of a two-layered \emph{multiplex} network data structure is assumed here to propose a new Generalized 3-Stage Least Squares (G3SLS) estimator that is consistent, asymptotically normally distributed, and also easy to implement using widely-used existing statistical software because of its closed-form definition. The underlying assumptions are general enough to accommodate common problems with observational data such as measurement error, simultaneity, and unobserved heterogeneity. Monte Carlo exercises confirm the good small sample performance of the proposed G3SLS estimator in these scenarios. An empirical application finds positive and significant peer effects in citations among research articles published in top general-interest journals in economics.
\end{abstract}
\noindent
{\it Keywords:} Instrumental Variables; Linear-in-Means Models; Multidimensional Networks; Multilayered Networks; Multiplex Networks
\\
\noindent
{\it JEL code:} A1, C21, C31, C51, I23, J24.
\vfill
\newpage
\spacingset{1.45}
\section{Introduction}\label{sec:intro}
In social science research, understanding the causal mechanisms behind individual outcomes is a central challenge. With network data, the outcome of a unit, such as the performance of a firm or the academic achievement of a student, is often influenced by its own characteristics and those of its peers (those to whom they are connected). The \emph{linear-in-means model} is a standard tool for quantifying these interactions, which are broadly categorized as direct ($\boldsymbol{\gamma}$), peer ($\beta$), and contextual effects ($\boldsymbol{\delta}$). The model is formalized as follows.
\begin{equation}
\label{intro}
y_{i}=\alpha+\beta\sum_{j \neq i}w_{i,j}y_{j}+\sum_{j \neq i}w_{i,j}\mathbf{x}_{j}^{\top}\boldsymbol{\delta}+ \mathbf{x}_{i}^{\top}\boldsymbol{\gamma}+v_{i}\text{,}
\end{equation}
\noindent where $i, j \in \{1, \dots, n\}$ represents the nodes within the network, $\mathbf{x}$ is a vector of attributes associated with the units, and $w_{i,j}$ takes a value of 1 if the nodes $i$ and $j$ are connected (indicating a link or edge) and 0 otherwise. The term $v_{i}$ accounts for an unobserved random error, and $n$ is the total number of nodes (or observations) in the network. The entire network structure is represented by an adjacency matrix $n \times n$, $\mathbf{W}$, where the element at position $(i,j)$ corresponds to $w_{i,j}$. The parameters in \eqref{intro} are straightforward to identify when the network structure is exogenous; that is, it is not correlated with the unobserved error term of the model. The error term $\mathbf{v} = [v_1, \dots, v_n]^{\top}$ has an expected value of zero conditional on the observed covariates and the network structure, i.e., $E[\mathbf{v}|\mathbf{X}, \mathbf{W}] = \mathbf{0}$, where $\mathbf{X} = [\mathbf{x}_1, \dots, \mathbf{x}_n]^{\top}$ is the matrix of covariates; see \cite{Paula2017} and references therein.
However, when the assumption of exogeneity does not hold—meaning that $E[\mathbf{v}|\mathbf{X}, \mathbf{W}] \neq \mathbf{0}$—the situation becomes more complex. This occurs when the network structure or the error term is correlated with the covariates, which introduces endogeneity. Endogeneity can arise from various sources, such as simultaneity bias (where outcomes and connections influence each other), errors in the measurement of network links, or homophily (the tendency of similar units to form connections). This challenge remains a focal area of study in social science research, as highlighted by works such as \cite{Johnsson2019} and \cite{Chan_et_al_social_effects}.
As in \cite{Chan_et_al_social_effects}, our paper uses a type of multilayered (multidimensional) network data structure known as \emph{multiplex} networks to consistently estimate and perform correct inference on the structural social parameters in \eqref{intro} with a potentially endogenous network structure $\mathbf{W}$. In particular, it is assumed that the researcher observes another set of social ties between the original observations $\{w_{0;i,j}\}_{j=1,j\neq i}^n$ in the form of an adjacency matrix $\mathbf{W}_0$ $n\times n$, which is exogenous in the usual sense, i.e., $E_{\mathbf{X},\mathbf{W}_0}[\mathbf{v}]=\mathbf{0}$, instead of leaving the correlation structure between the errors $\mathbf{v}$ and the original social structure $\mathbf{W}$ in \eqref{intro} unspecified. This type of data structure is becoming increasingly popular in the Social Sciences; see, e.g., \cite{Jackson_Multiplex} in Anthropology, \cite{manta2021} in Economics, \cite{Aldasoro_Alves_JFinStab} in Finance, \cite{An_et_al_2026} in Sociology, and \cite{PoliSci_Multiplex} in Political Science. See \cite{Boccaletti2014} and \cite{Kivela_multilayer_network_2014} for an up-to-date survey of the mathematical underpinning of this type of data.
Our proposed estimator of parameters in \eqref{intro} is performed in three steps; we shall hereafter refer to it as the Generalized Three-Stage least squares estimator (G3SLS). Using the linearity of the model and simple linear projection arguments in each step, the resulting estimation procedure is very simple to implement. It is readily available in \texttt{Stata}, i.e., \cite{netivreg_stata_journal}. Furthermore, the estimator is shown to be asymptotically normal at the standard root--$n$ rate of convergence. The estimation effect from the multistage procedure is fully characterized, and a consistent asymptotic variance-covariance estimator is proposed as well.
Our approach is related to a growing literature that uses instrumental variables for identification in linear-in-means models, notably \cite{Kelejian_et_al_2014}, \cite{Konig_et_al_2019}, and \cite{Lee_et_al_2021}. These papers generally address network endogeneity by constructing instruments derived from observable dyadic attributes or estimated link formation probabilities. In contrast, the G3SLS estimator proposed here leverages the multiplex data structure by assuming that the researcher observes a distinct exogenous layer $\mathbf{W}_0$ that correlates with the endogenous layer $\mathbf{W}$. Unlike \cite{Kelejian_et_al_2014}, who achieves identification by assuming that the elements of the endogenous weighting matrix are linear functions of observable exogenous dyadic variables (such as distance or border length) and then estimating these weights to construct the instrument matrix, our identification strategy utilizes a linear projection onto a distinct, observed exogenous network layer $\mathbf{W}_0$. Furthermore, unlike \cite{Konig_et_al_2019}, who addresses endogeneity by explicitly modeling the strategic network formation process to construct instruments based on predicted link probabilities, our method does not require a specific description of the network formation process, allowing it to remain robust to different sources of endogeneity, such as simultaneity and measurement error. Additionally, while \cite{Lee_et_al_2021} constructs instruments by weighting peers' attributes with predicted link probabilities derived from a logistic regression on exogenous dyadic characteristics, our method directly utilizes the observed connections in the exogenous network $\mathbf{W}_0$ for identification, thereby bypassing the need to estimate a link formation model.
The Monte Carlo exercise showcases the versatility of the proposed estimator by analyzing its performance in three separate data generating processes that consider different scenarios, such as omitted variable bias, measurement error, and unobserved homophily in network formation. The simulation results show that the proposed G3SLS estimator performs well in terms of bias and root mean squared error, with sample sizes as small as 50 observations. In addition to the simulated data, this paper also presents an application to real data on publication outcomes in Economics \citep{netivreg_g3sls}. The use of web scraping and existing data on authors' research fields, education, and employment history allows for the creation of two types of professional ties among scholars, namely co-authorship and alumni connections. These multiplex networks are then used to uncover positive and significant peer effects in terms of citations among articles published by these scholars, as well as significant positive effects of research teams that are gender diverse on the quality of a paper, measured in terms of citation outcomes after controlling for other articles' characteristics, such as the number of pages, number of bibliographic references, co-authoring with current and previous editors, and various network fixed effects.
The structure of the paper is as follows. Section \ref{background} defines multiplex networks with an example and provides various scenarios in which the required exogeneity requirement for one of the network layers can be naturally satisfied. Section \ref{indet} introduces the model and identification conditions for the parameters of interest. Section \ref{est} describes the multi-step estimation procedure, asymptotic properties, and valid asymptotic standard errors. Section \ref{mc} presents the small sample properties of the proposed estimator in various Monte Carlo exercises, while Section \ref{emp} discusses an empirical application. Finally, Section \ref{conclusion} concludes.
\ref{Appendix_A} and \ref{Appendix_B} contain, respectively, all mathematical derivations of the main results and the conditions used to establish the asymptotic properties of the proposed estimator, while intermediate steps are collected in \ref{Appendix_C}. The supplementary material reports further numerical experiments evaluating the proposed estimator across various scenarios, see \ref{Appendix_D}. Finally, \ref{Appendix_E} and the references therein provide details on the real data application and several robustness checks.
\section{Background\label{background}}
\subsection{Multiplex Networks}
The cornerstone of the identification of social parameters in \eqref{intro} in this paper is the complete observability of more than one type of social interaction between economic agents. These data structures are known as \emph{multi-dimensional} or \emph{multi-layered} networks; see, e.g., \cite{Boccaletti2014} and \cite{Kivela_multilayer_network_2014} for up-to-date comprehensive surveys and references therein. Following the definition of \citeauthor{Boccaletti2014} (\citeyear{Boccaletti2014}) a multilayer network is a pair $\mathcal{M}=(\mathcal{G},\mathcal{C})$, where $\mathcal{G}=\{g_{m};\quad m \in \{1,\dots,M\}\}$ is a family of $M$ graphs, and $\mathcal{C}$ is the set of interconnections between nodes of different layers $g_{\alpha}$ and $g_{\beta}$ with $\alpha\neq \beta$. When the same nodes are in each layer and there are no connections between different nodes in different layers except with themselves, these networks are called \emph{multiplex}; see, e.g., \cite{Jackson_Multiplex}. The latter is the specific type of data structure used in this paper for identification.
Figure \ref{F1}(a) presents an example of a two layer multiplex network, where there are two different types of connections between the same four nodes, i.e., the blue and red edges. The network $G_{0}$ depicted with the blue edges represents a (possibly predetermined) network, while the red edges in $G$ represent another (possibly endogenous) network. In this example, $M=2$, $g_{\alpha}=\{(1,2),(1,4),(3,4)\}$, and $g_{\beta}=\{(1,2),(1,3),(3,4)\}$ are described in the corresponding adjacency matrices $\mathbf{W}$ and $\mathbf{W}_0$ in panel (c). The relationship between the networks is given by the different types of edges that are shared by the same nodes; see, i.e., the flattened version of the multiplex network in Figure \ref{F1}(b). In this example, many of the red connections exist where the blue connections also exist, showing an important relationship between the networks. The matrix product $\widetilde{\mathbf{W}}\equiv\mathbf{W}_{0}\mathbf{W}$ in panel (c) encapsulates a measure of correlation between the networks that is relevant for the identification and estimation results presented in Sections \ref{indet} and \ref{est} below, respectively. The $(i,i)$th element of $\widetilde{\mathbf{W}}$ represents the number of nodes $j \neq i$ that are connected to node $i$ by both types of edges (blue and red), e.g., the position $(1,1)$ of $\widetilde{\mathbf{W}}$ equals one because only node $2$ is connected with node $1$ by both types of edges. Additionally, the $(i,j)$th element of $\widetilde{\mathbf{W}}$ contains the number of two length paths that start with a blue connection from $\mathbf{W}_{0}$ and are followed by a red connection from $\mathbf{W}$. These types of paths are called \textit{inter-layer intransitive triads} (`a friend's relative in the family layer is not a friend in the friendship layer') hereafter. For example, the position $(1,4)$ in $\widetilde{\mathbf{W}}$ is equal to one because there is an inter-layer intransitive triad that changes the edge colors connecting nodes 1 and 4 via node 3.
\begin{figure}[!htb]
\caption{Example of a Two-Layered Multiplex Network and Their Adjacency Matrices}\label{F1}
\centering
\begin{minipage}{0.33\textwidth}
\vspace{7mm}
\begin{tikzpicture}[multilayer=3d]
\Vertex[x=-0.1, y = 0.2, label = 1,layer=1, color = red!60]{1b}
\Vertex[x=0.5, y = 0.8, label = 2,layer=1, color = red!60]{2b}
\Vertex[x=1.7, y = 0.1, label = 3, layer=1, color = red!60]{3b}
\Vertex[x=1, y = -1, label = 4,layer=1, color = red!60]{4b}
\Edge[color = red](1b)(2b)
\Edge[color = red](3b)(4b)
\Edge[color = red](1b)(4b)
\Vertex[x=0, y = -0.4, label = 1,layer=2, color = blue!60]{1a}
\Vertex[x=0.5, y = 0.4, label = 2,layer=2, color = blue!60]{2a}
\Vertex[x=1.7, y = -0.2, label = 3, layer=2, color = blue!60]{3a}
\Vertex[x=1.1, y = -1.6, label = 4,layer=2, color = blue!60]{4a}
\Edge[color = blue](1a)(2a)
\Edge[color = blue](1a)(3a)
\Edge[color = blue](4a)(3a)
\SetLayerDistance{-3}
\Plane[x=-0.6,y=-1.7,width=3,height=3,color=red!40,layer=1]
\Plane[x=-1.7,y=-0.7,width=3,height=3, color=blue!40,layer=2]
\Edge[style=dashed](1a)(1b)
\Edge[style=dashed](2a)(2b)
\Edge[style=dashed](4a)(4b)
\Edge[style=dashed](3a)(3b)
\Text[x = -0, y = -2.2, layer = 1, color = red, opacity = 2, rotation=30]{$G$}
\Text[x = -0.3, y = -1.5, layer = 2, color = blue, opacity = 2, rotation=30]{$G_0$}
\end{tikzpicture}
\begin{center}
Panel (a)
\end{center}
\end{minipage}
\hfill
\begin{minipage}{0.33\textwidth}
\flushright
\vspace{16mm}
\begin{tikzpicture}
\Vertex[x=0, y=0, color = white, label = 1]{1}
\Vertex[x=1,y=1, color = white, label = 2]{2}
\Vertex[x=2, y=0, color = white, label = 3]{3}
\Vertex[x=1, y=-1, color = white, label = 4]{4}
\Edge[color = blue, bend = -30](1)(2)
\Edge[color = blue](1)(3)
\Edge[color = blue, bend = -30](3)(4)
\Edge[color = red, bend = 30](1)(2)
\Edge[color = red](1)(4)
\Edge[color = red, bend = -30](4)(3)
\Text[x = -0.5 , y = -0.2, layer = 1, color = red]{$\mathbf{G}$}
\Text[y = 2.3, x = 0.5, layer = 1, color = blue]{$\mathbf{G}_{0}$}
\end{tikzpicture}
\vspace{18mm}
\begin{center}
Panel (b)
\end{center}
\end{minipage}
\hfill
\begin{minipage}{0.30\textwidth}
\vspace{5mm}
\centering
\footnotesize
\[
\textcolor{red}{\mathbf{W}}=
\begin{bmatrix}
0 & \textcolor{red}{1} & 0 & \textcolor{red}{1} \\
\textcolor{red}{1} & 0 & 0 & 0 \\
0 & 0 & 0 & \textcolor{red}{1} \\
\textcolor{red}{1} & 0 & \textcolor{red}{1} & 0
\end{bmatrix}
\]
\[
\textcolor{blue}{\mathbf{W}_{0}}=
\begin{bmatrix}
0 & \textcolor{blue}{1} & \textcolor{blue}{1} & 0 \\
\textcolor{blue}{1} & 0 & 0 & 0 \\
\textcolor{blue}{1} & 0 & 0 & \textcolor{blue}{1} \\
0 & 0 & \textcolor{blue}{1} & 0
\end{bmatrix}
\]
\[
\textcolor{blue}{\mathbf{W}_{0}}\textcolor{red}{\mathbf{W}}=
\begin{bmatrix}
1 & 0 & 0 & 1 \\
0 & 1 & 0 & 1 \\
1 & 1 & 1 & 1 \\
0 & 0 & 0 & 1
\end{bmatrix}
\]
\normalsize
\vspace{5mm}
\begin{center}
Panel (c)
\end{center}
\end{minipage}
\vspace{0.3cm}
\begin{minipage}{1\textwidth}
\footnotesize
Note: Panel (a) displays an undirected two-layered \emph{multiplex} network structure among four nodes. The red edges showcase a possibly endogenous network, while the blue edges display exogenous connections among these agents -- with associate adjacency matrices in Panel (c). The figure in panel (b) displays the flattened version of the multiplex network. It facilitates the visualization of the inter-layer intransitive triads such as those between nodes 1 and 4 or nodes 2 and 4, i.e., matrix $\widetilde{\mathbf{W}}$.
\end{minipage}
\end{figure}
\subsection{Some Examples\label{some_examples}}
In economics, \cite{Goldsmith-Pinkham2013} and \cite{Chan_et_al_social_effects} have used multidimensional network data structures for causal inference. \cite{Goldsmith-Pinkham2013} used previously observed sets of connections of the same type as explanatory variables in their endogenous network formation model, while \cite{Chan_et_al_social_effects} used different types of peers, i.e., roommates, classmates, seatmates, studymates, and friends, to study peer effects on academic performance. The identification strategy put forward in this paper utilizes the latter type of multiplex networks by simply pointing out that some of these types of connections are exogenous in nature (blue edges in Figure \ref{F1}(b)), e.g., `roommates' \citep{Sacerdote_QJE}, `ethnic groups' \citep{Reza2019}, and `classmates' \citep{Ammermueller_Pischke_JLO}, since they are not directly chosen by the units of observation, while others are likely endogenous (red edges in Figure \ref{F1}(b)), e.g., `friends' \citep{Fruehwirth_RevStat} or `studymates' \citep{Chan_et_al_social_effects} due to potential homophily in unobserved characteristics.
Similarly, in applications about publication outcomes in Academia, see, e.g., \cite{Newman_2004a} in Physics or \cite{AoAS_coauthorship} in Statistics, co-authorship connections are potentially endogenous (red edges in Figure \ref{F1}(b)). Information relating to authors who received their Ph.D. from the same institution can be considered pre-determined (blue edges in Figure \ref{F1}(b)) and therefore exogenous. The empirical application in Section \ref{emp} below utilizes this observation in an application where the outcome variable in \eqref{intro} is the natural logarithm of citations of a peer-reviewed research article published in the top journals of general-interest in economics between 2002 and 2006.
Finally, it should be noted that the best way to guaranty the exogeneity of the network $G_{0}$ would be by (quasi) randomizing individuals into groups. This has been done before. For example, \cite{carrell2013} randomly assigned freshmen students at the United States Air Force Academy to peer groups (squadrons), while \cite{falk2006} randomly assigned workers to groups in the lab. However, as pointed out by \cite{carrell2013}, these types of experimental designs do not consider the endogenous process of individuals sorting into more granular peer groups such as friends or study partners based on the intrinsic characteristics of the individuals. If the network of interest (for example, the friendship network) is observed by the researcher and can be represented by $G$, the proposed method here can be used to find the causal social effects generated by the network of interest, $G$.
\section{The Model\label{indet}}
This section presents the basic model used to quantify the peer effects in a linear-in-means regression model in \eqref{intro}. As is often the case with observational data, a set of sufficient conditions on the joint data generating process is presented here that guaranties the social parameters are uniquely recovered (identification) from the estimating sample. Let $\mathbf{W}$ be an $n \times n$ stochastic adjacency matrix corresponding to a random network. The matrix $\mathbf{W}$ may or may not be row normalized; the identification results apply to both types of adjacency matrices. The structural model \eqref{intro} in matrix form is given by:
\begin{equation}
\label{E1}
\mathbf{y}=\alpha\boldsymbol{\iota}+\beta\mathbf{Wy}+\mathbf{WX}\boldsymbol{\delta}+\mathbf{X}\boldsymbol{\gamma}+\mathbf{v}
\text{,}
\end{equation}
\noindent where $\alpha$ and $\beta$ are scalar structural parameters, $\boldsymbol{\gamma}$ and $\boldsymbol{\delta}$ are $k \times 1$ vectors of structural parameters, $\mathbf{y}$ is an $n\times1$ vector of outcomes, $\boldsymbol{\iota}$ is a $n\times1$ vector of ones, and $\mathbf{X}$ is an $n\times k$ matrix of exogenous covariates. Under the assumption that $E_{\mathbf{X},\mathbf{W}}[\mathbf{v}]=\mathbf{0}$, regressors $\mathbf{X}$ and $\mathbf{WX}$ are exogenous, but $\mathbf{Wy}$ is not because of simultaneity. Since \eqref{E1} can be thought of as a spatial autoregressive model, \cite{Kelejian1998,Kelejian_Prucha_1999_ER}, \cite{Lee2003}, and \cite{Bramoulle2009}, as well as \cite{BL:2011}, proposed using $[\mathbf{X},\mathbf{WX},\mathbf{W}^{2}\mathbf{X}]$ as the matrix of instruments for the matrix of \emph{endogenous} regressors $[\mathbf{X},\mathbf{WX},\mathbf{Wy}]$ in a Generalized Two-Stage Least Squares (G2SLS) procedure. The use of powers of the adjacency matrix of the form $\mathbf{W}^{2}\mathbf{X}$ as instruments for $\mathbf{W}\mathbf{y}$ is common in the social science literature. The rationale behind the validity of those instruments is that the existence of intransitive triads (contained in $\mathbf{W}^{2}$) guaranties that if nodes $i$ and $j$ are connected, and nodes $j$ and $k$ (but not $i$ and $k$) are connected, then $\mathbf{x}_{k}$ affects $\mathbf{y}_{i}$ but only through its effect on $\mathbf{y}_{j}$. Notice that only one variable of this matrix is endogenous, i.e., $\mathbf{Wy}$, since $E_{\mathbf{X}}[\mathbf{v}]=\mathbf{0}$, and $E_{\mathbf{X},\mathbf{W}}[\mathbf{v}]=\mathbf{0}$ by the tower property of conditional expectations.
This paper assumes instead that $E_{\mathbf{X},\mathbf{W}_{0}}[\mathbf{v}]=\mathbf{0}$, where $\mathbf{W}_{0}$ is an exogenous undirected adjacency matrix that may or may not be row normalized as well. In this way, $\mathbf{W}$ is allowed to be \emph{endogenous}, i.e., $E_{\mathbf{X},\mathbf{W}}[\mathbf{v}]\neq\mathbf{0}$, which invalidates the instruments previously proposed, rendering the G2SLS estimator inconsistent. Notice that there is a set of $k+1$ endogenous regressors in equation \eqref{E1}, i.e., $\mathbf{WX}$ and $\mathbf{Wy}$, given that $E_{\mathbf{X}}[\mathbf{v}]=\mathbf{0}$ by the law of iterated expectations applied to $E_{\mathbf{X},\mathbf{W}_{0}}[\mathbf{v}]=\mathbf{0}$, but $E_{\mathbf{X},\mathbf{W}}[\mathbf{v}]\neq\mathbf{0}$ by the tower property of conditional expectations. Assumptions \ref{A1} and \ref{A2} below state the necessary conditions for $\mathbf{W}_{0}$ to be a valid instrument.
\begin{assumption}
There exists an $n \times n$ adjacency matrix $\mathbf{W}_{0}$ such that $E_{\mathbf{X},\mathbf{W}_{0}}[\mathbf{v}]=\mathbf{0}$.\label{A1}
\end{assumption}
From the original condition, notice that by applying the law of iterated expectations to $E_{\mathbf{X},\mathbf{W}_{0}}[\mathbf{v}]=\mathbf{0}$, the expression $E_{\mathbf{W}_{0}}[\mathbf{v}]=\mathbf{0}$ is obtained. Then, for all $i$ and $E_{\mathbf{W}_{0}}[v_i]=0$, this consequently implies that $E[\mathbf{w}_{0;i}v_i]=0$ by the conditioning theorem, where $\mathbf{w}_{0;i}$ represents the $i$th row of $\mathbf{W}_{0}$. Additionally, Assumption \ref{A1} also implies $E[\mathbf{x}_{i}v_i]=\mathbf{0}$ and $E[\mathbf{W}_{0}\mathbf{x}_{i}v_i]=\mathbf{0}$ by the same logic. The latter implication corresponds to Assumption 1 in \cite{Chan_et_al_social_effects}, and therefore Assumption \ref{A1} above is stronger. Note that equation \eqref{E1}, Assumption \ref{A1}, and the multiplex data structure together imply an exclusion restriction on the adjacency matrix $\mathbf{W}_{0}$. The exogenous adjacency matrix should only affect the outcome $y_{i}$ through its correlation with the adjacency matrix of interest, $\mathbf{W}$. A setup where individuals are (quasi-) randomized into groups is likely to guaranty the validity of the exclusion restriction on $\mathbf{W}_0$. In these settings, individuals are randomly selected into groups that do not necessarily have the potential to create network effects \citep[see, e.g.,][]{carrell2013}. On the other hand, in observational studies, a potential exclusion restriction argument could be based on predetermined networks with respect to the outcome. A network formed in the past is likely to affect outcomes only indirectly through its correlation with contemporaneous relevant networks, as in our real data application.
Consider the regressors formed with the endogenous matrix $\mathbf{W}$, i.e., $\mathbf{W}\mathbf{y}$ and $\mathbf{W}\mathbf{X}$. Let $\mathbf{S}$ be a $n \times (k+1)$ matrix given by $\mathbf{S}\equiv[\mathbf{y} \quad \mathbf{X}]$, and $\boldsymbol{\theta}\equiv(\beta,\boldsymbol{\delta}^{\top})^{\top}$ be a $(k+1) \times 1$ vector of parameters such that $\beta\mathbf{Wy}+\mathbf{WX}\boldsymbol{\delta}=\mathbf{WS}\boldsymbol{\theta}$. Therefore, equation \eqref{E1} can be written as
\begin{equation}
\label{E2}
\mathbf{y}=\alpha\boldsymbol{\iota}+\mathbf{WS}\boldsymbol{\theta}+\mathbf{X}\boldsymbol{\gamma}+\mathbf{v}
\text{,}
\end{equation}
\noindent where, given Assumption \ref{A1}, the endogenous adjacency matrix $\mathbf{W}$, with $i$th row denoted as $\mathbf{w}_{i}$, can be instrumented by $\mathbf{w}_{0;i}$. The above implies that instruments can be constructed by combining the predetermined instrumental matrix $\mathbf{W}_{0}$ and the regressor defined by $\mathbf{S}$ in equation \eqref{E2}, i.e., this is formalized in the following assumption.
\begin{assumption}
\label{A2}
Let $\mathbf{\Pi}$ be the $(k+1)\times(k+1)$ matrix of coefficients from the system regression
\begin{equation}
\label{E3}
\mathbf{WS} =\mathbf{W}_{0}\mathbf{S\Pi}+\mathbf{U}\text{,}
\end{equation}
\noindent where the $n \times(k+1)$ matrix of system errors $\mathbf{U}$ is such that $E_{\mathbf{W}_{0}\mathbf{y}, \mathbf{W}_{0},\mathbf{X}}[\mathbf{U}]=\mathbf{O}$ (a matrix of zeros), $E[\mathbf{S}^{\top}\mathbf{w}_{0;i}\mathbf{w}_{0;i}^{\top}\mathbf{S}]$ is positive definite, and \emph{rank}$(\mathbf{\Pi})=k+1$. Furthermore, the first row of $\mathbf{\Pi}$, $\boldsymbol{\pi}_{1}$, is such that $\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta}< 1/\lambda_{\textup{max}}$, where $\lambda_{\textup{max}}$ is the largest eigenvalue of $\mathbf{W}_{0}$ and $\boldsymbol{\theta}\equiv(\beta,\boldsymbol{\delta}^{\top})^{\top}$.
\end{assumption}
Note that Assumption \ref{A2} implies that $\mathbf{\Pi}$ is uniquely determined by the joint probability distribution of $(\mathbf{S}^{\top}\mathbf{w}_{i},\mathbf{S}^{\top}\mathbf{w}_{0;i})$. The rank condition is necessary for identification. Given that rank$(\mathbf{\Pi})\leq \min \{\text{rank}(E[\mathbf{S}^{\top}\mathbf{w}_{0;i}\mathbf{w}_{0;i}^{\top}\mathbf{S}]^{-1}),\text{rank}(E[\mathbf{S}^{\top}\mathbf{w}_{0;i}\mathbf{w}_{i}^{\top}\mathbf{S}])\}$, a necessary condition for rank$(\mathbf{\Pi})=k+1$ is that $\text{rank}(E[\mathbf{S}^{\top}\mathbf{w}_{0;i}\mathbf{w}_{i}^{\top}\mathbf{S}])=k+1$, which would be equivalent to the \emph{relevance} condition in the classical Instrumental Variable literature. Similarly, as is the case in this literature, equations \eqref{E2} and \eqref{E3} create a natural relationship between $\mathbf{\Pi}$ and structural parameters that is later used for estimation.
As in Section \ref{background}, the condition that $E[\mathbf{S}^{\top}\mathbf{w}_{0;i}\mathbf{w}_{0;i}^{\top}\mathbf{S}]$ is positive definite imposes restrictions on the product matrix $\widetilde{\mathbf{W}}\equiv\mathbf{W}_{0}\mathbf{W}$ for a large enough sample size. Note that if the matrix $\widetilde{\mathbf{W}}=\mathbf{O}$ does not have full rank, this rank condition fails. Thus, identification requires the networks $G$ and $G_{0}$ to be somewhat correlated so that the rank condition imposed by Assumption \ref{A2} holds in the population. This correlation condition effectively requires the existence of edges that overlap in both layers and inter-layer intransitive triads. The second part of Assumption \ref{A2} is likely to hold when $\mathbf{W}_{0}$ is right-stochastic since $0<\lambda_{\textup{max}}<1$ is satisfied in this case. Given the reduced form relation in \eqref{E3}, a reduced form for $\mathbf{y}$ can be constructed by substituting equations \eqref{E3} into \eqref{E2}.
\begin{equation}
\label{E4}
\mathbf{y}=\alpha\boldsymbol{\iota}+\mathbf{W}_{0}\mathbf{S}\mathbf{\Pi}\boldsymbol{\theta}+\mathbf{X}\boldsymbol{\gamma}+\mathbf{e}
\text{,}
\end{equation}
\noindent where $\mathbf{e}\equiv\mathbf{U}\boldsymbol{\theta}+\mathbf{v}$. Note that in \eqref{E4}, $E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{e}]\neq \mathbf{0}$ because of the simultaneity of $\mathbf{W}_{0}\mathbf{y}$ that still persists. To see that, note $E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{e}]=E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{U}\boldsymbol{\theta}]+E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{v}]=E[\mathbf{S}^{\top}\mathbf{W}_{0}E_{\mathbf{S},\mathbf{W}_{0}}[\mathbf{U}]\boldsymbol{\theta}]+E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{v}] =E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{v}]$, where $E[\mathbf{W}_{0}\mathbf{X}^{\top}\mathbf{v}]=\mathbf{0}$ by Assumption \ref{A1}; however, $E[\mathbf{W}_{0}\mathbf{y}^{\top}\mathbf{v}]\neq \mathbf{0}$, which implies $E[\mathbf{S}^{\top}\mathbf{W}_{0}\mathbf{e}]\neq\mathbf{0}$. Therefore, finding the reduced form for $\mathbf{y}$ requires decomposing the matrix $\mathbf{S}$ and the vector $\boldsymbol{\theta}$. Equation \eqref{E4} can be written as
\begin{equation}
\label{E5}
\mathbf{y}[\mathbf{I}-(\boldsymbol{\pi}_{1}^{\top}\boldsymbol{\theta})\mathbf{W}_{0}]=\alpha\boldsymbol{\iota}+[\gamma_1 \mathbf{I}+ (\boldsymbol{\pi}^{\top}_{2}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_1+\dots+[\gamma_k \mathbf{I}+ (\boldsymbol{\pi}^{\top}_{k+1}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_k+\mathbf{e}
\text{,}
\end{equation}
\noindent where $\boldsymbol{\pi}_{j}$ is a $(k+1) \times 1$ vector containing the $j$th row of the coefficient matrix $\mathbf{\Pi}$, such that $\boldsymbol{\pi}^{\top}_{j}\boldsymbol{\theta}$ is a scalar representing the inner product of the two vectors of parameters. The parameter $\gamma_j$ is the $j$th entry of the $k \times 1$ vector $\boldsymbol{\gamma}$, the matrix $\mathbf{I}$ represents the $n \times n$ identity matrix, and $\mathbf{x}_l$ is an $n \times 1$ vector containing the $j$th column of the matrix $\mathbf{X}$. The second part of Assumption \ref{A2} implies that $\mathbf{I}-(\boldsymbol{\pi}_{1}^{\top}\boldsymbol{\theta})\mathbf{W}_{0}$ is invertible. Then, the reduced form for $\mathbf{y} $ is given by:
\begin{equation}
\label{E6}
\mathbf{y}=[\mathbf{I}-(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})\mathbf{W}_{0}]^{-1}\{\alpha\boldsymbol{\iota}+[\gamma_1 \mathbf{I}+ (\boldsymbol{\pi}^{\top}_{2}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_1+\dots+[\gamma_k \mathbf{I}+ (\boldsymbol{\pi}^{\top}_{k+1}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_k+\mathbf{e}\}
\text{,}
\end{equation}
\noindent where now $E[\mathbf{x}_{k}^{\top}\mathbf{e}]=E[\mathbf{x}_{k}^{\top}\mathbf{U}\boldsymbol{\theta}]+E[\mathbf{x}_{k}^{\top}\mathbf{v}]=0$ and $E[\mathbf{x}_{k}^{\top}\mathbf{W}_{0}\mathbf{e}]=E[\mathbf{x}_{k}\mathbf{W}_{0}^{\top}\mathbf{U}\boldsymbol{\theta}]+E[\mathbf{x}_{k}^{\top}\mathbf{W}_{0}^{\top}\mathbf{v}]=0$ for all $k$. The second part of Assumption \ref{A2} also implies that \eqref{E6} can be expressed as
\begin{align}
\label{E7}
\mathbf{y}=&\sum_{r=0}^{\infty}(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})^{r}\mathbf{W}_{0}^{r}\left\{\alpha\boldsymbol{\iota}+[\gamma_1 \mathbf{I}+ (\boldsymbol{\pi}^{\top}_{2}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_1+\dots+[\gamma_k \mathbf{I}+ (\boldsymbol{\pi}^{\top}_{k+1}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_k+\mathbf{e}\right\} \\
\label{E8}
=& [\mathbf{I}-(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})\mathbf{W}_{0}]^{-1}\alpha \boldsymbol{\iota}+\gamma_1 \mathbf{x}_1+[\gamma_1(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})+\boldsymbol{\pi}^{\top}_{2}\boldsymbol{\theta}]\sum_{r=0}^{\infty}(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})^{r}\mathbf{W}_{0}^{r+1}\mathbf{x}_1+\dots+\gamma_k \mathbf{x}_k \nonumber \\
&+[\gamma_k(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})+\boldsymbol{\pi}^{\top}_{k+1}\boldsymbol{\theta}]\sum_{r=0}^{\infty}(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})^{r}\mathbf{W}_{0}^{r+1}\mathbf{x}_k+\sum_{r=0}^{\infty}(\boldsymbol{\pi}^{\top}_{1}\boldsymbol{\theta})^{r}\mathbf{W}_{0}^{r}\mathbf{e}\text{.}
\end{align}
When $k=1$, we have $\boldsymbol{\theta}=(\beta,\delta)^{\top}$, $\boldsymbol{\gamma}=\gamma$, and $\mathbf{\Pi}$ reduces to a $2 \times 2$ matrix. For that particular case, \eqref{E8} reduces to
\begin{align}
\label{E9}
\mathbf{y}=&\alpha/(1-\pi_{1,1}\beta-\pi_{1,2}\delta) \boldsymbol{\iota}+[\gamma(\pi_{1,1}\beta+\pi_{1,2}\delta)+\pi_{2,1}\beta+\pi_{2,2}\delta]\sum_{r=0}^{\infty}(\pi_{1,1}\beta+\pi_{1,2}\delta)^{r}\mathbf{W}_{0}^{r+1}\mathbf{x} \nonumber \\
&+\gamma \mathbf{x} +\sum_{r=0}^{\infty}(\pi_{1,1}\beta+\pi_{1,2}\delta)^{r}\mathbf{W}_{0}^{r}\mathbf{e}\text{.}
\end{align}
\noindent From \eqref{E9}, note that $\gamma (\pi_{1,1}\beta+\pi_{1,2}\delta)+\pi_{2,1}\beta+\pi_{2,2}\delta \neq 0$ is a relevant condition to ensure that $\mathbf{W}_{0}^{2}\mathbf{X}$ is a valid instrument for $\mathbf{W}_{0}\mathbf{y}$ in the case when $k=1$. For the general case of $k$ covariates, the condition is generalized to $\beta(\gamma_k\pi_{1,1}+\pi_{k,1})+\sum_{l=1}^{k}\delta_l(\gamma_l\pi_{1,l+1}+\pi_{k,l+1}) \neq 0$ for all $k$. This is now formalized in the following main identification result:
\begin{theorem}
\label{T1}
Let Assumptions \ref{A1}, \ref{A2} hold and $\beta(\gamma_k\pi_{1,1}+\pi_{k,1})+\sum_{l=1}^{k}\delta_l(\gamma_l\pi_{1,l+1}+\pi_{k,l+1}) \neq 0$ for all $k$. If the matrices $\mathbf{I}$, $\mathbf{W}_{0}$, and $\mathbf{W}_{0}^{2}$ are linearly independent, then the parameters $\alpha$, $\beta$, $\boldsymbol{\gamma}$, and $\boldsymbol{\delta}$ in \eqref{E1} are identified.
\end{theorem}
Assumptions \ref{A1} and \ref{A2} do not require the layers in the multiplex networks to be undirected or unweighted. Thus, the results in Theorem \ref{T1} apply to the general case of potentially directed and weighted networks. Also, Theorem \ref{T1} is based on the existence of a set of regressors for which Assumption \ref{A1} applies, in particular, it requires the regressors in $\mathbf{x}_{i}$ to be exogenous. Note that only one such regressor is necessary to identify the peer effects parameter $\beta$. It is possible to control for other observable characteristics that are not orthogonal to the errors, as long as they are uncorrelated with the set of exogenous regressors.
The condition $\beta(\gamma_k\pi_{1,1}+\pi_{k,1})+\sum_{l=1}^{k}\delta_l(\gamma_l\pi_{1,l+1}+\pi_{k,l+1}) \neq 0$ on the structural parameters has a straightforward interpretation. From equation \eqref{E8} note that this condition involves the composed coefficients associated with the variables $\mathbf{W}_{0}^{r+1}\mathbf{x}_{l}$ for $r \in [0, \infty)$ and $l \in\left\{1,\cdots,k\right\}$. This variable can be interpreted as the reduced form social effects generated by the $l$th observable characteristic. Thus, this restriction on the parameters can be interpreted as a condition guaranteeing that the peer and contextual effects in the network space generated by $\mathbf{W}_{0}$ do not cancel each other out. The linear independence condition guaranties enough exclusion restrictions in the set of simultaneous equations generated by the linear model in \eqref{E4}. This condition can be numerically verified from the observed adjacency matrix $\mathbf{W}_{0}$. Section \ref{corr_effects} in the supplemental materials provides an extension of model \eqref{E1} that weakens the requirement that $\mathbf{W}_0$ is completely predetermined.
\subsection{Network Exogeneity Failure}\label{corr_effects}
In situations where it is difficult to find an adjacency matrix that is completely predetermined and, therefore, not correlated with the unobserved characteristics in the outcome equation, consider the following extension of model \eqref{E1}.
\begin{align}
\label{E25}
&\mathbf{y}=\alpha\boldsymbol{\iota}+\beta\mathbf{Wy}+\mathbf{WX}\boldsymbol{\delta}+\mathbf{X}\boldsymbol{\gamma}+\mathbf{v},\\ \nonumber
&\mathbf{v}=\boldsymbol{\lambda}_{0}+\mathbf{v}^{+}
\text{,}
\end{align}
\noindent where $\boldsymbol{\lambda}_{0}$ represents the correlated effects in the exogenous network $\mathbf{W}_{0}$, the structure of $\boldsymbol{\lambda}_{0}$ is such that it allows for unobserved heterogeneity that is common for all individuals in the same group but varies across individuals in different groups. This specification explicitly accommodates scenarios where group assignment is non-random and correlated with unobservables—such as selection into schools or training programs based on latent ability—a challenge recently highlighted by \cite{Sheng_et_al_2025} in the context of endogenous group formation. The model in \eqref{E25} contains an important implicit assumption, i.e., the relevant characteristics that determine the structure of the exogenous network are common among individuals in the same group.
The relevant assumption for the identification of the transformed model in \eqref{E25} is given by condition $E_{\mathbf{X},\mathbf{W}_{0},\boldsymbol{\lambda}_{0}}[\mathbf{v}^{+}]=\mathbf{0}$. However, $E_{\mathbf{X},\mathbf{W}_{0},\mathbf{W}}[\lambda_0]$ is allowed to be any function of the conditioning random variables. Note that the error $\mathbf{v}^{+}$ is allowed to contain unobserved heterogeneity that is correlated with $\mathbf{W}$ in any form. However, the correlation of $\mathbf{v}^{+}$ with $\mathbf{W}_{0}$ is only allowed through $\boldsymbol{\lambda}_{0}$. Applying group differences to equation \eqref{E25} eliminates the network-specific unobserved heterogeneity for the exogenous network. After applying the transformation, the structural model becomes
\begin{equation}
\label{E26}
(\mathbf{I}-\mathbf{W}_{0})\mathbf{y}=\beta(\mathbf{I}-\mathbf{W}_{0})\mathbf{W}\mathbf{y}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{X}\boldsymbol{\gamma}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{W}\mathbf{X}\boldsymbol{\delta}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{v}^{+}.
\end{equation}
\noindent Analogously to the classic IV estimation in the panel data literature, equation \eqref{E3} is transformed following the same approach as in equation \eqref{E26} to obtain the expression
\begin{equation}
\label{E3s}
(\mathbf{I}-\mathbf{W}_{0})\mathbf{WS} =(\mathbf{I}-\mathbf{W}_{0})\mathbf{W}_{0}\mathbf{S\Pi}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{U}\text{.}
\end{equation}
\noindent Using the notation introduced earlier and considering the transformed model in \eqref{E3s}, the structural equation in \eqref{E26} can be re-written as
\begin{align}
\label{E27}
(\mathbf{I}-\mathbf{W}_{0})\mathbf{y}&=(\mathbf{I}-\mathbf{W}_{0})\mathbf{W}\mathbf{S}\boldsymbol{\theta}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{X}\boldsymbol{\gamma}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{v}^{+} \\
\label{E28}
&=(\mathbf{I}-\mathbf{W}_{0})\mathbf{W}_{0}\mathbf{S}\mathbf{\Pi}\boldsymbol{\theta}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{X}\boldsymbol{\gamma}+(\mathbf{I}-\mathbf{W}_{0})\mathbf{e}^{+},
\end{align}
\noindent where $\mathbf{e}^{+}=\mathbf{U}\boldsymbol{\theta}+\mathbf{v}^{+}$. The next result establishes the identification in this setting through an application of Theorem \ref{T1} to the model described in equation \eqref{E28}.
\begin{proposition}
\label{P1}
Let $\beta(\gamma_k\pi_{1,1}+\pi_{k,1})+\sum_{l=1}^{k}\delta_l(\gamma_l\pi_{1,l+1}+\pi_{k,l+1}) \neq 0$ for all $k$ in model \eqref{E28}. Then the parameters $\alpha$, $\beta$, $\boldsymbol{\gamma}$, and $\boldsymbol{\delta}$ are identified if and only if the matrices $\mathbf{I}$, $\mathbf{W}_{0}$, $\mathbf{W}_{0}^{2}$, and $\mathbf{W}_{0}^{3}$ are linearly independent.
\end{proposition}
\subsection{Identification Failure}
In practice, even in situations where one has access to an exogenous network $\mathbf{W}_{0}$, identification can still fail. Specifically, identification requires $\widetilde{\mathbf{W}}\equiv\mathbf{W}_{0}\mathbf{W}\neq\mathbf{O}$, and this is guaranteed if there is an overlap between the two network layers through the existence of inter-layer intransitive triads. However, if $\mathbf{W}_0$ is extremely sparse, the number of these triads may be insufficient, leading to a problem where the rank of $\mathbf{\Pi}$ in the assumption \ref{A2} is technically full, but $\mathbf{\Pi}$ is numerically close to the zero matrix. Conversely, if $\mathbf{W}_0$ is very dense (e.g., a complete graph or a block-diagonal matrix with identical group sizes), the powers $\mathbf{I}$, $\mathbf{W}_0$, and $\mathbf{W}_0^2$ may become linearly dependent, failing the requirements of Theorem \ref{T1}.
If $\mathbf{W}$ represents close friendship ties in a school, for example, while $\mathbf{W}_0$ represents associations through participation in unrelated, non-overlapping extracurricular activities, the potential lack of shared edges or paths will likely cause $\widetilde{\mathbf{W}}$ to be the zero matrix, violating Assumption \ref{A2}. On the other hand, if individuals are assigned to groups of identical size with no connections between groups, $\mathbf{W}_0$ becomes a block-diagonal matrix where $\mathbf{W}_0^2$ is a linear combination of $\mathbf{I}$ and $\mathbf{W}_0$, precluding identification \citep[see, e.g.,][]{Lee2007}. This highlights a trade-off when applying the transformation proposed in Section \ref{corr_effects}. While highly regular group-based structures naturally motivate the condition $(\mathbf{I}-\mathbf{W}_{0})\boldsymbol{\lambda}_{0}=\mathbf{0}$ by allowing $\mathbf{W}_{0}$ to act as a perfect within-group averaging operator, this extreme regularity simultaneously undermines the linear independence conditions required by Theorem \ref{T1} and Proposition \ref{P1}. Consequently, successful identification in the presence of exogenous correlated effects requires a structural balance: the network must exhibit a group-based architecture to effectively differentiate $\boldsymbol{\lambda}_{0}$, but it must also possess sufficient irregularity—such as variation in group sizes or varied within-group connectivity—to maintain the linear independence of the network powers.
\section{Generalized Three-Stage Least Squares Estimation\label{est}}
Given the point identification in Theorem \ref{T1}, this section describes a multi-step procedure to estimate the parameters of interest in \eqref{E1}. An important feature of this proposed estimator is that it is already implemented in \texttt{Stata} \citep[see, e.g.,][]{netivreg_stata_journal}. First, notice that \eqref{E4} can be re-written as
\begin{equation}
\label{En}
\mathbf{y}=\alpha\boldsymbol{\iota}+\mathbf{W}_{0}\mathbf{S}\boldsymbol{\theta}^{\ast}+\mathbf{X}\boldsymbol{\gamma}+\mathbf{e}
\text{,}
\end{equation}
\noindent where $\boldsymbol{\theta}^{\ast}=\mathbf{\Pi}\boldsymbol{\theta}$. Note that \eqref{En} cannot be directly estimated due to the simultaneity in $\mathbf{W}_{0}\mathbf{y}$. Therefore, the idea of the G2SLS estimator of \cite{Kelejian1998,Kelejian_Prucha_1999_ER} is extended here to propose a Generalized Three-Stage least squares estimator (G3SLS), similar to \cite{BD:2015}, for the structural parameters $(\alpha,\boldsymbol{\gamma}^{\top},\boldsymbol{\theta}^{\top})$ as follows:\\
\noindent\underline{Step One:}
\noindent One starts by estimating the projection coefficient in equation \eqref{E3} by (system) Ordinary Least Squares, i.e.,
\begin{equation}
\label{E17}
\widehat{\mathbf{\Pi}}=(\widehat{\boldsymbol{\pi}}_1,\widehat{\boldsymbol{\pi}}_2,\ldots,\widehat{\boldsymbol{\pi}}_{k+1})^{\top}=(\mathbf{S}^{\top}\mathbf{W}_{0}^{2}\mathbf{S})^{-1}\mathbf{S}^{\top}\widetilde{\mathbf{W}}\mathbf{S}\text{,}
\end{equation}
\noindent where $\widehat{\boldsymbol{\pi}}_1^{\top}=(\mathbf{S}^{\top}\mathbf{W}_{0}^{2}\mathbf{S})^{-1}\mathbf{S}^{\top}\widetilde{\mathbf{W}}\mathbf{y}$, and each $\widehat{\boldsymbol{\pi}}_j^{\top}=(\mathbf{S}^{\top}\mathbf{W}_{0}^{2}\mathbf{S})^{-1}\mathbf{S}^{\top}\widetilde{\mathbf{W}}\mathbf{x}_{j}$ for $j=1, \ldots ,k$. Lemma \ref{l1} in \ref{Appendix_C} states its consistency under the Assumptions listed in \ref{Appendix_B}. In this step, one saves the matrices of estimated coefficients, $\widehat{\mathbf{\Pi}}$, and residuals, $\widehat{\mathbf{U}}=\mathbf{WS}-\mathbf{W}_{0}\mathbf{S}\widehat{\mathbf{\Pi}}$. Matrices $\mathbf{W}_{0}^{2}$ and $\widetilde{\mathbf{W}}\equiv\mathbf{W}_{0}\mathbf{W}$ have a socioeconomic interpretation as explained in Section \ref{background}. The matrix $\mathbf{W}_{0}^{2}$ contains the number of connections that individual $i$ has on the main diagonal, and the number of individuals that separate individuals $i$ and $j$ on the off-diagonal elements. Similarly, $\widetilde{\mathbf{W}}$ contains the number of connections that individual $i$ shares in both networks on the main diagonal, while the off-diagonal elements contain individual $i$'s number of connections from $\mathbf{W}_{0}$ that are connected with individual $j$ in $\mathbf{W}$ for all $i\neq j$. Hence, the matrix $\widetilde{\mathbf{W}}$ represents a measure of association between the endogenous adjacency matrix $\mathbf{W}$ and the exogenous matrix $\mathbf{W}_{0}$.\\
\noindent\underline{Step Two:}
\noindent Rewrite \eqref{En} as
\begin{equation}
\label{En2}
\mathbf{y}=\mathbf{D}_{0}\boldsymbol{\psi}^{\ast}+\mathbf{e},
\end{equation}
\noindent where $\mathbf{D}_{0}=[\mathbf{\iota},\mathbf{X},\mathbf{W}_{0}\mathbf{y},\mathbf{W}_{0}\mathbf{X}]$ is a $n\times (2k+2)$ matrix, $\boldsymbol{\psi}^{\ast}=(\alpha,\boldsymbol{\gamma}^{\top},\boldsymbol{\theta}^{\ast\top})^{\top}$ is a $(2k+2) \times 1$ vector of parameters, and $\mathbf{e}=\mathbf{U}\boldsymbol{\theta}+\mathbf{v}$. From equation \eqref{En2}, the set of parameters $\boldsymbol{\psi}^{\ast}$ can be estimated using 2SLS. Let $\mathbf{Z}=[\mathbf{\iota},\mathbf{X},\mathbf{W}^{2}_{0}\mathbf{X},\mathbf{W}_{0}\mathbf{X}]$ be the matrix of instruments for the variables in $\mathbf{D}_{0}$. For the general case where $k>1$, the model is overidentified and the standard two-stage least squares (2SLS) estimator is given by:
\begin{equation}
\label{E19}
\widehat{\boldsymbol{\psi}}^{\ast}_{\text{2SLS}}=(\mathbf{D}_{0}^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\mathbf{D}_{0})^{-1} \mathbf{D}_{0}^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\mathbf{y}.
\end{equation}
\noindent Under Assumptions in \ref{Appendix_B}, Lemma \ref{l2} in \ref{Appendix_C} establishes its consistency. Taking into account the estimator $\widehat{\boldsymbol{\psi}}^{\ast}_{\text{2SLS}}=(\widehat{\alpha}_{\text{2SLS}},\widehat{\boldsymbol{\gamma}}_{\text{2SLS}}^{\top},\widehat{\boldsymbol{\theta}}^{\ast\top}_{\text{2SLS}})^{\top}$, an estimator for the structural parameters of the social effects can be recovered from the relation $\boldsymbol{\theta}^{\ast}=\mathbf{\Pi}\boldsymbol{\theta}$, i.e., using $\widehat{\boldsymbol{\Pi}}$ from the first step, a consistent estimator for the social effects in equation \eqref{E1} is given by $\widehat{\boldsymbol{\theta}}=\widehat{\mathbf{\Pi}}^{-1}\widehat{\boldsymbol{\theta}}^{\ast}_{\text{2SLS}}$ -- Lemma \ref{l3} in \ref{Appendix_C} proves its consistency under the assumptions listed in \ref{Appendix_B} and its performance in our Monte Carlo designs below is reported in \ref{Appendix_D} of the supplement. In this step, the vector of the estimated coefficients is stored $\widehat{\alpha}_{\text{2SLS}}$, $\widehat{\boldsymbol{\gamma}}_{\text{2SLS}}$, and $\widehat{\boldsymbol{\theta}}$.\\
\noindent\underline{Step Three:}
\noindent This last step consists of constructing the optimal instrument for the regressor $\mathbf{W}_{0}\mathbf{y}$ in \eqref{En}, i.e., the optimal instrument (in the classical sense) is given by $E_{\mathbf{X},\mathbf{W}_0}[\mathbf{W}_0\mathbf{y}]$. Note that from the reduced-form equation in \eqref{E5} one has:
\begin{align}
\label{E20}
&E_{\mathbf{X},\mathbf{W}_0}[\mathbf{W}_0\mathbf{y}](\boldsymbol{\psi},\mathbf{\Pi})\\
&=\mathbf{W}_0[\mathbf{I}-(\boldsymbol{\pi}_{1}^{\top}\boldsymbol{\theta})\mathbf{W}_{0}]^{-1}\left\{\alpha\boldsymbol{\iota}+[\gamma_1 \mathbf{I}+ (\boldsymbol{\pi}_{2}^{\top}\boldsymbol{\theta})\mathbf{W}_{0}] \mathbf{x}_1+\dots+[\gamma_k \mathbf{I}+ (\boldsymbol{\pi}_{k+1}^{\top}\boldsymbol{\theta})\mathbf{W}_{0}]\mathbf{x}_{k}\right\}.\nonumber
\end{align}
\noindent A valid estimator of \eqref{E20} is then given by $E_{\mathbf{X},\mathbf{W}_0}(\widehat{\boldsymbol{\psi}},\widehat{\mathbf{\Pi}})$, where $\widehat{\boldsymbol{\psi}}=(\widehat{\alpha}_{\text{2SLS}},\widehat{\boldsymbol{\gamma}}_{\text{2SLS}}^{\top},\widehat{\boldsymbol{\theta}}^{\top})^{\top}$ from the second step. Now rewrite \eqref{E1} as
\begin{equation}
\label{En3}
\mathbf{y}=\mathbf{D}\boldsymbol{\psi}+\mathbf{v}\text{,}
\end{equation}
\noindent where $\mathbf{D} = [\boldsymbol{\iota},\mathbf{X},\mathbf{W}\mathbf{S}]$ and $\boldsymbol{\psi}=(\alpha,\boldsymbol{\gamma}^{\top},\boldsymbol{\theta}^{\top})^{\top}$. Using $\widehat{\boldsymbol{\Pi}}$ from the first step, let the matrix $\widehat{\mathbf{D}}$ be defined as $\widehat{\mathbf{D}}=[\boldsymbol{\iota},\mathbf{X},\mathbf{W}_{0}\mathbf{S}\widehat{\mathbf{\Pi}}]=[\boldsymbol{\iota},\mathbf{X},\mathbf{W}_{0}\mathbf{S}]\widehat{\mathbf{\Gamma}}=[\mathbf{\iota},\mathbf{X},\mathbf{W}_{0}\mathbf{y},\mathbf{W}_{0}\mathbf{X}]\widehat{\mathbf{\Gamma}}=\mathbf{D}_{0}\widehat{\mathbf{\Gamma}}$, where
\begin{equation}
\label{E18}
\widehat{\mathbf{\Gamma}}= \begin{bmatrix}
\mathbf{I}_{k+1} & \mathbf{O}_{k+1} \\
\mathbf{O}_{k+1} & \widehat{\mathbf{\Pi}} \\
\end{bmatrix},
\end{equation}
\noindent is a $(2k+2) \times (2k+2)$ matrix, $\mathbf{O}_{k+1}$ is a $(k+1)\times(k+1)$ matrix of zeros, and $\mathbf{I}_{k+1}$ represents the identity matrix of order $k+1$. Let $\widetilde{\mathbf{Z}}^{\ast}$ be the matrix of optimal instruments for $\mathbf{D}_{0}\mathbf{\Gamma}$, where $\widetilde{\mathbf{Z}}^{\ast}=\mathbf{Z}^{\ast}\mathbf{\Gamma}$, and $\mathbf{Z}^{\ast}=[\boldsymbol{\iota},\mathbf{X},E_{\mathbf{X},\mathbf{W}_0}[\mathbf{W}_0\mathbf{y}](\boldsymbol{\psi},\mathbf{\Pi}),\mathbf{W}_{0}\mathbf{X}]$. Analogous to the previous transformation, by additionally considering the estimator of the conditional expectation of $\mathbf{W}_{0}\mathbf{y}$ in the second step, an estimator of the optimal matrix of instruments is then given by
\begin{equation}
\label{EW0y}
\widehat{\widetilde{\mathbf{Z}}}^{\ast}=\widehat{\mathbf{Z}}^{\ast}\widehat{\mathbf{\Gamma}}.
\end{equation}
\noindent where $\widehat{\mathbf{Z}}^{\ast}=[\boldsymbol{\iota},\mathbf{X},E_{\mathbf{X},\mathbf{W}_0}[\mathbf{W}_0\mathbf{y}](\widehat{\boldsymbol{\psi}},\widehat{\mathbf{\Pi}}),\mathbf{W}_{0}\mathbf{X}]$. Therefore, the G3SLS in this just-identified case is
\begin{equation}
\label{E21}
\widehat{\boldsymbol{\psi}}_{\text{G3SLS}}=(\widehat{\widetilde{\mathbf{Z}}}{}^{\ast\top}\widehat{\mathbf{D}})^{-1}\widehat{\widetilde{\mathbf{Z}}}{}^{\ast\top}\mathbf{y}.
\end{equation}
\noindent Finally, from this step, one saves the resulting vector of estimated coefficients $\widehat{\boldsymbol{\psi}}_{\text{G3SLS}}$ and the residuals $\widehat{\mathbf{v}}=\mathbf{y}-\mathbf{D}\widehat{\boldsymbol{\psi}}_{\text{G3SLS}}$.\\
After defining $\mathbf{M}_{\mathbf{W}_0}\equiv\mathbf{I}_{n}-\mathbf{W_{0}}\mathbf{S}(\mathbf{S}^{\top}\mathbf{W}_{0}^{2}\mathbf{S})^{-1}\mathbf{S}^{\top}\mathbf{W_{0}}$, the following theorem establishes the asymptotic normality of the proposed G3SLS estimator.
\begin{theorem}
\label{dist}
Let Assumptions \ref{A1}, \ref{A2}, and \ref{D1}--\ref{D8} in \ref{Appendix_B} hold, then $\widehat{\boldsymbol{\psi}}_{\text{\emph{G3SLS}}}=\boldsymbol{\psi}+o_p(1)$ and $n^{1/2}(\widehat{\boldsymbol{\psi}}_{\text{\emph{G3SLS}}}-\boldsymbol{\psi})\overset{d}{\longrightarrow}N(\mathbf{0},\mathbf{V}_{\boldsymbol{\psi}})$, where
\[
\mathbf{V}_{\boldsymbol{\psi}}=(\mathbf{\Gamma}^{\top}\mathbf{Q}_{\mathbf{Z}^{\ast}\mathbf{D}_{0}}\mathbf{\Gamma})^{-1}\boldsymbol{\Omega}(\mathbf{\Gamma}\mathbf{Q}^{\top}_{\mathbf{Z}^{\ast}\mathbf{D}_{0}}\mathbf{\Gamma}^{\top})^{-1},
\]
$\mathbf{Q}_{\mathbf{Z}^{\ast}\mathbf{D}_{0}}=\lim_{n\to\infty}n^{-1}\mathbf{Z}^{\ast\top}\mathbf{D}_{0}$, $\boldsymbol{\Omega}=\lim_{n\to\infty}n^{-1}\sum_{i=1}^{n}\widetilde{\mathbf{z}}_{i}^{\ast}\widetilde{\mathbf{z}}_{i}^{\ast\top}E_{\widetilde{\mathbf{Z}^{\ast}}}[e_{i}^{\ast2}]$, $\widetilde{\mathbf{z}}_{i}^{\ast}$ represents the $i$th row of $\widetilde{\mathbf{Z}}^{\ast}$, and $e_{i}^{\ast}$ denotes the $i$th element of the $n\times 1$ vector $\mathbf{e}^{\ast}=\mathbf{M}_{\mathbf{W}_0}\mathbf{U}\boldsymbol{\theta}+\mathbf{v}$.
\end{theorem}
\noindent The asymptotic variance-covariance matrix, $\mathbf{V}_{\boldsymbol{\psi}}$, takes into account the estimation effects of the various steps, i.e., the error $\mathbf{e}^{\ast}$ associated with the estimator for $\boldsymbol{\psi}$ is composed of structural errors $\mathbf{v}$ and also errors in the first step $\mathbf{U}$. Note that $\mathbf{V}_{\boldsymbol{\psi}}$ is a $(2k+2) \times (2k+2)$ matrix containing the variances and covariances of all the structural coefficients in equation \eqref{E1}. To obtain the variance-covariance matrix of the parameters of the social effects $\boldsymbol{\theta}$, we need to extract the appropriate sub-component $(k+1) \times (k+1)$ from $\mathbf{V_{\boldsymbol{\psi}}}$. Let $\mathbf{J}$ be a $(k+1) \times (2k+2)$ matrix such that $\mathbf{J}=[\mathbf{O}_{k+1},\mathbf{I}_{k+1}]$. Then, the variance-covariance matrix of the social parameters of interest is given by
\begin{equation}
\mathbf{V}_{\boldsymbol{\theta}}=\mathbf{J}\mathbf{V_{\boldsymbol{\psi}}}\mathbf{J}^{\top}.\label{V_theta}
\end{equation}
\subsubsection*{Standard Errors Calculation}
\noindent A consistent estimator for $\mathbf{V_{\boldsymbol{\psi}}}$ can be calculated as
\begin{equation}
\label{E23}
\widehat{\mathbf{V}}_{\boldsymbol{\psi}}=(n^{-1}\widehat{\widetilde{\mathbf{Z}}}{}^{\ast\top}\widehat{\mathbf{D}})^{-1} (n^{-1}\sum\nolimits_{i=1}^{n}\widehat{\widetilde{\mathbf{z}}}{}^{\ast}_{i}\widehat{\widetilde{\mathbf{z}}}{}^{\ast\top}_{i}\widehat{e}_{i}^{\ast 2}) (n^{-1}\widehat{\mathbf{D}}^{\top}\widehat{\widetilde{\mathbf{Z}}}{}^{\ast})^{-1},
\end{equation}
\noindent where $\widehat{\mathbf{e}}^{\ast}=\mathbf{M}_{\mathbf{W}_0}\widehat{\mathbf{U}}\widehat{\boldsymbol{\theta}}+\widehat{\mathbf{v}}$. The residuals $\widehat{\mathbf{U}}$ are obtained from the first step, $\widehat{\boldsymbol{\theta}}$, and the residuals $\widehat{\mathbf{v}}$ are taken from the last step above. Therefore, a consistent estimator of \eqref{V_theta} is simply given by $\widehat{\mathbf{V}}_{\boldsymbol{\theta}}=\mathbf{J}\widehat{\mathbf{V}}_{\boldsymbol{\psi}}\mathbf{J}^{\top}$. The standard errors for the coefficients of interest are then the squared root of the main diagonal elements of this matrix after dividing them by $n$.
Note that, as with standard instrumental variable methods, the finite-sample moments of the proposed G3SLS estimator here may not exist under certain conditions. As established in the simultaneous equations literature \citep[see, e.g.,][]{mariano2001simultaneous}, the existence of their finite moments typically depends on the degree of over-identification, specifically, the difference between the number of valid instruments used and the number of endogenous regressors. When this degree of over-identification is low, the estimator may lack finite moments of higher orders.
\section{Monte Carlo Experiments\label{mc}}
In this section, we showcase the proposed estimator's strong performance and versatility across three distinct data-generating processes (DGPs). The endogeneity in these DGPs is generated by a classic omitted variable (Design 1), measurement error in the connections (Design 2), and an unobserved homophily with the simultaneous determination of network formation and outcomes (Design 3). Our Design 3 corresponds to the unobserved characteristics with homophily scenario (Design 1) presented in \cite{Chan_et_al_social_effects}. Similarly, our Design 2 is a modified version of the misclassified links scenario (Design 2) in that same paper, originally from \cite{Lewbel_Qu_Tang}. We include a detailed description of these scenarios here for completeness and to highlight their relevance in social science research. A total of 1,500 data sets $\left\{y_i,x_i,\{w_{i,j}\}_{j=1,j\neq i}^n,\{w_{0;i,j}\}_{j=1,j\neq i}^n\right\}_{i=1}^n$ with $n\in\left\{50,100,200\right\}$ are generated from \eqref{E1} by setting $k=1$, $\beta=0.7$, $\alpha=\delta=\gamma=1$, and drawing $\left\{x_i\right\}_{i=1}^n$ as a random sample from a normal distribution with mean zero and variance 3. The other data components are constructed as follows:
\subsubsection*{Design 1: Unobserved Heterogeneity\label{d1}}
As in \cite{Johnsson2019}, the network formation process follows \cite{Graham2017}, i.e., links are formed according to the rule $w_{i,j}(\psi)\equiv\mathbb{I}[z_{i}z_{j} + \psi(a_{i}+a_{j}) - u_{i,j} \geq 0]$, where $\mathbb{I}(\cdot)$ denotes the indicator function that equals one if its argument is true and zero otherwise. The scalar parameter $\psi$ acts as a switch to activate or deactivate the individual degree heterogeneity in the network formation process, so that the endogenous adjacency matrix is $\mathbf{W} = [w_{i,j}(1)]$, while the exogenous matrix is given by $\mathbf{W}_{0} = [w_{i,j}(0)]$. The exogenous variable $z$ takes on values -1 and 1 with a probability of 0.5 (these values imply a strong taste for homophilic matching; see \citealt{Graham2017}), $u_{ij}$ is drawn from a logistic distribution with mean zero and scale parameter 1, and $a_{i}=\alpha_{L} \mathbb{I}[z_{i}=-1]+\alpha_{H} \mathbb{I}[z_{i}=1]+\xi_{i}$, where $\alpha_{L}=-3/2$ and $\alpha_{H}=1$ are both parameters controlling the extent to which the degree heterogeneity $a_i$ is correlated with the observable exogenous characteristic $z_i$. The error $\xi_i$ is drawn from a re-centered Beta($1/4$,$3/4$) distribution (implying the often found right-skewed degree distribution in empirical applications, \citealp[see, e.g., ][]{Johnsson2019}). On the other hand, individual outcomes are constructed as in \eqref{E1}, where $\mathbf{v}=m\times \boldsymbol{a}+\boldsymbol{\varepsilon}_{0}$ and $m\in\{10,12\}$ measure the importance of degree heterogeneity; i.e., the higher it is, the higher the level of endogeneity will be. The error vector $\boldsymbol{\varepsilon}_{0}$ is drawn from a multivariate standard normal distribution independently of everything else.
\subsubsection*{Design 2: Misclassified Links\label{d2}}
This is an adapted iteration of the Monte Carlo design originally presented by \cite{Lewbel_Qu_Tang}. As in \cite{Chan_et_al_social_effects}, the true DGP includes a hidden adjacency matrix $\mathbf{W}_{0}^{\ast}=[w_{0;i,j}^{\ast}]$, which is derived from a typical random network model by \citeauthor{Erdos1959}'s \citeyearpar{Erdos1959}, featuring a density of 0.05 across a network of size $n$. However, it is assumed that the researcher only has access to an adjacency matrix $\mathbf{W}=[w_{i,j}]$ where links are randomly misclassified. This misclassification is represented as $w_{i, j}=w_{0;i,j}^{\ast}e_{1;i,j}+(1-w_{0;i,j}^{\ast}) e_{2;i,j}$ for $i \neq j$. Additionally, there is an exogenous adjacency matrix $\mathbf{W}_{0}=[w_{0;i,j}]$ with $w_{0;i,j}=w_{0;i,j}^{*}b_{1;i,j}+(1-w_{0;i,j}^{*}) b_{2;i,j}$ for $i \neq j$. The variables $e_{1;i,j}$, $e_{2;i,j}$, $b_{1;i,j}$, and $b_{2;i,j}$ are independently drawn Bernoulli random variables from one another for all $i \neq j$, with probabilities set at 0.75, 0, $1-\tau$, and 0.002, respectively. The design parameter $\tau\in\left\{0.01,0.05\right\}$ dictates the likelihood of misclassification in $\mathbf{W}_0$. It is crucial to note, consistent with \cite{Lewbel_Qu_Tang}, that non-existent links in $\mathbf{W}$ are never misidentified, whereas such misclassification is allowed in $\mathbf{W}_0$ but at a meager chance of 0.2\%. Nevertheless, this setup directly links the vector of individual outcomes to the proportion of misclassification in $\mathbf{W}$ for each $i$. The $n\times 1$ outcome vector $\mathbf{y}$ is constructed according to equation \eqref{E1}, where $\mathbf{v}=\boldsymbol{\varepsilon}_{1}+\boldsymbol{\varepsilon}_{2}$, $\varepsilon_{1,i}=\frac{1}{n}\sum_{j=1}^{n}w_{0;i,j}^{*}e_{1;i,j}$, and $\boldsymbol{\varepsilon}_{2}$ originate from a multivariate standard normal distribution that is independent of all other elements.
\subsubsection*{Design 3: Unobserved Characteristics with Homophily\label{d3}}
Within this framework, the outcome variable for individual $i$, denoted as $y_i$, alongside connections $\left\{w_{i,j}\right\}_{j=1,j\neq i}^n$, are simultaneously determined by a shared idiosyncratic unobserved feature associated with homophily, $\varepsilon_{3;i}^{\ast}$. Initially, an exogenous adjacency matrix $\mathbf{W}_{0}=[w_{0;i,j}]$ is constructed from a \citeauthor{Erdos1959}'s \citeyearpar{Erdos1959} random graph featuring a density of 0.01, along with an $n\times 1$ vector $\boldsymbol{\varepsilon}_{3}^{\ast}=[\varepsilon_{3;1}^{\ast},\ldots,\varepsilon_{3;n}^{\ast}]^\top$ drawn from a multivariate standard normal distribution. Subsequently, the elements of the endogenous adjacency matrix $\mathbf{W}=[w_{i,j}]$ are computed as
\begin{equation*}
w_{ij}=
\begin{cases}
\mathbb{I}[|\varepsilon_{3;i}^{\ast}-\varepsilon_{3;j}^{\ast}|<\widehat{F}_{\varepsilon_3^\ast}^{-1}(0.99)]\times (1-w_{0;i,j}) + w_{0;i,j} & \text{; if $\varepsilon_{3;i}^{\ast}>\Phi^{-1}(0.99)$,} \\
\mathbb{I}[|\varepsilon_{3;i}^{\ast}-\varepsilon_{3;j}^{\ast}|<\widehat{F}_{\varepsilon_3^\ast}^{-1}(0.99)] \times w_{0;i,j} & \text{; if $\varepsilon_{3;i}^{\ast}<\Phi^{-1}(0.01)$,} \\
w_{0;i,j} & \text{; otherwise},
\end{cases}
\end{equation*}
\noindent where $\widehat{F}_{\varepsilon_3^{\ast}}^{-1}(0.99)$ denotes the 99\% empirical quantile of the components of the vector $\boldsymbol{\varepsilon}_3^{\ast}$, where $\varepsilon_{3;k}^{\ast}$ signifies its $k$th component, and $\Phi^{-1}(\cdot)$ is the inverse operation of the cumulative distribution function for a standard normal variable. This framework encapsulates the concept of homophily, suggesting that agents with higher values of $\varepsilon_{3}$ are inclined to forge or sustain relationships with others who also possess high $\varepsilon_{3}$ values, while disconnecting from those with lower values of this unique, unobservable trait. The $n\times 1$ outcome vector, $\boldsymbol{y}$, is derived from \eqref{E1} by defining $\mathbf{v}=m \times \boldsymbol{\varepsilon}_{3}+\boldsymbol{\varepsilon_{4}}$, where $m$ is selected from $\left\{1,3\right\}$, $\boldsymbol{\varepsilon_{4}}$ is sampled from a multivariate standard normal distribution, and the elements of $\boldsymbol{\varepsilon}_{3}$ are specified as
\begin{equation*}
\varepsilon_{3;i}=
\begin{cases}
\varepsilon_{3;i}^{\ast} & \text{; if $\varepsilon_{3;i}^{\ast}<\Phi^{-1}(0.01)$ or $\varepsilon_{3;i}^{\ast}>\Phi^{-1}(0.99)$,} \\
0 & \text{; otherwise}.
\end{cases}
\end{equation*}
\subsubsection*{Results\label{MC_results}}
\noindent Tables \ref{tab:design1_heterogeneity_main}--\ref{tab:design3_homophily_main} report the Monte Carlo results for the three estimators across the alternative network designs. In addition to the proposed G3SLS estimator in \eqref{E21}, we evaluate the standard Ordinary Least Squares (OLS) estimator -- which ignores the network endogeneity -- and the Generalized Two-Stage Least Squares (G2SLS) estimator -- which sets $\mathbf{W}=\mathbf{W}_0$ in our proposed estimator. All adjacency matrices are row-normalized prior to estimation, following \cite{liu2014}.
Tables \ref{tab:design1_heterogeneity_main}--\ref{tab:design3_homophily_main} summarize the finite-sample performance of the estimators under each design. For the peer ($\beta$), contextual ($\delta$), and direct ($\gamma$) effects in \eqref{E1}, we report bias, standard deviation (SD), root mean squared error (RMSE), and inter-quartile range (IQR) across alternative values of the design parameters $m$ or $\tau$ and different sample sizes. Several regularities emerge from the Monte Carlo evidence when focusing on bias and RMSE. Across the three structural parameters — peer, contextual, and direct effects — the proposed G3SLS estimator systematically delivers lower bias relative to OLS and G2SLS, particularly as the sample size increases. The improvements in bias directly translate into lower RMSE. For all three parameters, the RMSE of G3SLS declines steadily with larger samples and remains uniformly below that of OLS and G2SLS across designs. By contrast, the higher bias observed under naive OLS and conventional G2SLS contributes to persistently larger RMSE values, a pattern consistent with the link misclassification mechanisms discussed in \cite{Chandrasekhar_unpub_2016}. Overall, the Monte Carlo evidence indicates that G3SLS has larger Monte Carlo variance than both OLS and G2SLS in DGPs 1 across all sample sizes and parameters, and for the peer effect parameter in DGPs 2 and 3 across all sample sizes. This might be attributed to the potential lack of finite higher moments in these designs, rather than the lack of consistency; i.e., the G3SLS achieves a more accurate estimation of all social interaction parameters, as reflected in both the lower bias and the systematically smaller RMSE.
Finally, to assess the practical value of the third estimation step in Section \ref{est}, we compare the finite-sample properties of the efficient G3SLS estimator against the intermediate 2SLS estimator defined in equation \eqref{E19} in Step 2. Detailed results are provided in Appendix \ref{2SLS_MC} of the Supplementary Material (Tables \ref{tab:design1_heterogeneity_2sls}-\ref{tab:design3_homophily_2sls}). The analysis confirms that while the 2SLS estimator is consistent, the G3SLS estimator yields substantial efficiency gains, particularly in designs with significant unobserved heterogeneity, thereby justifying the additional estimation step.
\begin{landscape}
\begin{table}[!htbp]
\centering
\small
\begin{threeparttable}
\caption{Estimator Performance under Unobserved Degree Heterogeneity (Design 1)}
\label{tab:design1_heterogeneity_main}
\begin{tabular}{ccc c cccc c cccc c cccc}
\toprule
& & & &\multicolumn{4}{c}{Peer effects} & &\multicolumn{4}{c}{Contextual effects} & &\multicolumn{4}{c}{Direct effects} \\
\cmidrule(lr){5-8} \cmidrule(lr){10-13} \cmidrule(lr){15-18}
$m$ & Estimator & $n$ & & Bias & SD & RMSE & IQR & & Bias & SD & RMSE & IQR & & Bias & SD & RMSE & IQR \\
\midrule
\multirow{9}{*}{10} & \multirow{3}{*}{OLS} & 50 & & 0.359 & 0.046 & 0.361 & 0.072 & & -1.119 & 0.875 & 1.420 & 1.357 & & -0.449 & 0.616 & 0.762 & 0.983 \\
& & 100 & & 0.352 & 0.031 & 0.353 & 0.049 & & -1.104 & 0.576 & 1.245 & 0.897 & & -0.441 & 0.417 & 0.607 & 0.643 \\
& & 200 & & 0.348 & 0.022 & 0.349 & 0.033 & & -1.110 & 0.404 & 1.181 & 0.633 & & -0.472 & 0.287 & 0.552 & 0.442 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G2SLS} & 50 & & 0.691 & 0.317 & 0.760 & 0.306 & & -1.612 & 2.233 & 2.754 & 2.917 & & -0.588 & 1.087 & 1.235 & 1.502 \\
& & 100 & & 0.742 & 0.268 & 0.789 & 0.278 & & -2.143 & 1.786 & 2.789 & 2.500 & & -0.816 & 0.811 & 1.150 & 1.158 \\
& & 200 & & 0.810 & 0.272 & 0.854 & 0.289 & & -2.532 & 1.551 & 2.969 & 2.203 & & -1.034 & 0.714 & 1.256 & 0.957 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G3SLS} & 50 & & 0.074 & 0.986 & 0.989 & 1.283 & & 0.250 & 6.687 & 6.689 & 8.261 & & 0.500 & 1.594 & 1.670 & 2.372 \\
& & 100 & & 0.027 & 0.773 & 0.773 & 1.032 & & 0.214 & 4.136 & 4.140 & 5.281 & & 0.562 & 1.079 & 1.217 & 1.629 \\
& & 200 & & -0.005 & 0.549 & 0.549 & 0.791 & & 0.245 & 2.486 & 2.497 & 3.500 & & 0.562 & 0.718 & 0.912 & 1.124 \\
\midrule
\multirow{9}{*}{12} & \multirow{3}{*}{OLS} & 50 & & 0.361 & 0.047 & 0.364 & 0.074 & & -1.132 & 1.043 & 1.539 & 1.628 & & -0.454 & 0.736 & 0.865 & 1.163 \\
& & 100 & & 0.354 & 0.031 & 0.355 & 0.050 & & -1.110 & 0.687 & 1.305 & 1.078 & & -0.441 & 0.501 & 0.667 & 0.776 \\
& & 200 & & 0.350 & 0.022 & 0.351 & 0.034 & & -1.118 & 0.481 & 1.217 & 0.749 & & -0.476 & 0.344 & 0.588 & 0.533 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G2SLS} & 50 & & 0.732 & 0.304 & 0.793 & 0.305 & & -1.770 & 2.619 & 3.160 & 3.341 & & -0.627 & 1.266 & 1.413 & 1.758 \\
& & 100 & & 0.795 & 0.236 & 0.829 & 0.266 & & -2.199 & 1.952 & 2.940 & 2.632 & & -0.840 & 0.914 & 1.241 & 1.371 \\
& & 200 & & 0.827 & 0.204 & 0.852 & 0.240 & & -2.543 & 1.541 & 2.974 & 2.249 & & -1.044 & 0.732 & 1.275 & 0.997 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G3SLS} & 50 & & 0.114 & 1.004 & 1.011 & 1.295 & & -0.051 & 7.860 & 7.857 & 9.720 & & 0.473 & 1.902 & 1.959 & 2.851 \\
& & 100 & & 0.115 & 0.824 & 0.831 & 1.088 & & -0.069 & 4.997 & 4.995 & 6.546 & & 0.523 & 1.267 & 1.371 & 1.941 \\
& & 200 & & 0.015 & 0.608 & 0.608 & 0.832 & & 0.148 & 2.927 & 2.930 & 4.032 & & 0.550 & 0.857 & 1.018 & 1.328 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
\footnotesize
\setstretch{1}
\item \textit{Notes:} This table reports Monte Carlo results under unobserved degree heterogeneity. The parameter $m$ controls the intensity of heterogeneity in node degrees and is common across individuals. Reported statistics include bias, standard deviation (SD), root mean squared error (RMSE), and inter-quantile range (IQR) for the peer effect ($\beta$=0.7), contextual ($\delta=1$), and direct ($\gamma=1$) Effects.
\end{tablenotes}
\end{threeparttable}
\normalsize
\end{table}
\end{landscape}
\begin{landscape}
\begin{table}[!htbp]
\centering
\small
\begin{threeparttable}
\caption{Estimator Robustness to Network Link Misclassification (Design 2)}
\label{tab:design2_misclassification_main}
\begin{tabular}{ccc c cccc c cccc c cccc}
\toprule
& & & &\multicolumn{4}{c}{Peer effects} & &\multicolumn{4}{c}{Contextual effects} & &\multicolumn{4}{c}{Direct effects} \\
\cmidrule(lr){5-8} \cmidrule(lr){10-13} \cmidrule(lr){15-18}
$\tau$ & Estimator & $n$ & & Bias & SD & RMSE & IQR & & Bias & SD & RMSE & IQR & & Bias & SD & RMSE & IQR \\
\midrule
\multirow{9}{*}{0.01} & \multirow{3}{*}{OLS} & 50 & & -0.387 & 0.129 & 0.408 & 0.201 & & -0.425 & 0.314 & 0.528 & 0.497 & & 0.411 & 0.220 & 0.466 & 0.346 \\
& & 100 & & -0.525 & 0.087 & 0.533 & 0.137 & & -0.628 & 0.153 & 0.646 & 0.244 & & 0.203 & 0.098 & 0.225 & 0.154 \\
& & 200 & & -0.555 & 0.062 & 0.559 & 0.098 & & -0.778 & 0.091 & 0.783 & 0.149 & & 0.092 & 0.047 & 0.103 & 0.073 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G2SLS} & 50 & & -0.538 & 0.222 & 0.582 & 0.305 & & -0.152 & 0.469 & 0.493 & 0.718 & & 0.475 & 0.251 & 0.538 & 0.387 \\
& & 100 & & -0.649 & 0.133 & 0.663 & 0.184 & & -0.457 & 0.215 & 0.505 & 0.326 & & 0.237 & 0.112 & 0.262 & 0.179 \\
& & 200 & & -0.688 & 0.092 & 0.695 & 0.140 & & -0.623 & 0.122 & 0.635 & 0.191 & & 0.114 & 0.051 & 0.125 & 0.080 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G3SLS} & 50 & & 0.739 & 0.466 & 0.873 & 0.668 & & 2.406 & 1.961 & 3.103 & 2.783 & & 0.051 & 0.092 & 0.105 & 0.138 \\
& & 100 & & 0.249 & 0.188 & 0.312 & 0.294 & & 0.995 & 0.648 & 1.187 & 0.976 & & 0.023 & 0.059 & 0.063 & 0.092 \\
& & 200 & & 0.043 & 0.175 & 0.180 & 0.274 & & 0.235 & 0.337 & 0.411 & 0.535 & & 0.007 & 0.038 & 0.039 & 0.059 \\
\midrule
\multirow{9}{*}{0.05} & \multirow{3}{*}{OLS} & 50 & & -0.387 & 0.129 & 0.408 & 0.201 & & -0.425 & 0.314 & 0.528 & 0.497 & & 0.411 & 0.220 & 0.466 & 0.346 \\
& & 100 & & -0.525 & 0.087 & 0.533 & 0.137 & & -0.628 & 0.153 & 0.646 & 0.244 & & 0.203 & 0.098 & 0.225 & 0.154 \\
& & 200 & & -0.555 & 0.062 & 0.559 & 0.098 & & -0.778 & 0.091 & 0.783 & 0.149 & & 0.092 & 0.047 & 0.103 & 0.073 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G2SLS} & 50 & & -0.538 & 0.222 & 0.582 & 0.305 & & -0.152 & 0.469 & 0.493 & 0.718 & & 0.475 & 0.251 & 0.538 & 0.387 \\
& & 100 & & -0.649 & 0.133 & 0.663 & 0.184 & & -0.457 & 0.215 & 0.505 & 0.326 & & 0.237 & 0.112 & 0.262 & 0.179 \\
& & 200 & & -0.688 & 0.092 & 0.695 & 0.140 & & -0.623 & 0.122 & 0.635 & 0.191 & & 0.114 & 0.051 & 0.125 & 0.080 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G3SLS} & 50 & & 0.639 & 0.491 & 0.805 & 0.715 & & 2.604 & 2.007 & 3.287 & 2.918 & & 0.083 & 0.114 & 0.141 & 0.174 \\
& & 100 & & 0.187 & 0.206 & 0.279 & 0.317 & & 1.077 & 0.685 & 1.276 & 1.047 & & 0.038 & 0.067 & 0.077 & 0.106 \\
& & 200 & & 0.006 & 0.183 & 0.183 & 0.296 & & 0.283 & 0.350 & 0.450 & 0.545 & & 0.013 & 0.040 & 0.042 & 0.062 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
\footnotesize
\setstretch{1}
\item \textit{Notes:} This table reports Monte Carlo results under network link misclassification. The parameter $\tau$ controls the intensity of misclassification and is common across individuals. Reported statistics include bias, standard deviation (SD), root mean squared error (RMSE), and inter-quantile range (IQR) for the peer effect ($\beta$=0.7), contextual ($\delta=1$), and direct ($\gamma=1$) Effects.
\end{tablenotes}
\end{threeparttable}
\normalsize
\end{table}
\end{landscape}
\begin{landscape}
\begin{table}[!htbp]
\centering
\small
\begin{threeparttable}
\caption{Estimator Performance under Unobserved Homophily (Design 3)}
\label{tab:design3_homophily_main}
\begin{tabular}{ccc c cccc c cccc c cccc}
\toprule
& & & &\multicolumn{4}{c}{Peer effects} & &\multicolumn{4}{c}{Contextual effects} & &\multicolumn{4}{c}{Direct effects} \\
\cmidrule(lr){5-8} \cmidrule(lr){10-13} \cmidrule(lr){15-18}
$m$ & Estimator & $n$ & & Bias & SD & RMSE & IQR & & Bias & SD & RMSE & IQR & & Bias & SD & RMSE & IQR \\
\midrule
\multirow{9}{*}{1} & \multirow{3}{*}{OLS} & 50 & & 0.077 & 0.033 & 0.084 & 0.051 & & -0.201 & 0.119 & 0.234 & 0.189 & & -0.109 & 0.082 & 0.137 & 0.129 \\
& & 100 & & 0.071 & 0.022 & 0.074 & 0.035 & & -0.188 & 0.079 & 0.204 & 0.127 & & -0.107 & 0.056 & 0.121 & 0.091 \\
& & 200 & & 0.080 & 0.019 & 0.082 & 0.032 & & -0.181 & 0.062 & 0.191 & 0.099 & & -0.088 & 0.037 & 0.096 & 0.060 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G2SLS} & 50 & & 0.019 & 0.038 & 0.043 & 0.060 & & -0.051 & 0.127 & 0.137 & 0.195 & & -0.024 & 0.088 & 0.091 & 0.141 \\
& & 100 & & 0.010 & 0.026 & 0.028 & 0.040 & & -0.026 & 0.084 & 0.088 & 0.139 & & -0.013 & 0.060 & 0.061 & 0.094 \\
& & 200 & & 0.015 & 0.022 & 0.027 & 0.036 & & -0.034 & 0.065 & 0.073 & 0.104 & & -0.018 & 0.039 & 0.043 & 0.062 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G3SLS} & 50 & & -0.013 & 0.049 & 0.051 & 0.074 & & 0.035 & 0.149 & 0.153 & 0.231 & & 0.014 & 0.100 & 0.101 & 0.160 \\
& & 100 & & -0.006 & 0.030 & 0.030 & 0.046 & & 0.018 & 0.095 & 0.096 & 0.149 & & 0.006 & 0.066 & 0.066 & 0.103 \\
& & 200 & & -0.006 & 0.026 & 0.027 & 0.043 & & 0.014 & 0.072 & 0.074 & 0.113 & & 0.002 & 0.043 & 0.043 & 0.067 \\
\midrule
\multirow{9}{*}{3} & \multirow{3}{*}{OLS} & 50 & & 0.118 & 0.040 & 0.124 & 0.064 & & -0.304 & 0.152 & 0.340 & 0.236 & & -0.168 & 0.124 & 0.209 & 0.198 \\
& & 100 & & 0.098 & 0.026 & 0.101 & 0.040 & & -0.258 & 0.095 & 0.275 & 0.149 & & -0.149 & 0.077 & 0.168 & 0.118 \\
& & 200 & & 0.113 & 0.023 & 0.115 & 0.037 & & -0.256 & 0.075 & 0.267 & 0.122 & & -0.124 & 0.052 & 0.135 & 0.081 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G2SLS} & 50 & & 0.056 & 0.050 & 0.075 & 0.078 & & -0.143 & 0.174 & 0.225 & 0.265 & & -0.078 & 0.137 & 0.157 & 0.209 \\
& & 100 & & 0.032 & 0.034 & 0.046 & 0.055 & & -0.082 & 0.110 & 0.137 & 0.171 & & -0.047 & 0.086 & 0.098 & 0.135 \\
& & 200 & & 0.048 & 0.027 & 0.055 & 0.041 & & -0.109 & 0.080 & 0.135 & 0.123 & & -0.053 & 0.056 & 0.077 & 0.091 \\
\noalign{\vskip 0.5em}
& \multirow{3}{*}{G3SLS} & 50 & & -0.019 & 0.072 & 0.074 & 0.103 & & 0.051 & 0.215 & 0.221 & 0.324 & & 0.026 & 0.141 & 0.143 & 0.222 \\
& & 100 & & -0.007 & 0.039 & 0.040 & 0.062 & & 0.024 & 0.128 & 0.131 & 0.202 & & 0.007 & 0.083 & 0.083 & 0.130 \\
& & 200 & & -0.007 & 0.034 & 0.035 & 0.054 & & 0.017 & 0.094 & 0.096 & 0.142 & & 0.003 & 0.057 & 0.057 & 0.091 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
\footnotesize
\setstretch{1}
\item \textit{Notes:} This table reports Monte Carlo results under unobserved homophily. The parameter $m$ controls the intensity of homophilous link formation and is common across individuals. Reported statistics include bias, standard deviation (SD), root mean squared error (RMSE), and inter-quantile range (IQR) for the peer effect ($\beta$=0.7), contextual ($\delta=1$), and direct ($\gamma=1$) Effects.
\end{tablenotes}
\end{threeparttable}
\normalsize
\end{table}
\end{landscape}
\section{Application to Publication Outcomes in Economics\label{emp}}
\subsection{Background}
The availability of online scientific research repositories has resulted in a stable source of data to uncover scholars' professional connections. Additionally, when linked with scholars' biographical public information, other types of professional connections beyond observed Co-authorship can be uncovered \citep{Colussi2018}. However, a major challenge when performing causal inference is that co-authorship connections are inherently correlated with publication outcomes, such as citation counts. However, where scholars completed their graduate training, although related to the co-authorship network, does not a priori directly affect publication outcomes. This setting, therefore, directly lends itself to the use of the proposed estimator to calculate social effects in research citations.
\subsection{Empirical Example: Peer-effects in Publishing}
\noindent To illustrate the proposed method, a data set of 1,628 peer-reviewed articles published between 2000-2006 in the top-four general-interest journals in Economics is used here to fit the following specific case of the estimating equation in \eqref{intro}:
\begin{equation}
\label{est_eq}
y_{i,r,t}=\alpha+\beta\sum_{j \neq i}w_{i,j,t}y_{j,r,t}+\sum_{j \neq i}w_{i,j, t}\widetilde{\mathbf{x}}_{j, r,t}^{\top}\boldsymbol{\delta}+ \mathbf{x}_{i, r, t}^{\top}\boldsymbol{\gamma}+ \lambda_{r} + \lambda_{t} + \lambda_{0} + v_{i,r,t}\text{,}
\end{equation}
\noindent The dependent variable $y_{i,r,t}$ denotes the natural logarithm of the total citations received by the article $i$ in the journal $r$ within eight years after its publication. The term $w_{i,j,t}$ corresponds to the $(i,j)$ element of the adjacency matrix $\mathbf{W}$ that represents the \emph{Co-authorship} network at time $t$. A detailed description of this data set and the corresponding summary statistics can be found in \cite{netivreg_g3sls} and \cite{EstradaDingRosales2025}.
The vector of controls $\mathbf{x}_{i,r,t}$ includes indicators for whether current or former editors of the journal $r$ at time $t$ appear among the authors of the article $i$ (\texttt{Editor}) and for whether the co-authors of article $i$ are of different genders (\texttt{Different Gender}), which is set to zero for single-author papers. Additional article-level covariates comprise the total number of pages (\texttt{Number of Pages}), authors (\texttt{Number of Authors}), and bibliographic references (\texttt{Number of References}), as well as a dummy variable identifying articles that are isolated within the network (\texttt{Isolated}). Contextual or peer effects are computed only for the \texttt{Editor} and \texttt{Different Gender} variables, summarized in $\widetilde{\mathbf{x}}_{j,r,t}$. The specification also includes fixed effects for the journal ($\lambda_{r}$) and year ($\lambda_{t}$).
The structural error term $v_{i,r,t}$ satisfies $E_{\mathbf{X},\mathbf{W}}[\mathbf{v}] \neq 0$, reflecting potential endogeneity in the \emph{Co-authorship} network. To address this, the \emph{Alumni} network is employed as the exogenous matrix $\mathbf{W}_{0}$ in the identification theorem (\ref{T1}). Additionally, to mitigate possible deviations from the exogeneity condition, the model incorporates the corresponding \emph{Alumni-based} components ($\lambda_{0}$), as specified in equation (\ref{E25}).
The model \eqref{est_eq} is estimated in a rolling-regression setting for $t=2002$, $2003$, $2004$, $2005$, and $2006$; i.e., the estimation sample each year includes those from previous years. The results for 2000 and 2001 are not included because they suffer from degrees-of-freedom problems given the specification \eqref{est_eq}. Table \ref{empt3} summarizes the results using the proposed G3SLS estimator, while Tables \ref{empt4} and \ref{empt5} in \ref{Appendix_E} report the OLS and G2SLS results, respectively. In all cases, the asymptotic standard errors are clustered at the corresponding network component level. This is a natural way of clustering when utilizing network data because each component corresponds to a portion of the network that is disconnected from the others, allowing for articles within each component to be correlated but not between disconnected components.
The results show that the estimated peer effect ($\widehat{\beta}$) is consistently positive and becomes statistically significant at the 5\% level from 2004 onward, peaking at 0.676 in 2005. This confirms the presence of significant positive citation spillovers among articles connected within the \emph{Co-authorship} network. Turning to direct effect estimates ($\widehat{\boldsymbol{\gamma}}$), the indicator of gender-diverse research teams exhibits a stable, positive, and statistically significant impact on citation counts in all accumulated samples. This result demonstrates that gender diversity within research teams positively influences the academic impact of a paper, an outcome that aligns with recent evidence \cite{RHJS:2024} showing that gender-diverse teams also produce more readable and accessible research. Furthermore, standard article characteristics such as the number of pages and bibliographic references remain robust, positive predictors of citations, while isolated articles consistently incur a significant citation penalty. In contrast, the estimated contextual effects ($\widehat{\boldsymbol{\delta}}$) related to the presence of the editorial and the diversity of the gender are statistically indistinguishable from zero in all specifications. Finally, empirical evidence supporting the underlying exclusion restriction and relevance conditions necessary for identification is provided in \ref{Appendix_E} in the Supplementary Material.
\begin{table}
\centering
\small
\caption{Estimation Results for Social and Direct Effects}
\begin{tabular}{lccccc}
\hline
& \multicolumn{5}{c}{\emph{Co-author Network}} \\
\cline{2-6}
& 2002 & 2003 & 2004 & 2005 & 2006 \\
\hline
Peer Effects ($\widehat{\beta}$) & 0.520 & 0.476 & 0.542** & 0.676** & 0.570** \\
& (0.362) & (0.316) & (0.261) & (0.283) & (0.242) \\ \hline
Contextual Effects ($\boldsymbol{\widehat{\delta}}$) & & & & & \\
\hspace{0.2cm} \texttt{Editor} & 1.790 & -0.971 & -2.910 & -5.355 & -5.044 \\
& (5.363) & (3.587) & (3.111) & (4.135) & (4.188) \\
\hspace{0.2cm} \texttt{Different Gender} & -2.096 & -1.661 & -1.922 & -2.967 & -1.559 \\
& (2.829) & (1.793) & (1.838) & (2.65) & (1.737) \\ \hline
Direct Effects ($\boldsymbol{\widehat{\gamma}}$) & & & & \\
\hspace{0.2cm} \texttt{Editor} & 0.173 & 0.027 & -0.037 & 0.025 & 0.057 \\
& (0.116) & (0.124) & (0.131) & (0.116) & (0.130) \\
\hspace{0.2cm} \texttt{Different Gender} & 0.219* & 0.205* & 0.221** & 0.166* & 0.143* \\
& (0.131) & (0.108) & (0.092) & (0.085) & (0.080) \\
\hspace{0.2cm} \texttt{Number of Pages} & 0.029*** & 0.027*** & 0.023*** & 0.019*** & 0.018*** \\
& (0.004) & (0.004) & (0.004) & (0.004) & (0.004) \\
\hspace{0.2cm} \texttt{Number of Authors} & 0.072 & 0.087* & 0.072 & 0.097*** & 0.076** \\
& (0.060) & (0.050) & (0.044) & (0.038) & (0.031) \\
\hspace{0.2cm} \texttt{Number of References} & 0.012*** & 0.012*** & 0.011*** & 0.011*** & 0.012*** \\
& (0.003) & (0.002) & (0.002) & (0.002) & (0.001) \\
\hspace{0.2cm} \texttt{Isolated} & -0.223** & -0.236*** & -0.353*** & -0.399*** & -0.407*** \\
& (0.132) & (0.110) & (0.103) & (0.092) & (0.088) \\
\hline
$n$ & 729 & 961 & 1187 & 1412 & 1628 \\
$R^{2}$ & 0.172 & 0.199 & 0.178 & 0.125 & 0.148 \\
\hline
\end{tabular}
\vspace{0.4cm}
\begin{minipage}{1\textwidth}
\footnotesize
Note: Standard errors are in parentheses and are clustered at the specific network's components. Stars follow the key: * $p$ $<$ 0.10, ** $p$ $<$ 0.05, and *** $p$ $<$ 0.01, where $p$ stands for $p$-values. $R^2$ are calculated as the squared of the sample correlation coefficients between the observed outcomes and their fitted values. All specifications include indicator variables for Journal, Year and Alumni Network Components.
\end{minipage}
\label{empt3}
\end{table}
\section{Discussion\label{conclusion}}
\noindent In this paper, we propose a novel and computationally simple way to identify and consistently estimate social parameters in a linear-in-means model with endogenous network formation. Identification can be achieved by the inclusion of the multiplex network data structure, where at least one of the layers can be assumed to be exogenous (potentially pre-determined), and the multiplex structure is such that the layers correlate with each other. Our research shows that peer and contextual effects can be uniquely recovered from an estimating sample. Unlike current alternatives that require smoothing techniques and/or Bayesian methods, the resulting estimator is simple to compute, and it is already implemented in a user-written \texttt{Stata} command, see \cite{netivreg_stata_journal}. The asymptotic normality of the proposed multi-step estimator is established, and a consistent estimator of its asymptotic variance-covariance matrix is proposed for performing inference.
The type of endogeneity allowed in our framework is general enough to encompass settings with measurement error, sample selection, and correlation between the unobservables driving network formation and outcomes in a linear-in-means model. It is argued that the full observability of a multiplex data structure, with at least one of the layers being exogenous, is not a limiting data requirement that can easily be constructed in some cases based purely on characteristics the researcher usually observes. With this in mind, an empirical application is presented where the tools of web scraping and text mining are used to construct a data set consisting of all peer-reviewed research articles published in 4 of economics' top general-interest journals between 2000 and 2006. Using publicly available information on where the authors of these publications obtained their Ph.D. degrees from, an \emph{Alumni} network is constructed and argued to be pre-determined yet correlated with the observed \emph{Co-authorship} ties among these scholars. The results show the existence of positive peer effects in terms of citations among peer-reviewed research articles connected through co-authorship connections of their authors, as well as significant positive effects of research teams that are gender diverse on the quality of a paper measured in terms of citation outcomes similar to \cite{RHJS:2024}.
The results in our paper have the potential to be extended in different directions. From a theoretical perspective, this paper introduces the concept of multidimensional networks and uses, in particular, the structure of multiplex networks. Although the former is an active research topic in other fields, this paper proposes its usage for causal inference in the Social Sciences. An interesting extension would be to study the dynamics of peer effects using the setting proposed in this paper. Similarly, as with all IV-based estimation procedures, issues pertaining to weak or invalid instruments, many-instrument problems, or small-sample performance of the proposed estimator are valid concerns, but they are beyond the scope of this paper and therefore left for future research.
\section{Acknowledgements}
We are grateful to the Co-Editor, the Associate Editor, and two anonymous referees for their constructive comments and suggestions that significantly improved this manuscript. We also thank the discussants and participants at various conferences, workshops, and seminars where earlier versions of this research were presented for their valuable feedback. Portions of this paper are based on results from Juan Estrada's doctoral dissertation, \emph{Causal Inference in Multilayered Networks}, at Emory University. Kim P. Huynh dedicates this paper in the memory of Tony S. Wirjanto.