2.1 ERGMs and TERGMs

The TERGM has an Exponential Random Graph Model (ERGM) at its core. As an exponential family model, the ERGM comes with well-known properties (Wasserman & Pattison Reference Wasserman and Pattison1996; Robins et al. Reference Robins, Snijders, Wang, Handcock and Pattison2007; Lusher et al. Reference Lusher, Koskinen and Robins2013). The ERGM defines the probability to observe a network based on endogenous network terms and exogenous covariates. The most commonly used statistics are counts of substructures, for example, the number of reciprocated ties or the number of in-stars of some order, which may be combined with nodal or dyadic attributes. The statistics that are typically used are based on principled assumptions about the dependencies among the ties (Frank & Strauss Reference Frank and Strauss1986; Snijders et al. Reference Snijders, Pattison, Robins and Handcock2006). The expected prevalence of these statistics is determined by a statistical parameter. The probability to observe the realization x of a network is given by

\begin{align*}p_{{\rm ERGM}}\left(X=x\right)={\kappa }^{-1}{\rm exp}\left(\sum_k{{\theta }_ks_k(x)}\right),\end{align*}

where X is the random network state, $\theta $ is a statistical parameter, and s(x) is a vector of statistics describing the network; $\kappa $ is a normalizing constant.

For later use, we elaborate the difference between endogenous and exogenous effects. We represent the network by its adjacency matrix x = ( ${x}_{ij}$ ), where ${x}_{ij}$ is the indicator of the tie from node i to node j. An effect is exogenous if it is a linear function of the ${x}_{ij}$ , and endogenous if it is not. Dependence considerations for networks (see the cited ERGM literature) lead to effects depending on subgraph counts, which are equivalent to products of tie indicators. Examples are counts of two-stars, representing degree variability, and counts of triangles or other triadic configurations, representing transitivity. Endogenous effects imply emergence in the evolution of networks and imply dependence between the tie indicator variables, which is essential for representing a network by a statistical model. In contrast, the impact of nodal and dyadic covariates is usually represented by exogenous effects.

The TERGM is a model for network panel data where the t’th network observation is modeled by an ERGM in which functions of the preceding network observations can enter as exogenous variables. Various ways to specify this have been proposed (Robins & Pattison Reference Robins and Pattison2001; Hanneke et al. Reference Hanneke, Fu and Xing2010; Desmarais & Cranmer Reference Desmarais and Cranmer2012; Krivitsky & Handcock Reference Krivitsky and Handcock2014). The TERGM for two waves at times t-1 and t can be represented by the conditional probability function

(1) \begin{align}p_{{\rm TERGM}}(X(t)&=x(t)|X(t-1)=x(t-1)) \nonumber\\&={\kappa }^{-1}{\rm exp}\left(\sum_k{{\theta }_ks_k(x(t))}+\sum_h{{\theta }_hz_h}(x(t),x(t-1))\right)\end{align}

where x(t) and x(t-1) are the realizations, X(t) and X(t-1) are the random networks, s(x(t)) is a vector of statistics for network x(t) like above, and z(x(t), x(t-1)) is a vector of statistics of both networks. For the parameter $\theta $ , we denote by $\theta {}_{k}$ those pertaining to s(x(t)) and by $\theta {}_{h}$ those pertaining to z(x(t), x(t-1)). The difference with the standard ERGM lies in the extra statistics ${Z}_{h}$ (x(t),x(t-1)) that are memory terms depending on the networks at time t as well as time t–1. These are usually exogenous effects, that is, linear functions of the tie statistics ${x}_{ij}$ (t). A basic example of such a statistic representing the match between the two consecutive observations is the number of identical tie variables,

\begin{align*}z_h(x(t),x(t-1))=\sum_{ij}{\left(1-\left|x_{ij}(t)-x_{ij}(t-1)\right|\right)},\end{align*}

called the dyadic stability term by Leifeld & Cranmer (2019 a); it can be rewritten as a linear function of $x_{ij}(t)$ , modeling inertia.

2.2 The model specification in Leifeld & Cranmer’s (2019a) TERGM

