Social Network Analysis R Tutorial Part 10: Longitudinal Network Analysis

R
SNA
Author

Shuhan (Alice) Ai

Published

December 15, 2025

So far in this tutorial series, we’ve explored various aspects of social network analysis using cross-sectional data—networks captured at a single point in time. However, many social networks are dynamic, evolving over time as relationships form, dissolve, and transform. In this tutorial, we’ll explore longitudinal network analysis, focusing on modeling networks that change across multiple time periods.

After this tutorial, we can answer questions like:

1. Introduction to Longitudinal Networks

1.1 What are Longitudinal Networks?

Longitudinal networks (also called temporal networks or dynamic networks) are networks observed at multiple time points. Instead of a single snapshot at time t, we observe a sequence of snapshots: \(N_1, N_2, ..., N_T\), where each \(N_t\) represents the network at time t.

1.2 Why Model Temporal Dependencies?

A principle in longitudinal network modeling is that the probability of observing the network at time t depends on the network structure at time t-1. This makes intuitive sense:

  • Existing friendships are more likely to persist than new friendships to form randomly
  • Reciprocity takes time to develop (if A likes B at t-1, B is more likely to like A at t)
  • Network structures like transitivity (friends of friends becoming friends) unfold over time

This conditional dependence is expressed mathematically as:

\[P(N_t | N_{t-1}, N_{t-2}, ..., N_{t-K}, \theta)\]

where \(K\) represents the memory or lag of the model—how many previous time steps influence the current network.

1.3 Temporal ERGMs (TERGMs)

The Temporal ERGM (TERGM) extends the cross-sectional ERGM framework to longitudinal data. For a single time point t, the TERGM is:

\[P(N_t | N_{t-K}, ..., N_{t-1}, \theta) = \frac{\exp(\theta^T h(N_t, N_{t-1}, ..., N_{t-K}))}{\sum_{n \in \mathcal{N}} \exp(\theta^T h(n, N_{t-1}, ..., N_{t-K}))}\]

For multiple time periods (from K+1 to T), we model the joint probability as:

\[P(N_{K+1}, ..., N_T | N_1, ..., N_K, \theta) = \prod_{t=K+1}^{T} P(N_t | N_{t-K}, ..., N_{t-1}, \theta)\]

This assumes that conditional on the previous K networks, each time point is independent.

2. Memory Terms in TERGMs

Memory terms capture temporal processes and are essential for understanding network dynamics.

2.1 Common Memory Terms

Positive Autoregression (Edge Stability)
Counts edges that persist from t-1 to t:
\[h_a = \sum_{ij} N_t^{ij} N_{t-1}^{ij}\]

Interpretation: “Do existing ties tend to remain?”

Dyadic Stability
Counts both persistent edges AND persistent non-edges:
\[h_s = \sum_{ij} [N_t^{ij} N_{t-1}^{ij} + (1-N_t^{ij})(1-N_{t-1}^{ij})]\]

Interpretation: “Do dyadic relationships (present or absent) remain stable?”

Delayed Reciprocity
Counts reciprocated edges with a time lag:
\[h_r = \sum_{ij} N_t^{ij} N_{t-1}^{ji}\]

Interpretation: “If A nominated B at t-1, does B nominate A at t?”

3. The btergm Package

3.1 Installation and Loading

Code 
install.packages("xergm")
install.packages("btergm")
install.packages("statnet")
install.packages("texreg")
Code 
library(statnet)
library(xergm)
library(btergm)
library(texreg)
library(dplyr)

set.seed(12345)

3.2 Why MPLE with Bootstrap?

The btergm package uses Maximum Pseudolikelihood Estimation (MPLE) with bootstrap confidence intervals:

  • Speed: No simulation required for estimation
  • Consistency: MPLE is consistent as either n (nodes) or T (time periods) increases
  • Scalability: Works well with large networks

The bootstrap procedure corrects for MPLE’s underestimation of standard errors by:

  1. Sampling with replacement from the T-K observed networks
  2. Re-estimating the model for each bootstrap sample
  3. Using the distribution of bootstrap estimates for confidence intervals

The following two examples were extracted from this paper.

Leifeld, P., Cranmer, S. J., & Desmarais, B. A. (2018). Temporal exponential random graph models with btergm: Estimation and bootstrap confidence intervals. Journal of Statistical Software, 83, 1-36.

