Regression of binary network data with exchangeable latent errors

Abstract

Undirected, binary network data consist of indicators of symmetric relations between pairs of actors. Regression models of such data allow for the estimation of effects of exogenous covariates on the network and for prediction of unobserved data. Ideally, estimators of the regression parameters should account for the inherent dependencies among relations in the network that involve the same actor. To account for such dependencies, researchers have developed a host of latent variable network models; however, estimation of many latent variable network models is computationally onerous and which model is best to base inference upon may not be clear. We propose the probit exchangeable (PX) model for undirected binary network data that is based on an assumption of exchangeability, which is common to many of the latent variable network models in the literature. The PX model can represent the first two moments of any exchangeable network model. We leverage the EM algorithm to obtain an approximate maximum likelihood estimator of the PX model that is extremely computationally efficient. Using simulation studies, we demonstrate the improvement in estimation of regression coefficients of the proposed model over existing latent variable network models. In an analysis of purchases of politically aligned books, we demonstrate political polarization in purchase behavior and show that the proposed estimator significantly reduces runtime relative to estimators of latent variable network models, while maintaining predictive performance.

1. Introduction

Undirected binary network data measure the presence or absence of a relationship between pairs of actors and have recently become extremely common in the social and biological sciences. Some examples of data that are naturally represented as undirected binary networks are international relations among countries (Fagiolo et al., 2008), gene co-expression (Zhang & Horvath, 2005), and interactions among students (Han et al., 2016). We focus on an example of politically aligned books, where a relation exists between two books if they were frequently purchased by the same person on Amazon.com. Our motivations are estimation of the effects of exogenous covariates, such as the effect of alignment of political ideologies of pairs of books on the propensity for books to be purchased by the same consumer, and the related problem of predicting unobserved relations using book ideological information. For example, predictions of relations between new books and old books could be used to recommend new books to potential purchasers.

A binary, undirected network $\left \{ y_{ij} \in \{0,1 \}: i,j\in \{1,\ldots,n\}, i \lt j \right \}$ , which we abbreviate $\{y_{ij} \}_{ij}$ , may be represented as an $n \times n$ symmetric adjacency matrix which describes the presence or absence of relationships between unordered pairs of $n$ actors. The diagonal elements of the matrix $\{y_{ii}\;:\; i \in \{1,\ldots,n\}\}$ are assumed to be undefined, as we do not consider actor relations with him/herself. We use $\textbf{y}$ to refer to the $\binom{n}{2}$ vector of network relations formed by a columnwise vectorization of the upper triangle of the matrix corresponding to $\{y_{ij} \}_{ij}$ .

A regression model for the probability of observing a binary outcome is the probit model, which can be expressed

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

where $\epsilon _{ij}$ is a mean-zero normal random error, ${\textbf{x}}_{ij}$ is a fixed vector of covariates corresponding to relation $ij$ , and $\boldsymbol{\beta }$ is a vector of coefficients to be estimated. When each entry in the error network $\{ \epsilon _{ij} \}_{ij}$ is independent of the others, estimation of the probit regression model in (1) is straightforward and proceeds via standard gradient methods for maximum likelihood estimation of generalized linear models (Greene, 2003). The assumption of independence of $\{ \epsilon _{ij} \}_{ij}$ may be appropriate when the mean $\big\{{\textbf{x}}_{ij}^T{\boldsymbol{\beta }} \big\}_{ij}$ represents nearly all of the dependence in the network $\{ y_{ij}\}_{ij}$ . However, network data naturally contain excess dependence beyond the mean: the errors $\epsilon _{ij}$ and $\epsilon _{ik}$ both concern actor $i$ (see Faust & Wasserman, 1994, e.g., for further discussion of dependencies in network data). In the context of the political books data set, the propensity of “Who’s Looking Out For You?” by Bill O’Reilly to be purchased by the same reader as “Deliver Us from Evil” by Sean Hannity may be similar to the propensity of “Who’s Looking Out For You?” and “My Life” by Bill Clinton to be co-purchased simply because “Who’s Looking Out For You?” is a popular book. Or, in a student friendship network, the friendship that Julie makes with Steven may be related to the friendship that Julie makes with Asa due to Julie’s gregariousness. Unlike the case of typical linear regression, the estimator that maximizes the likelihood of the generalized linear regression model in (1), when assuming independence of each entry in the error network $\{ \epsilon _{ij} \}_{ij}$ , is not unbiased for $\boldsymbol{\beta }$ . Ignoring the excess dependence in $\{\epsilon _{ij} \}_{ij}$ can thus be expected to result in poor estimation of $\boldsymbol{\beta }$ and poor out-of-sample predictive performance. We observe this phenomenon in the simulation studies and analysis of the political books network (see Sections 7 and 8, respectively). Thus, estimators of $\boldsymbol{\beta }$ and ${\mathbb{P}}(y_{ij} = 1)$ in (1) for the network $\{ y_{ij} \}_{ij}$ should ideally account for the excess dependence of network data. A host of regression models exist in the literature that do just this; we briefly review these here.

