Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,117 characters · 12 sections · 68 citation commands
A Structural Model of Business Card Exchange Networks
\abstract{Social and professional networks affect labor market dynamics, knowledge diffusion and new business creation. To understand the determinants of how these networks are formed in the first place, we analyze a unique dataset of business card exchanges among a sample of over 240,000 users of the multi-platform contact management and professional social networking tool for individuals Eight. We develop a structural model of network formation with strategic interactions, and we estimate users' payoffs that depend on the composition of business relationships, as well as indirect business interactions. We allow heterogeneity of users in both observable and unobservable characteristics to affect how relationships form and are maintained. The model's stationary equilibrium delivers a likelihood that is a mixture of exponential random graph models that we can characterize in closed-form. We overcome several econometric and computational challenges in estimation, by exploiting a two-step estimation procedure, variational approximations and minorization-maximization methods. Our algorithm is scalable, highly parallelizable and makes efficient use of computer memory to allow estimation in massive networks. We show that users payoffs display homophily in several dimensions, e.g. location; furthermore, users unobservable characteristics also display homophily. }
Encounters are a seed of business. Thomas Edison and Henry Ford became friends after a chat at the convention of the Association of Edison Illuminating Companies in New York. Steve Jobs and Steve Wozniak met through a mutual friend, Bill Fernandez. Needless to say, these encounters eventually turned into great businesses, path-breaking innovations, and new products. Social and business networks that ultimately stem from such encounters create interesting economic phenomena, which has attracted many scholars to explore the subject. Indeed, even if they do not cause the birth of new businesses, professional networks play an important role in various economic activities at both individual and firm levels.\footnote{Professional networks provide information about vacancies and quality of job applicants through referrals in labor markets IoannidesLoury2004,CalvoArmengol2004, Calvo-ArmengolJackson2004, GaleottiMerlino2014, Galenianos2014 as documented in several empirical studies Beaman2012, BayerRossTopa2008. Firm-to-firm networks are suggested to catalyze aggregate fluctuations Acemoglu2012network, performance heterogeneity among firms Bernard2019production, agglomeration Miyauchi2021, and knowledge creation and diffusion Konig2018endogenous.} However, most of the existing studies presume that there are already more-or-less established relationships among firms or persons, and there are no many empirical studies about how these professional networks are formed in the first place.
In this paper, we employ a unique dataset with over $240,000$ individuals and the $670,000$ business connections among them to estimate a business network formation model, accounting for observable and unobservable individual characteristics that affect the willingness to form professional relationships. The data is a subset of the social network formed by users of Eight, a multi-platform contact management and professional social networking tool for individuals, provided by the Japanese company Sansan, Inc. In this social network, users connect with each other by exchanging business cards, which allows us to analyze a network of mostly face-to-face connections across a diverse spectrum of industries, occupations and locations in the whole Japan, on a scale that has not been used in previous work in network economics.
The Japanese labor market is an ideal setting to study the very beginning of business networks. Indeed, business cards are extensively used in Japan as a way of self-introduction, information sharing and establishing business relationships. One nature of business card exchanges is that when two persons exchange business cards, it is very likely that they are meeting each other for the first time. This aspect is different from other social networks such as Facebook and LinkedIn, where people are often acquainted with each other before they become connected on the platform. The fact that business card exchanges in many cases take place at the first meeting allows us to answer the question of how business networks emerge to begin with.
Our approach overcomes many estimation and empirical challenges posed by the scale of the network by carefully using the model's equilibrium implications as well as new and improved algorithms for estimation. We provide a theoretical framework for understanding face-to-face professional networking, where agents have observable and unobservable characteristics that affect their willingness to form professional connections. Their payoffs are also affected by link externalities such as popularity or common business connections. The equilibrium of the model provides the likelihood of observing a particular network of professional relationships at a particular point in time, that we use as likelihood of the data Mele2017, Mele2020, MeleZhu2021, BoucherMourifie2017. Unobserved heterogeneity is modeled as grouped random effects, thus providing a mixture model of network formation in equilibrium. We exploit a specification of the model with local externalities, inducing a likelihood that factorizes in between- and within-blocks contributions. This in turn reduces the computational challenge because links across unobserved types/blocks are conditionally independent.
To estimate the structural model with the massive Eight dataset, we develop a scalable two-step estimation algorithm, improving methods from VuEtAl2013 and BabkinEtAl2020 and including observable covariates. In the first step, we recover the unobservable heterogeneity by approximating the likelihood of the model with a stochastic blockmodel, thus abstracting from externalities within blocks. As shown in BabkinEtAl2020, this approximation works as long as the network is large and the number of unobservable agents' types (the size of the support for the random effect) is relatively large. We derive mean-field variational approximations for the likelihood of the model, and use an expectation-maximization algorithm to estimate the block structure WainwrightJordan2008, Bishop2006, BickelEtAl2013, BabkinEtAl2020. Furthermore, we use a minorization-maximization algorithm, which speeds up computation by several magnitudes with respect to standard maximization BabkinEtAl2020, VuEtAl2013. Our improved algorithm makes extensive use of the sparsity of the network and efficient sparse matrix algebra routines, as well as a scalable initialization algorithm Rosvall_2009, in order to further reduce the memory and time requirements for estimation. In the second step of the algorithm, we estimate the structural payoff parameters using a flexible pseudolikelihood estimator BoucherMourifie2017, conditioning on the estimated unobserved heterogeneity in the first step. This two-step procedure allows us to obtain reliable estimates in such a complex model using a large sample.
In the empirical implementation we allow parameters to be a function of between- and within-blocks memberships. We control for homophily in the location, industry and occupation of the users. Our results show that there is homophily in observable characteristics and users prefer to network with users in the same industry-occupation and location, other things being equal. We also find that users respond to popularity and tend to form and maintain links to popular users as well as users that have common connections.
Our estimated model can be used to improve the quality of recommendation systems and for counterfactual policy simulations of business networks, e.g. to assess the impact of new services or exogenous events on the networking in equilibrium Mele2020AEJPol. \\
We contribute to the network economics literature in three complementary ways.\footnote{For an extensive literature review please refer to Jackson2008, GrahamDePaula2020, DePaula2017, Chandrasekhar2016 .} First, we use a tractable structural model of network formation to understand individual and aggregate networking on the job, where agents payoffs depend on observable and unobservable characteristics, as well as linking externalities Mele2017, Mele2020, MeleZhu2021, BoucherMourifie2017, GrahamDePaula2020. The network formation process converges to a stationary equilibrium, corresponding to a mixture of exponential random graphs SchweinbergerHandcock2015, Mele2020, BabkinEtAl2020.
Second, we use unique data from a business card exchange platform to estimate the model, in particular preferences for networking that depend on observables, unobservables and endogenous equilibrium network features (equilibrium externalities). Our data contain digitized information about in-person business interactions, while most of the literature relies on in-person data collected through surveys,\footnote{A very popular dataset is Add Health, containing a survey of high school friendship networks. Many authors use this dataset, e.g. Mele2020,BoucherMourifie2017. Similar datasets are collected in development economics, e.g. BanerjeeEtAl2013.} or data from online platforms, where interactions occur exclusively in the online media.
Third, we propose a scalable estimation algorithm for this class of models, by mixing several approximation and estimation methods. Most of the literature on structural network formation models has relied on small networks to estimate the models, because of the complexities of computing equilibria and the presence of externalities in the model of network formation. In the current work, we also include unobserved heterogeneity in the model, thus increasing the computational complexity even further. Previous work has approached estimation in different ways. Some authors do not include externalities Graham2014, some use equilibrium properties and subnetworks to reduce the computational burden Sheng2020, DePaulaEtAl2014, others exploit pseudolikehood methods BoucherMourifie2017 or incomplete information Leung2015. We exploit the fact that in equilibrium our model generates networks in the class of hierarchical exponential random graph models SchweinbergerHandcock2015, Mele2020, a mixture model of network formation with complex dependencies among links. Estimation of such models via Bayesian methods is intractable for large networks SchweinbergerHandcock2015, Mele2020, SchweinbergerStewart2020, BabkinEtAl2020. Maximum likelihood estimation also does not scale well with the size of the network. We use variational approximations WainwrightJordan2008,Bishop2006, BickelEtAl2013 and efficient algorithms to estimate the unobserved heterogeneity; this is achieved by using state-of-the-art computational algorithms VuEtAl2013, BabkinEtAl2020, reducing the memory usage and with appropriate initialization of the maximization Rosvall_2009. The second step of the algorithm uses pseudolikelihood methods, conditional on the estimated latent block structure to estimate the structural payoff parameters BoucherMourifie2017. Our two-step procedure is similar to ideas proposed in empirical industrial organization or recent work by BonhommeLamadonManresa2019 for bipartite networks.
The remainder of the paper is organized as follows. In section (ref) we briefly describe the data from Eight. Section (ref) develops and analyzes the theoretical model. Section (ref) describes our two-steps estimation algorithm, and results are shown in Section (ref). Section (ref) concludes. Additional details about the computations are provided in appendix (ref).
We employ data from Eight,\footnote{\url{https://8card.net/en}} a multi-platform contact management and professional social networking tool for individuals launched in 2012 and provided by Sansan, Inc. Founded in 2007, Sansan, Inc. is a Japanese company that offers business card-based services for corporations and individuals. It is the largest provider in the Japanese market, with its corporate service holding over 80% of the market share.\footnote{\url{https://ir.corp-sansan.com/en/ir/news/news391044820607604029/main/0/link/Presentation%20Material%20for%20FY2020%20Q2%20(EN)_revision.pdf}} With more than 2.8 million registered users, Eight is the leading professional social network in Japan. Centered on business cards, it acts as a contact manager, as well as a networking tool, and offers functionality such as a home feed, user profiles, and instant messaging. Users connect within the context of the Eight network by scanning each other’s business cards or by sending online friendship requests. Other services by Eight include paid premium plans for individuals and companies, and the direct recruiting platform Eight Career Design.
Eight allows users to scan business cards with a smartphone’s camera, and to set the date the encounter happened\footnote{Some information in English about the business card database can be found in \url{https://datalp.sansan-dsoc.com/}}. OCR algorithms extract information from the business card image, including the name and company of the person and the office address, among other items. This makes it possible to identify individuals and organizations.\footnote{A gentle introduction in English to the digitization process can be found in \url{https://en.sansan-dsoc.com/data/imagerecognition/}} In order to register for the service, users need to scan their own business cards (hereafter called profile card). In case of changes, users can update their profile card information by a simple scan.
The data from Eight that is used in this research includes only anonymized information on connections formed between January and December of 2019 among users that have agreed with Eight's Terms of Service. Nodes represent Eight users who have uploaded a profile card at least once by the end of 2019. We keep only nodes for which all covariates used in the analysis have non-missing values and that belong to the largest connected component of the resulting network. An edge is formed between user A and user B when either A scans B’s business card or vice versa. We exclude purely digital connections, and concentrate only on face-to-face encounters. We assume that a business card exchange is bilateral, and therefore the network is undirected. We also consider only the first contact between a pair of users, so that each link has a weight of 1. Self-loops are excluded.
We obtain node attributes from the latest profile card uploaded by the user in order to cover three sources of homophily: geographic proximity, job type similarity and industrial similarity. We extract the user's job category from the job description in its profile card, and assign users an industrial category based on the user's place of employment. We construct an industry-occupation covariate as the interaction between the industrial category of the user's company and the user's occupation code, which can take any of 8,006 unique values.
For measuring geographic proximity we employ the Hexagonal Hierarchical Spatial Index (hereafter H3) created and open sourced by Uber\footnote{\url{https://eng.uber.com/h3/}}. It is an indexing system that projects the sphere of the Earth into an icosahedron and constructs a grid of nested hexagons, each one of which is assigned a unique identifier or index. The H3 indexing system supports 15 resolutions, where hexagons at higher resolutions have a smaller average area\footnote{A table with the mean area per hexagon at each resolution can be found at \url{https://h3geo.org/docs/core-library/restable/}}, and hexagons within the same resolution do not overlap.
In order to assign an H3 index to a node we perform geocoding on the office address in the user's profile card and obtain its latitude and longitude coordinates. The geocoding mechanism is based on data obtained through the Location Reference Information Download Service \footnote{Location Reference Information Download Service (Geospatial Information Authority of Japan) \url{https://nlftp.mlit.go.jp/index.html}}, and address geolocation data \footnote{Address Geolocation Data (Geospatial Information Authority of Japan) \url{https://www.gsi.go.jp/kihonjohochousa/jukyo_jusho.html}}, both provided by the Geospatial Information Authority of Japan. Addresses within Japanese cities are represented using three nested subdivisions: chome, ban and gou (building frontage), in order of granularity. The quality of geocoding varies depending on the availability of data at each region and the level of granularity. Among the addresses in the dataset, 52% could be matched up to the gou level (highest level of accuracy), and 43% to the ban level. The remaining 4.6% is matched at the chome level. We assign each node the H3 index that contains its coordinates at the resolution of 8. The user's location is thus represented by a hexagon of roughly 0.74 $\text{km}^2$. The choice of the resolution represents a trade-off between the capability to model homophily, and memory requirements of its usage for block recovery. Choosing a resolution that is too low or too high would prevent us from capturing the effect of spatial homophily, as too few/many nodes would be located in the same tile, and too low a resolution imposes high memory requirements to the matrix representation of homophily in this dimension.
One alternative would be to perform matching at the zip code level; however, regions sharing the same zip code can differ greatly in area, and important differences may arise between large cities and the countryside. At resolution 8, H3 index similarity captures more business connections than zip code similarity while still being sparse enough for keeping the computation manageable. Although H3 index similarity does not provide a measure of spatial homophily at large distances, the area of tiles at resolution 8 is large enough to contain several high rise office buildings and commercial areas, and therefore captures the cost of business connections at a local level. The H3 index covariate in the dataset has 28,392 unique values.
The resulting network has $242,223$ nodes and $682,920$ edges. The network is very sparse, with a density of roughly \num{2.3e-5}. The network contains $27,289$ triangles and $9,661,321$ 2-stars. All the 47 prefectures and 1,650 cities (roughly 96% of the total Japanese cities), are represented in the sample. Nodes based in Japan’s largest urban areas in Tokyo and Osaka account for 50% of the nodes. The data is highly geographically concentrated. The most common H3 index is shared by 2,603 nodes, and 90% of the tiles contain 10 or less nodes.
34.3% of the sample is composed by persons in Sales-related occupations, followed by company directors (13.1%). 61.5% of the sample holds the ranks of staff. Nodes are more evenly distributed across industrial categories, with the largest category, IT companies, accounting for only 3.3% of the sample. Same industry-occupation connections represent a 2.3% of the total connections, while same H3 index connections represent roughly 1.8%.
The degree distribution resembles a power law, just like many other large social networks, as shown in Figure (ref).
We model the decision of users to create professional relationships through in-person interactions or via the business card exchange platform. The set of users is $\mathcal{I}=\lbrace 1, 2, ..., n \rbrace$ and each user $i\in \mathcal{I}$ is characterized by a vector of observable characteristics $\bm{x}_i$, such as gender, location, etc. Additionally, each user is characterized by a $K$-dimensional vector $\bm{z}_i$ that is unobservable to the researcher, but it is observed by other users. The vector $\bm{z}_i = (z_{i1}, ..., z_{iK}) $ is interpreted as an assignment to one of $K$ types; we say that user $i$ belongs to type $k$ if $z_{ik}=1$ and $z_{i\ell} = 0$ for all $\ell\neq k$.
Business card exchanges are recorded in the adjacency matrix $\bm{g}$, whose generic element $g_{ij}=1$ if users $i$ and $j$ have exchanged a business card, and $g_{ij}=0$ otherwise.
We model users' objective function as a function of the network $\bm{g}$, observable characteristics $\bm{x}$, unobservable types $\bm{z}$ and parameter vector $\bm{\theta} = ( \bm{\alpha}, \bm{\beta},\bm{\psi},\bm{\gamma})$
The payoff of direct interactions $u_{ij}(\bm{\alpha},\bm{\beta}):=u(\bm{x}_i,\bm{x}_j,\bm{z}_i,\bm{z}_j;\bm{\alpha},\bm{\beta})$ includes both costs and benefits of interacting. The payoff is a function of observable characteristics $(\bm{x}_i, \bm{x}_j)$ as well as unobservable types $(\bm{z}_i, \bm{z}_j)$. User $i$ receives a net benefit $u_{ij}(\bm{\alpha},\bm{\beta})$ for a link to user $j$. The second part of the payoff is the effect of popularity $w_{ijr}(\bm{\psi}):=w(\bm{x}_i,\bm{x}_j,\bm{x}_r,\bm{z}_i,\bm{z}_j,\bm{z}_r;\bm{\psi})$. If user $i$ forms a link to $j$, she receives an indirect payoff $w_{ijr}(\bm{\psi})$ from each link formed by $j$. Therefore we can interpret the second term in the utility function as a weighted payoff from popularity, where the weights are functions of observable and unobservable characteristics. Finally, the third term in the payoff is the effect of transitivity $v_{ijr}(\bm{\gamma}):=v(\bm{x}_i,\bm{x}_j,\bm{x}_r,\bm{z}_i,\bm{z}_j,\bm{z}_r;\bm{\gamma})$, or the payoff from common connections. Each user $i$ receives a payoff $v_{ijr}(\bm{\gamma})$ from each user $r$ that is connected to both $i$ and $j$. In the standard strategic network formation literature, direct connections are assumed costly, while indirect connections are free. In this model we do not need to assume that, as the payoff structure allows for costly indirect benefits as well, in principle.
We conceptualize the network formation process as a sequential game, where users form links over time; however, the researcher only observes the network at a particular point in time.\footnote{There is a growing literature in network econometrics using sequential network formation to improve tractability and obtain an equilibrium selection rule. See DePaula2017, ChristakisEtAl2010, Graham2020, GrahamDePaula2020, Chandrasekhar2016, Mele2017,Mele2020,MeleZhu2021, Jackson2008, JacksonWatts2001 for examples. } We thus focus on the analysis of stationary equilibria of the model.
In each period two users meet and decide whether to exchange their business cards, by maximizing the surplus generated by the interaction. Before making a decision the users also observe the matching quality of their link, which is also unobserved by the researcher. Formally, in period $t=0$, each user is randomly assigned to a unobservable type $z_i$, drawn from a Multinomial Distribution
Conditional on the realization of the types' assignment $\bm{z}$, the network is formed over time according to the following sequence:
This process of network formation generates a Markov Chain of networks, where each network only depends on the previous period's network. In each period, only one link is updated and only two users are actively playing, best-responding to the previous period link decisions of the other players. To characterize the long-run behavior of the network, we make the following formal assumptions.
The first assumption, imposes that any pair has a positive probability of meeting, however small. This implies that in the long-run two users have (possibly infinitely) many opportunities to form and delete a link. The second assumption restricts the preferences so that we can identify the externality effects. The final assumption about the matching shock is standard in random utility models and is common in the network econometrics literature GrahamDePaula2020, DePaula2017, Chandrasekhar2016, Graham2020.
Under these assumptions we can show that the game of network formation is a potential game and it converges to a unique stationary equilibrium distribution over networks. Therefore, given the block structure $\bm{z}$, in the long-run we expect to see a network $\bm{g}$ with probability $\pi(\bm{g},\bm{x},\bm{z};\bm{\theta})$, as shown in the next proposition Mele2017, Mele2020, MeleZhu2021.
The proof can be found in Mele2020. Proposition (ref) shows that -- after conditioning on the realized unobservable heterogeneity $\bm{z}$ -- the network formation game admits a representation as a potential game MondererShapley2006, where all the deterministic incentives of the users to form links are captured by the potential function (ref). Indeed, we can show that
for all user pairs $i,j\in \mathcal{I}$.\footnote{See Mele2017, MeleZhu2021, Mele2020 for a formal proof of this statement.} This means that all the incentives of each pair of users to create or delete a link, net of the matching quality, are described by the aggregate potential function.
The potential game characterization in the proposition implies that the equilibrium pairwise stable networks (with transfers) can be obtained by finding the (local) maxima of the potential $Q(\bm{g},\bm{x},\bm{z};\bm{\theta})$. Therefore, in the long-run we expect to see the pairwise stable networks with high probability, according to the stationary distribution $\pi(\bm{g},\bm{x},\bm{z};\bm{\theta})$. In such equilibria, the surplus generated by each link is not necessarily split equally between the users involved in the relationship. This allows us to model the fact that some networking relationships are asymmetric or players have different bargaining power Jackson2008.
Estimation of this model is challenging because the likelihood depends on the normalizing constant (ref) that is hard to evaluate even with modern supercomputers Snijders2002.\footnote{The constant is the sum of the exponential of potential functions over all possible network configurations. This sum thus includes $2^{(n(n-1)/2}$ terms. Even considering parallelization of the computations, the exact computation of the normalizing constant is either impractical or infeasible for most network sizes. In our data we have around $n=240,000$. Additionally, the standard MCMC methods used in the literature to estimate ERGMs converge too slowly for our data, as in the best case scenario the algorithms converve in $n^2 \log(n)$ steps BhamidiEtAl2011, Mele2017. } This is especially true for the size of our dataset.
To get around some of these challenges, we exploit a particular specification of the model to obtain a likelihood that can be factorized in between- and within-blocks components, after conditioning on the unobserved block structure. This factorization crucially decreases the complexity of computations.
To estimate the block structure, we approximate the model using a stochastic blockmodel. Given the estimated block structure we estimate the full model using maximum pseudolikelihood estimators.
These methods bypass the need to compute the likelihood and the normalizing constant, thus allowing estimation in massive networks. In this section we provide several details about the specification of the payoff functions, the estimation of the unobserved heterogeneity (the block structure) and the estimation of the structural payoff parameters.
The model specification is crucial for a tractable estimation procedure, thus we adopt the following specification with local externalities Mele2020, Schweinberger2020, SchweinbergerHandcock2015, BabkinEtAl2020. The parameter $\bm{\alpha}$ only depends on the unobservable types, taking value $\alpha_w$ if the users belong to the same type; otherwise it is $\alpha_b$. The parameters $\bm{\beta}$ contain the net marginal benefits of observable characteristics and we assume that they vary within-types ($\bm{\beta}_w$) and between-types ($\bm{\beta}_b$). Finally, we specify externalities as local, that is we assume that the externality is part of the payoff only if all the users involved in the relationship belong to the same unobserved type.
The specification differs from other papers using the HERGM framework. In fact, we allow the parameters for the observable covariates to vary within and between blocks, while most papers assume homogeneity for the entire network SchweinbergerHandcock2015, BabkinEtAl2020. This allows us more flexibility in estimation.
The specification with local transitivity and local popularity is convenient for estimation and computation. Indeed, we can show that the potential function ((ref)) can be decomposed in the sum of within- and between-community potentials. Let $\bm{g}_{k,l}$ denote the sub-network among users in blocks $\mathcal{C}_{k}$ and $\mathcal{C}_{l}$. Let $\bm{x}^{(k)}$ denote the observable covariates of users in community $\mathcal{C}_k$. Let' define the functions:
Then the potential function can be re-written as
From a practical standpoint, this decomposition implies that the likelihood of the network factorizes as product of within- and between-community likelihoods, facilitating estimation and identification.
where the within-community and between-communities normalizing constants are, respectively
Notice that the between-community potential $Q_{k,l}(\bm{g}_{k,l},\bm{x}^{(k)},\bm{x}^{(l)}, \bm{z};\bm{\theta}) $ does not include the link externalities (transitivity and popularity). Therefore, the second part of likelihood ((ref)) is the product of conditionally independent links,
To summarize, Assumption (ref) guarantees independence of between-communities links; on the other hand, within-community links may have strong dependence. In aggregate, our model maintains the complex correlation structure of exponential family random graphs (ERGMs) locally, but allows for weaker dependence among links globally.
The likelihood factorization described in the previous section attenuates some of the computational issues in estimation. However, most applications to date have focused on networks of few hundred nodes when using a Bayesian estimation strategy SchweinbergerHandcock2015,Mele2020; and networks with few thousands nodes when using an approximate maximum likelihood strategy BabkinEtAl2020. Our data contain hundreds of thousands nodes and therefore we have to use alternative methods and computational strategies to obtain a computationally tractable estimation method.
Our approximate algorithm consists of two steps, as suggested in BabkinEtAl2020. In step 1 we estimate the block structure $\hat{\bm{z}}$, approximating the likelihood of the model as the one of a stochastic blockmodel. We then use a variational approximation to obtain a tractable lower bound of the log-likelihood and accelerate the estimation using a minorization-maximization algorithm suggested in VuEtAl2013. In step 2, given the estimated block structure $\hat{\bm{z}}$, we estimate the parameters of the model $(\bm{\alpha}, \bm{\beta}, \bm{\psi}, \bm{\gamma})$ using maximum pseudolikelihood estimators.
This procedure is based on two considerations. First, the likelihood of a stochastic block-model imposes the same probability of the original likelihood on between-block links. Therefore, the approximation only involves the within-block sub-networks. As long as the network is large, most of the probability mass is on the between-block links, and therefore this approximation works well.\footnote{Formal statements are contained in BabkinEtAl2020.}
Second, while the likelihood of a stochastic block-model is intractable, there exist variational methods of inference to recover its parameters. Variational methods maximize a lower bound to the likelihood, recovering an estimated block structure. The asymptotic analysis shows that variational estimates are consistent and asymptotically normal BickelEtAl2013, DaudinEtAl2008. Computations can be sped up by using minorization-maximization techniques VuEtAl2013.\footnote{Alternatively, spectral methods hold promise in dealing with massive network data AthreyaEtAl2018a, MeleEtAl2021. } However, the implementation in VuEtAl2013 does not take into account the observable covariates, which are crucial in our application. Therefore, we extend their algorithm to include (discrete) covariates.
We present the two steps in the following subsections, while providing more technical details in Appendix (ref).
In step 1, we approximate the log-likelihood of the model, as if there are no link externalities, i.e. we rewrite the likelihood as if $(\psi, \gamma)= (0,0)$. This approximation works as long as we have many blocks, that is when $K$ is relatively high compared with the size of the network $n$.
The full likelihood of our model can be written as follows
Conditional on the community structure $\bm{z}$, the probability that we observe network $\bm{g}$ is given by $\pi(\bm{g},\bm{x},\bm{z};\bm{\theta})$: this corresponds to the probability of observing the network in the long run, that is
The complete likelihood ((ref)) is obtained by multiplying the likelihood ((ref)) by the probability of firm types/communities $\bm{z}$, that is $P_{\bm{\eta}}\left(\bm{Z}=\bm{z}\right)$, given by a multinomial distribution
and summing over all possible community structures $\bm{z}\in\mathcal{Z}$.\\
Our estimation method is based on the observation that if the externalities are not present in the model, $\psi=0$ and $\gamma=0$, the likelihood is the same as the one of a standard K-block stochastic blockmodel with nodal covariates. Therefore we consider the approximation
To estimate the block-structure $\bm{z}$ we use a variational approximation for stochastic blockmodels, and compute the lower bound of the log-likelihood WainwrightJordan2008, Bishop2006, BabkinEtAl2020. Let $q(\bm{z})$ be an approximating distribution over blocks $\bm{z}$. Then the lower bound $\ell_B (\bm{g},\bm{x};\bm{\alpha},\bm{\beta}, \bm{\eta})$ is obtained via an application of Jensen's inequality
The variational method finds the best approximating distribution $q(\bm{z})$, by finding the best lower bound. Because this variational problem is also intractable and cannot be solved in closed-form, we choose $q(\bm{z})$ from a tractable family of distributions, as suggested in the literature WainwrightJordan2008. For our model, it is natural to choose a multinomial distribution $q_{\bm{\xi}}(\bm{z})$
that can be optimized with respect to the variational parameters $\bm{\xi}_i$'s to obtain a tractable bound $\ell_B (\bm{g},\bm{x},\bm{\alpha}, \bm{\beta}, \bm{\eta}; \bm{\xi})$
where $\log \pi_{ij,kl}(g_{ij}, \bm{x},\bm{z})$ is the log-likelihood of a link between nodes in blocks $k$ and $l$
and $u_{ij,kl}(\bm{\alpha},\bm{\beta})=u(\bm{x}_i,\bm{x}_j,z_{ik}=z_{jl}=1,\bm{z};\bm{\alpha},\bm{\beta})$ is the direct net benefit payoff of user $i$ in block $k$ from forming a link with user $j$ in block $l$. \\
Minorization-Maximization. This estimation framework for stochastic blockmodels is relatively standard in the literature and it enjoys several asymptotic properties and guarantees BickelEtAl2013. In particular, the estimates are consistent and asymptotically normal. The EM algorithm consists of iteratively updating the parameters via an expectation and a maximization step, whose updates are available in closed-form DaudinEtAl2008, BickelEtAl2013.
However, maximizing the lower bound $\ell_B (\bm{g},\bm{x},\bm{\alpha}, \bm{\beta}, \bm{\eta}; \bm{\xi})$ with respect to $\bm{\xi}$ can still be impractical in very large datasets, because the iterative update for each $\xi_{ik}$ depends on $(n-1)K$ other terms $\xi_{jl}$. These updates are time consuming. Additionally, the iterative algorithm used for computing the lower bound approximation is a local algorithm and may get stuck in a local maximum.
To alleviate these computational problems, we extend the Minorization-Maximization methods of VuEtAl2013 in two complementary directions. First we allow the algorithm to incorporate (discrete) covariates. The original algorithm is designed for stochastic blockmodels without any observable covariates, so this is a significant improvement in terms of applicability of the method. Second, we provide an efficient computational algorithm that exploits matrix algebra rather than nested loops in computation, to speed up computations by a factor of 14000. This allows estimation in massive networks. Our implementation takes advantage of the sparsity of the graph by making use of sparse matrices where possible in order to make efficient use of the memory, which is a problem when dealing with massive datasets.
The idea of minorization algorithms is to find a function that approximates the lower bound $\ell_B (\bm{g},\bm{x},\bm{\alpha}, \bm{\beta}, \bm{\eta}; \bm{\xi})$, while being simpler to maximize. In practice, a function $M\left(\bm{\xi}; \bm{g},\bm{x},\bm{\alpha},\bm{\beta},\bm{\eta}, \bm{\xi}^{(s)} \right)$ minorizes the likelihood lower bound $\ell_B (\bm{g},\bm{x},\bm{\alpha}, \bm{\beta}, \bm{\eta}; \bm{\xi})$ at parameter $\bm{\xi}^{(s)}$ and iteration $s$ if
where $\bm{\alpha}, \bm{\beta}, \bm{\eta} $ and $\bm{\xi}^{(s)}$ are fixed.
We follow VuEtAl2013 and use the following function for a stochastic block model lower bound
The main difference from VuEtAl2013 is that our model includes observable (discrete) covariates. Therefore, the updates of the maximization are slightly different.
As in a standard Variational EM algorithm, we can write down the parameter updates in closed-form. The update rules for $\bm{\xi}$, $\bm{\eta}$, and $\pi_{ij;kl}(g_{ij}, \bm{x},\bm{z})$ follow
and
for $ k, l = 1, \ldots, K$ and $d, \chi_{1}, \ldots, \chi_{p} \in \{0, 1\} $, respectively. In the formula for $\pi_{ij;kl}^{(s+1)}(d, \chi_{1}, \ldots, \chi_{p},\bm{z}))$ we have used the notation $\chi_{p,ij}$ to denote and indicator variable equal to 1 if the (discrete) nodal covariate $p$ of $i$ and $j$ are the same, i.e. $\chi_{p,ij} = \bm{1}\lbrace x_{ip}=x_{jp} \rbrace $. Generalizations of this specification are allowed.
The estimated block structure $\widehat{\bm{z}}$ is obtained by choosing the modal block assignment, that is $\widehat{z}_{ik}=1$ if $\widehat{\xi}_{ik} \geq \widehat{\xi}_{i\ell} $ for all $\ell \neq k$ and $\widehat{z}_{il}=0$ for all $l\neq k$.
Conditioning on the estimate of $\widehat{\bm{z}}$, we estimate the structural parameters $\bm{\theta} = \left(\bm{\alpha}, \bm{\beta},\psi,\gamma \right)$ by maximum pseudolikelihood (MPLE) methods BoucherMourifie2017, Snijders2002, BabkinEtAl2020. This amounts to maximize the product of the conditional link probabilities.
Formally, given the estimated $\widehat{\bm{z}}$, we compute the conditional probability of a link
where $\Lambda(u) = e^u /(1+e^u)$ is the logistic function. The pseudolikelihood function is
and the estimator is the maximizer of the log-pseudolikelihood
The asymptotic framework for the maximum pseudolikelihood estimator was recently analyzed in BoucherMourifie2017. It can be shown that the estimates are consistent and asymptotically normal under some regularity conditions.
We estimated our model for the network presented in Section (ref) with a maximum of 1,500 blocks, using 250 iterations of the EM algorithm. After recovering the estimated block structure $\widehat{\bm{z}}$ and controlling for the node covariates described in Section (ref), we estimate the structural parameters using the maximum pseudolikelihood method and accounting for node covariates on the block recovery step. For comparison, we also perform the estimation without taking node covariates into account for block recovery. Our implementation of the algorithm used to obtain these results is available in the lighthergm R package, which can be found at \url{https://github.com/sansan-inc/lighthergm}. We present the results of each step in detail in the subsections below.
The block structure estimation step is by far the most computationally intensive part of the estimation. For this application we employed an Ubuntu Linux machine with 128 GB of memory and 64 processor cores. We set the maximum number of types/blocks to 1,500. The computation is performed with about 35 GB of memory for the block recovery step accounting for node covariates, although it can be performed with well under 32 GB of memory when node covariates are not employed. All processor cores are in use during most of the calculation time.
First, we initialize the blocks by using the Infomap algorithm by Rosvall_2009. Infomap presents several advantages over other clustering algorithms for our particular use case. First, Infomap's time complexity is linear in the number of edges, which makes it a good choice for initializing the block memberships on very sparse networks. Yang2016 show that Infomap performs better than other algorithms with similar time complexities at the same values of the mixing parameter. Furthermore, Infomap's performance at recovering the true communities is independent of the network size. In comparison, the default initialization algorithm on the original hergm R package version 4.1-7 is Walktrap Pascal2005, which, despite having properties that make it a good candidate, has a space complexity that is quadratic in the number of nodes, making it an expensive choice for clustering large networks.\\
Starting at the initial block structure, we apply 250 iterations of the fast EM algorithm. Each EM algorithm iteration takes approximately 14 minutes, for a total of 38.3 hours to complete the whole EM iteration part of the block structure estimation step when accounting for node covariates. In comparison, VuEtAl2013 employ a similar variational approach on a network of $131,000$ nodes with only $20$ blocks and without nodal characteristics with 100 EM iterations taking a total of 24 hours.
Figure (ref) shows the improvement in the target function's Lower Bound at each iteration. The improvements are monotonic, and converge after close to 100 iterations, although the improvement keeps being positive after 250 iterations. In order to understand how much the block structure changes with the number of iterations, we compute the Yule's coefficient with respect to the initial block structure obtained by Infomap. The Yule's coefficient measures the similarity between two block structures regardless of the labels. It takes values between 0 and 1, where higher values mean a higher similarity. We find that the final blocks differ considerably from the initial structure. After 100 iterations, the Yule's coefficient is roughly 0.74, and after 250 iterations it goes down to 0.05.
Figure (ref) shows the sizes of all the blocks in ascending order of the number of affiliated nodes. The dotted line marks the median block size of 59 nodes. We observe that a few blocks contain a large number of nodes, and the largest block contains 37,615 nodes, representing a 15.5% of the total nodes. Looking at the distribution as a whole, blocks are in general quite homogeneous in size. The block size distribution is strongly concentrated around the median, and has an interquartile range of only 22 nodes. Figure (ref) displays a subset of the network and the estimated block affiliation of the nodes, and shows that the recovered block structure in fact represents areas of the graph with denser connectivity.
Figure (ref) shows the relationship between the number of nodes and the share of the five largest blocks by prefecture. Nodes in hub prefectures such as Tokyo and Osaka tend to be distributed across many blocks, with no single dominating cluster. On the other hand, smaller prefectures tend to be dominated by a few large blocks. Okinawa is a clear outlier, with the largest block accounting for over 40% of its nodes.
Finally, when covariates are not employed at the block recovery step, a more skewed distribution of block sizes is obtained. The size of the largest block in this case is 85,512, and the median block size is 72. The Yule's coefficient between this partition and the one obtained when employing node covariates is 0.37, suggesting that in fact, accounting for homophily on observable characteristics can have an impact on the structure of the resulting partition, holding the number of EM iterations and the initial block structure constant.
We take the community structure obtained in the previous step and proceed to estimate the within-block and between-block ERGM parameters. Given the assumptions in our model, we separate the between-block connections from the within-block ones, and estimate the parameters for each set independently employing maximum pseudolikelihood estimation. This step requires fewer resources and processing time compared to the previous step. The estimates for both sets of parameters are shown in Table (ref).
Standard errors correspond to the ones obtained by each separate maximum pseudolikelihood estimation and thus do not consider the error in the block structure recovery step. Regarding the between-block model results, estimates are similar regardless of whether node covariates were employed to inform the block recovery step. This is to be expected, given that most of the possible connections are across blocks, despite the differences in the final block structure. Coefficients are all significant at the 1% level. The coefficient for the edges term reflects the sparsity of the network, while geospatial and industrial-occupational homophily are significant factors explaining business connections across blocks.
Results for the within-block connection model show that random connections within the same community are slightly more likely than across blocks, suggesting that persons prefer forming business connections among peers within the same communities, everything else constant, although a formal statistical test is required. We observe a significant preference for transitivity and popularity, which highlights the importance of externalities on the network formation process within communities. Similar to connections across communities, homophily in location, industry and occupation is an important component of the utility of business connections within communities. Coefficients for the edges and externalities terms do not differ greatly depending on whether node covariates are employed for block recovery; however, the importance of homophily when explaining business connections within blocks is higher when the block recovery step does not account for covariates.
Networking on the job is an important determinant of mobility and career advancement in many labor markets. In this paper we have studied a network of business relationships using the digital trace of business card exchanges from Eight, a platform for the digitization of business cards containing data from the whole Japan. Our sample contains about 240,000 users of the platform.
Our empirical analysis is guided by a theoretical equilibrium model of network formation where users form relationships based on their preferences for observables, unobservables, and endogenous network features. The stationary equilibrium characterizes the likelihood of observing a network in the data, and we estimate the parameters using approximate maximum likelihood methods. Crucially, the unobserved heterogeneity is discrete, and the equilibrium is a mixture of exponential random graphs SchweinbergerHandcock2015, Mele2020.
We rely on a two-step approach to estimation, first developed in BabkinEtAl2020. The first step involves an approximate clustering of the nodes, to estimate the unobserved (discrete) type distribution. The second step estimates the structural payoff parameters using a pseudolikelihood estimatior BoucherMourifie2017, Snijders2002.
We propose several algorithmic improvements to the model-based clustering algorithm in VuEtAl2013, to include discrete nodal covariates and to speed-up computations through a mix of variational approximations, fast sparse matrix algebra routines and minorization-maximization methods. These improvements allow estimation of the structural model using a massive dataset with about 240,000 users, controlling for (discrete) observable characteristics.
Our analysis shows that this massive business network contains a large number of small (unobserved) business communities. A standard exponential random graph model is unable to capture this feature. This confirms that including unobserved heterogeneity in the network formation model is crucial to understand the business networking patterns in this data.
Our scalable method will allow network researchers to estimate complex models using massive datasets. Previous work on the econometric analysis of large networks has been limited by the complexity of estimation algorithms, and for most studies the definition of a large network has been mostly limited to a few thousands of nodes. Our algorithmic improvements makes it possible to analyze networks with hundreds of thousands of nodes, while using relatively few resources.
Furthermore, additional improvements in computational speed and scalability can be obtained, e.g by using GPUs. Larger networks can be handled at a lower cost by means of distributed computing. Crucially, the space complexity of our implementation depends heavily on the size of the matrix of variational parameters, which is a dense matrix and grows with the number of clusters. Since it is reasonable to expect that the number of unobservable blocks grows with the size of the network, the dimensions of this matrix may impose a limitation to the size of networks that can be feasibly analyzed.
We also acknowledge that our implementation of the clustering algorithm makes use of the fact that the network is sparse and the discrete covariates follow the same sparse pattern. In particular, we employ feature adjacency matrices to facilitate matrix algebra. Our current implementation requires these matrices to be sufficiently sparse to fit into memory. This complication arises because the algorithm requires the creation of a number of sparse matrices that grows with the square of the number of discrete covariates. Both of these issues impose limitations to the type and number of covariates that can be employed. We believe that it is possible to further improve the speed of the block structure recovery step and at the same time break the dependency on the sparsity of the discrete features, which should make it possible to employ more and better covariates, and to run more iterations of the EM algorithm, thus improving its block recovery capabilities. We expect to extend our algorithm to include these improvements in future versions of this research, and that solving some of this additional complications may contribute to popularize the industrial use of exponential random graph models for the analysis of large social networks and their simulation.