2. Latent variable network models

In this section, we briefly summarize a number of latent variable network models in the literature that are used to capture excess dependence in network observations. All latent variable network models we consider here may be written in the common form

(2) \begin{align} &{\mathbb{P}}(y_{ij} = 1 ) ={\mathbb{P}} \!\left ( \mu _{ij} + f_{\boldsymbol{\theta }} (\textbf{v}_i, \textbf{v}_j) +{\textbf{x}}i_{ij} \gt 0 \right ), \quad \\ &\textbf{v}_i \stackrel{iid}{\sim } (\textbf{0}, \boldsymbol{\Sigma }_v), \qquad{\textbf{x}}i_{ij} \stackrel{iid}{\sim } \textrm{N}(0, \sigma ^2),\nonumber \end{align}

where $\textbf{v}_i \in{\mathbb{R}}^K$ with mean $\textbf{0}$ and covariance matrix $\boldsymbol{\Sigma }_v$ , and $\mu _{ij}$ is fixed. We avoid specifying a distribution for the latent vectors $\{ \textbf{v}_{i} \}_{i=1}^n$ , although they are often taken to be multivariate Gaussian. We set the total variance of the latent variable representation to be $1 = \sigma ^2 + \textrm{var}[f_{\boldsymbol{\theta }} (\textbf{v}_i, \textbf{v}_j)]$ , since it is not identifiable. The function of the latent variables $f_{\boldsymbol{\theta }} \;:\;{\mathbb{R}}^K \times{\mathbb{R}}^K \rightarrow{\mathbb{R}}$ , parametrized by $\boldsymbol{\theta }$ , serves to distinguish the latent variable network models discussed below. Regression latent variable network models are formed when the latent mean is represented as a linear function of exogenous covariates ${\textbf{x}}_{ij} \in{\mathbb{R}}^p$ , such that $\mu _{ij} ={\textbf{x}}_{ij}^T{\boldsymbol{\beta }}$ . The latent nodal random vectors $\{ \textbf{v}_{i} \}_{i=1}^n$ represent excess network dependence—beyond the mean $\mu _{ij}$ . Since relations $y_{ij}$ and $y_{ik}$ share latent vector $\textbf{v}_i$ corresponding to shared actor $i$ , and thus, $y_{ij}$ and $y_{ik}$ have related distributions through the latent function $f_{\boldsymbol{\theta }} (\textbf{v}_i, \textbf{v}_j)$ . Many popular models for network data may be represented as in (2), such as the social relations model, the latent position model, and the latent eigenmodel.

2.1 Social relations model

The social relations model was first developed for continuous, directed network data (Warner et al., Reference Warner, Kenny and Stoto1979; Wong, Reference Wong1982; Snijders & Kenny, Reference Snijders and Kenny1999). In the social relations model for binary network data (Hoff, Reference Hoff2005), $f_{\boldsymbol{\theta }}(\textbf{v}_i, \textbf{v}_j) = \textbf{v}_i + \textbf{v}_j$ and $\textbf{v}_i = a_i \in \mathbb{R}$ for each actor $i$ , such that

(3) \begin{align} &{\mathbb{P}}(y_{ij} = 1 ) ={\mathbb{P}} \!\left ({\textbf{x}}_{ij}^T{\boldsymbol{\beta }} + a_i + a_j +{\textbf{x}}i_{ij} \gt 0 \right ), \quad \\ &a_i \stackrel{iid}{\sim } \!\left(0, \sigma ^2_a\right), \qquad{\textbf{x}}i_{ij} \stackrel{iid}{\sim } \textrm{N}(0, \sigma ^2). \nonumber \end{align}

Each actor’s latent variable $\{ a_i \}_{i=1}^n$ may be thought of as the actor’s sociability: large values of $a_i$ correspond to actors with a higher propensity to form relations in the network. The random $\{ a_i \}_{i=1}^n$ in (3) also accounts for the excess correlation in network data; any two relations that share an actor, for example, $y_{ij}$ and $y_{ik}$ , are marginally correlated.

2.2 Latent position model

A more complex model for representing excess dependence in social network data is the latent position model (Hoff et al., Reference Hoff, Raftery and Handcock2002). The latent position model extends the idea of the social relations model by giving each actor $i$ a latent position $\textbf{u}_i$ in a Euclidean latent space, for example ${\mathbb{R}}^{K}$ . Then, actors whose latent positions are closer together in Euclidean distance are more likely to share a relation:

(4) \begin{align} &{\mathbb{P}}(y_{ij} = 1 ) ={\mathbb{P}} \!\left ({\textbf{x}}_{ij}^T{\boldsymbol{\beta }} + a_i + a_j - || \textbf{u}_i - \textbf{u}_j ||_2 +{\textbf{x}}i_{ij} \gt 0\right ), \quad \\ &a_i \stackrel{iid}{\sim } \!\left(0, \sigma ^2_a\right), \quad \textbf{u}_i \stackrel{iid}{\sim } (0, \boldsymbol{\Sigma }_u), \quad{\textbf{x}}i_{ij} \stackrel{iid}{\sim } \textrm{N}(0, \sigma ^2). \nonumber \end{align}

In the form of (2), the latent position model contains latent random vector $\textbf{v}_{i} = [a_i, \textbf{u}_i]^T \in{\mathbb{R}}^{K+1}$ , and $f_{\boldsymbol{\theta }}(\textbf{v}_{i}, \textbf{v}_j ) = a_i + a_j - || \textbf{u}_i - \textbf{u}_j ||_2$ . Hoff et al. (Reference Hoff, Raftery and Handcock2002) show that the latent position model is capable of representing transitivity, that is, when $y_{ij} = 1$ and $y_{jk} = 1$ , it is more likely that $y_{ik} = 1$ . Models that are transitive often display a pattern observed in social network data: a friend of my friend is also my friend (Wasserman & Faust, Reference Wasserman and Faust1994).

2.3 Latent eigenmodel

The latent eigenmodel also associates each actor with a latent position $\textbf{u}_i$ in a latent Euclidean space; however, the inner product between latent positions (weighted by symmetric parameter matrix $\Lambda$ ) measures the propensity of actors $i$ and $j$ to form a relation, rather than the distance between positions (Hoff, Reference Hoff2008):

(5) \begin{align} &{\mathbb{P}}(y_{ij} = 1 ) ={\mathbb{P}} \!\left ({\textbf{x}}_{ij}^T{\boldsymbol{\beta }} + a_i + a_j + \textbf{u}_i^T \Lambda \textbf{u}_j +{\textbf{x}}i_{ij} \gt 0 \right ), \quad \\ &a_i \stackrel{iid}{\sim } \!\left(0, \sigma ^2_a\right), \quad \textbf{u}_i \stackrel{iid}{\sim } (0, \boldsymbol{\Sigma }_u), \quad{\textbf{x}}i_{ij} \stackrel{iid}{\sim } \textrm{N}(0, \sigma ^2). \nonumber \end{align}

In the context of (2), the function $f_{\boldsymbol{\theta }}(\textbf{v}_{i}, \textbf{v}_{j}) = a_i + a_j + \textbf{u}_i^T \Lambda \textbf{u}_j$ for the latent eigenmodel, where the parameters $\boldsymbol{\theta }$ are the entries in $\Lambda$ and $\textbf{v}_{i} = [a_i, \textbf{u}_i]^T \in{\mathbb{R}}^{K+1}$ . Hoff (Reference Hoff2008) shows that the latent eigenmodel is capable of representing transitivity and that the latent eigenmodel generalizes the latent position model given sufficiently large dimension of the latent vectors $K$ .