A method used to account for excess dependence in regression of binary network data is the estimation of generalized linear mixed models, which were first introduced for repeated measures studies (Stiratelli et al., 1984; Breslow & Clayton, 1993). In these models, a random effect, that is, latent variable, is estimated for each individual in the study, to account for possible individual variation. Warner et al. (1979) used latent variables to account for excess network dependence when analyzing data with continuous measurements of relationships between actors, and Holland & Leinhardt (1981) extended their approach to networks consisting of binary observations. Hoff et al. (2002) further extended this approach to include nonlinear functions of latent variables, and since then, many variations have been proposed (Handcock et al., 2007; Hoff, 2008; Sewell & Chen, 2015). We refer to parametric network models wherein the observations are independent conditional on random latent variables as “latent variable network models,” which we discuss in detail in Section 2. Separate latent variable approaches may lead to vastly different estimates of $\boldsymbol{\beta }$ , and it may not be clear which model’s estimate of $\boldsymbol{\beta }$ , or prediction, to choose. Goodness-of-fit checks are the primary method of assessing latent variable network model fit (Hunter et al., 2008b); however, selecting informative statistics is a well-known challenge. Finally, latent variable network models are typically computationally burdensome to estimate, often relying on Markov chain Monte Carlo methods.

Another approach to estimating covariate effects on network outcomes is the estimation of exponential random graph models, known as ERGMs. ERGMs represent the probability of relation formation using a generalized exponential family distribution, ${\mathbb{P}}(y_{ij} = 1) \propto \textrm{exp}(\textbf{t}(\textbf{y}_{ij},{\textbf{x}}_{ij})^T \theta )$ , where $\theta$ is a vector of parameters to be estimated. In this flexible formulation, the effects of the exogenous covariates are included in the network statistics $\textbf{t}(\textbf{y}_{ij},{\textbf{x}}_{ij})$ . ERGMs also account for excess network dependence using the network statistics $\textbf{t}(\textbf{y}_{ij},{\textbf{x}}_{ij})$ , such as counts of the number of observed triangles or the number of “2-stars”—pairs of indicated relations that share an actor. ERGMs were developed by Frank & Strauss (1986) and Snijders et al. (2006) and are typically estimated using Markov chain Monte Carlo (MCMC) approximations to posterior distributions (Snijders, 2002; Handcock et al., 2019; Hunter et al., 2008a). ERGMs have been shown to be prone to place unrealistic quantities of probability mass on networks consisting of all “1”s or all “0”s (Handcock et al., 2003; Schweinberger, 2011), and the estimation procedures may be slow to complete (Caimo & Friel, 2011). Further, parameter estimates typically cannot be generalized to populations outside the observed network (Shalizi & Rinaldo, 2013).

A final approach to account for excess network dependence is to explicitly model the correlation among network observations. This is the approach we take in this paper. In this approach, an unobserved normal random variable, $z_{ij}$ , is proposed to underlie each data point, such that $y_{ij} = \textbf{1}[z_{ij} \gt 0]$ for $\textbf{z} \sim \textrm{N}(\textbf{X}{\boldsymbol{\beta }}, \boldsymbol{\Omega }(\boldsymbol{\theta }))$ . In this formulation, excess dependence due to the network is accounted for in $\boldsymbol{\Omega }$ . The parameters $\boldsymbol{\beta }$ and $\boldsymbol{\theta }$ of the distribution of the unobserved normal random variables $\{ z_{ij} \}_{ij}$ may be estimated using likelihood methods. For example, Ashford & Sowden (1970) propose likelihood ratio hypothesis tests and Ochi & Prentice (1984) give closed-form parameter estimators for studies of repeated observations on the same individual, such that $\boldsymbol{\Omega }(\boldsymbol{\theta })$ is block diagonal. In more general scenarios, such as unrestricted correlation structures, methods such as semi-parametrics (Connolly & Liang, 1988), pseudo-likelihoods (Le Cessie & Van Houwelingen, 1994), and MCMC approximations to EM algorithms (Chib & Greenberg, 1998; Li & Schafer, 2008) are employed for estimation.