4. Simple Example: International Alliances

Before diving into complex temporal models with changing actor sets, let’s start with a simpler example: the network of international defense alliances.

4.1 Load Alliance Data

Code 
# Load the alliance dataset
data("alliances", package = "xergm.common")

# Examine the data structure
cat("Number of time periods:", length(allyNet), "\n")
Number of time periods: 20
Code 
cat("Number of countries:", network.size(allyNet[[1]]), "\n")
Number of countries: 164

This dataset contains:

  • allyNet: List of 20 network objects (1981-2000), each with 164 countries
  • contigMat: Binary matrix indicating shared borders (constant over time)
  • LSP: Time-varying matrix of shared partners from previous year
  • warNet: Time-varying matrices of militarized interstate disputes
  • Node attributes: cinc (military capabilities), polity (regime authority)

4.2 Visualize Alliance Networks (last three)

Code 
par(mfrow = c(1, 3), mar = c(0, 0, 1, 0))
for (i in (length(allyNet) - 2):length(allyNet)) {
  plot(allyNet[[i]], main = paste("t =", i),
       vertex.cex = 0.7, edge.col = "grey50")
}

4.3 Model 1a: Pooled TERGM (No Temporal Dependencies)

Let’s start with a pooled TERGM that treats all 20 time periods as independent.

“In a first simple model, a TERGM without any temporal dependencies is estimated. This corresponds to a pooled ERGM – pooling across T networks that are assumed to be independent of each other – where the estimates reflect the average effects across the 20 time points. The (somewhat unrealistic) assumption here is that the consecutive networks are independent from each other.”

Code 
model_1a <- btergm(
  allyNet ~ edges + 
    gwesp(0, fixed = TRUE) +            # Geometrically weighted ESP
    edgecov(LSP) +                      # Shared partners (from t-1)
    edgecov(warNet) +                   # Militarized disputes
    nodecov("polity") +                 # Regime authority
    nodecov("cinc") +                   # Military capabilities
    absdiff("polity") +                 # Absolute difference in polity
    absdiff("cinc") +                   # Absolute difference in capabilities
    edgecov(contigMat),                 # Shared border
  R = 50,                               # Bootstrap replications
  parallel = "no",
  verbose = FALSE
)

summary(model_1a)
                         Estimate  Boot mean    2.5%   97.5%
edges                  -5.0498497 -5.0433724 -5.1925 -4.9231
gwesp.fixed.0           1.6242771  1.6191803  1.5700  1.6835
edgecov.LSP[[i]]        2.1065906  2.1180989  1.8509  2.3401
edgecov.warNet[[i]]     0.2866346  0.2827274  0.1897  0.3585
nodecov.polity         -0.0050764 -0.0052387 -0.0122  0.0012
nodecov.cinc            8.3308516  8.1533635  3.6648 12.4279
absdiff.polity         -0.1257948 -0.1266169 -0.1422 -0.1132
absdiff.cinc           -4.8494093 -4.6991871 -8.7845 -0.7636
edgecov.contigMat[[i]]  3.2145679  3.2165368  3.0952  3.3212
  • edges: Negative coefficient indicates sparse network baseline
  • gwesp.fixed.0: Positive coefficient indicates tendency toward triadic closure
  • edgecov.LSP: Positive and significant coefficient indicates shared partners matter
  • edgecov.warNet: Positive coefficient indicates countries in disputes form alliances
  • absdiff.polity: Negative coefficient indicates similar regime types more likely to ally
  • edgecov.contigMat: Positive and significant coefficient indicates neighbors form alliances

4.4 Goodness-of-Fit for Model 1a

Code 
gof_1a <- gof(model_1a, nsim = 50, 
              statistics = c(esp, geodesic, deg))
                t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 t=11 t=12 t=13 t=14
allyNet (row)   164 164 164 164 164 164 164 164 164  164  164  164  164  164
allyNet (col)   164 164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (row)       164 164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (col)       164 164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (row)    164 164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (col)    164 164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (row) 164 164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (col) 164 164 164 164 164 164 164 164 164  164  164  164  164  164
                t=15 t=16 t=17 t=18 t=19 t=20
allyNet (row)    164  164  164  164  164  164
allyNet (col)    164  164  164  164  164  164
LSP (row)        164  164  164  164  164  164
LSP (col)        164  164  164  164  164  164
warNet (row)     164  164  164  164  164  164
warNet (col)     164  164  164  164  164  164
contigMat (row)  164  164  164  164  164  164
contigMat (col)  164  164  164  164  164  164
                t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 t=11 t=12 t=13 t=14