In addition to transitivity, a second phenomenon observed in social networks is structural equivalence, wherein different groups of actors in the network form relations in a similar manner to others in their group. One form of structural equivalence is associative community structure, where the social network may be divided into groups of nodes that share many relations within group, but relatively few relations across groups. Such behavior is common when cliques are formed in high school social networks or around subgroups in online social networks. A form of structural equivalence is when actors in a given group are more likely to form relations with actors in other groups than with actors in their own group, for example, in networks of high-functioning brain regions when performing cognitively demanding tasks (Betzel et al., Reference Betzel, Bertolero and Bassett2018). Two models that are aimed at identifying subgroups of nodes that are structurally equivalent are the latent class model of Nowicki & Snijders (Reference Nowicki and Snijders2001) and the mixed membership stochastic blockmodel (Airoldi et al., Reference Airoldi, Blei, Fienberg and Xing2008). Hoff (Reference Hoff2008) shows that the latent eigenmodel is capable of representing stochastic equivalence in addition to transitivity and that the latent eigenmodel generalizes latent class models given sufficiently large dimension of the latent vectors $K$ . For this reason, we focus on the latent eigenmodel, and the simpler social relations model, as reference models in this paper.

3. Exchangeable network models

To motivate the formulation of the proposed model, we briefly discuss the theory of exchangeable random network models and their relationship to latent variable network models. A random network model for $\{\epsilon _{ij} \}_{ij}$ is exchangeable if the distribution of $\{\epsilon _{ij} \}_{ij}$ is invariant to permutations of the actor labels, that is, if

(6) \begin{align}{\mathbb{P}}\!\left ( \{\epsilon _{ij} \}_{ij} \right ) ={\mathbb{P}} \!\left ( \{\epsilon _{\pi (i) \pi (j) } \}_{ij} \right ), \quad \end{align}

for any permutation $\pi (.)$ . There is a rich theory of exchangeable network models, dating back to work on exchangeable random matrices (Hoover, Reference Hoover1979; Aldous, Reference Aldous1981), upon which we draw in this section.

All the latent variable network models discussed in Section 2 have latent error networks $\{\epsilon _{ij} \}_{ij}$ that are exchangeable, where we define $\epsilon _{ij} = f_{\boldsymbol{\theta }}(\textbf{v}_{i}, \textbf{v}_{j}) +{\textbf{x}}i_{ij}$ from (2), the random portion of a general latent variable network model. Further, under constant mean $\mu _{ij} = \mu$ , all the latent variable network models for the observed network $\{ y_{ij} \}_{ij}$ in Section 2 are exchangeable. In fact, any exchangeable network model may be represented by a latent variable network model. Specifically, the theory of exchangeable network models states that every exchangeable random network model may be represented in the following form (see, for example, Lovász & Szegedy, Reference Lovász and Szegedy2006; Kallenberg, Reference Kallenberg2006):

(7) \begin{align} &{\mathbb{P}}(\epsilon _{ij} = 1 ) ={\mathbb{P}} \!\left ( \mu + h (u_i, u_j) +{\textbf{x}}i_{ij} \gt 0 \right ), \quad \\ &u_i \stackrel{iid}{\sim } \textrm{Uniform}(0, 1), \quad{\textbf{x}}i_{ij} \stackrel{iid}{\sim } \rm{N}(0, \sigma ^2), \nonumber \end{align}

where the function $h \;:\; [0,1] \times [0,1] \rightarrow{\mathbb{R}}$ has finite integral $\int _{[0,1] \times [0,1]} h(u,v) du dv \lt \infty$ and serves to distinguish the various exchangeable network models. It can be shown that (7) is equivalent to the graphon representation of exchangeable random network models, where the graphon is the canonical probabilistic object of exchangeable random network models (Lovász & Szegedy, Reference Lovász and Szegedy2006; Borgs et al., Reference Borgs, Chayes, Cohn and Zhao2014). Noting that we may always map the random scalar $u_i$ to some random vector $\textbf{v}_i$ , the expression in (7) illustrates how every exchangeable random network model may be represented by a latent variable network model in the sense of (2).

3.1 Covariance matrices of exchangeable network models

The expression in (7) shows that any exchangeable network model for binary network data must correspond to a latent random network $\{\epsilon _{ij} \}_{ij}$ that is continuous and exchangeable. The covariance matrix of any undirected exchangeable network model has the same form and contains at most two unique nonzero values. (Marrs et al. (Reference Marrs, Fosdick and McCormick2017) show that directed exchangeable network models with continuous values all have covariance matrices of the same form with at most five unique nonzero terms). This fact can be seen by simply considering the ways that any pair of relations can share an actor. In addition to a variance, the remaining covariances are between relations that do and do not share an actor:

(8) \begin{align} \textrm{var}[\epsilon _{ij}] = \sigma ^2_\epsilon, \quad \textrm{cov}[\epsilon _{ij}, \epsilon _{ik}] \;:\!=\; \rho, \quad \textrm{cov}[\epsilon _{ij}, \epsilon _{kl}] = 0, \quad \end{align}

where the indices $i,j,k,$ and $l$ are unique. It is easy to see the second equality holds for any pair of relations that share an actor by the exchangeability property, that is, by permuting the actor labels. The third equality results from the fact that the only random elements in (7) are the actor random variables $u_i$ , $u_j$ , and the random error ${\textbf{x}}i_{ij}$ . When the random variables corresponding to two relations $\epsilon _{ij}$ and $\epsilon _{kl}$ share no actor, the pair of relations are independent by the generating process. Finally, we note that exchangeable network models have relations that are marginally identically distributed, and thus, relations therein have the same expectation and variance. That said, in the generalized linear regression case of (2), the means $\mu _{ij} ={\textbf{x}}^T_{ij}{\boldsymbol{\beta }}$ are non-constant, and thus, the observations $\{y_{ij} \}_{ij}$ are not exchangeable; only the latent error network $\{ \epsilon _{ij} \}_{ij}$ is exchangeable in the generalized linear regression case. In the proposed model, rather than put forth a particular parametric model for the latent network $\{ \epsilon _{ij} \}_{ij}$ , we simply model the covariance structure outlined in (8), which is sufficient to represent the covariance structure of any exchangeable network model for the errors.

4. The probit exchangeable (PX) model

In this section, we propose the probit exchangeable network regression model, which we abbreviate as the “PX” model. In the PX model, the vectorized mean of the network is characterized by a linear combination of covariates, $\textbf{X}{\boldsymbol{\beta }}$ , where $\boldsymbol{\beta }$ is a $p$ -length vector of coefficients that are the subject of inference and $\textbf{X}$ is a $\binom{n}{2} \times p$ matrix of covariates. The excess network dependence beyond that captured in $\textbf{X}{\boldsymbol{\beta }}$ is represented by an unobservable mean-zero error vector $\boldsymbol{\epsilon }$ , a vectorization of $\{ \epsilon _{ij} \}_{ij}$ , that is exchangeable in the sense of (6). The PX model is

(9) \begin{align} &{\mathbb{P}}(y_{ij} = 1) ={\mathbb{P}} \!\left ({\textbf{x}}_{ij}^T{\boldsymbol{\beta }} + \epsilon _{ij} \gt 0 \right ), \quad \\ &\boldsymbol{\epsilon } \sim \textrm{N}( \textbf{0}, \boldsymbol{\Omega }), \nonumber \end{align}