We now evaluate Leifeld & Cranmer’s (2019a) model specification based on the replication script empirical.R (Leifeld & Cranmer Reference Leifeld and Cranmer2019b). Leifeld & Cranmer (Reference Leifeld and Cranmer2019a) use the btergm package for all analyses, which uses functionality for the established “ergm” and “RSiena” package in R (Hunter et al. 2008; Ripley et al. Reference Ripley, Snijders, Boda, Vörös and Preciado2021). Leifeld & Cranmer (Reference Leifeld and Cranmer2019a) state that they use “the same specification” (p. 42) in their TERGM as in the SAOM and give the statistics used in Equations (15)–(21) in their article. The ERGM equivalents of the SAOM terms “in-degree popularity,” “out-degree popularity,” and “out-degree activity” (Equations (19)–(21)) are “in-2-stars,” “mixed 2-paths,” and “out-2-stars” (Frank & Strauss Reference Frank and Strauss1986). The replication script allows to identify how they have translated Equations (19)–(21) used as endogenous effects in the SAOM to statistics used in the TERGM.

The specification can be found on lines 243–247 of the empirical script:

The three terms <monospace>nodeicov(“idegsqrt”) + nodeicov(“odegsqrt”) + nodeocov (“odegsqrt”)</monospace> are meant to correspond to the endogenous effects “in-degree popularity,” “out-degree popularity,” and “out-degree activity”, in the SAOM. However, rather than endogenous network terms, “<monospace>idegsqrt</monospace>” and “<monospace>odegsqrt</monospace>” refer to exogenous vertex attributes representing a square root transformation of nodes’ indegree and outdegree in the dependent variable, respectively, as defined in lines 222–230 of the script:

This means that the terms <monospace>nodeicov(“idegsqrt”) + nodeicov(“odegsqrt”) + nodeocov(“odegsqrt”)</monospace> model (i) the tendency of nodes with high observed indegree to receive ties, (ii) the tendency of nodes with high observed outdegree to receive ties, and (iii) the tendency of nodes with high observed outdegree to send ties, respectively. The first and third of these express, respectively, that the indegrees predicted in the model are similar to observed indegrees, and that outdegrees predicted in the model are similar to observed outdegrees. These are tautological terms, by definition receiving high parameter estimates in networks. The second term expresses that modeled indegrees are similar to observed outdegrees. This is a circular term at the level of networks, but not a direct tautology, and the resulting parameter estimate is small. An intuitive way of conceiving of what is happening is that a crucial summary—vertex degrees—of the observed dependent network is used to predict this same feature of the dependent network.

This specification means that the model is not of the form (1), but of the form

(2) $$\matrix{ {{p_{{\rm{TERGM}};\;{\rm{LC}}}}(X(t)} \hfill & { = x(t)|X(t - 1) = x(t - 1);\;{x_{{\rm{obs}}}})} \hfill \cr {} \hfill & { = {\kappa ^{ - 1}}{\rm{exp}}\left( {\sum\limits_k {{\theta _k}{s_k}(x(t))} + \sum\limits_h {{\theta _h}{z_h}} (x(t),x(t - 1))} \right.} \hfill \cr {} \hfill & {\quad \left. { + \sum\limits_l {{\phi _l}} {u_l}(x(t),{x_{{\rm{obs}}}}(t))} \right).} \hfill \cr } $$

Added to the model (1) is a further set of statistics, here denoted $u_l(x(t),x_{{\rm obs}}(t))$ , that is dependent on a particular realization x ${}_{obs}$ (t) of the networks. A positive parameter associated with these statistics $u_l(x(t),\,x_{{\rm obs}}(t))$ , expressing an aspect of the similarity between x(t) and x ${}_{obs}$ (t), implies that any realization of the network that is closer to the actual observation according to this match has a higher probability to be observed. Thus, the appropriate endogenous specification (1) is turned into circularity. This gives no indication about the salient features of a network, but only artificially improves the fit of the model to this particular dataset without uncovering any dependencies. Leifeld & Cranmer’s (2019a) analysis uses three such terms, depending on the observed indegree and outdegree sequences. These are exogenous terms treating the observed network as a covariate for itself. Taken to the extreme, in this spirit we could include the entire observed network as a dyadic covariate among the predictors of the network—this would lead to perfect fit of the model, but would have no explanatory value whatsoever.