allyNet (row)   164 164 164 164 164 164 164 164 164  164  164  164  164  164
allyNet (col)   164 164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (row)       164 164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (col)       164 164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (row)    164 164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (col)    164 164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (row) 164 164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (col) 164 164 164 164 164 164 164 164 164  164  164  164  164  164
                t=15 t=16 t=17 t=18 t=19 t=20
allyNet (row)    164  164  164  164  164  164
allyNet (col)    164  164  164  164  164  164
LSP (row)        164  164  164  164  164  164
LSP (col)        164  164  164  164  164  164
warNet (row)     164  164  164  164  164  164
warNet (col)     164  164  164  164  164  164
contigMat (row)  164  164  164  164  164  164
contigMat (col)  164  164  164  164  164  164
Code 
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
plot(gof_1a)

Poor fit: The model systematically misses the observed network structure, particularly in:

  • Edge-wise shared partners
  • Geodesic distances
  • Degree distribution

This suggests we need temporal dependencies.

4.5 Model 1b: TERGM with Temporal Dependencies

Now add temporal terms to capture network dynamics:

Code 
model_1b <- btergm(
  allyNet ~ edges + 
    gwesp(0, fixed = TRUE) +
    edgecov(LSP) +
    edgecov(warNet) +
    nodecov("polity") +
    nodecov("cinc") +
    absdiff("polity") +
    absdiff("cinc") +
    edgecov(contigMat) +
    memory(type = "stability") +                # Dyadic stability (NEW)
    timecov(transform = function(t) t) +        # Linear time trend (NEW)
    timecov(warNet, transform = function(t) t), # War × time interaction (NEW)
  R = 50,
  parallel = "no",
  verbose = FALSE
)

summary(model_1b)
                              Estimate Boot mean     2.5%    97.5%
edges                        -3.139220 -3.167361  -4.5438  -2.4270
gwesp.fixed.0                 1.790187  1.774662   1.3700   1.9787
edgecov.LSP[[i]]              0.113878  0.141360   0.0590   0.5084
edgecov.warNet[[i]]          -1.137333 -1.005569  -2.4126   0.5632
nodecov.polity               -0.015363 -0.014300  -0.0860   0.0622
nodecov.cinc                 -0.163806  0.761174 -18.4340  16.6924
absdiff.polity               -0.068772 -0.077040  -0.1753  -0.0233
absdiff.cinc                  3.274009  1.526167 -12.4441  14.7695
edgecov.contigMat[[i]]        1.597682  1.515217   0.4001   2.4837
edgecov.memory[[i]]           5.097553  5.407778   4.6339   7.5021
edgecov.timecov1[[i]]         0.072603  0.084207   0.0079   0.2789
edgecov.timecov2.warNet[[i]]  0.054838  0.047244  -0.1070   0.1320

4.6 New Temporal Terms Explained

memory(type = "stability"):
Captures whether dyadic relationships (alliances or lack thereof) persist over time. Positive coefficient indicates high stability—existing alliances tend to continue.

timecov(transform = function(t) t):
Tests for a linear time trend in alliance density. Positive coefficient means the network gets denser over time.

timecov(warNet, transform = function(t) t):
Tests whether the effect of militarized disputes on alliances changes over time. Non-significant here, suggesting the war-alliance relationship is stable across years.

4.7 Compare Alliance Models

Code 
screenreg(list(model_1a, model_1b),
          custom.model.names = c("Model 1a", "Model 1b"),
          single.row = TRUE)

=====================================================================================
                              Model 1a                    Model 1b                   