where we note that the variance of $\epsilon _{ij}$ is not identifiable, and thus, we choose $\textrm{var}[\epsilon _{ij}] = 1$ without loss of generality. We focus on normally distributed unobserved errors $\boldsymbol{\epsilon }$ in this paper; however, other common distributions, such as the logistic distribution, could be used. We note that the normal distribution assumption implies that (9) is a typical probit regression model, but with correlation among the observations due to network structure.

As discussed in Section 3, under the exchangeability assumption, the covariance matrix of the latent error network $\textrm{var}[\boldsymbol{\epsilon }] = \boldsymbol{\Omega }$ has at most two unique nonzero parameters. Taking $\textrm{var}[\epsilon _{ij}] = 1$ , the covariance matrix of $\boldsymbol{\epsilon }$ has a single parameter $\rho = \textrm{cov}[\epsilon _{ij}, \epsilon _{ik}]$ . We may thus write

(10) \begin{align} \boldsymbol{\Omega }(\rho ) = \mathcal{S}_1 + \rho \mathcal{S}_2, \quad \end{align}

where we define the binary matrices $\{ \mathcal{S}_{i} \}_{i=1}^3$ indicating unique entries in $\boldsymbol{\Omega }$ . The matrix $\S _1$ is a diagonal matrix indicating the locations of the variance in $\boldsymbol{\Omega }$ , and $\S _2$ and $\S _3$ indicate the locations in $\boldsymbol{\Omega }$ corresponding to the covariances $\textrm{cov}[\epsilon _{ij}, \epsilon _{ik}]$ , and $\textrm{cov}[\epsilon _{ij}, \textrm{cov}[\epsilon _{ij}, \epsilon _{kl}]$ , respectively, where the indices $i,j,k,$ and $l$ are unique.

The PX model unifies many of the latent variable network models discussed in Sections 2 and 3. Similar to (7), the PX model may be seen as representing the covariance structure of the latent variables $\{ f_{\boldsymbol{\theta }}(\textbf{v}_i, \textbf{v}_j) +{\textbf{x}}i_{ij} \}_{ij}$ with $\{ \epsilon _{ij} \}_{ij}$ , the unobservable error network of the PX model in (9). As both networks $\{ f_{\boldsymbol{\theta }}(\textbf{v}_i, \textbf{v}_j) +{\textbf{x}}i_{ij} \}_{ij}$ and $\{ \epsilon _{ij} \}_{ij}$ are exchangeable, they have covariance matrices of the same form (see discussion in Section 3). As every exchangeable random network model may be represented by a latent variable network model, the PX model may represent the latent correlation structure of any exchangeable network model, yet without specifying a particular exchangeable model. Further, we now show that the PX model is equivalent to the social relations model under certain conditions.

Exact distributional equivalence between the PX model and latent variable models other than the social relations model will typically not hold. For example, the latent eigenmodel in (5) includes non-Gaussian random variables, so that exact distributional equivalence is impossible. Similarly, it appears likely that the general latent variable model in (2) may generate non-Gaussian random variables through the function $f_{\boldsymbol{\theta }} (\textbf{v}_i, \textbf{v}_j)$ . Importantly however, there does exist $\rho$ such that the covariance of the latent errors of every pair of relations, $\textrm{cov}[\epsilon _{ij}, \, \epsilon _{kl}]$ , is equal to the covariance of the latent errors in any exchangeable latent variable model, $\textrm{cov}[f_{\boldsymbol{\theta }} (\textbf{v}_i, \textbf{v}_j) +{\textbf{x}}i_{ij},\, f_{\boldsymbol{\theta }} (\textbf{v}_k, \textbf{v}_l) +{\textbf{x}}i_{kl}]$ . Hence, the PX model may be seen as a generalized exchangeable latent variable model that focuses all modeling effort on the first two moments of the data.

Proposition 4.1 states that the PX model and social relations model are equivalent under normality of their latent error networks. In principle, the social relations model is simply a generalized linear mixed model; however, existing software packages, such as lme4 in R (Bates et al., Reference Bates, Mächler, Bolker and Walker2015), do not appear to accommodate the random effects specification of the social relations model in (3) since the indices $i$ and $j$ pertain to random effects $a_i$ and $a_j$ from the same set (as opposed to $a_i$ and $b_j$ in a random crossed design). Nevertheless, the estimation scheme proposed in Section 5 employs the same strategies as those commonly used to estimate generalized linear mixed models (Littell et al., Reference Littell, Milliken, Stroup, Wolfinger and Oliver2006; Gelman & Hill, Reference Gelman and Hill2006). In the estimation algorithm in lme4, the marginal likelihood of the data is approximated and then maximized using numerical approximations with respect to $\boldsymbol{\beta }$ and random effects variance, for example $\sigma ^2_a$ in the social relations model. Rather than an approximate likelihood, we propose maximizing the true likelihood with respect to $\boldsymbol{\beta }$ and $\rho$ , yet also use numerical approximations to accomplish this maximization.

It is important to note that, although the latent errors $\{ \epsilon _{ij} \}_{ij}$ in the PX model form an exchangeable random network, the random network $y_{ij}$ represented by the PX model is almost certainly not exchangeable. For example, each $y_{ij}$ may have a different marginal expectation $\Phi \big({\textbf{x}}_{ij}^T{\boldsymbol{\beta }}\big)$ . Then, the relations in the network are not marginally identically distributed, which is a necessary condition for exchangeability. Further, the covariances between pairs of relations, say $y_{ij}$ and $y_{ik}$ , depend on the marginal expectations:

\begin{align*} \textrm{cov}\big[y_{ij}, y_{ik} \big] &= E \!\left [ y_{ij} y_{ik} \right ] - E \!\left [ y_{ij} \right ] E \!\left [ y_{ik} \right ] = \int _{-{\textbf{x}}^T_{ij}{\boldsymbol{\beta }}}^\infty \int _{-{\textbf{x}}^T_{ik}{\boldsymbol{\beta }}}^\infty dF_\rho - \Phi \!\left({\textbf{x}}^T_{ij}{\boldsymbol{\beta }}\right) \Phi \!\left({\textbf{x}}^T_{ik}{\boldsymbol{\beta }}\right). \end{align*}

Here, $dF_\rho$ is the bivariate standard normal distribution with correlation $\rho$ . Since the covariance $\textrm{cov}[y_{ij}, y_{ik} ]$ depends on the latent means ${\textbf{x}}^T_{ij}{\boldsymbol{\beta }}$ and ${\textbf{x}}^T_{ik}{\boldsymbol{\beta }}$ , $\textrm{cov}[y_{ij}, y_{ik} ]$ is only equal to $\textrm{cov}[y_{ab}, y_{ac} ]$ when the latent means are equal. As a result, although the covariance matrix of the unobserved errors $\boldsymbol{\Omega }$ is of a simple form with entries $\{1, \rho, 0 \}$ , the covariances between elements of the vector of observed relations $\textbf{y}$ are heterogeneous (in general) and depend on $\rho$ in a generally more complicated way.

5. Estimation

In this section, we propose an estimator of $\{{\boldsymbol{\beta }}, \rho \}$ in the PX model that approximates the maximum likelihood estimator (MLE). The algorithm we propose is based on the EM algorithm (Dempster et al., Reference Dempster, Laird and Rubin1977). Although the covariance matrix for the PX model is highly structured, as in (10), a closed-form expression for the MLE does not appear available. While we explored pseduo-likelihood pairwise approximations (also called “composite likelihoods” in some literature) to the complete PX likelihood (Heagerty & Lele, Reference Heagerty and Lele1998), we found no substantial advantage—neither in performance nor runtime—over the proposed estimation scheme in this paper.

The proposed estimation algorithm consists of alternating computation of the expected complete likelihood with maximization with respect to $\rho$ and $\boldsymbol{\beta }$ , iterating until convergence. Since the algorithm iterates expectation and two maximization steps, we term it the EMM algorithm. To improve algorithm efficiency, we initialize $\boldsymbol{\beta }$ at the ordinary probit regression estimator (assuming independence of the latent errors), and initialize $\rho$ with a mixture estimator based on possible values of $\rho$ such that $\boldsymbol{\Omega }$ is positive definite, as detailed in Section A.1. The complete EMM algorithm is presented in Algorithm 1. In the following text, we detail the EMM algorithm, beginning with maximization with respect to $\rho$ , and then proceeding to maximization with respect to $\boldsymbol{\beta }$ . We define $\gamma _i = E\!\left[\boldsymbol{\epsilon }^T \mathcal{S}_i \boldsymbol{\epsilon } \mid \textbf{y},{\widehat{\rho }}^{(\nu )}, \widehat{{\boldsymbol{\beta }}}^{(\nu )} \right]/ |\Theta _i|$ , where $\Theta _i$ is the set of relation pairs indicated by binary matrices $\mathcal{S}_i$ . By default, we typically set $\tau = 10^{-2}$ and $\delta = 10^{-1}$ .

Algorithm 1. EMM estimation of the PX model

5.1 Expectation

Consider the log-likelihood, $\ell _{\textbf{z}}$ , of the latent continuous random vector $\textbf{z}$ . Taking the expectation of $\ell _{\textbf{z}}$ conditional on $\textbf{y}$ , the expectation step for a given iteration $\nu$ of the EM algorithm is

(11) \begin{align} &E\!\left[\ell _{\textbf{z}} \, | \, \textbf{y}, \rho ={\widehat{\rho }}^{(\nu )},{\boldsymbol{\beta }}=\widehat{{\boldsymbol{\beta }}}^{(\nu )} \right] \nonumber \\ & \hspace{.25in} = -\frac{1}{2}\textrm{log } 2 \pi |\boldsymbol{\Omega }| - \frac{1}{2} E\!\left [ (\textbf{z} - \textbf{X}{\boldsymbol{\beta }})^T \boldsymbol{\Omega }^{-1} (\textbf{z} - \textbf{X}{\boldsymbol{\beta }}) \, | \, \textbf{y}, \rho ={\widehat{\rho }}^{(\nu )},{\boldsymbol{\beta }}=\widehat{{\boldsymbol{\beta }}}^{(\nu )} \right ], \quad \end{align}

where ${\widehat{\rho }}^{(\nu )}$ and $\widehat{{\boldsymbol{\beta }}}^{(\nu )}$ are the estimators of $\rho$ and $\boldsymbol{\beta }$ at iteration $\nu$ . In discussing the maximization step for $\rho$ , we will show that that the $\rho$ update depends on the data through the expectations denoted by $\gamma _i$ for $i \in \{1,2,3 \}$ . In discussing the maximization step for $\beta$ , we will show that that the $\beta$ update depends on the data only through the expectation $E\!\left[\boldsymbol{\epsilon } \, | \, \textbf{y},{\widehat{\rho }}^{(\nu )}, \widehat{{\boldsymbol{\beta }}}^{(\nu )} \right]$ .

5.1.1 Approximations

The computation of $E\!\left[\boldsymbol{\epsilon } \, | \, \textbf{y},{\widehat{\rho }}^{(\nu )}, \widehat{{\boldsymbol{\beta }}}^{(\nu )}\right]$ in (14) is nontrivial, as it is a $\binom{n}{2}$ -dimensional truncated multivariate normal integral. We exploit the structure of $\boldsymbol{\Omega }$ to compute $E\!\left[\boldsymbol{\epsilon } \, | \, \textbf{y},{\widehat{\rho }}^{(\nu )}, \widehat{{\boldsymbol{\beta }}}^{(\nu )}\right]$ using the law of total expectation. A Newton-Raphson algorithm, along with an approximate matrix inverse, is employed to compute an approximation of $E\!\left[\boldsymbol{\epsilon } \, | \, \textbf{y},{\widehat{\rho }}^{(\nu )}, \widehat{{\boldsymbol{\beta }}}^{(\nu )}\right]$ . Details of the implementation of the approximations are given in Section A.2.

The expectations $\{ \gamma _i \}_{i=1}^3$ require the computation of $\binom{n}{2}$ -dimensional truncated multivariate normal integrals, which are onerous for even small networks. Thus, we make two approximations to $\{ \gamma _i \}_{i=1}^3$ to reduce the runtime of the EMM algorithm. First, we compute the expectations conditioning only on the entries in $\textbf{y}$ that correspond to the entries in $\boldsymbol{\epsilon }$ being integrated, for example, instead of computing $E[\epsilon _{jk} \epsilon _{lm} \, | \, \textbf{y}]$ , we compute $E[\epsilon _{jk} \epsilon _{lm} \, | \, y_{jk}, y_{lm}]$ . This first approximation is most appropriate when $\rho$ is small, since $y_{lm}$ is maximally informative for $\epsilon _{jk}$ when $\rho$ is large (for $l,m,j,$ and $k$ distinct). Second, we find empirically that $\gamma _2 = E\!\left[\boldsymbol{\epsilon }^T \mathcal{S}_2 \boldsymbol{\epsilon } \, | \, \textbf{y} \right]/ |\Theta _2|$ is approximately linear in $\rho$ , since this sample mean of conditional expectations concentrates around a linear function of $\rho$ . Thus, we compute $\gamma _2$ for $\rho =0$ and $\rho = 1$ and use a line connecting these two values to compute $\gamma _2$ for arbitrary values of $\rho$ (see evidence of linearity of $\gamma _2$ for the political books network in Appendix E). The details of the approximations to $\{ \gamma _i \}_{i=1}^3$ are given in Section A.3.

5.2 Maximization with respect to $\boldsymbol\rho$

To derive the maximization step for $\rho$ , we use the method of Lagrange multipliers, since differentiating (11) directly with respect to $\rho$ gives complex nonlinear equations that are not easily solvable. We first define the set of parameters $\{\phi _i \}_{i=1}^3$ , representing the variance and two possible covariances in $\boldsymbol{\Omega }$ ,

\begin{align*} \textrm{var}\!\left[\epsilon _{ij}\right] = \phi _1, \quad \textrm{cov}\!\left[\epsilon _{ij}, \epsilon _{ik}\right] = \phi _2 = \rho, \quad \textrm{cov}\!\left[\epsilon _{ij}, \epsilon _{kl}\right] = \phi _3, \end{align*}

where the indices $i,j,k,$ and $l$ are distinct. In addition, we let $\textbf{p} = [p_1, p_2, p_3]$ parametrize the precision matrix $\boldsymbol{\Omega }^{-1} = \sum _{i=1}^3 p_i \mathcal{S}_i$ , which has the same form as the covariance matrix $\boldsymbol{\Omega }$ (see Marrs et al. (Reference Marrs, Fosdick and McCormick2017) for a similar result when $\{ \epsilon _{ij} \}_{ij}$ forms a directed network). The objective function, incorporating the restrictions that $\phi _1 = 1$ and $\phi _3 = 0$ , is

\begin{align*} Q_{\textbf{y}}(\boldsymbol{\phi }) \;:\!=\; E[\ell _{\textbf{z}} \, | \, \textbf{y} ] + \frac{1}{2}\lambda _1 (\phi _1 - 1) + \frac{1}{2}\lambda _3 \phi _3, \end{align*}

where $\boldsymbol{\phi } = [\phi _1, \phi _2, \phi _3]$ and the “ $\frac{1}{2}$ ” factors are included to simplify algebra. Then, differentiating $Q_{\textbf{y}}$ with respect to $\textbf{p}$ , $\lambda _1$ , and $\lambda _3$ , the estimators for $\rho$ , $\{ \lambda _1, \lambda _3 \}$ are

(12) \begin{align}{\widehat{\rho }} &= \gamma _2 - \frac{1}{|\Theta _2|} \begin{bmatrix} \frac{\partial \phi _1}{\partial p_2} & \frac{\partial \phi _3}{\partial p_2} \end{bmatrix}^T \begin{bmatrix} \lambda _1 \\ \lambda _3 \end{bmatrix} \quad \end{align}

(13) \begin{align} \begin{bmatrix} \widehat{\lambda }_1 \\ \widehat{\lambda }_3 \end{bmatrix} &= \begin{bmatrix} \frac{\partial \phi _1}{\partial p_1} & \frac{\partial \phi _3}{\partial p_1} \\ \frac{\partial \phi _1}{\partial p_3} & \frac{\partial \phi _3}{\partial p_3} \end{bmatrix}^{-1} \begin{bmatrix} |\Theta _1| & 0 \\ 0 & |\Theta _3| \end{bmatrix} \begin{bmatrix} \gamma _1 - 1 \\ \gamma _3 \end{bmatrix}, \quad \end{align}

where again $ \Theta _i$ is the set of pairs of relations $(jk, lm)$ that share an actor in the $i^{\textrm{th}}$ manner, for $i \in \{1,2,3 \}$ . For instance, $\Theta _2$ consists of pairs of relations of the form $(jk, jl)$ , where $j,k,$ and $l$ are distinct indices. In (12) and (13), the partial derivatives $\left \{ \partial \phi _i/ \partial p_j \right \}$ are available in closed form and are easily computable in $O(1)$ time using the forms of $\boldsymbol{\Omega }$ and $\boldsymbol{\Omega }^{-1}$ . See Appendix B for details.

Alternation of the estimators for $\rho$ and $\{ \lambda _1, \lambda _3 \}$ in (12) and (13) constitutes a block coordinate descent for $\rho = \phi _2$ subject to the constraints $\phi _1 = 1$ and $\phi _3 = 0$ . This block coordinate descent makes up the maximization step of the EMM algorithm for $\rho$ .

5.3 Maximization with respect to $\beta$

The maximization step with respect to $\boldsymbol{\beta }$ in the EMM algorithm can be obtained directly. Setting the derivative of (11) with respect to $\boldsymbol{\beta }$ equal to zero, the maximization step for $\boldsymbol{\beta }$ is

(14) \begin{align} \widehat{{\boldsymbol{\beta }}}^{(\nu +1)} & = \widehat{{\boldsymbol{\beta }}}^{(\nu )} + \left ( \textbf{X}^T \boldsymbol{\Omega }^{-1} \textbf{X} \right )^{-1} \textbf{X}^T \boldsymbol{\Omega }^{-1} E\!\left[\boldsymbol{\epsilon } \, | \, \textbf{y},{\widehat{\rho }}^{(\nu )}, \widehat{{\boldsymbol{\beta }}}^{(\nu )} \right], \quad \end{align}

where we use the identity $\boldsymbol{\epsilon } = \textbf{z} - \textbf{X}{\boldsymbol{\beta }}$ . In Appendix B, we show that the leading terms of the unique entries in $\Omega ^{-1}$ , $\textbf{p}$ , depend only on $\rho$ through a multiplicative factor,

\begin{align*} \textbf{p} \approx f(\rho )[g_1(n),\, g_2(n), \,g_3(n)]^T. \end{align*}

Thus, we may factor $f(\rho )$ out of (14), and the $\boldsymbol{\beta }$ maximization is asymptotically $\rho$ -free (except for the expectation term). Similarly, the maximization with respect to $\rho$ in Section 5.2 is free from $\boldsymbol{\beta }$ except for the expectation term. Hence, only a single maximization step with respect to each $\boldsymbol{\beta }$ and $\rho$ is required for each expectation.

5.4 Consistency of the EMM estimator

The complete multivariate normal likelihood for $\textbf{z}$ is a non-curved, identifiable likelihood. Then, it is known that each expectation and maximization step in an EM algorithm increases the current likelihood value (Wu, Reference Wu1983). Whenever there is a unique, single local maximum, the EM algorithm yields consistent, and efficient, estimators. We make a series of approximations to the expectations in the EMM algorithm to reduce computational demands, so that the theory of EM estimator convergence may not be directly applicable. Yet, we find that the EMM estimators, $\widehat{\beta }_{\textrm{EMM}}, \widehat{\rho }_{\textrm{EMM}} \}$ , maintain consistency. Taking the leading terms of ${\widehat{\rho }}_{\textrm{EMM}}$ ,