Moreover, as the degree sequence is used in the estimation of parameters, it also impacts the parameter estimates of other statistics directly through the strong correlations between statistics in network models. Thus, not only are the parameters $\phi_l$ associated with statistics $u_l(x(t),x_{{\rm obs}}(t))$ tautological, but also all other model parameters ${\theta }_k$ and ${\theta }_h$ will be systematically distorted. In the analysis by Leifeld & Cranmer (Reference Leifeld and Cranmer2019a), this can be seen in the astonishing finding that there is no support for transitivity in friendship dynamics. As we show in the following subsection, we find clear support for transitivity when replacing the circular exogenous terms with conventional endogenous terms, in line with the very vast literature on friendship networks.

2.3 Replication with correct ERGM terms

In this section, we replicate the analysis by Leifeld & Cranmer (2019 a) using conventional ERGM terms instead of the artificially exogenous variables ^{Footnote 3} . A model that would be the direct analogue of the SAOM specification has the three star-statistics, transitive and cyclic triples, all of which are sufficient statistics in a Markov graph derived from first principles and assumptions about dependence among the ties (Frank & Strauss Reference Frank and Strauss1986). This Markov model cannot be estimated, something which is to be expected for most networks ^{Footnote 4} . While the pseudo-maximum likelihood estimation (MLE) can be determined—that is, pseudo-MLE will produce parameter estimates—these estimates correspond to a degenerate model ^{Footnote 5} .

Following standard practice (Snijders et al. Reference Snijders, Pattison, Robins and Handcock2006; Hunter & Handcock Reference Hunter and Handcock2006; Lusher et al. Reference Lusher, Koskinen and Robins2013), we have to replace the star statistics and the triadic statistics with non-degenerate ones, in particular their equivalent, geometrically weighted statistics, derived from another set of principled dependence assumptions (Robins et al. Reference Robins, Pattison and Wang2009), in order to specify a model that can be estimated. We substitute the terms <monospace>nodeicov(“idegsqrt”)</monospace>, <monospace>nodeicov(“odegsqrt”)</monospace>, and <monospace>nodeocov(“odegsqrt”)</monospace> with the ERGM terms that model degree dispersion, in particular geometrically weighted in-stars, two-paths, and geometrically weighted out-stars (statnet terms <monospace>gwidegree, twopath</monospace>, and <monospace>gwodegree</monospace>). Further, we substitute the triadic terms transitive triplets (<monospace>ttriple</monospace> and <monospace>ctriple</monospace>) with their geometrically weighted versions (<monospace>dgwesp(type = “OTP”)</monospace> and <monospace>dgwesp(type = “ITP”)</monospace>).

In the SAOM analysis, we also substitute the terms for transitive and cyclic triplets with the geometrically weighted versions. This is done for two reasons. First, even though model terms in ERGM and SAOM analyses are necessarily different in their dependence assumptions (Block et al. Reference Block, Stadtfeld and Snijders2019), using geometrically weighted versions in both makes the results more comparable. Second, the specification in the tutorial article from 2010 is not “canonical” as claimed by Leifeld & Cranmer (2019 a) (p. 42), and so a more contemporary specification should be used for both models.

The results from both models, presented in Table 1, are remarkably similar. The “Transitive ties” parameter and the “GWESP Transitive” parameter need to be interpreted together, as both model the same tendency of transitive closure, but with different functional forms. The combination of parameters shows that both models find a strong tendency toward friendships being transitive. No endogenous sorting of degrees is found, while the impact of the exogenous covariates “Same primary class” and sex is in the same direction. The main difference is that the “GWESP cyclic” term is significant in the TERGM analysis but not in the SAOM analysis. This will be related to the different formulations of the various GWESP and transitivity terms in this ERGM-type and SAOM-type model as outlined in Block et al. (Reference Block, Stadtfeld and Snijders2019), but some differences are to be expected because the models are different (see literature cited in footnote 1). We conclude that the substantive insights we can draw from either model are not very different, but that understanding how these differences come about is a more complex task than attributing this to simple differences in model fit.

