Triadic Network Formation

Triadic Network Formation

\, and Cavit Pakel\footnote{University of Oxford. Email: [email removed]}.} }

Introduction

Networks are prevalent in economics, and they appear in many contexts: firms link to suppliers, countries trade in products, and researchers collaborate on projects. Increasingly detailed datasets capture such interactions across multiple dimensions, revealing structures that go beyond simple dyads. Triadic relationships in particular increasingly feature in empirical work; e.g., in production decisions (BDM22), import-export link formation of multinational corporations (CLMT24), and also in investigations of the effect of economic reliance on a trade partner on political alignment (KLR24). At the same time, sparsity appears to be a pervasive feature of network data.\footnote{A typical example is firm-level production network data, which is the focus of a well-established literature in economics (CTS19). These networks cover the universe of firms in a country and are typically characterised by the overwhelming majority of the firms having very few connections (BMUM18 and BMY19).} In modelling such networks, it is standard to account for unobserved heterogeneity using fixed effects.

In this paper we focus on the formation of sparse triadic networks, and develop econometric tools for estimation in a nonlinear binary choice logit model of link formation with dyad-level heterogeneity. Crucially, in such models it is not the links themselves but specific configurations of links---wirings---that carry information about the structural parameters. Our key finding is that under dyad-level heterogeneity, once sparsity exceeds a threshold the network may generate infinitely many links yet fail to produce enough informative wirings for inference. This is a fundamental departure from the behaviour under node-level heterogeneity: indeed we show that under node-level heterogeneity informative wirings accumulate automatically with link formation. We characterise precisely how sparse is `too sparse', deriving thresholds that determine when consistency holds and when asymptotic normality can be attained.

Our findings carry important practical implications. Dyad-level effects (e.g. importer--exporter, importer--product, exporter--product) offer a richer and more realistic representation of heterogeneity than classical node-level effects (e.g. importer, exporter, product). Yet our results reveal that this flexibility requires sparsity to remain within certain limits, which may not hold in more granular or disaggregate data. When these conditions are not met, one option is to revert to node-level heterogeneity, though such a modelling choice may not always be appropriate given the context.

To formalise our theoretical framework and contributions, our link formation model is given by

where $i,j,k$ are the indices for the three network dimensions, $X_{ijk}$ is the vector of observable covariates with the corresponding parameter vector $\beta_0$, and $\varepsilon_{ijk}$ is the triad-specific logistic random shock. $A_{ij}$, $B_{jk}$ and $C_{ik}$ are the dyad-level fixed effects. Our interest is in estimating $\beta_0$. Under sparse network asymptotics, fixed effects are generally not point-identified, which in turn leads to failure to identify $\beta_0$. We use a conditional likelihood approach based on sufficient statistics for the fixed effects in order to point-identify $\beta_0$. For dyadic link formation, sufficient statistics have been proposed in the form of informative wirings between groups of four nodes, called tetrads (Charbonneau17, Graham17, Jochmans18). In the triadic setting we propose a conditional likelihood estimator based on hexad subnetworks---the hexad logit estimator---and establish the conditions under which it is consistent and asymptotically normal.

To provide some more precise discussion on the collection of informative wirings, let $N$ be the number of nodes in each of the three parts of the network and define $\rho_N = \mathrm P (Y_{ijk}=1)$, the unconditional probability of link formation. The rate at which $\rho_N$ tends to zero measures the level of sparsity. A necessary condition for statistical analysis is that the average node degree in the network remains bounded from below by zero as $N \to \infty$, which ensures asymptotic non-emptiness of the network---but still allows the network to be sparse. Consequently, as $N\to\infty$, the network accumulates links. Under node-level heterogeneity, this asymptotic accumulation of links leads to a corresponding accumulation of informative wirings. This holds both in dyadic link formation models, as previously shown, and in triadic link formation, as we establish in this paper.

We show that under dyad-level heterogeneity accumulation of links does not necessarily translate into accumulation of informative wirings: specifically, if sparsity exceeds a certain threshold such that $\rho_N$ tends to zero faster than $N^{-3/2}$, the network will generate infinitely many links and yet produce an asymptotically vanishing number of informative wirings. The reason for this phenomenon lies in the nature of informative wirings. The richer heterogeneity structure translates into wirings that are more intricate, depending on a larger number of links and arranged in highly specific patterns. As a result, they are relatively less likely to appear in data than their counterparts under node-level heterogeneity. Put differently, noticing that the average degree per node in our setting is $O(N^2 \rho_N )$, collecting sufficiently many informative wirings requires networks that are less sparse than just non-empty.

Collecting informative wirings is equivalent to accumulating information on $\beta_0$. The situation identified here therefore has important consequences on the asymptotic behaviour of the hexad logit estimator $\widehat \beta$. First, we show that $\rho_N = O(N^{-\delta})$ with $0\leq \delta<3/2$ is indeed a necessary condition for $\widehat \beta \to_p \beta_0$. Furthermore, we also show that asymptotic normality does not necessarily follow unless $\rho_N = O(N^{-\delta})$ with $0\leq \delta<1$. As will be discussed in more detail, this condition is necessary to ensure that informative wirings accumulate faster than links as the network grows. We show that this is automatically satisfied for the dyadic and triadic models with node-level heterogeneity, but not for the triadic model under dyad-level heterogeneity. This requirement to accumulate informative wirings at a faster rate relates to the conditional likelihood score's interpretation as a U-statistic. In particular, unless $\rho_N = O(N^{-\delta})$ with $0\leq \delta<1$, the score behaves analogous to a highly degenerate U-statistic. Interestingly, this degeneracy occurs in a discrete fashion in the sense that as soon as the U-statistic becomes degenerate this degeneracy is of the highest possible order.\footnote{We note that this degeneracy is different from the one identified by Graham17 for the score of the tetrad logit estimator where he can still obtain asymptotic normality. That appears to be a typical feature of the score function in network data, and also appears in our analysis. However, the degeneracy spelled out here is of a different nature.} The resulting asymptotic distribution will be Gaussian chaos.

Literature review. Our work contributes to three strands of literature:

First, it builds on the work on binary logit link formation in sparse dyadic networks. Charbonneau17, Graham17, and Jochmans18 identify sufficient statistics in the form of tetrad wirings, laying the foundation for conditional likelihood estimation. Recent extensions include distribution regressions (Szini25) and ordered choice models in dyadic settings (MPZ25). Our contribution is to establish the corresponding theory for triadic link formation in a tripartite setting with a full set of dyad-level fixed effects, where the sufficient statistics take the form of hexads. For surveys of the broader literature on econometrics of networks, see dePaula17, Graham17, and GdP20.

Second, our approach is connected to the literature on sufficient statistics and the incidental parameter problem. The idea of eliminating fixed effects through conditioning goes back to Rasch60, Andersen70, and Chamberlain80. In panel data, identification issues arise in short panels (ArellanoBonhomme11); in networks, the analogue is sparsity. DBW25 provide a recent review and also a general result for obtaining sufficient statistics in the binary logit setting. See also BonhommeDano24 who consider the functional differencing approach.

Third, our work relates to the emerging literature on nonlinear three-way models. While linear three-way specifications with heterogeneity have been widely studied (BMW17), nonlinear counterparts are more recent. Recent examples are due to Stammann23, who analyses dense three-way nonlinear panels with dyad-specific effects and proposes bias correction methods, and WeidnerZylkin21 and YZ23 who consider the gravity model.

Model

We model the formation of connections among three nodes from three different partite sets (also called parts). Let $i \in \{1,\ldots,N_1\}$, $j \in \{1,\ldots,N_2\}$, and $k \in \{1,\ldots,N_3\}$ denote nodes from these three sets, respectively. We assume $N=N_1=N_2=N_3$ for notational simplicity. The key unit of analysis is a triad $(i,j,k)$ which is an ordered 3-tuple with node $i$ belonging to the first part, $j$ to the second part, and $k$ to the third part. The set of all $N^3$ such triads is given by $\mathbb T_N = \{1,\ldots,N\}^3$.

The binary variable $Y_{ijk}$ denotes formation of a hyperedge joining the nodes of the triad $(i,j,k)$. The hyperedge that connects the nodes of the triad $(i,j,k)$ is denoted by $[i,j,k]$. The set of nodes and the associated hyperedges form a hypergraph that we call a network. Formally, we consider $3$-partite $3$-uniform networks in the sense that hyperedges connect exactly three nodes, and no two nodes of a hyperedge belong to the same part. Figure (ref) visualises one such network with $N=4$ nodes in each of the three parts. In this example we observe the hyperedges $[i_1,j_2,k_1]$, $[i_2,j_1,k_2]$ and $[i_3,j_4,k_1]$.

We model a triad's decision to form a hyperedge as a binary choice logit model with additive fixed effects. Let $X_{ijk} = h(X_i,X_j,X_k) \in \mathbb R^P$ be a triad-specific covariate (a known transformation of the node-specific observable characteristics) and $\beta_0$ be the associated $P \times 1$ parameter vector. Let also $\varepsilon_{ijk}$ be a triad-specific logistic error with the CDF $\Lambda(u) = \exp(u)/(1+\exp(u))$. We model the link formation decision as

where $A_{ij},B_{jk},C_{ik}$ are dyad-specific fixed effects. For instance, in a setting where firm $i$ imports input $j$ from the exporting country $k$, these correspond to firm-input, input-exporter and firm-exporter effects, respectively. This overlapping heterogeneity structure provides the most general additive fixed effect structure in the given framework. Finally, we assume that the random shocks $\varepsilon_{ijk}$ are independent across all triads, conditional on covariates and fixed effects.

In what follows, we use bold notation when collecting objects across the network. In particular, $\mathbf Y = \left(Y_{ijk}\right)_{\mathbb T_N}$ and $\mathbf X = \left(X_{ijk}\right)_{\mathbb T_N}$ collect the link formation information and observable characteristics across all triads, respectively. $\mathbf{Y}=\mathbf{y}$ denotes a particular realisation of the network. The fixed effects for the triad $(i,j,k)$ are denoted $\mathrm{F}_{ijk} = (A_{ij},B_{jk},C_{ik})$ and fixed effects across all triads are gathered in $\mathbf{F} = \left(\mathrm{F}_{ijk}\right)_{\mathbb T_N}$.

The next assumption formalises the discussion made so far.

Identification

In sparse networks, the fixed effects will generally not be point-identified, as connections involving different node-pairs are typically observed in very few instances, if at all. This is analogous to the identification failure in short nonlinear panels (see, e.g., Chamberlain10 and ArellanoBonhomme11), and is an instance of the incidental parameter issue. As in the panel data literature, for models with logistic random shocks it is usually possible to overcome this problem by a conditional likelihood approach. In what follows, we develop such a conditional likelihood approach that eliminates the fixed effects through carefully constructed sufficient statistics, thereby achieving identification of the common parameter $\beta_0$. Our sufficient statistics will be based on hexad subnetworks.

Hexad subnetworks

A conditional likelihood function that applies to the full network can in principle be constructed by conditioning on a sufficient statistic for the unobserved heterogeneity of the entire network. However, the resulting conditional likelihood would be computationally intractable, even for moderately large networks. Focussing on lower dimensional subnetworks, on the other hand, maintains the statistical properties needed for consistent estimation of $\beta_0$ while remaining computationally practical. This principle has been successfully applied in dyadic network models using tetrad subnetworks Charbonneau17, Graham17, Jochmans18. To handle triadic networks, we use hexad subnetworks.

A hexad is a collection of six nodes. Our analysis is based on a specific type of hexad, consisting of two nodes from each of the three parts of the original full hypergraph. Let $\sigma=(i_1,j_1,k_1,i_2,j_2,k_2)$ be a generic hexad consisting of nodes $(i_1,i_2)$ from part 1, $(j_1,j_2)$ form part 2, and $(k_1,k_2)$ form part 3. Let $\Sigma$ be the set of all possible hexads consisting of two nodes from each of the three parts. There are $m_N = |\Sigma| = N^3(N-1)^3$ hexads in this structure. Notice that each hexad subnetwork contains $2^3$ distinct triads (and therefore anywhere from 0 to $2^3$ hyperedges). Let $X_\sigma$ and $\mathrm{F}_\sigma$ collect the covariates and fixed effects for all triads in the subnetwork generated by the hexad $\sigma$. Hexad subnetworks provide enough structure to find sufficient statistics for the fixed effects $\mathrm{F}_\sigma$, as we will show in the next section. The resulting hexad-specific conditional likelihood functions can then be combined to obtain a (composite) conditional likelihood function for estimation of $\beta_0$.

Let $d_\sigma(\cdot)$ return the degree of a chosen node in $\sigma$; that is, the number of hyperedges that this node belongs to in the subnetwork generated by $\sigma$. The degree sequence for $\sigma=(i_1,j_1,k_1,i_2,j_2,k_2)$ is then given by $D_{\sigma} = (d_\sigma(i_1),d_\sigma(j_1),d_\sigma(k_1),d_\sigma(i_2),d_\sigma(j_2),d_\sigma(k_2)).$ The degree sequence is a key quantity in obtaining a conditional likelihood function; see Graham17.

As an illustration of the concepts introduced thus far, consider the hexads in Figure (ref). Hyperedges are illustrated as line or curve segments of the same colour. In Wiring 1 there are three hyperedges: $[i_1,j_1,k_2]$, $[i_1,j_1,k_1]$ and $[i_2,j_2,k_1]$. Wiring 2, on the other hand, contains four hyperedges: $[i_1,j_1,k_1]$, $[i_1,j_2,k_2]$, $[i_2,j_1,k_1]$ and $[i_2,j_1,k_2]$. The degree sequences of Wirings 1 and 2 are $(2,2,2,1,1,1)$ and $(2,3,2,2,1,2)$, respectively.

Informative wirings and sufficiency

Finding a sufficient statistic for $\mathrm{F}_\sigma$ requires finding degree sequences that admit at least two different wirings that contain the same fixed effects. Consequently, such degree sequences can be used to isolate variation attributable to $\beta_0$ rather than to the fixed effects, thereby obtaining identification of $\beta_0$. This intuition has been used successfully in Charbonneau17, Graham17 and Jochmans18 in dyadic link formation.

We are specifically interested in using the minimal degree sequence that permits identification---by minimal we mean that there exists no other identifying degree sequence where at least one node has a lower degree (and no nodes have a higher degree). In Appendix (ref) we show that in the given setting this is achieved by the degree sequence $(2,2,2,2,2,2)$, which yields two wirings that can be used to construct a conditional likelihood function. We call these two wirings informative wirings or identifying wirings. For a generic hexad $\sigma=(i_1,j_1,k_1,i_2,j_2,k_2)$ these two informative wirings correspond to the following collections of hyperedges:

These are illustrated in Figure (ref): the first informative wiring corresponds to the left panel whereas the second is given in the right panel. Notice that the two informative wirings do not share any common hyperedges despite having the same degree sequence. However, and crucially, they both contain all the fixed effects in $\mathrm{F}_\sigma$. This can also be confirmed by noticing that both wirings contain all possible dyadic links: since our fixed effects are dyad-specific, this confirms that both wirings contain all the fixed effects in $\mathrm{F}_\sigma$. Consequently, the differences between the wirings should be attributable to covariates, enabling point-identification of $\beta_0$.

Letting $\overline{Y}_{abc} = 1-Y_{abc}$, the indicator functions for these wirings are given by

In other words $S_{\sigma,1}$ indicates whether the hexad $\sigma$ admits the first informative wiring or not (and similarly for $S_{\sigma,2}$). Notice that a hexad cannot admit both informative wirings at once. Therefore, the binary variable

indicates whether the hexad $\sigma$ is informative or not.

Theorem (ref) confirms that conditional on $\sigma$ admitting one of the two identifying wirings, the probability of observing the first informative wiring is independent of $\mathrm{F}_\sigma$ and has the logistic form. It directly follows that the conditional probability of observing the second informative wiring is also independent of $\mathrm{F}_\sigma$.

The hexad logit estimator

In this section we propose a conditional likelihood estimator for $\beta_0$, the hexad logit estimator. For a generic hexad $\sigma=(i_1,j_1,k_1,i_2,j_2,k_2)$ define

noting that $W_\sigma = W_{\sigma,1} - W_{\sigma,2}$ as defined in Theorem (ref). Define, furthermore, the shorthand notation $p_{\sigma, c}(\beta) = \mathrm{P}(S_{\sigma,c} = 1 | S_\sigma = 1, X_\sigma)$ where $c\in\{1,2\}$ indexes the two informative wirings defined in the previous section. Therefore, remembering Theorem (ref), we have

The conditional likelihood function that combines information across all informative hexads is given by

We note that this is a composite log-likelihood function due to dependence across hexads that share common triads. Our proposed hexad logit estimator is given by

where $\mathcal B$ is a compact parameter space containing the true parameter $\beta_0$. Defining

the corresponding score function is given by

Also, the Hessian contribution for each hexad $\sigma$ is equal to

Derivations for the score and Hessian are provided in Appendix (ref).

Large sample theory

This section develops the large sample theory for the hexad logit estimator, deriving conditions on network sparsity under which it is consistent and asymptotically normal.

A quantity that plays an important role in the asymptotic properties of the hexad logit estimator is the unconditional probability of link formation, defined as

Asymptotic sparsity of the network is conceptualised by $\rho_N$ tending to zero as $N$ increases. Importantly, $\rho_N$ cannot go to zero faster than $O(1/N^2)$ or the average node degree $N^2 \rho_N$ will also converge to zero, rendering the network asymptotically empty.\footnote{See, for example, Assumption 4.(ii) of Graham17 for a corresponding requirement in dyadic link formation.}

A second key quantity is

the unconditional probability that a hexad $\sigma$ is informative. Recalling that informative wirings consist of four hyperedges, we have $p_N = O(\rho_N^4)$. This quantity directly controls the rate at which informative wirings accumulate and is therefore essential to the asymptotic behaviour of the hexad logit estimator.

We next introduce the assumptions that will be used in obtaining consistency and asymptotic normality.

Assumptions (ref), (ref) and (ref) are analogous to those made by Graham17 and Jochmans18 in dyadic link formation models. Compactness of regressors in Assumption (ref) is not essential and can be exchanged with appropriate assumptions on the moments of regressors. Assumption (ref) is standard. Assumption (ref) is a full-rank condition on the Hessian and guarantees that the objective function has a unique optimiser in the limit.

The dyad-specific fixed effects structure admits more complex forms of heterogeneity. This leads to fundamental differences with both dyadic and triadic models that involve node-specific fixed effects only. In particular, under the dyad-specific heterogeneity structure, even when the network is asymptotically non-empty, it can still be too sparse in the sense that in the limit one can accumulate infinitely many links but asymptotically zero informative wirings. This is a novel phenomenon, absent from models with dyad-level fixed effects. We formally establish and discuss these points in detail below, and provide a precise measure of when sparse becomes too sparse in terms of link formation probability $\rho_N$.

Consistency

This subsection presents the first key result of our large-sample analysis: the consistency of the hexad logit estimator under specific conditions on network sparsity.

This theorem formally establishes that for consistency to hold, link formation probability must tend to zero slower than $1/N^{3/2}$. In other words, in the given framework one needs something more than asymptotic non-emptiness of the network, which holds under the weaker condition that $N^2\rho_N$ remains above zero for $N$ sufficiently large. This is a fundamental deviation from the dyadic (and also triadic) model with node-level heterogeneity, as we discuss below.

The stronger condition on $\delta$ in Theorem (ref) is essentially caused by the more complex heterogeneity structure admitted by the dyad-specific effects. This translates into more complex informative wirings, involving a greater number of hyperedges and therefore forming with lower probability. Consequently, under the dyad-specific heterogeneity structure one requires more than asymptotic non-emptiness to collect sufficiently many informative wirings, leading to the requirement $0\leq \delta < 3/2$.

To formalise this intuition, we first investigate the structure of informative wirings for the dyadic and triadic models with node-specific heterogeneity. The dyadic link formation model is given by

which has already been considered in the literature. In particular, we know that in the bipartite case informative wirings consist of two edges only (Charbonneau17 and Jochmans18). The triadic link formation model with node-specific effects is given by

We develop the theory for this variant in Appendix (ref) and find that there are four informative wirings for this model, each of which consists of two hyperedges (see Section (ref) and Figure (ref)).

As opposed to both these models, we have already seen that in the triadic model with dyad-effects, informative wirings consist of four hyperedges. Wirings with a greater number of hyperedges form with (much) lower probability as the probability of wiring formation decreases exponentially with the number of hyperedges involved in the wiring. This is why the current model requires more than an asymptotically non-empty network: informative wirings accumulate relatively slowly so something denser than a barely non-empty network is required.

The reason for the specific threshold value of $3/2$ for $\delta$ is a novel pathological scenario under the triadic model with dyad-level heterogeneity: unless $\delta < 3/2$, as $N\to\infty$ the network will not produce infinitely many informative wirings in the limit even though there will be asymptotically infinitely many hyperedges. In fact, under $\delta >3/2$ there will be asymptotically zero informative wirings. This clearly is a problematic scenario, as collection of informative wirings is essential to convergence. Furthermore, this is a peculiarity of the model in (ref): as it turns out, the dyadic and triadic models with node-level heterogeneity automatically avoid such a scenario.

To establish these points formally, let $\widecheck{\rho}_N$ be the edge formation probability under the dyadic model (ref) and $\widetilde{\rho}_N$ be the hyperedge formation probability under the triadic model (ref). As defined earlier in (ref), $\rho_N$ is the corresponding probability for our main model in (ref). In what follows, a link refers to an edge in the dyadic model and a hyperedge in the triadic models. The expected number of links in the network is given by the product of the total number of dyads/triads and the relevant probability of link formation. This is $N^2\widecheck{\rho}_N$ for the dyadic model, $N^3\widetilde{\rho}_N$ for the triadic model with node-specific effects, and $N^3{\rho}_N$ for the triadic model with dyad-specific effects. The expected number of informative wirings, on the other hand, is proportional to the product of the total number of tetrads/hexads in the network and the probability of a generic tetrad/hexad being informative. As explained previously, under node-level heterogeneity an informative wiring consists of two edges/hyperedges. Then, the probability of observing an informative wiring is $O(\widecheck \rho _N^2)$ and $O(\widetilde \rho _N^2)$ under dyadic and triadic models with node-specific effects, respectively. In the triadic model with dyad-level effects, on the other hand, informative wirings form with probability $O(\rho_N^4)$. Noting that a dyadic model has $O(N^4)$ tetrads whereas a triadic model contains $O(N^6)$ hexads finally yields the asymptotic rates $O(N^4\widecheck \rho_N^2)$, $O(N^6\widetilde \rho_N^2)$, and $O(N^6 \rho_N^4)$. These calculations are summarised in Table (ref). We also state the average node degree for each setting.

Our calculations now clearly reveal that in the models with node-specific effects, the expected number of informative wirings is proportional to the square of the expected number of links. Therefore, informative wirings accumulate automatically with links: by design, one simply cannot have infinitely many links without also accumulating infinitely many informative wirings in the limit. Under the triadic link formation model with dyad-level effects, however, it is now clear that unless $\delta<3/2$, one will end up in the limit with infinitely many links but zero (or finite) informative wirings, and fail to achieve consistency.\footnote{The importance of having both infinitely many links and infinitely many informative wirings is directly reflected in the proof of Theorem (ref). See in particular equation (ref) and the surrounding discussion. For consistency to hold, the right hand side of equation (ref) must go to zero, and this happens only if the number of both the links and the informative wirings tend to infinity. The corresponding theoretical arguments for the dyadic model and the triadic model with node-level effects can be found in equations (ref) and (ref), respectively.} We also confirm that in the models with node-level effects, asymptotic non-emptiness (i.e. average node degree being above 0 for sufficiently large $N$) automatically guarantees accumulation of infinitely many informative wirings as $N\to\infty$. However, for the triadic model with dyad-level effects non-emptiness does not guarantee this.

Asymptotic normality

We now turn to the asymptotic distribution of the hexad logit estimator, establishing the conditions for its normality.

Theorem (ref) shows that the hexad logit estimator has the $\sqrt{N^3\rho_N}$ convergence rate. This is in line with the extant results in the literature, where the convergence rate is the root of the expected number of links across the network.

What is different is the stronger requirement $0\leq \delta <1$ on link formation probability. To understand why, first recall that the score is given by

This score function is analogous to a U-statistic, as previously noted in the literature for the score functions of dyadic link formation models. This analogy is helpful in obtaining the asymptotic normality of the score function by using results from the U-statistic literature. Importantly, however, under $\delta\geq 1$ the score $Z_N$ becomes akin to a highly degenerate U-statistic---and the degeneracy of a U-statistic typically leads to asymptotic non-normality.

To establish these points more formally, first consider the variance decomposition of the score, which is central to its asymptotic behaviour:

where

and ${\sum \sum}_{[\sigma,\sigma']_q}$ denotes summation over all hexad pairs $(\sigma,\sigma')$ with $q$ common nodes. This decomposition is derived and analysed in detail in Section (ref) in the Appendix.

That $\overline {\mathcal{C}}_{0,N}=\overline{\mathcal{C}}_{1,N}=\overline{\mathcal{C}}_{2,N}=0$ follows from the conditional independence of triad-specific shocks: hexad-pairs that share less than three common nodes cannot share a common triad and are therefore conditionally independent. More importantly, when $\overline{\mathcal{C}}_{3,N}$ is the leading term in (ref), a valid (triad-level) H\'{a}jek projection exists that is asymptotically equivalent to the score, and has a limiting Normal distribution. This proof approach is typical for U-statistics and is key to obtaining asymptotic normality of $\widehat \beta$.

However, our analysis also reveals that $\overline{\mathcal{C}}_{3,N}$ is not necessarily the leading term of the decomposition. Indeed, by equation (ref) in the Appendix we have,

It is straightforward to confirm that, as long as $\delta <2$, $\overline{\mathcal{C}}_{4,N}$ and $\overline{\mathcal{C}}_{5,N}$ are always dominated by $\overline{\mathcal{C}}_{3,N}$. However, $\overline{\mathcal{C}}_{6,N}$ dominates $\overline{\mathcal{C}}_{3,N}$ if $\delta > 1$. This would render $Z_N$ analogous to a degenerate U-statistic.\footnote{Typically, a (textbook) U-statistic would be called degenerate if the decomposition component for the covariance between terms with one common index ($\overline{\mathcal{C}}_{1,N}=0$ in our case) is zero. That we are able to obtain asymptotic normality under $\overline{\mathcal{C}}_{1,N}=\overline{\mathcal{C}}_{2,N}=0$ is a consequence of the random shock $\varepsilon_{ijk}$ being conditionally independent across triads. This triadic independence is what obtains asymptotic normality and therefore triadic information corresponds to what Serfling calls the “basic information”. As discussed in detail in Appendix (ref), the appropriate H\'{a}jek projection is therefore on the triadic information and the term $\overline{\mathcal{C}}_{3,N}$ captures the variance of this projection. Therefore, degeneracy in our case corresponds to $\overline{\mathcal{C}}_{3,N}=0$.} This usually leads to non-normal limiting distributions known as Gaussian chaos. While it is possible to pin down the resulting distribution, it is practically difficult to utilise, as a simple formula usually does not exist; see, e.g., Chapter 5.5.2 of Serfling and Chapter 12.3 of vdV. Consequently, we use $0\leq \delta<1$ as a sufficient condition for asymptotic normality.

The tension between the number of links and the number of identifying wirings that we discussed in the previous section is also at work here, and rather subtly. To see this, notice one can also write

which reveals the number of links and the number of informative wirings in the denominators of $\overline{\mathcal{C}}_{3,N}$ and $\overline{\mathcal{C}}_{6,N}$. For consistency, it was enough that both the number of links, $O(N^3 \rho_N)$, and number of informative wirings $O(N^6\rho_N^4)$, go to infinity. The above display reveals that for non-degeneracy we additionally need the informative wirings to accumulate faster than the number of links; otherwise, $\overline{\mathcal{C}}_{3,N}$ will not be the leading term. Intuitively, if the network is not too sparse (that is if $0\leq \delta <1$), the hyperedges we accumulate will be enough to construct an eventually even larger number of informative wirings. However, under $1\leq \delta <3/2$, the network is too sparse for sufficiently fast accumulation of informative hexads.

The dyadic and triadic models with node-level effects have no such degeneracy issue; see, in particular, Sections (ref) and (ref) for derivations of variance decompositions analogous to (ref) for these two models. From our discussion of Table (ref) we already know that these models do not admit a scenario where one can have many links but slow accumulation of informative wirings.

Simulation study

In this section, we study the finite sample performance of the hexad logit estimator. The data generating process follows the model in equation (ref), specified for a scalar covariate $X_{ijk}$ and parameter $\beta_0$:

The errors $\varepsilon_{ijk}$ follow a standard logistic distribution, the covariate $X_{ijk} \sim \mathcal{N}(0,1)$ is drawn independently across triads and independently of the errors and fixed effects, and the true parameter is set to $\beta_0 = 1$.

We mimic the simulation design in Jochmans18 and adapt it to the triadic setting. To induce degree heterogeneity and dependence driven by dyad-level fixed effects, each node $r \in \{1,\ldots,N\}$ is assigned a tendency parameter

so, lower indices imply higher tendency values. Sparsity is governed by a common scale parameter $c_N$ that depends on the sample size. We consider three regimes:

Dyad-specific fixed effects are then constructed by taking the negative of the average of node tendencies and applying the sparsity parameter:

This specification has three key features. First, the tendency profile induces degree heterogeneity: lower-indexed nodes have higher $\phi_r$ and therefore receive more negative fixed effects. Consequently, the lower a node is indexed, the less likely it is to form links. Second, larger values of the sparsity parameter $c_N$ yield sparser networks by reducing the overall probability of link formation. Third, setting $c_N = 0$ removes the fixed effects entirely, so link probabilities depend only on $X_{ijk}\beta_0$ and the logistic shock, producing the densest networks in our study.

We examine $N \in \{20, 25, 30, 40, 50\}$ and the three sparsity regimes defined above. Each design cell uses $K = 1000$ replications. Inference is based on a triad-clustered sandwich variance estimator that sums hexad scores over shared triads $(i,j,k)$. Convergence tolerance is set at $10^{-8}$ with a maximum of 1000 Newton-Raphson iterations.

Table (ref) reports, for each $(N, c_N)$ design cell, the following: the Monte Carlo mean and standard deviation of $\widehat{\beta}$; the Monte Carlo average of the estimated standard error $\text{se}(\widehat{\beta})$; the root mean squared error relative to $\beta_0 = 1$; the calibration ratio $\overline{\text{se}}/\text{SD}$ comparing the average estimated standard error to the empirical standard deviation; coverage at 90% and 95% based on normal approximations; and the rejection probability of $H_0: \beta = 0$ at the 5% level (power under $\beta_0 = 1$). We also provide Q-Q plots to compare the sample distribution of $\widehat \beta$ to the Normal distribution with mean and variance given by the Monte Carlo average and variance of $\widehat \beta$ (Figures (ref) and (ref)). In these Q-Q plots we also include the results for $c_N = \log N$, which is sparser than the settings introduced already.

Several observations stand out: First, especially across dense and log-log regimes, sampling distributions are tightly centered at $\beta_0$ with small standard deviations. As would be expected, dispersion increases as informativeness declines; that is, as $N$ decreases and/or the network becomes sparser. This is most pronounced at smaller sample sizes but attenuates as $N$ reaches 50. Second, the Q-Q plots confirm that the estimator is generally approximately normally distributed. One exception is the sparsest setting of $c_N=\ln N$, especially at the low sample size of $N=30$. However, we observe that the approximation improves quickly as $N$ increases to 50. This again highlights that information accumulation depends on both the sample size and the level of sparsity. Third, estimated standard errors are conservative in finite samples, especially when $N$ is small and the network is more sparse, yielding coverage at or above nominal levels. However, even under the more sparse settings, $\overline{\text{se}}/\text{SD}$ and the coverage rates improve quickly with $N$. Fourth, the test of $H_0: \beta = 0$ at the 5% level has practically unit rejection probability across all reported designs, reflecting strong power against the fixed alternative $\beta_0 = 1$.

All in all, the hexad logit estimator performs well even at moderately large $N$. Although confidence bands are relatively conservative, we note that they become tighter as $N$ increases. When $N$ is in the hundreds, the network will generate a large number of informative wirings (even under more sparse settings) which will lead to tighter confidence bands.

Conclusion

In this paper we introduced a triadic link formation model with dyad-specific effects and developed a likelihood estimator based on hexad subnetworks---the hexad logit estimator. Our theoretical analysis shows that the estimator is consistent when the unconditional link probability satisfies $\rho_N = O(N^{-\delta})$ with $0 \le \delta < 3/2$. Furthermore, asymptotic normality is obtained with rate $\sqrt{N^3 \rho_N}$ when $0 \le \delta < 1$. These requirements on $\rho_N$ go beyond asymptotic non-emptiness because, under dyad-level heterogeneity, accumulation of links need not lead to asymptotic accumulation of informative wirings. The stated sparsity thresholds provide precise conditions under which regular estimation and inference are attainable. These are novel features specific to the model at hand, which we do not observe in link formation models with node-level heterogeneity. We conjecture that models for higher dimensional graphs---which allow for more complicated forms of heterogeneity---will be subject to analogous and more severe sparsity constraints. A more general investigation of this is the subject of ongoing work.

Our findings also carry implications for applied research. While dyad-level effects offer a richer and often more realistic way to capture heterogeneity, we see that this flexibility imposes a limit on the level of sparsity. This can be an important consideration in granular or disaggregate datasets which are becoming increasingly common in economic analysis. Aggregating data along an appropriate dimension could be a practical first step in such situations.