-------------------------------------------------------------------------------------
edges                             -5.05 [-5.19; -4.92] *      -3.14 [ -4.54; -2.43] *
gwesp.fixed.0                      1.62 [ 1.57;  1.68] *       1.79 [  1.37;  1.98] *
edgecov.LSP[[i]]                   2.11 [ 1.85;  2.34] *       0.11 [  0.06;  0.51] *
edgecov.warNet[[i]]                0.29 [ 0.19;  0.36] *      -1.14 [ -2.41;  0.56]  
nodecov.polity                    -0.01 [-0.01;  0.00]        -0.02 [ -0.09;  0.06]  
nodecov.cinc                       8.33 [ 3.66; 12.43] *      -0.16 [-18.43; 16.69]  
absdiff.polity                    -0.13 [-0.14; -0.11] *      -0.07 [ -0.18; -0.02] *
absdiff.cinc                      -4.85 [-8.78; -0.76] *       3.27 [-12.44; 14.77]  
edgecov.contigMat[[i]]             3.21 [ 3.10;  3.32] *       1.60 [  0.40;  2.48] *
edgecov.memory[[i]]                                            5.10 [  4.63;  7.50] *
edgecov.timecov1[[i]]                                          0.07 [  0.01;  0.28] *
edgecov.timecov2.warNet[[i]]                                   0.05 [ -0.11;  0.13]  
-------------------------------------------------------------------------------------
Num. obs.                     534640                      507908                     
=====================================================================================
* Null hypothesis value outside the confidence interval.

4.8 Key Findings from Model Comparison

  1. Memory dominates: memory(stability) has the strongest effect—alliance relationships are highly persistent

  2. Effect changes:

    • edgecov.LSP: Coefficient drops dramatically (from 2.11 to 0.11) once we control for stability
    • edgecov.warNet: Becomes non-significant in TERGM
    • edgecov.contigMat: Still significant but weaker (from 3.21 to 1.60)

  3. Time trend: Positive and significant—international alliance network becomes slightly denser over time

  4. Interpretation: Without controlling for temporal stability, we overestimate the importance of shared partners and conflicts. The real driver is inertia—past alliances strongly predict future alliances.

4.9 Goodness-of-Fit for Model 1b

Code 
gof_1b <- gof(model_1b, nsim = 50,
              statistics = c(esp, geodesic, deg))
 t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 t=11 t=12 t=13 t=14 t=15 t=16 t=17
   1   2   3   4   5   6   7   8   9   10   11   12   13   14   15   16   17
 t=18 t=19 t=20
   18   19   20
 t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 t=11 t=12 t=13 t=14 t=15 t=16 t=17
   1   2   3   4   5   6   7   8   9   10   11   12   13   14   15   16   17
 t=18 t=19 t=20
   18   19   20
                      t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 t=11 t=12 t=13 t=14
allyNet (row)         164 164 164 164 164 164 164 164  164  164  164  164  164
allyNet (col)         164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (row)             164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (col)             164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (row)          164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (col)          164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (row)       164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (col)       164 164 164 164 164 164 164 164  164  164  164  164  164
memory (row)          164 164 164 164 164 164 164 164  164  164  164  164  164
memory (col)          164 164 164 164 164 164 164 164  164  164  164  164  164
timecov1 (row)        164 164 164 164 164 164 164 164  164  164  164  164  164
timecov1 (col)        164 164 164 164 164 164 164 164  164  164  164  164  164
timecov2.warNet (row) 164 164 164 164 164 164 164 164  164  164  164  164  164
timecov2.warNet (col) 164 164 164 164 164 164 164 164  164  164  164  164  164
                      t=15 t=16 t=17 t=18 t=19 t=20
allyNet (row)          164  164  164  164  164  164
allyNet (col)          164  164  164  164  164  164
LSP (row)              164  164  164  164  164  164
LSP (col)              164  164  164  164  164  164
warNet (row)           164  164  164  164  164  164
warNet (col)           164  164  164  164  164  164
contigMat (row)        164  164  164  164  164  164
contigMat (col)        164  164  164  164  164  164
memory (row)           164  164  164  164  164  164
memory (col)           164  164  164  164  164  164
timecov1 (row)         164  164  164  164  164  164
timecov1 (col)         164  164  164  164  164  164
timecov2.warNet (row)  164  164  164  164  164  164
timecov2.warNet (col)  164  164  164  164  164  164
                      t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9 t=10 t=11 t=12 t=13 t=14