Table 1. Results of TERGM and SAOM analysis using standard endogenous ERGM terms

Notes: All analyses performed using standard best practises. Significance levels: * = 0.05; ** = 0.01; *** = 0.001; ${}^\circ$ = not tested.

Notes: All analyses performed using standard best practises. Significance levels: * = 0.05; ** = 0.01; *** = 0.001; ${}^\circ$ = not tested.

2.4 Out-of-sample prediction

Out-of-sample analysis plays a major role in Leifeld & Cranmer’s (2019 a) results and final recommendations. In out-of-sample analysis, a dataset is split into training data and test data. Leifeld & Cranmer (Reference Leifeld and Cranmer2019a) choose the first three waves of the Knecht network data as training data and the test data is the fourth wave. Out-of-sample analysis estimates a model using the training data, after which data beyond the training data are simulated based on those estimates and compared to the test data. Importantly, the model should use no information from the test data to generate the predictions. However, the “out-of-sample” prediction using Leifeld & Cranmer’s (2019a) TERGM specification violates this principle.

We saw above that in the estimation of parameters for each wave the empirically observed degree sequences for waves at the end of the period were extracted and used as exogenous nodal covariates in the model of change for the period. The same was done to generate likely outcomes of wave 4 (out-of-sample prediction).

One of the wrapper functions of the btergm package is the “<monospace>gof</monospace>” function that facilitates simulating out-of-sample predictions for a specified model. This “<monospace>gof</monospace>” function is prominently applied in section 4 of this article for the out-of-sample comparison. The earlier error is repeated here in lines 254–260:

The crucial part here is “<monospace>formula = friendship[3:4] </monospace>”, which tells the algorithm to use the “covariates” stored in the network object <monospace>friendship[[4]]</monospace> in the simulations. However, these covariates include square root transformed information about observed indegrees and outdegrees at wave 4. The “<monospace>gof</monospace>” function as used here thus simulates a model using future information (the future in- and out-degree of nodes). The same circularity as before is even more evident here; part of the test data is used to predict the test data. This mistake appears to be encoded in the <monospace>btergm</monospace> package that is used in Leifeld & Cranmer’s (2019a) analysis.

This analysis cannot then be regarded as out-of-sample testing. This circularity has major consequences for the article results. The comparison case, that is, the out-of-sample predictions of the SAOM, makes no use of any information from the fourth wave, biasing the comparison toward the TERGM.

2.5 Replication with correctly specified ERGM terms

It would in principle be possible that a properly specified TERGM could still perform substantially better in out-of-sample prediction, even if the substantive conclusions do not differ much. To test this, we used the TERGM estimates obtained as outlined above.

The results using the identical metrics for comparison as employed by Leifeld & Cranmer (Reference Leifeld and Cranmer2019a) are presented in Figure 1. The GOFs for the auxiliary statistics “Edge-wise shared partners,” “Dyad-wise shared partners,” “Indegree distribution,” and “Geodesic distances” show no clear trend favoring either model. The ROC and PR curves are, for what they are worth, very similar between the SAOM and the re-specified TERGM. In sum, these analyses indicate no gain in predictive value between a correctly specified TERGM and the SAOM corroborating earlier findings and theoretical understanding of such processes (Block et al. Reference Block, Koskinen, Hollway, Steglich and Stadtfeld2018). ^{Footnote 6}

Figure 1. Replication of Figures 4 and 5 from Leifeld & Cranmer’s (2019a) article with ERGM terms in the model that do not use future information. This shows that the claim that there is a substantial margin in performance differences is false.

In sum, it was on the basis of a supposed performance margin that Leifeld & Cranmer (Reference Leifeld and Cranmer2019a) recommend researchers estimate both models and use the wrapper functions provided in btergm to identify the one with the better predictive fit. However, we have seen that this comparison based on the btergm functionality and thus this conclusion is fundamentally undermined by their repeated error.