\begin{align*}{\widehat{\rho }}_{\textrm{EMM}} &= \frac{1}{2} + \frac{1}{n^3} \sum _{jk, lm \in \Theta _2} E\!\left[\epsilon _{jk} \epsilon _{lm} \mid y_{jk}, y_{lm} \right] - \frac{1}{n^2} \sum _{jk} E\!\left[\epsilon _{jk}^2 \mid y_{jk} \right]\ldots \\ & \ldots - \frac{2}{n^4} \sum _{jk, lm \in \Theta _3} E\!\left[\epsilon _{jk} \mid y_{jk}\right] E\!\left[\epsilon _{lm} \mid y_{lm}\right] + O(n^{-1}). \end{align*}

which has expectation $E[{\widehat{\rho }}_{\textrm{EMM}}] = \rho + O(n^{-1})$ . Then, consistency can be established by showing that the variance of ${\widehat{\rho }}_{\textrm{EMM}}$ tends to zero. We provide details in Appendix C. The estimator $\widehat{\beta }_{\textrm{EMM}}$ is particularly difficult to analyze, as $E[\epsilon _{jk} \mid \textbf{y}]$ depends on every entry in $\textbf{y}$ , and because we approximate this expectation. We provide a sketch for a proof of consistency in Appendix C by bounding the distance between the EMM estimator for $\beta$ and the true MLE, $\left|\!\left|\widehat{\beta }_{\textrm{EMM}} - \widehat{\beta }_{\textrm{MLE}} \right|\right|_2^2$ , using an easier-to-analyze estimator for $\beta$ which replaces $E[\epsilon _{jk} \mid \textbf{y}]$ with $E[\epsilon _{jk} \mid y_{jk}]$ . As in the argument for consistency of ${\widehat{\rho }}_{\textrm{EMM}}$ , we establish consistency of the bounding estimator by showing the expectation is asymptotically equal to the true value of $\beta$ and that the variance of the bounding estimator tends to zero. We also discuss performance of $\widehat{\beta }_{\textrm{EMM}}$ under model misspecification, showing that it maintains consistency even under violation of the normality and exchangeability assumptions.