In this paper, we propose the probit exchangeable (PX) model, a parsimonious regression model for undirected binary network data based on an assumption of exchangeability of the unobserved normal random variables $\{ z_{ij} \}_{ij}$ . The assumption of exchangeability is pervasive in random network models and, in fact, underlies many of the latent variable network models (see Section 3 for a detailed discussion of exchangeability). ¹ We show that, under exchangeability, the excess network dependence in $\{ z_{ij} \}_{ij}$ may be quantified using a single parameter $\rho$ such that $\boldsymbol{\Omega }(\boldsymbol{\theta }) = \boldsymbol{\Omega }(\rho )$ . This fact remains regardless of the particular exchangeable generating model, and thus, our approach can be seen as subsuming exchangeable latent network variable models, at least up to the second moment of their latent distributions. The proposed model may be rapidly estimated using an expectation-maximization (EM) algorithm to attain a numerical approximation to the maximum likelihood estimator, where we make approximations in the expectation step for runtime considerations. The estimation scheme we employ is similar to those used to estimate generalized linear mixed models in the literature (Littell et al., 2006; Gelman & Hill, 2006).

This paper is organized as follows. As latent variable network models are strongly related to our work, we review them in detail in Section 2. We provide supporting theory for exchangeable random network models and their connections to latent variable network models in Section 3. In Section 4, we define the PX model and then the estimation thereof in Section 5. In Section 6, we give a method for making predictions on unobserved relations. We provide simulation studies demonstrating consistency of the proposed estimation algorithm and demonstrating the improvement with the proposed model over latent variable network models in estimating $\boldsymbol{\beta }$ in Section 7. We analyze a network of political books in Section 8, demonstrating the reduction in runtime when PX model, and compare its out-of-sample performance to existing latent variable network models. A discussion with an eye toward future work is provided in Section 9.

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., 1979; Wong, 1982; Snijders & Kenny, 1999). In the social relations model for binary network data (Hoff, 2005), $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., 2002). 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. (2002) 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, 1994).

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, 2008):

(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 (2008) 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., 2018). Two models that are aimed at identifying subgroups of nodes that are structurally equivalent are the latent class model of Nowicki & Snijders (2001) and the mixed membership stochastic blockmodel (Airoldi et al., 2008). Hoff (2008) 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.

2.4 Drawbacks

The latent variable network models discussed in this section were developed based on the patterns often observed in real-world social networks. Latent variable network models contain different terms to represent the social phenomena underlying these patterns, and thus, different models may lead to substantially different estimates of $\boldsymbol{\beta }$ . It may not be clear which model’s estimate of $\boldsymbol{\beta }$ , or which model’s prediction of $\{ y_{ij} \}_{ij}$ , is best. Generally, latent variable network models are evaluated using goodness-of-fit checks (Hunter et al., 2008b), rather than rigorous tests, and it is well-known that selecting informative statistics for the goodness-of-fit checks is challenging. The latent variable network models described in this section are typically estimated using a Bayesian Markov chain Monte Carlo (MCMC) approach, which may be slow, especially for large data sets. Some recent advances do directly attempt to maximize the likelihood of network models with latent spaces (Ma et al., 2020; Zhang et al., 2022); however, public software implementations of these methods do not appear available, and they require certain covariate types (relation-level and actor-level, respectively) and certain latent space structures, such as the Euclidean distance latent space.

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, 1979; Aldous, 1981), 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, 2006; Kallenberg, 2006):

(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, 2006; Borgs et al., 2014). 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. (2017) 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.

Proposition 4.1. Suppose that the random effects $\{a_i \}_{i=1}^n$ for the social relations model in (3) are normally distributed. Then, there exists $\rho \in [0, 1/2]$ such that $\{y_{ij} \}_{ij}$ in the PX model in (9) is equal in distribution to $\{y_{ij} \}_{ij}$ as specified by the social relations model in (3).