allyNet (row)         164 164 164 164 164 164 164 164  164  164  164  164  164
allyNet (col)         164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (row)             164 164 164 164 164 164 164 164  164  164  164  164  164
LSP (col)             164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (row)          164 164 164 164 164 164 164 164  164  164  164  164  164
warNet (col)          164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (row)       164 164 164 164 164 164 164 164  164  164  164  164  164
contigMat (col)       164 164 164 164 164 164 164 164  164  164  164  164  164
memory (row)          164 164 164 164 164 164 164 164  164  164  164  164  164
memory (col)          164 164 164 164 164 164 164 164  164  164  164  164  164
timecov1 (row)        164 164 164 164 164 164 164 164  164  164  164  164  164
timecov1 (col)        164 164 164 164 164 164 164 164  164  164  164  164  164
timecov2.warNet (row) 164 164 164 164 164 164 164 164  164  164  164  164  164
timecov2.warNet (col) 164 164 164 164 164 164 164 164  164  164  164  164  164
                      t=15 t=16 t=17 t=18 t=19 t=20
allyNet (row)          164  164  164  164  164  164
allyNet (col)          164  164  164  164  164  164
LSP (row)              164  164  164  164  164  164
LSP (col)              164  164  164  164  164  164
warNet (row)           164  164  164  164  164  164
warNet (col)           164  164  164  164  164  164
contigMat (row)        164  164  164  164  164  164
contigMat (col)        164  164  164  164  164  164
memory (row)           164  164  164  164  164  164
memory (col)           164  164  164  164  164  164
timecov1 (row)         164  164  164  164  164  164
timecov1 (col)         164  164  164  164  164  164
timecov2.warNet (row)  164  164  164  164  164  164
timecov2.warNet (col)  164  164  164  164  164  164
Code 
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
plot(gof_1b)

Much better fit: The TERGM with temporal dependencies captures the observed network structure more accurately across all diagnostics.

5. Advanced Example: Friendship Networks in a Dutch School

Now let’s tackle a more complex example with changing actor composition.

5.1 Load and Explore Data

We’ll use the Knecht friendship network dataset, which tracks friendship nominations among 26 students across four waves.

Code 
data("knecht", package = "xergm.common")

cat("Number of time periods:", length(friendship), "\n")
Number of time periods: 4
Code 
cat("Network dimensions at each wave:\n")
Network dimensions at each wave:
Code 
for (i in 1:length(friendship)) {
  cat("  Wave", i, ":", dim(friendship[[i]]), "\n")
}
  Wave 1 : 26 26 
  Wave 2 : 26 26 
  Wave 3 : 26 26 
  Wave 4 : 26 26

The data includes:

  • friendship: List of 4 adjacency matrices (one per wave)
  • demographics: Data frame with student attributes (sex, etc.)
  • primary: Matrix indicating shared primary school attendance

5.2 Data Preparation

Code 
#label rows and columns with node IDs for tracking
for (i in 1:length(friendship)) {
  rownames(friendship[[i]]) <- 1:nrow(friendship[[i]])
  colnames(friendship[[i]]) <- 1:ncol(friendship[[i]])
}

#label the primary school covariate
rownames(primary) <- rownames(friendship[[1]])
colnames(primary) <- colnames(friendship[[1]])

#extract and label sex variable
sex <- demographics$sex
names(sex) <- 1:length(sex)

Handle Missing Data:

Code 
#remove nodes that dropped out (coded as 10)
friendship <- handleMissings(friendship, na = 10, method = "remove")

#impute missing values (unit non-response) with modal value (0)
friendship <- handleMissings(friendship, na = NA, method = "fillmode")

Node 21 dropped out after wave 2, so waves 3-4 have only 25 students.

Convert to Network Objects:

Code 
for (i in 1:length(friendship)) {
  #adjust sex variable to match current network dimensions
  s <- adjust(sex, friendship[[i]])
  
  #convert matrix to network object
  friendship[[i]] <- network(friendship[[i]])
  
  #set vertex attributes
  friendship[[i]] <- set.vertex.attribute(friendship[[i]], "sex", s)
  
  #calculate degree-based covariates
  idegsqrt <- sqrt(degree(friendship[[i]], cmode = "indegree"))
  friendship[[i]] <- set.vertex.attribute(friendship[[i]], "idegsqrt", idegsqrt)
  
  odegsqrt <- sqrt(degree(friendship[[i]], cmode = "outdegree"))
  friendship[[i]] <- set.vertex.attribute(friendship[[i]], "odegsqrt", odegsqrt)
}

#check network sizes
sapply(friendship, network.size)
t1 t2 t3 t4 
26 26 25 25

5.3 Visualize Networks