6. Prediction

In this section, we describe how to use the PX model, and the approximations in service of Algorithm 1, to make predictions for an unobserved network relation without undue computational cost. The predicted value we seek is the probability of observing $y_{jk} = 1$ given all the other values $\textbf{y}_{-jk}$ , where $\textbf{y}_{-jk}$ is the vector of observations $\textbf{y}$ excluding the single relation $jk$ . As in estimation, the desired probability is again equal to a $\binom{n}{2}$ -dimensional multivariate truncated normal integral, which is computationally burdensome. Thus, we approximate the desired prediction probability

(15) \begin{align} \mathbb{P}( y_{jk} = 1 \, | \, \textbf{y}_{-jk} ) &= E \!\left [E \!\left [ \textbf{1}\!\left [\epsilon _{jk} \gt -{\textbf{x}}^T_{jk}{\boldsymbol{\beta }}\right ] \, | \,{\boldsymbol{\epsilon }}_{-jk} \right ] \, | \, \textbf{y}_{-jk} \right ], \quad \\ &\approx \Phi \!\left (\frac{ E[{\epsilon }_{jk} \, | \, \textbf{y} ] +{\textbf{x}}^T_{jk}{\boldsymbol{\beta }}}{ \sigma _n} \right ).\nonumber \end{align}