Proof. As the PX and social relations models are probit regression models with the same mean structure, given by $\textbf{X}{\boldsymbol{\beta }}$ , it is sufficient to show that their latent covariance matrices are equivalent, that is, that $\textrm{var}[ \{a_i + a_j +{\textbf{x}}i_{ij} \}_{ij} ] = \textrm{var}[ \{ \epsilon _{ij} \}_{ij} ]$ . By exchangeability, the latent covariance matrices of the PX and social relations models have the same form and by assumption have variance 1. It is easy to see that, given $\sigma ^2_a \le 1$ (a necessary condition for $\textrm{var}[\epsilon _{ij}]=1$ ), we may take $\rho = \sigma ^2_a/2$ for the PX model, which establishes equality in the model distributions.

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., 2015), 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., 2006; Gelman & Hill, 2006). 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., 1977). 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, 1998), 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. (2017) 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, 1983). 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., 2017). 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. ² 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, 2004), 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, 2005) and climate science (Shi et al., 2017).

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., 2017), 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.

9. Discussion

In this paper, we present the PX model, a probit regression model for undirected, binary networks. The PX model adds a single parameter—latent correlation $\rho$ —to the ordinary probit regression model that assumes independence of observations. Our focus in this paper is estimation of the effects of exogenous covariates on the observed network, $\boldsymbol{\beta }$ , and prediction of unobserved network relations. Thus, we do not present uncertainty estimators for $\widehat{{\boldsymbol{\beta }}}$ or $\widehat{\rho }$ . However, practitioners estimating the PX model may require uncertainty estimators to perform inference. Development and evaluation of estimators of the uncertainty in estimators of network data is nontrivial; indeed, entire papers are dedicated to this task for the simpler linear regression case (see, for example, Aronow et al., 2015; Marrs et al., 2017). Future development of uncertainty estimators for the PX model may draw upon existing literature for uncertainty in EM estimators (Louis, 1982) and the numerical approximations in this paper.

A popular notion in the analysis of network data is the presence of higher-order dependencies, meaning beyond second order (Hoff, 2005). The representation of triadic closure, a form of transitivity—the friend of my friend is likely to also be my friend— is one motivation for the latent eigenmodel (Hoff, 2008). The PX model does represent triadic closure to a degree. One can show that, given two edges of a triangle relation exist, $y_{ij} = y_{jk} = 1$ , the probability that the third edge exists, ${\mathbb{P}}(y_{ik} = 1)$ , increases as $\rho$ increases. However, the increase in probability describing triadic closure under the PX model is fixed based on the estimated value of $\rho$ , which is informed only by the first two moments of the data when using the EMM estimator. It may be desirable to develop a test for whether the PX model sufficiently represents the level of triadic closure as suggested by the data. One such test might compute the empirical probability that ${\mathbb{P}}(y_{ik} = 1 \, | \, y_{ij} = y_{jk} = 1)$ and compare this statistic to its distribution under the null that the PX model is the true model with correlation parameter $\rho ={\widehat{\rho }}$ . Future work consists of theoretical development of the distributions of the test statistic(s) of choice under the null. Statistics of interest will likely be related to various clustering coefficients in the network literature (Wasserman & Faust, 1994; Watts & Strogatz, 1998).

We focus on the probit model in this paper. However, we find that this choice may limit the degree of covariance in the observed network $\{ y_{ij} \}_{ij}$ that the PX model can represent. For constant mean ${\textbf{x}}_{ij}^T{\boldsymbol{\beta }} = \mu$ , the maximum covariance the PX model can represent is bounded by

(21) \begin{align} \textrm{cov}[y_{ij}, y_{ik}] \, \le \, \lim _{\rho \rightarrow 1/2} \int _{-\mu }^\infty \int _{-\mu }^\infty dF_\rho - \Phi (\mu )^2, \quad \end{align}

where $dF_\rho$ is the bivariate standard normal distribution with correlation $\rho$ . The use of different latent distributions for $\boldsymbol{\epsilon }$ other than normal may allow a model analogous to the PX model to represent a larger range of observed covariances $\textrm{cov}[y_{ij}, y_{ik}]$ . Future work may consider a logistic distribution for $\boldsymbol{\epsilon }$ , as some researchers prefer to make inference with logistic regression models for binary data due to the ease of interpretation.