Code 
par(mfrow = c(2, 2), mar = c(0, 0, 2, 0))
for (i in 1:length(friendship)) {
  plot(friendship[[i]], 
       main = paste("Wave", i),
       usearrows = TRUE,
       edge.col = "grey50",
       vertex.col = c("lightpink", "lightblue")[get.vertex.attribute(friendship[[i]], "sex")])
}

5.4. Model 1: Pooled TERGM (No Temporal Dependencies)

As a baseline, we estimate a pooled TERGM that treats all time periods as independent:

Code 
model_2a <- btergm(
  friendship ~ edges + 
    mutual + 
    ttriple + 
    transitiveties + 
    ctriple + 
    nodeicov("idegsqrt") + 
    nodeicov("odegsqrt") + 
    nodeocov("odegsqrt") + 
    nodeofactor("sex") + 
    nodeifactor("sex") + 
    nodematch("sex") + 
    edgecov(primary),
  R = 100,
  parallel = "no",
  verbose = FALSE
)

summary(model_2a)
                     Estimate Boot mean     2.5%    97.5%
edges                -9.17591  -9.35366 -10.1717  -8.7292
mutual                2.96172   2.92790   2.3336   3.4775
ttriple               0.21316   0.20708   0.1102   0.3012
transitiveties        0.36794   0.34223   0.2347   0.4245
ctriple              -0.67677  -0.66110  -0.8700  -0.4783
nodeicov.idegsqrt     1.17509   1.23593   1.0359   1.5003
nodeicov.odegsqrt    -0.28432  -0.29985  -0.5536  -0.1675
nodeocov.odegsqrt     1.19846   1.24079   1.1665   1.3927
nodeofactor.sex.2     0.60720   0.61228   0.4327   0.7825
nodeifactor.sex.2     0.22578   0.23272   0.1075   0.3651
nodematch.sex         1.76805   1.78158   1.5765   2.0602
edgecov.primary[[i]]  1.05123   0.99961   0.7622   1.3889

5.5 Model 2: TERGM with Temporal Dependencies

Now we add temporal dependencies:

Code 
model_2b <- btergm(
  friendship ~ edges + 
    mutual + 
    ttriple + 
    transitiveties + 
    ctriple + 
    nodeicov("idegsqrt") + 
    nodeicov("odegsqrt") + 
    nodeocov("odegsqrt") + 
    nodeofactor("sex") + 
    nodeifactor("sex") + 
    nodematch("sex") + 
    edgecov(primary) +
    delrecip +                        # Delayed reciprocity
    memory(type = "stability"),       # Dyadic stability
  R = 100,
  parallel = "no",
  verbose = FALSE
)

summary(model_2b)
                      Estimate Boot mean     2.5%    97.5%
edges                 -9.54228  -9.90972 -12.4210  -8.7203
mutual                 2.17294   2.18165   1.8416   2.8489
ttriple                0.13423   0.12416   0.0345   0.2411
transitiveties         0.31509   0.33219   0.2884   0.3945
ctriple               -0.55434  -0.55656  -0.8214  -0.4185
nodeicov.idegsqrt      1.28064   1.34390   1.1152   1.6191
nodeicov.odegsqrt     -0.13179  -0.14224  -0.3148  -0.0238
nodeocov.odegsqrt      1.49624   1.57413   1.3781   2.0224
nodeofactor.sex.2      0.52990   0.57497   0.3846   0.8317
nodeifactor.sex.2      0.28694   0.31416  -0.1334   0.6005
nodematch.sex          1.49245   1.54311   1.3044   2.0275
edgecov.primary[[i]]   0.42395   0.43225  -0.2828   0.9817
edgecov.delrecip[[i]]  0.67436   0.74217   0.3309   1.5025
edgecov.memory[[i]]    0.77922   0.80008   0.6761   1.0279

5.6 New Temporal Terms

delrecip (Delayed Reciprocity):
Tests whether a tie from j→i at t-1 leads to i→j at t. Positive coefficient indicates reciprocation develops over time.

memory(type = "stability"):
Captures dyadic stability—both persistent ties AND persistent non-ties. Positive coefficient means relationships (present or absent) tend to be stable.

Important: When using temporal terms with lag 1, only waves 2-4 are used as outcomes (wave 1 provides the initial state).

5.7 Model Comparison

Code 
screenreg(list(model_2a, model_2b),
          custom.model.names = c("Model 2a", "Model 2b"),
          single.row = TRUE)

===========================================================================
                       Model 2a                   Model 2b                 