The approximation in (15) is based on the fact that $[\epsilon _{jk} \, | \,{\boldsymbol{\epsilon }}_{-jk}]$ is normally distributed:

(16) \begin{align} &\epsilon _{jk} \, | \,{\boldsymbol{\epsilon }}_{-jk} \sim \textrm{N}\!\left (m_{jk}, \sigma ^2_n\right ), \nonumber \\ m_{jk} &= -{\sigma _n^2} \textbf{1}_{jk}^T \!\left ( p_2 \mathcal{S}_2 + p_3 \mathcal{S}_3 \right )\widetilde{\boldsymbol{\epsilon }}_{-jk},\quad \sigma ^2_n = \frac{1}{p_1 }, \end{align}

where $\textbf{1}_{jk}$ is the vector of all zeros with a one in the position corresponding to relation $jk$ and, for notational simplicity, we define $\widetilde{\boldsymbol{\epsilon }}_{-jk}$ is the vector $\boldsymbol{\epsilon }$ with a zero in the entry corresponding to relation $jk$ . We note that the diagonal of the matrix $p_2 \mathcal{S}_2 + p_3 \mathcal{S}_3$ consists of all zeros so that $m_{jk}$ is free of $\epsilon _{jk}$ . Then, the inner expectation in (15) is

(17) \begin{align} E \!\left [ \textbf{1}[\epsilon _{jk} \gt -{\textbf{x}}^T_{jk}{\boldsymbol{\beta }} ] \, | \,{\boldsymbol{\epsilon }}_{-jk} \right ] = \Phi \!\left ( \frac{m_{jk} +{\textbf{x}}^T_{jk}{\boldsymbol{\beta }} }{\sigma _n}\right ). \quad \end{align}

Of course, $m_{jk}$ depends on ${\boldsymbol{\epsilon }}_{-jk}$ which is unknown, and thus, we replace $m_{jk}$ with its conditional expectation $E[m_{jk} \, | \, \textbf{y}_{-jk}] = E[\epsilon _{jk} \, | \, \textbf{y}_{-jk}]$ .

Computing $E[\epsilon _{jk} \, | \, \textbf{y}_{-jk}]$ is extremely difficult; however, computing $E[\epsilon _{jk} \, | \, \textbf{y}]$ proves feasible if we exploit the structure of $\boldsymbol{\Omega }$ . Thus, we approximate the desired expectation by imputing $y_{jk}$ with the mode of the observed data:

(18) \begin{align} E \!\left [ \epsilon _{jk} \, | \, \textbf{y}_{-jk} \right ] \approx E \!\left [ \epsilon _{jk} \, | \, \textbf{y}_{-jk}, y_{jk} = y^* \right ] = E \!\left [ \epsilon _{jk} \, | \, \textbf{y} \right ], \quad \end{align}

where $y^*$ is the mode of $\textbf{y}_{-jk}$ . The error due to this approximation is small and shrinks as $n$ grows. Substituting (18) for $m_{jk}$ in (17) gives the final expression in (15).

7. Simulation studies

In this section, we describe three simulation studies. The first verifies that the performance of the EMM estimator in Algorithm 1 provides improvement over standard probit regression. The second simulation study verifies consistency of the EMM estimators of $\boldsymbol{\beta }$ , and compares the performance of these estimators to the estimators of $\boldsymbol{\beta }$ from the social relations model and the latent eigenmodel. The third simulation study evaluates the robustness of the PX model, and EMM algorithm, to the assumption that the latent random variables are normally distributed.

For both simulation studies, we generated data with mean consisting of three covariates and an intercept:

(19) \begin{align} y_{ij} = \textbf{1} \Big [\beta _0 + \beta _1 \textbf{1}[x_{1i} \in C] \textbf{1}[x_{1j} \in C]+\beta _2|x_{2i}-x_{2j}|+\beta _3 x_{3ij} + \epsilon _{ij} \gt 0 \Big ]. \quad \end{align}

In the model in (19), $\beta _0$ is an intercept; $\beta _1$ is a coefficient on a binary indicator of whether individuals $i$ and $j$ both belong to a pre-specified class $C$ ; $\beta _2$ is a coefficient on the absolute difference of a continuous, actor-specific covariate $x_{2i}$ ; and $\beta _3$ is that for a pair-specific continuous covariate $x_{3ij}$ . We fixed ${\boldsymbol{\beta }} = [\beta _0, \beta _1, \beta _2, \beta _3]^T$ at a single set of values. Since the accuracy of estimators of $\boldsymbol{\beta }$ may depend on $\textbf{X}$ , we generated $20$ random design matrices $\textbf{X}$ for each sample size of $n \in \{ 20, 40, 80 \}$ actors. We emphasize that, although these may appear to be only moderately sized networks, each consists of $\binom{n}{2} \in \{ 190, \, 780, \, 3160\}$ observations. For each design matrix, we simulated 100 error realizations of $\{ \epsilon _{ij} \}_{ij}$ , with distribution that depended on the generating model. When generating from the PX model, half of the total variance in $\epsilon _{ij}$ was due to correlation $\rho = 1/4$ , and the remaining half was due to the unit variance of $\epsilon _{ij}$ . When generating from the latent eigenmodel in (5), one-third of the variance in $\epsilon _{ij}$ was due to each term $a_i + a_j$ , $\textbf{u}_i^T \Lambda \textbf{u}_j$ , and ${\textbf{x}}i_{ij}$ , respectively. For additional details of the simulation study procedures, see Section D.1.

7.1 Evaluation of approximations in Algorithm 1

To evaluate the efficacy of the approximations described in the estimation procedure in Algorithm 1, we simulated from (19) for a single $\textbf{X}$ with $n=40$ (larger $n$ caused multivariate normal integral failures in R). We simulated 100 networks from the PX model in (9) using this $\textbf{X}$ , for each value of $\rho \in \{0.1, 0.25, 0.4 \}$ (we note that we require $\rho \lt 1/2$ for the error covariance matrix $\Omega$ to be positive definite). For each realization, we estimated $\beta$ in the PX model using EMM in Algorithm 1. To estimate $\beta$ in the standard probit model, we used the function glm in R. To compute the MLE, we numerically optimized the data log-likelihood using the Broyden-Fletcher-Goldfarb-Shanno algorithm as implemented in the optim function in R, initializing at the true values of $\{ \beta, \rho \}$ .

