Welcome to Part 7 of our Social Network Analysis R Tutorial series! In our previous tutorial (Part 6), we explored the theoretical foundations of ERGMs, understanding the mathematical framework behind these powerful models. We covered the exponential family form, sufficient statistics, and the challenges of model degeneracy.
In this post, we will advance into sophisticated ERGM modeling by incorporating complex structural terms that capture network dependencies beyond simple dyadic relationships. We’ll explore geometrically weighted statistics, perform comprehensive MCMC diagnostics, and learn to interpret results in meaningful sociological contexts.
After completing this tutorial, you will be able to:
rm(list = ls())
library(network) # Network object manipulation
library(ergm) # ERGM fitting
library(networkdata) # Network datasets
library(sna) # Social network analysis tools
library(parallel) # Speed up MCMC
library(ergm.tapered) # Tapered ERGM for handling degeneracy
library(amen)
# Install ergm.tapered if needed
# devtools::install_github("statnet/ergm", ref="tapered")Why ergm.tapered?
The ergm.tapered package provides an alternative estimation method that better handles models with strong dependence between terms. When you include multiple structural terms (like gwesp and gwdsp), standard ERGM estimation can struggle with convergence. Tapered ERGM adds a “tapering” mechanism that helps stabilize estimation.
Software available on GitHub as ergm.tapered: https://github.com/statnet/ergm.tapered
You can use these following code to install ergm or ergm.tapered:
install.packages("ergm", repos="http://www.stat.ucla.edu/~handcock")
install.packages("ergm.tapered", repos="http://www.stat.ucla.edu/~handcock")or
install.packages("devtools")
devtools::install_github("statnet/ergm.tapered")Before diving into applications, let’s understand the sophisticated structural terms available in ERGMs. These terms capture complex network dependencies while avoiding the degeneracy issues that plague simpler specifications.
Why Geometrically Weighted Terms? Traditional count statistics (like triangle counts) can lead to model degeneracy because they treat all configurations equally. Geometrically weighted terms apply diminishing returns - the first few shared partners matter more than additional ones, making the model more stable and interpretable.
Here’s a detailed reference table for advanced ERGM structural terms:
| Function | Syntax | Description | Interpretation | When to Use |
|---|---|---|---|---|
| istar(k) | istar(2) |
k-in-stars: counts nodes with exactly k incoming edges | Popularity effects in directed networks | Testing for preferential attachment to popular nodes |
| ostar(k) | ostar(2) |
k-out-stars: counts nodes with exactly k outgoing edges | Activity effects in directed networks | Testing for variation in node activity levels |
| gwodegree(decay) | gwodegree(0.5, fixed=TRUE) |
Geometrically weighted out-degree distribution | Controls out-degree heterogeneity with diminishing returns | Modeling sender activity patterns without degeneracy |
| gwidegree(decay) | gwidegree(0.5, fixed=TRUE) |
Geometrically weighted in-degree distribution | Controls in-degree heterogeneity with diminishing returns | Modeling popularity/receiver effects without degeneracy |
| gwesp(decay) | gwesp(0.5, fixed=TRUE) |
Geometrically weighted edgewise shared partners | Transitivity: “friend of friend is friend” with diminishing returns | Capturing clustering while avoiding degeneracy |
| gwdsp(decay) | gwdsp(0.5, fixed=TRUE) |
Geometrically weighted dyadwise shared partners | Shared partners between any dyad (connected or not) | Modeling local clustering and “small world” effects |
| ctriple | ctriple or ctriad |
Cyclic triples: A→B→C→A patterns | Non-hierarchical, cyclic relationships | Testing for cyclic vs hierarchical structures |
| ttriple | ttriple or ttriad |
Transitive triples: A→B→C with A→C | Hierarchical, transitive closure | Testing for transitivity and balance theory |
Check the model terms using help("ergm-terms")
Choosing Decay Parameters The decay parameter (α) in geometrically weighted terms controls how quickly the weight diminishes:
Start with α = 0.5 as a default and adjust based on model diagnostics.
One of the most frustrating challenges in ERGM modeling is model degeneracy (or near-degeneracy). This occurs when a model places almost all its probability mass on extreme network configurations—like empty graphs, complete graphs, or other boundary cases.
A classic example is the triangle model:
\[P_\eta(Y = y) = \frac{\exp\{\eta_1\text{edge}(y) + \eta_2\text{triangle}(y)\}}{c(\eta_1, \eta_2)}\]
This model is near-degenerate for most values of \(\eta_2 > 0\). When we try to include theoretically meaningful terms like triangles, GWESP, or other clustering measures, the model often becomes unstable and produces unreliable results.
Near-degeneracy has constrained practical ERGM modeling:
This is why we’ve had to be so careful about model specification in previous tutorials—we’ve been working around the degeneracy problem rather than solving it directly.
The Tapered ERGM has the form:
\[q(y|\theta, \tau) = \frac{1}{Z(\theta, \tau)} \exp\left\{\sum_k \theta_k g_k(y) - \sum_k \tau_k(\mu_k(\theta, \tau) - g_k(y))^2\right\}\]
where:
The tapering term acts as a soft constraint that penalizes network configurations far from the mean \(\mu\). Think of it as:
The tapering term has a beautiful theoretical interpretation. The Tapered ERGM arises from a constrained maximum entropy problem:
\[\begin{align} \text{maximize}_q &\quad -\sum_y q(y)\log(q(y)) \\ \text{subject to} &\quad E_q[g_i(Y,X)] = \mu_i \\ &\quad E_q[(\mu_i - g_i(Y,X))^2] \leq \sigma^2_i, \quad \forall i \end{align}\]
In other words, we’re finding the distribution that: 1. Has maximum entropy (is least constrained) 2. Matches expected network statistics 3. Also constrains the variance of those statistics to be at most \(\sigma^2_i\)
This variance constraint is what prevents degeneracy.
When Tapering is Zero (\(\tau = 0\))
The model reduces to standard ERGM, and \(\theta_k\) represents conditional log-odds as usual.
When Tapering is Non-Zero (\(\tau > 0\))
The interpretation changes slightly. The conditional log-odds of a tie becomes:
\[\log\left(\frac{P(Y_{ij}=1|Y^c_{ij})}{P(Y_{ij}=0|Y^c_{ij})}\right) = \sum_k \Delta t_k(Y_{ij})[\theta_k + \tau_k \delta_{kij}]\]
The term \(\tau_k \delta_{kij}\) represents a bias introduced by tapering. However, this bias:
How do we determine \(\tau\)? The lecture materials suggest using kurtosis as a measure of potential degeneracy:
\[\text{Kurt}[g(Y)] = \frac{\mu_4}{\mu_2^2} = E[Z^4] \text{ where } Z \text{ is standardized}\]
Key facts about kurtosis: - \(\text{Kurt}[g(Y)] \geq 1\) always - Gaussian: \(\text{Kurt} = 3\) - Uniform: \(\text{Kurt} = 1.8\) - Higher kurtosis suggests bimodality/degeneracy
The tapering is determined by constraining variance:
\[V(\theta, \tau) = Var_q(g(Y)) = \sigma^2\]
The ergm.tapered package automatically estimates appropriate \(\tau\) values based on: 1. The observed network statistics 2. Desired variance constraints 3. Ensuring non-degeneracy
Automatic Selection: in practice, the package ergm.tapered() will
Use Tapered ERGM when (recommended):
Advantages:
ergm.tapered package)Trade-offs:
The Tapered ERGM represents a major breakthrough in network modeling, finally allowing us to use the full richness of ERGM specifications without worrying too much about the degeneracy.
We’ll work with a network collected by Charles Kadushin (1995) of 28 members of the French financial elite. This network captures friendship ties and includes rich vertex attributes that allow us to explore multiple dimensions of homophily.
Data Source:
Kadushin, C. 1995. Friendship among the French financial elite. American Sociological Review 60:202-221.
data(ffef)
help(ffef)
network.size(ffef)[1] 28network.edgecount(ffef)[1] 66gden(ffef)[1] 0.1746032list.vertex.attributes(ffef) [1] "birthdate" "birthplace" "cabinet" "clubs" "eliteprom"
[6] "elitevote" "ena" "enayear" "fathers-lev" "finance-min"
[11] "igyear" "inspec-gen" "masons" "na" "normal-sch"
[16] "party" "polytech" "polyyear" "prestige" "religion"
[21] "sciencepoly" "socialreg" "topboards" "university" "vertex.names"
[26] "zipcode" The network includes several important vertex covariates:
ena_status <- ffef %v% "ena"
vertex_colors <- ifelse(ena_status == 1, "coral", "skyblue")
plot(ffef, vertex.col = vertex_colors,
main = "French Financial Elite Network Colored by ENA Status",
vertex.cex = 1.5, edge.col = "gray70")
Observations:
#using ergm.tapered() for better handling of strong dependencies
fit_elite <- ergm.tapered(
ffef ~
edges + # Baseline density
nodematch("ena") + # ENA homophily
nodematch("prestige") + # Prestige homophily
nodematch("party") + # Political party homophily
gwesp(0.5, fixed = TRUE) + # Transitive closure (clustering)
gwdsp(0.5, fixed = TRUE), # Shared partners
verbose = FALSE
)
summary(fit_elite) Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -4.74578 0.52795 0 -8.989 < 1e-04 ***
nodematch.ena 1.55048 0.37570 0 4.127 < 1e-04 ***
nodematch.prestige 0.66878 0.26532 0 2.521 0.01171 *
nodematch.party 1.22964 0.37519 1 3.277 0.00105 **
gwesp.fixed.0.5 0.46892 0.20191 0 2.322 0.02021 *
gwdsp.fixed.0.5 0.19187 0.06753 0 2.841 0.00449 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The estimated tapering scaling factor is 2.
Null Deviance: 524.0 on 378 degrees of freedom
Residual Deviance: 301.4 on 372 degrees of freedom
AIC: 313.4 BIC: 337 (Smaller is better. MC Std. Err. = 0)coef_table <- data.frame(
Term = names(coef(fit_elite)),
Coefficient = round(coef(fit_elite), 3),
SE = round(summary(fit_elite)$coefficients[, "Std. Error"], 3),
OR = round(exp(coef(fit_elite)), 3),
p_value = round(summary(fit_elite)$coefficients[, "Pr(>|z|)"], 4)
)
print(coef_table) Term Coefficient SE OR p_value
edges edges -4.746 0.528 0.009 0.0000
nodematch.ena nodematch.ena 1.550 0.376 4.714 0.0000
nodematch.prestige nodematch.prestige 0.669 0.265 1.952 0.0117
nodematch.party nodematch.party 1.230 0.375 3.420 0.0010
gwesp.fixed.0.5 gwesp.fixed.0.5 0.469 0.202 1.598 0.0202
gwdsp.fixed.0.5 gwdsp.fixed.0.5 0.192 0.068 1.212 0.0045baseline_prob <- plogis(coef(fit_elite)["edges"])
ena_match_prob <- plogis(coef(fit_elite)["edges"] + coef(fit_elite)["nodematch.ena"])
prestige_match_prob <- plogis(coef(fit_elite)["edges"] + coef(fit_elite)["nodematch.prestige"])
party_match_prob <- plogis(coef(fit_elite)["edges"] + coef(fit_elite)["nodematch.party"])
cat("Baseline (no matches):", round(baseline_prob, 4), "\n")Baseline (no matches): 0.0086 cat("Same ENA status:", round(ena_match_prob, 4), "\n")Same ENA status: 0.0393 cat("Same prestige level:", round(prestige_match_prob, 4), "\n")Same prestige level: 0.0167 cat("Same party affiliation:", round(party_match_prob, 4), "\n")Same party affiliation: 0.0289 Edges: \(\beta = -4.63\) represent the baseline log-odds of a friendship tie existing between any two randomly selected individuals in the networks, while controlling other predictors. The strong negative coefficient indicates this is a very sparse network. The \(OR = 0.010\) means the baseline odds of a tie is about 1.0%. The baseline robability of friendship is approximately 0.96%.
nodematch.ena: \(\beta = 1.57\) measures homophily on ENA attendance. The \(OR = 4.79\) means that two individuals who share ENA status (both graduate from ENA or both did not) have 4.79 times higher odds of being friends compared to two individuals with different ENA status. The coefficient is strong and highly significant. And the probability of a friendship tie between two individuals who share the same ENA status is 4.45%.
nodematch.prestige: \(\beta = 0.68\) measures homophily on prestige level (having particule and/or social register listing).The coefficient is moderate and significant (p = 0.015). The \(OR = 1.97\) mean that two individuals with the same prestige level have 1.97 times higher odds of being friends compared to those with different prestige levels. The probability of a friendship tie between two individuals who share the same prestige level is 1.88%.
nodematch.party: \(\beta = 1.25\) measures homophily on political party affiliation. The coefficient is strong and highly significant (p < 0.001).The \(OR = 3.49\) means that two individuals from the same political party have 3.49 times higher odds of being friends compared to those from different parties. The probability of a friendship tie between two individuals who share the same party affiliation is 3.28%.
GWESP(0.5): \(\beta = 0.43\) measures transitivity and clustering through Geometrically Weighted Edgewise Shared Partners. This term captures the “friend of a friend is a friend” phenomenon. The coefficient is moderate and significant (p = 0.013). The \(OR = 1.53\) indicates that networks with more shared partners (triangles) are more probable. The French financial elite network exhibits strong transitivity. If person A is friends with person B, and person B is friends with person C, then A and C are significantly more likely to also be friends.
GWDSP(0.5): \(\beta = 0.19\) measures local network structure through Geometrically Weighted Dyadwise Shared Partners. This term captures the tendency for dyads (pairs) to share partners, even if they’re not directly connected. The coefficient is modest but significant (p = 0.003). The \(OR = 1.21\) indicates preference for configurations where individuals share common friends. Beyond direct transitivity, there’s a tendency for individuals to be embedded in overlapping social circles. This creates a “small world” structure where even if two people aren’t direct friends, they’re likely to share mutual acquaintances.
plogis(coef["edges"] + coef["gwesp.fixed.0.5"])This does NOT work because:
Summary:
| Term Type | Example Terms | Can use plogis(edges + term)? |
Why/Why not? |
|---|---|---|---|
| Node attributes | nodematch(), nodecov(), absdiff() |
YES | Binary 0/1 or known numeric covariate values |
| Dyad attributes | edgecov(), dyadcov() |
YES (with care) | Known covariate value for specific dyad |
| Structural terms | gwesp(), gwdsp(), triangle, ttriad, ctriad |
NO | Value depends on entire network configuration |
| Degree terms | gwdegree(), degree(), idegree(), odegree() |
NO | Depends on node’s position in network |
mcmc_diag <- mcmc.diagnostics(fit_elite, vars.per.page = 6)Sample statistics summary:
Iterations = 7168:131072
Thinning interval = 512
Number of chains = 1
Sample size per chain = 243
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
edges 0.2551 3.396 0.2178 0.2178
nodematch.ena -0.1481 3.732 0.2394 0.2394
nodematch.prestige 0.2675 3.286 0.2108 0.2368
nodematch.party -0.2551 2.165 0.1389 0.1539
gwesp.fixed.0.5 -0.2306 6.289 0.4034 0.3635
gwdsp.fixed.0.5 0.7706 17.451 1.1195 1.1195
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
edges -6.00 -2.000 0.0000 2.000 6.00
nodematch.ena -7.00 -2.500 0.0000 2.000 7.00
nodematch.prestige -5.00 -2.000 0.0000 2.000 7.00
nodematch.party -4.00 -2.000 0.0000 1.000 4.00
gwesp.fixed.0.5 -12.36 -4.163 -0.2767 3.834 11.94
gwdsp.fixed.0.5 -31.39 -11.157 0.5425 12.347 32.89
Are sample statistics significantly different from observed?
edges nodematch.ena nodematch.prestige nodematch.party
diff. 0.2551440 -0.1481481 0.2674897 -0.25514403
test stat. 1.1713321 -0.6187983 1.1295006 -1.65758049
P-val. 0.2414653 0.5360492 0.2586867 0.09740217
gwesp.fixed.0.5 gwdsp.fixed.0.5 Overall (Chi^2)
diff. -0.2306391 0.7705798 NA
test stat. -0.6345827 0.6883491 8.7038475
P-val. 0.5257006 0.4912330 0.2113753
Sample statistics cross-correlations:
edges nodematch.ena nodematch.prestige nodematch.party
edges 1.0000000 0.65124476 0.3056581 0.100520838
nodematch.ena 0.6512448 1.00000000 0.2225766 0.069973569
nodematch.prestige 0.3056581 0.22257659 1.0000000 0.107209446
nodematch.party 0.1005208 0.06997357 0.1072094 1.000000000
gwesp.fixed.0.5 0.6453890 0.54044787 0.1958230 -0.024200392
gwdsp.fixed.0.5 0.6351486 0.21440212 0.1562678 0.008799192
gwesp.fixed.0.5 gwdsp.fixed.0.5
edges 0.64538902 0.635148626
nodematch.ena 0.54044787 0.214402119
nodematch.prestige 0.19582297 0.156267804
nodematch.party -0.02420039 0.008799192
gwesp.fixed.0.5 1.00000000 0.569487104
gwdsp.fixed.0.5 0.56948710 1.000000000
Sample statistics auto-correlation:
Chain 1
edges nodematch.ena nodematch.prestige nodematch.party
Lag 0 1.0000000000 1.0000000000 1.00000000 1.00000000
Lag 512 0.0003085203 0.0619330983 0.11374068 0.26791210
Lag 1024 0.0045786107 0.0146450292 -0.04058172 0.03309929
Lag 1536 -0.0373848575 0.0465258243 -0.04305592 -0.12742631
Lag 2048 0.0058741033 0.0217632475 -0.06955911 -0.05032647
Lag 2560 0.0544176158 -0.0006588757 0.07730979 0.03956701
gwesp.fixed.0.5 gwdsp.fixed.0.5
Lag 0 1.00000000 1.00000000
Lag 512 -0.10606986 0.01751455
Lag 1024 -0.02840634 -0.07995494
Lag 1536 -0.04271817 -0.05239432
Lag 2048 0.06281154 0.07101157
Lag 2560 -0.01374112 0.04090792
Sample statistics burn-in diagnostic (Geweke):
Chain 1
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
edges nodematch.ena nodematch.prestige nodematch.party
-1.2426 -0.3939 1.4785 3.3660
gwesp.fixed.0.5 gwdsp.fixed.0.5
-0.7608 0.3016
Individual P-values (lower = worse):
edges nodematch.ena nodematch.prestige nodematch.party
0.2140208588 0.6936776354 0.1392719443 0.0007625596
gwesp.fixed.0.5 gwdsp.fixed.0.5
0.4467701978 0.7629851534
Joint P-value (lower = worse): 0.0518683 .
MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).MCMC Diagnostic Components Explained:
Let’s extract some key diagnostic statistics:
gof_elite <- gof(fit_elite, GOF = ~degree + esp + dsp + triadcensus,
verbose = FALSE)
plot(gof_elite)
From the sample statistics vs. observed values, the statistics show good agreement between simulated and observed values. All terms have p-value > 0.05, The overall chi-squared test p = 0.022, indicating the model successfully reproduces the observed network statistics. The trace plots show relatively good mixing, fluctuating randomly around zero, which means MCMC chains are exploring the parameter space well and not getting trapped in any particular region. All the density plots are approximately centered near zeor, and no extreme skewness and multinodality (multiple peaks), indicating the sampling stable parameter estimation. Moreover, for the Geweke Diagnostic (Burn-in Test), all p-values are well above 0.05, indicating no evidence of non-convergence. The joint p-value of 0.820 indicates good overall convergence. Overall, the model successfully converged and the parameter estimates are reliable for inference.
Now let’s examine balance theory using directed friendship nominations in a sixth-grade classroom, we’ve used this dataset in the previous tutorial (Part 5):
Heider’s Balance Theory In a balanced network:
Balanced networks have positive transitive triad coefficients and negative cyclic triad coefficients.
library(networkdata)
data(hansell)
help(hansell)
#summary(hansell)
table(hansell %v% "sex")
female male
14 13 triad.census(hansell) 003 012 102 021D 021U 021C 111D 111U 030T 030C 201 120D 120U 120C 210 300
[1,] 771 982 248 234 115 150 68 97 126 2 13 32 32 20 28 7#transitive triads
triads <- triad.census(hansell)
triads[9] # 030T[1] 126#cyclic triads
triads[10] # 030C [1] 2fit_balance <- ergm.tapered(
hansell ~
edges + # Baseline density
mutual + # Reciprocity
nodematch("sex", diff = TRUE) + # Sex homophily (differentiated)
ttriple + # Transitive triples
ctriple, # Cyclic triples
verbose = FALSE
)
summary(fit_balance) Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -3.05448 0.12122 0 -25.198 <1e-04 ***
mutual 0.27994 0.36194 0 0.773 0.439
nodematch.sex.female 0.87380 0.13343 0 6.549 <1e-04 ***
nodematch.sex.male 1.05969 0.16671 0 6.356 <1e-04 ***
ttriple 0.32832 0.01771 0 18.544 <1e-04 ***
ctriple -0.40282 0.07707 0 -5.226 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The estimated tapering scaling factor is 2.
Null Deviance: 973.2 on 702 degrees of freedom
Residual Deviance: 633.1 on 696 degrees of freedom
AIC: 645.1 BIC: 672.5 (Smaller is better. MC Std. Err. = 0)balance_coef <- data.frame(
Term = c("Edges", "Mutual", "Same Sex (Female)", "Same Sex (Male)",
"Transitive Triples", "Cyclic Triples"),
Coefficient = round(coef(fit_balance), 3),
OR = round(exp(coef(fit_balance)), 3))
print(balance_coef, row.names = FALSE) Term Coefficient OR
Edges -3.054 0.047
Mutual 0.280 1.323
Same Sex (Female) 0.874 2.396
Same Sex (Male) 1.060 2.885
Transitive Triples 0.328 1.389
Cyclic Triples -0.403 0.668Evidence for Balance model
The model provides strong evidence for balance theory:
This pattern suggests the classroom network organizes into balanced, sex-segregated friendship groups with clear hierarchical structures rather than cyclic patterns.
Transitive ties (ttriple): \(\beta = 0.33\), the coefficient measures the tendency for transitive triadic closure in the network, it’s significant. The \(OR = 1.39\) means that for each additional transitive triple configuration, the odds of the network configuration increased by 39%. When student A nominates student B as a friend, and B nominates C, there is a significantly higher-than-random probability that A will also nominate C. This creates cohesive friendship groups with closed triadic structures.
Cyclic ties (ctriple): \(\beta = -0.40\), the coefficient measures the tendency for cyclic triadic structures in the network, where: \(A→B→C→A\) forms a cycle, but without the full transitive closure. The negative and highly significant coefficient (p < 0.001) indicates a strong avoidance of cyclic triadic structures. \(OR = 0.67\) means that cyclic triple configurations decrease the odds of the network configuration by 33%. The negative coefficient shows that cyclic patterns are less common than would be expected in a random network.
# Fit alternative model for French elite network
fit_elite_extended <- ergm.tapered(
ffef ~
edges +
nodematch("ena") +
nodematch("prestige") +
nodematch("party") +
nodematch("masons") + # Add Masonic membership
nodematch("religion") + # Add religion
nodecov("topboards") + # Board memberships as covariate
absdiff("topboards") + # Difference in board memberships
gwesp(0.5, fixed = TRUE) +
gwdsp(0.5, fixed = TRUE),
verbose = FALSE
)
# Model comparison
anova_result <- anova(fit_elite, fit_elite_extended)
print(anova_result)Analysis of Variance Table
Model 1: ffef ~ Taper(~edges + nodematch("ena") + nodematch("prestige") +
nodematch("party") + gwesp(0.5, fixed = TRUE) + gwdsp(0.5,
fixed = TRUE), coef = .taper.coef, m = .taper.center)
Model 2: ffef ~ Taper(~edges + nodematch("ena") + nodematch("prestige") +
nodematch("party") + nodematch("masons") + nodematch("religion") +
nodecov("topboards") + absdiff("topboards") + gwesp(0.5,
fixed = TRUE) + gwdsp(0.5, fixed = TRUE), coef = .taper.coef,
m = .taper.center)
Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
NULL 378 524.02
1 6 222.63 372 301.39 <2e-16 ***
2 4 0.56 368 300.83 0.9674
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC comparison
aic_comparison <- data.frame(
Model = c("Base", "Extended"),
AIC = c(AIC(fit_elite), AIC(fit_elite_extended)),
BIC = c(BIC(fit_elite), BIC(fit_elite_extended)),
Parameters = c(length(coef(fit_elite)), length(coef(fit_elite_extended)))
)
print(aic_comparison) Model AIC BIC Parameters
1 Base 313.3867 336.9960 6
2 Extended 320.8267 360.1756 10Model Selection Guidelines
In our case, the simpler model is preferred: lower AIC, similar fit, better convergence.
Butland et al (2005) “Interaction network containing conserved and essential protein complexes in Escherichia coli” reported a network of protein-protein interactions (bindings) that we obtained from http://pil.phys.uniroma1.it/~gcalda/cosinsite/extra/data/proteins/.
This data is available in the networkdata package.
Convert the edgelist to a directed network (The el2sm function may be helpful).
Fit various tapered ERGM models to the network using ergm.tapered. Consider terms documented under ergm-terms. Good candidates include istar, ostar, gwodegree, gwidegree, dgwest, dgwdsp, ctriple, ttriple. Check the MCMC diagnostics with mcmc.diagnostics. Overall, how does the goodness-of-fit look?
Model 1 is a basic model, only include edges + mutual terms, The mutual term has an estimate of -Inf, may indicate the network don’t have many mutual ties.
Model 2 is a degree distribution model, includes edges + gwidegree + gwodegree, MCMC diagnostics show relatively good convergence, sample statistics not significantly different from observed (overall p = 0.548) and the trace plot show good mixing.
Model 3 added dgwesp + dgwdsp + ctriple + ttriple based on the model 2, the ctriple term has estimate -14.20 with NA standard error, which may suggest there are zero or nearly zero cyclic triads in the network. MCMC diagnostics show reasonable convergence, with overall chi-squared test p = 0.530. Lower AIC (6770 vs 7083) suggests better fit than Model 2.
library(networkdata)
library(amen)
library(network)
data(butland_ppi)
help(butland_ppi)
head(butland_ppi) V1 V2
[1,] 3 2
[2,] 4 1
[3,] 5 1
[4,] 5 2
[5,] 5 3
[6,] 5 4#get all unique vertices and create mapping
all_vertices <- sort(unique(c(butland_ppi[, 1], butland_ppi[, 2])))
n_vertices <- length(all_vertices)
vertex_map <- data.frame(original_id = all_vertices, new_id = 1:n_vertices)
#remap edgelist to sequential IDs
edgelist_remapped <- butland_ppi
edgelist_remapped[, 1] <- match(butland_ppi[, 1], vertex_map$original_id)
edgelist_remapped[, 2] <- match(butland_ppi[, 2], vertex_map$original_id)
#convert edgelist to sociomatrix
sm <- el2sm(edgelist_remapped)
#create network object
ppi_net <- network(sm, directed = TRUE)
network.size(ppi_net)[1] 270network.edgecount(ppi_net)[1] 716network.edgecount(ppi_net)[1] 716library(parallel)
library(ergm.tapered)
set.seed(123)
#model 1: basic model with edges and reciprocity
#fit_ppi1 <- ergm.tapered(ppi_net ~ edges + mutual, r= 1)
fit_ppi1 <- ergm(ppi_net ~ edges + mutual)
summary(fit_ppi1)Call:
ergm(formula = ppi_net ~ edges + mutual)
Monte Carlo Maximum Likelihood Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -4.59802 0.04037 1 -113.9 <1e-04 ***
mutual -Inf 0.00000 0 -Inf <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning: The following terms have infinite coefficient estimates:
mutual #mcmc.diagnostics(fit_ppi1)
#model 2: add in-degree / outdegree distribution effects
fit_ppi2 <- ergm.tapered(ppi_net ~ edges + gwidegree(0.5, fixed = TRUE) + gwodegree(0.5, fixed = TRUE))
summary(fit_ppi2) Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -3.23637 0.04762 0 -67.959 <1e-04 ***
gwideg.fixed.0.5 -4.94911 0.16593 0 -29.827 <1e-04 ***
gwodeg.fixed.0.5 -1.48985 0.16165 0 -9.216 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The estimated tapering scaling factor is 2.
Null Deviance: 100687 on 72630 degrees of freedom
Residual Deviance: 7077 on 72627 degrees of freedom
AIC: 7083 BIC: 7110 (Smaller is better. MC Std. Err. = 0)#mcmc.diagnostics(fit_ppi2)
#model 3: add GW edgewise shared partners + GW Dyadwise Shared Partners
fit_ppi3 <- ergm.tapered(ppi_net ~ edges + gwidegree(0.5, fixed = TRUE) + gwodegree(0.5, fixed = TRUE) +
dgwesp(0.5, fixed = TRUE) + dgwdsp(0.5, fixed = TRUE) + ctriple + ttriple)
summary(fit_ppi3) Results:
Estimate Std. Error MCMC % z value Pr(>|z|)
edges -3.68349 0.07278 0 -50.609 <1e-04 ***
gwideg.fixed.0.5 -3.62306 0.15319 0 -23.650 <1e-04 ***
gwodeg.fixed.0.5 -0.07758 0.17652 0 -0.440 0.66
gwesp.OTP.fixed.0.5 1.25150 0.10950 0 11.429 <1e-04 ***
gwdsp.OTP.fixed.0.5 -0.13459 0.01107 0 -12.156 <1e-04 ***
ctriple -14.20162 NA NA NA NA
ttriple 0.25770 0.03041 1 8.473 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The estimated tapering scaling factor is 2.
Null Deviance: 100687 on 72630 degrees of freedom
Residual Deviance: 6756 on 72623 degrees of freedom
AIC: 6770 BIC: 6835 (Smaller is better. MC Std. Err. = 0)mcmc.diagnostics(fit_ppi3)Sample statistics summary:
Iterations = 516096:10240000
Thinning interval = 8192
Number of chains = 1
Sample size per chain = 1188
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
edges 0.7273 20.073 0.5824 1.0869
gwideg.fixed.0.5 0.3315 7.002 0.2032 0.2438
gwodeg.fixed.0.5 0.4258 7.554 0.2192 0.4120
gwesp.OTP.fixed.0.5 -1.1183 12.920 0.3748 0.6775
gwdsp.OTP.fixed.0.5 -0.9406 49.632 1.4400 1.8440
ctriple 0.0000 0.000 0.0000 0.0000
ttriple -1.4705 32.450 0.9415 3.4166
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
edges -38.00 -13.000 0.0000 15.000 41.00
gwideg.fixed.0.5 -13.45 -4.267 0.2103 4.910 14.42
gwodeg.fixed.0.5 -14.44 -4.514 0.3375 5.519 15.11
gwesp.OTP.fixed.0.5 -27.19 -10.008 -1.2674 8.071 22.80
gwdsp.OTP.fixed.0.5 -97.18 -34.265 -1.1433 33.403 96.41
ctriple 0.00 0.000 0.0000 0.000 0.00
ttriple -63.00 -25.000 -3.0000 20.000 63.33
Are sample statistics significantly different from observed? edges gwideg.fixed.0.5 gwodeg.fixed.0.5 gwesp.OTP.fixed.0.5
diff. 0.7272727 0.3315362 0.4257717 -1.11833353
test stat. 0.6691541 1.3596334 1.0333857 -1.65064489
P-val. 0.5033972 0.1739459 0.3014234 0.09881111
gwdsp.OTP.fixed.0.5 ctriple ttriple Overall (Chi^2)
diff. -0.9405732 0 -1.4705387 NA
test stat. -0.5100722 NaN -0.4304121 5.1678864
P-val. 0.6100009 NaN 0.6668959 0.5297571
Sample statistics cross-correlations: edges gwideg.fixed.0.5 gwodeg.fixed.0.5
edges 1.00000000 0.4269642 0.65956349
gwideg.fixed.0.5 0.42696417 1.0000000 0.30435499
gwodeg.fixed.0.5 0.65956349 0.3043550 1.00000000
gwesp.OTP.fixed.0.5 0.22779159 -0.1895341 -0.06397202
gwdsp.OTP.fixed.0.5 0.36242725 0.1921583 0.20738634
ctriple NA NA NA
ttriple -0.01901124 -0.1382923 -0.14616591
gwesp.OTP.fixed.0.5 gwdsp.OTP.fixed.0.5 ctriple ttriple
edges 0.22779159 0.36242725 NA -0.01901124
gwideg.fixed.0.5 -0.18953414 0.19215830 NA -0.13829233
gwodeg.fixed.0.5 -0.06397202 0.20738634 NA -0.14616591
gwesp.OTP.fixed.0.5 1.00000000 0.23252966 NA 0.64993306
gwdsp.OTP.fixed.0.5 0.23252966 1.00000000 NA -0.04517776
ctriple NA NA 1 NA
ttriple 0.64993306 -0.04517776 NA 1.00000000
Sample statistics auto-correlation:
Chain 1
edges gwideg.fixed.0.5 gwodeg.fixed.0.5 gwesp.OTP.fixed.0.5
Lag 0 1.00000000 1.000000000 1.00000000 1.00000000
Lag 8192 0.21494244 0.180153729 0.23245618 0.42194314
Lag 16384 0.09749634 0.042029652 0.15218526 0.25848757
Lag 24576 0.03196039 0.001603261 0.11438570 0.17421111
Lag 32768 0.01160266 0.023297475 0.07232346 0.12904988
Lag 40960 0.04270366 0.040651136 0.07850349 0.09324514
gwdsp.OTP.fixed.0.5 ctriple ttriple
Lag 0 1.000000000 NaN 1.0000000
Lag 8192 0.102487501 NaN 0.7013075
Lag 16384 0.079645513 NaN 0.5554517
Lag 24576 0.029373624 NaN 0.4579118
Lag 32768 0.066712902 NaN 0.3884092
Lag 40960 0.009440932 NaN 0.3525965
Sample statistics burn-in diagnostic (Geweke):Chain 1
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
edges gwideg.fixed.0.5 gwodeg.fixed.0.5 gwesp.OTP.fixed.0.5
1.1512 0.6234 2.6941 -0.3126
gwdsp.OTP.fixed.0.5 ctriple ttriple
-2.0375 NaN 0.3523
Individual P-values (lower = worse):
edges gwideg.fixed.0.5 gwodeg.fixed.0.5 gwesp.OTP.fixed.0.5
0.249633746 0.533037272 0.007058556 0.754567964
gwdsp.OTP.fixed.0.5 ctriple ttriple
0.041598024 NaN 0.724576484
Joint P-value (lower = worse): 0.1153546 .
MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).Hansell, S. (1984). Cooperative groups, weak ties, and the integration of peer friendships. Social Psychology Quarterly, 47(4), 316-328.
Kadushin, C. (1995). Friendship among the French financial elite. American Sociological Review, 60, 202-221.
Snijders, T.A.B., Pattison, P.E., Robins, G.L., & Handcock, M.S. (2006). New specifications for exponential random graph models. Sociological Methodology, 36(1), 99-153.
Hunter, D.R., & Handcock, M.S. (2006). Inference in curved exponential family models for networks. Journal of Computational and Graphical Statistics, 15(3), 565-583.
Blackburn, B., & Handcock, M. S. (2023). Practical network modeling via tapered exponential-family random graph models. Journal of Computational and Graphical Statistics, 32(2), 388-401.
This tutorial is based on UCLA STATS 218: Social Network Analysis taught by Prof Mark Handcock.
4. Social Balance and Transitivity in Network Analysis
4.1 Definition of Social Balance
Social balance and transitivity are concepts in network analysis that help us understand how relationships form and evolve in social networks. These concepts, rooted in social psychology and graph theory, provide insights into the micro-level mechanisms that shape macro-level network structures.
Heider’s Balance Theory
Fritz Heider’s balance theory examines relationships among three actors (P-O-X triads) where each dyadic relationship can be positive (+) or negative (-). According to Heider:
The key principle: multiply the signs of all three edges. If the product is positive, the triad is balanced; if negative, it’s unbalanced.
The following plot shows the four foundamental triad configurations in balance theory. Balanced triads (top row) are psychologically stable, while unbalanced triads (bottom row) create cognitive tension that motivates relationship change.
4.2 From Undirected to Directed Networks
While Heider’s original theory focused on undirected relationships, real social networks often involve directed ties (e.g., “A nominates B as a friend” doesn’t guarantee reciprocity). For directed networks, we use:
Transitivity in Directed Networks
Transitivity captures the “friend of a friend” phenomenon:
4.3 Network-Level Implications
Balance theory has profound implications for understanding macro-level network structure:
Clustered Communities: If balance holds throughout a network, we expect clear group boundaries with positive ties within groups and negative ties between groups
Structural Balance Theorem: A completely balanced signed network must have one of two structures:
Weak Balance: Real networks often show “weak balance” where the all-negative triad (-)(-)(-) is also allowed, permitting multiple clusters rather than just two factions
4.4 Measuring Balance in Real Networks
We can empirically test balance theory using:
ttriad(transitive triads count)ctriad(cyclic triads count)