---------------------------------------------------------------------------
edges                    -9.18 [-10.17; -8.73] *    -9.54 [-12.42; -8.72] *
mutual                    2.96 [  2.33;  3.48] *     2.17 [  1.84;  2.85] *
ttriple                   0.21 [  0.11;  0.30] *     0.13 [  0.03;  0.24] *
transitiveties            0.37 [  0.23;  0.42] *     0.32 [  0.29;  0.39] *
ctriple                  -0.68 [ -0.87; -0.48] *    -0.55 [ -0.82; -0.42] *
nodeicov.idegsqrt         1.18 [  1.04;  1.50] *     1.28 [  1.12;  1.62] *
nodeicov.odegsqrt        -0.28 [ -0.55; -0.17] *    -0.13 [ -0.31; -0.02] *
nodeocov.odegsqrt         1.20 [  1.17;  1.39] *     1.50 [  1.38;  2.02] *
nodeofactor.sex.2         0.61 [  0.43;  0.78] *     0.53 [  0.38;  0.83] *
nodeifactor.sex.2         0.23 [  0.11;  0.37] *     0.29 [ -0.13;  0.60]  
nodematch.sex             1.77 [  1.58;  2.06] *     1.49 [  1.30;  2.03] *
edgecov.primary[[i]]      1.05 [  0.76;  1.39] *     0.42 [ -0.28;  0.98]  
edgecov.delrecip[[i]]                                0.67 [  0.33;  1.50] *
edgecov.memory[[i]]                                  0.78 [  0.68;  1.03] *
---------------------------------------------------------------------------
Num. obs.              2500                       1850                     
===========================================================================
* Null hypothesis value outside the confidence interval.

5.8 Goodness-of-Fit Assessment

Code 
gof_temporal <- gof(model_2b, 
                    nsim = 50,
                    statistics = c(esp, dsp, geodesic, deg, triad.undirected))
                 t=2 t=3 t=4
friendship (row)  26  25  25
friendship (col)  26  25  25
primary (row)     26  26  26
primary (col)     26  26  26
delrecip (row)    26  26  25
delrecip (col)    26  26  25
memory (row)      26  26  25
memory (col)      26  26  25
  label time   object where
1    21    3 networks   row
2    21    3 networks   col
3    21    4 networks   row
4    21    4 networks   col
                            t=2 t=3 t=4
maximum deleted nodes (row)   0   1   1
maximum deleted nodes (col)   0   1   1
remaining rows               26  25  25
remaining columns            26  25  25
                 t=2 t=3 t=4
friendship (row)  26  25  25
friendship (col)  26  25  25
primary (row)     26  25  25
primary (col)     26  25  25
delrecip (row)    26  25  25
delrecip (col)    26  25  25
memory (row)      26  25  25
memory (col)      26  25  25
Code 
par(mfrow = c(2, 3), mar = c(4, 4, 2, 1))
plot(gof_temporal)

5.9 Network Evolution Over Time

Code 
#calculate network statistics for each wave
net_stats <- data.frame(
  wave = 1:4,
  nodes = sapply(friendship, network.size),
  edges = sapply(friendship, network.edgecount),
  density = sapply(friendship, network.density),
  reciprocity = sapply(friendship, grecip, measure = "edgewise"),
  transitivity = sapply(friendship, gtrans)
)

print(net_stats)
   wave nodes edges   density reciprocity transitivity
t1    1    26    91 0.1400000   0.6153846    0.5102740
t2    2    26   117 0.1800000   0.5641026    0.5066163
t3    3    25   133 0.2216667   0.5413534    0.5153129
t4    4    25   119 0.1983333   0.5546218    0.4016667
Code 
#visualize trends
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))

plot(net_stats$wave, net_stats$density, 
     type = "b", pch = 16, col = "steelblue",
     xlab = "Wave", ylab = "Density",
     main = "Network Density")

plot(net_stats$wave, net_stats$reciprocity, 
     type = "b", pch = 16, col = "darkgreen",
     xlab = "Wave", ylab = "Reciprocity",
     main = "Edgewise Reciprocity")

plot(net_stats$wave, net_stats$transitivity, 
     type = "b", pch = 16, col = "darkred",
     xlab = "Wave", ylab = "Transitivity",
     main = "Transitivity")

This tutorial is based on UCLA STATS 218: Social Network Analysis taught by Prof Mark Handcock.