In the left panel of Figure 1, we evaluate the performance of the EMM estimator by comparing the root mean square error (RMSE) between the EMM coefficient estimate, $\widehat{\beta }_{\textrm{EMM}}$ , and the MLE obtained by the optimization procedure $\widehat{\beta }_{\textrm{MLE}}$ . As a baseline, we compute the RMSE between $\widehat{\beta }_{\textrm{MLE}}$ and the true value $\beta$ . If the approximations in the EMM algorithm are small, we expect the RMSE between $\widehat{\beta }_{\textrm{EMM}}$ and $\widehat{\beta }_{\textrm{MLE}}$ to be much smaller than the RMSE between $\widehat{\beta }_{\textrm{MLE}}$ and $\beta$ . Generally, the RMSE between $\widehat{\beta }_{\textrm{EMM}}$ and $\widehat{\beta }_{\textrm{MLE}}$ is smaller than the RMSE between $\widehat{\beta }_{\textrm{MLE}}$ and $\beta$ . However, the discrepancy between the two RMSEs decreases as the true $\rho$ grows. As a reference, the MSE between $\widehat{\beta }_{\textrm{Std. probit}}$ and $\widehat{\beta }_{\textrm{MLE}}$ is also shown in the left panel of Figure 1; the EMM estimator is closer to $\widehat{\beta }_{\textrm{MLE}}$ than the standard probit estimator is to $\widehat{\beta }_{\textrm{MLE}}$ for all values of $\rho$ . Raw RMSE values between the estimators and the truth, shown in Figure 2, confirm that the EMM algorithm does perform better than standard probit in RMSE with respect to estimation of $\beta$ . The results of this simulation study suggest that the EMM algorithm improves estimation of $\beta$ over the standard probit estimator for $\rho \gt 0$ and that the EMM estimator is reasonably close to the MLE, signifying the approximations in the EMM algorithm are reasonable.

Figure 1. The left panel depicts performance in estimating $\beta$ : RMSE between the EMM estimator and the MLE (RMSE $(\widehat{\beta }_{\textrm{MLE}} - \widehat{\beta }_{\textrm{EMM}})$ ), between the MLE and the truth (RMSE $(\widehat{\beta }_{\textrm{MLE}} -{\beta })$ ), and between the MLE and the standard probit estimator (RMSE $(\widehat{\beta }_{\textrm{MLE}} - \widehat{\beta }_{\textrm{Std. probit}})$ ). The right panel depicts performance in estimating $\rho$ : RMSE between the MLE and the EMM estimator (RMSE $({\widehat{\rho }}_{\textrm{MLE}} -{\widehat{\rho }}_{\textrm{EMM}})$ ) and between the MLE and the truth (RMSE $({\widehat{\rho }}_{\textrm{MLE}} - \rho )$ ). The RMSEs are plotted as a function of the true values of $\rho$ , and solid vertical lines denote Monte Carlo error bars. Some points obscure their Monte Carlo error bars.

Figure 2. The left panel depicts the RMSE in estimating $\beta$ using the EMM algorithm, MLE, and standard probit regression. The right panel depicts the same for $\rho$ . The MSEs are plotted as a function of the true values of $\rho$ , and solid vertical lines denote Monte Carlo error bars.

In the right panel of Figure 1, the EMM estimator of $\rho$ is closer to the MLE, ${\widehat{\rho }}_{\textrm{MLE}}$ , than the MLE is close to the true value of $\rho$ for all values of $\rho$ examined. This fact suggests that the approximation error in estimating $\rho$ in the EMM algorithm is small. Further, the raw RMSE values shown in Figure 2 illustrate that ${\widehat{\rho }}_{\textrm{EMM}}$ may be as good an estimator of $\rho$ as is ${\widehat{\rho }}_{\textrm{MLE}}$ . The approximations in ${\widehat{\rho }}_{\textrm{EMM}}$ appear to be stable over the range of $\rho$ values examined. Overall, since the degradation in performance of the EMM algorithm is most pronounced in estimation of $\beta$ , we postulate that the degradation may be due to the approximations in computing $E[\epsilon _{jk} \mid \textbf{y}]$ (see Section A.2).

7.2 Performance in estimation of $\beta$

To evaluate the performance of the PX estimator in estimating linear coefficients $\boldsymbol{\beta }$ , we compared estimates of $\boldsymbol{\beta }$ by the EMM algorithm to estimators of the social relations and latent eigenmodels on data generated from the PX model and data generated from the latent eigenmodel. We used the amen package in R to estimate the social relations model and latent eigenmodel (Hoff et al., Reference Hoff, Fosdick, Volfovsky and He2017). We again compared these estimators to the standard probit regression model assuming independence as a baseline, which we estimated using the function glm in R. We focused on the value of $\rho =0.25$ , in the center of the range of possible $\rho$ values.

In Figure 3, we plot the RMSE (scaled by $n^{1/2}$ ) of the $\beta$ coefficients estimated for the PX model, standard probit model, social relations model, and latent eigenmodels. We see that the EMM estimator for the PX model has a downward trend in $n^{1/2}$ RMSE with $n$ , and a reducing spread of $n^{1/2}$ RMSE with $n$ , for both the PX and latent eigenmodel generating models. These facts suggest that the PX estimator is consistent for $\boldsymbol{\beta }$ , at a rate $n^{1/2}$ or better, for both the PX and latent eigenmodel generating models, confirming the claims in Section 5.4. Further, the EMM estimator has the lowest median $n^{1/2}$ RMSE of any of the estimators for all entries in $\boldsymbol{\beta }$ , where $n^{1/2}$ RMSE is evaluated for each $\textbf{X}$ realization (across the error realizations), and the median is computed across the 20 $\textbf{X}$ realizations. We observe similar patterns for the correlation parameter $\rho$ ; see Section D.1. Interestingly, the superiority of the PX estimator holds whether we generate from the PX or latent eigenmodel, which suggests that any benefit in correctly specifying the latent eigenmodel is lost in the estimating routine. The larger $n^{1/2}$ RMSEs of the amen estimator of the social relations and latent eigenmodels are a result of bias; see Section D.1 for bias-variance decomposition of the MSEs.

Figure 3. Performance ( $n^{1/2}$ RMSE) of estimators of $\beta$ , for a given $\textbf{X}$ , when generating from the PX model (top row) and the latent eigenmodel (LE; bottom row). Variability captured by the boxplots reflects variation in RMSE with $\textbf{X}$ .

7.3 Runtimes

We evaluated the average runtimes of the algorithms used to estimate the simulated data. The average runtimes are plotted in Figure 4. The improvement in runtime offered by the EMM estimation scheme over SRM and LE MCMC estimation is several orders of magnitude. Interestingly, the runtime cost of EMM appears to grow faster than the MCMC routines and faster than standard probit regression. A contributing factor is the sum over $O(n^3)$ terms in the maximization of $\rho$ in the EMM algorithm. We have experimented with using only a random subset of $O(n^2)$ relation pairs in the maximization step, which results in gains in runtime with small cost in estimation performance. Such a tradeoff may become attractive for networks of sufficient size $n$ .

Figure 4. Average runtimes of various algorithms used on simulated data.

7.4 Evaluation of latent normality assumption

To evaluate the performance of the PX model under violation of the normality assumption on the latent errors ${\{ \epsilon _{ij}\}}_{ij}$ , we repeated the simulation study with $t$ -distributed latent random variables. Specifically, we simulated from (19), but replaced the latent error vector $\boldsymbol{\epsilon }$ with $\sigma ^{-1}$ $\Omega ^{1/2}{\textbf{u}}$ , where $\textbf{u}$ consists of independently and identically distributed $t$ random variables with 5 degrees of freedom. The scaling factor $\sigma = \sqrt{5/3}$ ensures that $\textbf{u}$ has unit population variance, for consistency with the Gaussian case, and $\Omega ^{1/2}$ is the matrix square root of $\Omega$ , with $\rho =0.25$ . This model thus has the same latent mean and covariance matrix as in the Gaussian case, but the latent errors have substantially heavier tails.

The left panel of Figure 5 shows the performance of the EMM algorithm in estimating ${\boldsymbol{\beta }}_1$ , compared to the standard probit regression estimates. As in the Gaussian case, the EMM algorithm produces estimates with $n^{1/2}$ RMSE tending to zero as $n$ grows. Also as in the Gaussian case, EMM estimation of the PX model improves estimation of $\boldsymbol{\beta }$ over standard probit regression. We observed the same results in estimation of the remaining coefficients (see Section D.2). Unlike the Gaussian case, $n^{1/2}$ RMSE in estimating $\rho$ did not appear to tend towards zero. However, in the right panel of Figure 5, the error in estimating $\rho$ scaled by $n^{1/4}$ , $n^{1/4}$ RMSE, does tend towards zero. This study confirms the claim in Section 5.4 that the EMM algorithm produces consistent estimators $\widehat{{\boldsymbol{\beta }}},{\widehat{\rho }} \}$ , even under violation of the normality assumption of the PX model.

Figure 5. The left panel depicts $n^{1/2}$ RMSE in estimating $\beta _1$ , using the EMM algorithm and standard probit regression, under $t$ distribution of the errors. The right panel depicts $n^{1/4}$ RMSE in estimating $\rho$ using the EMM algorithm in the same simulation. Variability captured by the boxplots reflects variation in RMSE with $\textbf{X}$ .

8. Analysis of a network of political books

We live in a time of political polarization. We investigate this phenomenon by analyzing a network of $n=105$ books on American politics published around the time of the 2004 presidential election.Footnote ² These data were compiled by Dr. Valdis Krebs using the “customers who bought this book also bought these books” list on Amazon.com. At the time, when browsing a particular book, Amazon listed the books that were bought by individuals who also bought the book in question. Thus, a relation between two books in the network indicates that they were frequently purchased by the same buyer on Amazon. Political books on the best-seller list of The New York Times were used as actors in the network. Finally, the books were labeled as conservative, liberal, or neutral based on each book’s description (Figure 6). Work by Dr. Krebs on a similar network was in a 2004 New York Times article (Eakin, Reference Eakin2004), where it was shown that there were many relations between books with similar ideologies yet relatively few across ideologies. The work by Dr. Krebs has inspired similar analyses of book purchasing networks in the fields of nanotechnology (Schummer, Reference Schummer2005) and climate science (Shi et al., Reference Shi, Shi, Dokshin, Evans and Macy2017).

Figure 6. Krebs’ political books network (left) and out-of-sample performance in 10-fold cross-validation, as measured by area under the precision-recall curve (PRAUC, right), plotted against mean runtime in the cross-validation. The estimators are standard probit assuming independent observations (Std. probit), the PX model as estimated by the EMM algorithm (PX), the social relations model estimator (SRM), and the latent eigenmodel estimator (LE).

To confirm previous work by Dr. Krebs, we develop a model that assigns a different probability of edge formation between books $i$ and $j$ depending on whether the books are ideologically aligned. By examining the network in Figure 6, we observe that neutral books appear to have ties than books that are labeled conservative or liberal. Thus, we add a nodal effect indicating whether either book in a relation is labeled neutral. The regression model specified is

(20) \begin{align}{\mathbb{P}}(y_{ij}=1) &={\mathbb{P}} (\beta _0 + \beta _1 \textbf{1}[c(i) = c(j)] \nonumber \\ &\hspace{.25in}+ \beta _2\textbf{1} \!\left [ \!\left \{c(i) = \textrm{neutral} \right \} \, \cup \, \!\left \{ c(j) = \textrm{neutral} \right \} \right ] + \boldsymbol{\epsilon }_{ij} \gt 0 ), \quad \\ \nonumber \boldsymbol{\epsilon } &\sim (\textbf{0}, \boldsymbol{\Sigma }), \end{align}

where $c(i)$ represents the class of book $i$ (neutral, conservative, or liberal), and the distribution and covariance matrix of $\boldsymbol{\epsilon }$ are determined by the particular model being estimated. In this section, we estimate the PX model (PX), the equivalent social relations model (SRM), the latent eigenmodel (LE), and, as a baseline, the standard probit regression model assuming independence of observations (which we label “std. probit”).

We used a 10-fold cross-validation to compare the out-of-sample predictive performance of the estimators and the runtimes of the algorithms for the models in question. We used the proposed EMM algorithm to estimate the PX model, the amen package in R to estimate the social relations model and latent eigenmodel (Hoff et al., Reference Hoff, Fosdick, Volfovsky and He2017), and the glm(.) command in the R package stats to estimate the standard probit model. We randomly divided the $\binom{105}{2}$ relations into 10 disjoint sets, termed “folds,” of roughly the same size. Then, for each fold, we estimated the models on the remaining nine folds and made predictions for the data in the fold that was not used for estimation (for details of estimation of the PX model with missing data, see Section A.4). Repeating this operation for each complete data set of out-of-sample predictions for each estimating model. The procedure to make marginal predictions from the PX model is described in Section 6. To compare with the PX model, we make marginal predictions from the social relations model and the latent eigenmodel, that is, by integrating over the random effect space. The predictions from the social relations model and the latent eigenmodel are automatically output from amen in the presence of missing data. The predictions from the standard probit model are marginal by default as there is no correlation structure.

Table 1. Results of fitting the Krebs political books data using the EMM estimator for the PX model and the amen estimator for the social relations and latent eigenmodels (SRM and LE, respectively). Point estimates for the coefficients are given to the left of the vertical bar, and runtimes (in seconds) and minimum effective sample sizes across the coefficient estimates are given to the right

We use area under the precision-recall curve (PRAUC) to measure performance of the predictions relative to the observed data, although using area under the receiver operating characteristic (ROC) yields the same conclusions (see Appendix E). In Figure 6, the proposed EMM estimator produces an improvement in PRAUC over standard probit prediction that is roughly equivalent to the improvement of the social relations model over standard probit, yet with an average runtime that is 45 times faster (about a minute compared with an hour). The latent eigenmodel produces an improvement in PRAUC over the proposed EMM algorithm and the social relations model, however, at the expense of significant increase in average runtime, that of about 3,000 times slower than EMM and taking almost three days to complete. Note that we selected the number of MCMC iterations for the social relations and latent eigenmodels that resulted in sets of samples from the posterior distributions (after burn-in) that had a effective sample sizes roughly equal to 100 independent samples of the $\boldsymbol{\beta }$ parameters. Increasing the number of iterations, which may be desirable, would result in even longer runtimes for the estimators of the social relations and latent eigenmodels. Taken together, the results of the cross-validation study suggest that the PX model accounts for a large portion of the correlation in network data with estimation runtime that, depending upon stopping criterion, is orders of magnitude faster the runtime than existing approaches.

To estimate the complete data set under the mean model in (20), we used the EMM algorithm for the PX model and the amen package for the social relations model (SRM) and latent eigenmodel (LE), which we ran for $1\times 10^6$ iterations after a burn-in of $5 \times 10^4$ iterations (with runtimes of roughly two hours for SRM and 17 hours for LE). The coefficient estimates in Table 1 suggest that books that share the same ideology are more likely to be frequently purchased together, as all $\widehat{\beta }_1 \gt 0$ . This positive coefficient estimate demonstrates political polarization in the network: conservative books are more likely to be purchased with other conservative books rather than with liberal books. The second coefficient estimate, $\widehat{\beta }_2 \gt 0$ , suggests that, relative to a random pair of ideologically misaligned books, pairs of books where at least one of the books is neutral are more likely to be purchased together. Neutral books are thus generally more likely to be purchased with books of disparate ideologies and have a unifying effect in the book network. Returning briefly to Table 1, the runtimes highlight that EMM reduces computational burden by order(s) of magnitude over existing approaches.