指数随机图模型完整代码：从网络结构到统计推断的落地路径

简介：这份资源围绕指数随机图模型（ERGM）展开，面向社会网络分析、复杂网络建模方向的学习者与研究者，帮助读者理解如何从关系视角刻画整体网络的形成机制。内容涵盖简单随机图模型、二元独立性模型、二元依赖性模型与高序依赖性模型四类建模思路，可用于解释并预测网络关系的生成，适合具备一定统计与网络分析基础的中高级读者参考。资源以zip压缩包形式提供，整体约11.02MB，文件总数与具体类型明细上游暂未提供，下载后可按模型类别自行整理学习。目前已有1427人学习下载，说明该主题在相关领域具有一定关注度。通过这份完整代码，读者可对照四类模型理解其结构差异与适用场景，把握网络成员属性、全局结构特征等因素如何融入建模过程，为自身的社会网络实证研究或复杂系统分析提供可复用的实现参考与思路借鉴。

1. 指数随机图模型完整代码：从网络结构到统计推断的落地路径

社交网络、蛋白质互作网络、贸易网络，这些看似毫不相关的系统，背后都有一个共同的分析需求：判断观测到的网络结构特征——比如互惠性、传递性、度分布的异质性——究竟是随机涌现的，还是存在统计上显著的生成机制。指数随机图模型（Exponential Random Graph Model, ERGM）就是回答这类问题的核心工具。它把网络看作一个随机变量，用局部结构统计量的线性组合来定义概率分布，再通过最大似然估计推断各统计量的权重。你手里如果有一份网络数据，想量化“互惠关系有多强”“三角形闭合是否超出随机预期”，ERGM 就是绕不开的方法。这篇笔记围绕“指数随机图模型完整代码”这个目标，把模型设定、参数估计、拟合诊断到结果解读的完整链路拆开，给出可直接复现的 Python 和 R 代码，同时把实操中容易翻车的地方标出来。适合已经掌握基础网络分析、需要把 ERGM 跑通并拿到可解释结果的从业者。

2. ERGM 的统计逻辑与代码实现选型

2.1 从网络统计量到概率分布：ERGM 到底在算什么

ERGM 的核心思想可以用一句话概括：一个网络出现的概率，取决于它包含哪些局部结构。形式化地，对于网络 $y$，其概率为：

$$P(Y=y|\theta) = \frac{\exp(\theta^T g(y))}{c(\theta)}$$

其中 $g(y)$ 是网络统计量向量，比如边数、互惠边数、三角形数、星形结构数等；$\theta$ 是对应的参数向量；$c(\theta)$ 是归一化常数，对所有可能网络求和。这个归一化常数是 ERGM 计算上最棘手的地方——一个含 $n$ 个节点的网络有 $2^{n(n-1)/2}$ 种可能的有向图，穷举不现实。所以实际估计不靠解析解，而是用 MCMC 从分布中采样，再通过随机近似来最大化似然。

理解这一点很关键：ERGM 不是把网络统计量当自变量去拟合某个因变量，它是在给网络本身建模。参数 $\theta_k$ 的含义是——在其他统计量保持不变的情况下，统计量 $g_k$ 每增加一个单位，网络出现概率的对数几率变化。正值表示该结构在网络中比随机预期更频繁，负值则相反。

常见的统计量项包括：

选哪些项进入模型，取决于你的研究假设。不是越多越好——项之间高度相关会导致估计不稳定，后面避坑章节会展开。

2.2 Python 还是 R：两条路线的取舍

ERGM 的代码实现主要有两个生态：R 的ergm包和 Python 的pyERGM/networkx+ 自定义 MCMC。选哪个取决于你的工作流。

R 的ergm是 Statnet 项目的一部分，经过近二十年迭代，统计量项最全，MCMC 采样和退化诊断最成熟。如果你做的是标准 ERGM 分析，R 是首选。Python 生态里pyERGM相对年轻，统计量项覆盖不如 R 全面，但如果你后续要接深度学习或自定义似然，Python 的灵活性更好。

我一般建议：先用 R 的ergm跑通标准分析，确认模型设定和结果合理后，如果需要嵌入更大的 Python 流水线，再用rpy2调用或迁移到pyERGM。这样避免一上来就在 Python 里手写 MCMC 采样器，调试成本太高。

下面给出两条路线的完整代码。先看 R 版本，因为它是“完整代码”最直接的落地形态。

2.3 R 版完整代码：从网络对象到 MCMC 估计

假设你有一份边列表数据edges.csv，包含from和to两列，以及节点属性node_attr.csv，包含id和group两列。完整流程如下：

# 安装依赖（首次运行） # install.packages(c("statnet", "ergm", "network", "coda")) library(statnet) library(ergm) library(network) library(coda) # 1. 读入数据并构建 network 对象 edges <- read.csv("edges.csv", stringsAsFactors = FALSE) node_attr <- read.csv("node_attr.csv", stringsAsFactors = FALSE) net <- network(edges[, c("from", "to")], directed = TRUE, vertices = data.frame(id = node_attr$id, group = node_attr$group)) # 2. 检查网络基本特征 summary(net) # 输出节点数、边数、密度、互惠边比例等 # 3. 设定 ERGM 模型 # edges: 基础密度 # mutual: 互惠性 # gwesp(0.5, fixed=TRUE): 几何加权边共享伙伴，衰减参数 0.5 # nodematch("group"): 同组偏好 model_formula <- net ~ edges + mutual + gwesp(0.5, fixed = TRUE) + nodematch("group") + gwdegree(0.5, fixed = TRUE) # 4. MCMC 最大似然估计 set.seed(42) ergm_fit <- ergm(model_formula, control = control.ergm( MCMC.samplesize = 20000, MCMC.burnin = 10000, MCMC.interval = 1000, seed = 42 )) # 5. 查看估计结果 summary(ergm_fit) # 6. 拟合优度诊断 gof_result <- gof(ergm_fit, GOF = ~ model + degree + espartners + distance) plot(gof_result) # 7. 提取系数和置信区间 coef_df <- data.frame( term = names(coef(ergm_fit)), estimate = coef(ergm_fit), se = summary(ergm_fit)$coefficients[, "Std. Error"], pvalue = summary(ergm_fit)$coefficients[, "Pr(>|z|)"] ) print(coef_df)

这段代码的逻辑链条是：先构建network对象（ERGM 的标准输入格式），再声明模型公式，然后调用ergm()执行 MCMC 估计，最后用gof()做拟合优度检验。

参数说明几个关键点。MCMC.samplesize控制每次迭代的采样量，太小会导致估计方差大，太大会拖慢速度，20000 是中等规模网络（100-500 节点）的常用起点。MCMC.burnin是预热期，让马尔可夫链从初始状态走到稳态，一般设为 samplesize 的一半到等量。MCMC.interval是采样间隔，用来降低样本自相关，1000 是保守值，如果诊断显示自相关低可以降到 100-200。

gwesp(0.5, fixed=TRUE)里的 0.5 是衰减参数，控制长路径三角形相对于短路径的权重衰减速度。固定为 TRUE 表示不估计这个参数，直接用设定值。实践中 0.5 是常用默认，但如果你的网络聚类很强，可以试 0.2 到 0.8 之间看拟合变化。

nodematch("group")假设同组节点之间的边有额外概率，对应同质性机制。如果你的属性有多个类别，它会为每个类别估计一个参数（或者用diff=TRUE只估计差异）。

2.4 Python 版完整代码：pyERGM 与自定义 MCMC

Python 路线用pyERGM包，安装方式是pip install pyergm。如果你的环境没有这个包，也可以用networkx配合自定义 Metropolis-Hastings 采样器，但代码量会大很多。这里给出pyERGM的标准用法：

import numpy as np import pandas as pd import networkx as nx from pyergm import ERGM from pyergm.statistics import edges, mutual, gwesp, nodematch # 1. 读入边列表和节点属性 edges_df = pd.read_csv("edges.csv") node_attr = pd.read_csv("node_attr.csv") # 2. 构建 networkx 图对象 G = nx.from_pandas_edgelist(edges_df, source="from", target="to", create_using=nx.DiGraph()) # 添加节点属性 attr_dict = dict(zip(node_attr["id"], node_attr["group"])) nx.set_node_attributes(G, attr_dict, "group") # 3. 定义 ERGM 统计量 # 每个统计量是一个函数，输入图对象，输出标量值 def edges_stat(g): return g.number_of_edges() def mutual_stat(g): return sum(1 for u, v in g.edges() if g.has_edge(v, u))) // 2 def gwesp_stat(g, alpha=0.5): # 几何加权边共享伙伴 total = 0.0 for u, v in g.edges(): # 计算 u 和 v 的共同邻居 common = set(g.successors(u)) & set(g.predecessors(v)) k = len(common) if k > 0: total += (1 - (1 - alpha) ** k) return total def nodematch_stat(g, attr="group"): count = 0 for u, v in g.edges(): if g.nodes[u].get(attr) == g.nodes[v].get(attr): count += 1 return count # 4. 设定模型并估计 model = ERGM( graph=G, statistics=[edges_stat, mutual_stat, gwesp_stat, nodematch_stat], stat_names=["edges", "mutual", "gwesp", "nodematch"] ) # 5. MCMC 估计 model.fit( mcmc_samples=20000, burnin=10000, step_interval=1000, seed=42 ) # 6. 输出结果 print(model.summary()) # 7. 拟合优度：模拟网络与观测网络对比 sim_graphs = model.simulate(n=100) obs_degree = [d for _, d in G.degree()] sim_degrees = [[d for _, d in sg.degree()] for sg in sim_graphs] # 后续可以用 matplotlib 画 degree 分布对比图

Python 版本的核心差异在于统计量需要自己定义为函数。pyERGM的ERGM类接受统计量函数列表，内部用 MCMC 采样估计参数。gwesp_stat的实现里，alpha控制衰减，(1 - (1 - alpha) ** k)是几何加权项，k 是共同邻居数。这个实现和 R 的gwesp在数学上等价，但数值上可能有细微差异，因为 R 内部做了优化。

model.fit()的参数含义和 R 版本一致：mcmc_samples对应MCMC.samplesize，burnin对应MCMC.burnin，step_interval对应MCMC.interval。seed保证可复现。

如果你用的pyERGM版本没有simulate方法，可以用model.sample()或手动从估计的分布中采样。具体看包的文档。

2.5 统计量项的选择逻辑与常见组合

选统计量项不是拍脑袋，要对应你的研究问题。几个常见组合：

基础模型：edges + mutual。适合只想控制密度和互惠性，看其他结构是否显著。

传递性模型：edges + mutual + gwesp。这是社交网络分析的标配，gwesp 捕捉“朋友的朋友也是朋友”的聚类倾向。

度异质性模型：edges + gwdegree。gwdegree 控制度分布的集中程度，正值表示存在枢纽节点。

同质性模型：edges + nodematch。检验属性相似的节点是否更容易连边。

完整模型：edges + mutual + gwesp + gwdegree + nodematch。适合探索性分析，但要注意项之间的共线性。

一个实用建议：从简单模型开始，逐步加项，每次加项后看 GOF 是否改善、MCMC 诊断是否正常。不要一上来就堆五六个项，退化风险很高。

3. 参数调优与 MCMC 诊断：让估计结果可信

3.1 MCMC 采样参数怎么设：burnin、samplesize、interval

MCMC 估计的可靠性取决于马尔可夫链是否收敛到目标分布。三个核心参数控制这个过程：

burnin（预热期）：链从初始网络出发，需要足够步数才能“忘记”初始状态，进入稳态分布。太小会导致估计偏差，太大浪费计算。经验规则：burnin 至少是 samplesize 的 50%，对于复杂模型（含 gwesp、gwdegree）建议等量或更多。我一般先设 burnin=10000，看诊断图再调整。

samplesize（采样量）：从稳态分布中采集的样本数。样本越多，估计方差越小，但计算时间线性增长。10000-50000 是常见范围。对于 100 节点以下的网络，10000 通常够用；500 节点以上建议 50000 起步。

interval（采样间隔）：相邻样本之间的 MCMC 步数。目的是降低自相关——如果每步都采样，相邻样本高度相关，有效样本量远小于名义样本量。interval=1000 意味着每 1000 步取一个样本。诊断显示自相关低时可以降到 100-200。

这三个参数的组合决定了有效样本量（ESS）。ergm的summary()会输出 ESS，一般要求每个参数的 ESS > 100，理想情况 > 500。如果 ESS 太低，优先增加 samplesize 或 interval。

3.2 退化诊断：什么时候模型“炸了”

ERGM 最让人头疼的问题是模型退化（degeneracy）。现象是：MCMC 链要么跑到全空图，要么跑到全连接图，或者在不同极端之间跳变，导致参数估计不收敛或标准误巨大。

退化的根源通常是模型设定不合理。比如只用edges + triangle（三角形计数）而不加 gwesp，当 triangle 参数为正时，模型会倾向于生成大量三角形，最终导致全连接图。因为 triangle 计数没有衰减，每增加一个三角形都线性增加概率，正反馈失控。

诊断退化的方法：

# 检查 MCMC 轨迹 mcmc_diag <- mcmc.diagnostics(ergm_fit) # 输出包含： # - 各统计量的轨迹图（trace plot） # - 自相关函数（ACF） # - Gelman-Rubin 诊断（如果跑了多条链）

轨迹图应该围绕均值平稳波动，没有趋势或跳变。如果看到统计量单调上升或下降，或者在某些值之间剧烈跳变，就是退化信号。

解决退化的常见手段：

1. 换统计量：把triangle换成gwesp，把degree换成gwdegree。几何加权项有衰减机制，避免正反馈失控。
2. 加约束：用offset或constraints限制网络空间。比如~ edges + gwesp + offset(edges)固定边数。
3. 调初始值：用control.ergm(init = ...)给参数一个合理起点，避免从极端值开始。
4. 降模型复杂度：去掉相关性高的项，减少参数数量。

3.3 拟合优度检验：GOF 图怎么看

gof()函数模拟一批网络（默认 100 个），比较模拟网络和观测网络在多个统计量上的分布。输出是一组箱线图，观测值用黑点标出。

看 GOF 图的原则：观测值应该落在模拟分布的中间区域，不能是极端离群值。如果观测值在箱线图之外，说明模型没有捕捉到该维度的结构特征。

常用的 GOF 维度：

• model：模型内的统计量，应该拟合得好（否则模型设定有问题）
• degree：度分布，检验是否捕捉到度异质性
• espartners：边共享伙伴分布，检验聚类模式
• distance：测地距离分布，检验网络连通性

如果degree的 GOF 差，说明需要加gwdegree项。如果espartners差，说明gwesp的衰减参数需要调。如果distance差，可能需要加gwdegree或考虑几何加权项的组合。

一个实操细节：GOF 模拟的样本量默认 100，对于复杂模型可能不够稳定。可以设gof(ergm_fit, GOF = ~ model + degree, control = control.gof.ergm(nsim = 500))增加模拟次数。

3.4 参数解释：从系数到实质含义

ERGM 系数是 log-odds 尺度，解释时需要转换。对于edges项，系数 $\theta_{edges}$ 表示在控制其他统计量后，增加一条边对 log-odds 的贡献。但这不是直接的概率变化，因为归一化常数也变。

更直观的解释用优势比（odds ratio）：$\exp(\theta_k)$ 表示统计量 $g_k$ 每增加一个单位，网络出现的优势乘以这个因子。比如mutual系数为 2.0，则 $\exp(2.0) \approx 7.39$，表示互惠边对的优势是非互惠边对的 7.39 倍（在控制其他项后）。

对于gwesp，系数为正表示三角形闭合的倾向强于随机预期。但 gwesp 的衰减参数影响解释——衰减越快，长路径三角形的权重越低。

nodematch系数为正表示同组节点之间的边比跨组边更可能出现。$\exp(\theta_{nodematch})$ 是同组边相对于跨组边的优势比。

标准误和 p 值用于判断显著性。但要注意：ERGM 的 p 值是 Wald 检验，基于渐近正态假设。对于小样本或退化模型，p 值可能不可靠。更稳健的方法是看置信区间是否包含 0。

4. 避坑与排查：ERGM 实操中的五个血泪教训

4.1 坑一：模型退化导致估计不收敛

现象：ergm()运行时间极长，输出警告 “MCMC did not converge” 或 “Model may be degenerate”，参数估计值极大（比如 |θ| > 10），标准误爆炸。

原因：统计量项组合导致概率分布集中在极端网络（全空或全连接）。最常见的是用triangle而不加衰减，或者gwesp的衰减参数设得太大（接近 1），导致长路径三角形权重过高。

解决：把triangle换成gwesp(0.5, fixed=TRUE)或gwesp(0.3, fixed=TRUE)。如果已经用了 gwesp 还退化，降低衰减参数到 0.2-0.3。检查是否有degree项，换成gwdegree(0.5, fixed=TRUE)。如果还不行，用offset(edges)固定边数，减少模型自由度。

4.2 坑二：MCMC 自相关过高导致 ESS 不足

现象：summary(ergm_fit)显示 ESS 远小于 samplesize，比如 samplesize=20000 但 ESS=50。参数标准误很大，置信区间宽。

原因：MCMC 采样间隔太小，相邻样本高度相关。或者模型本身混合速度慢（统计量项之间相关性高）。

解决：增大MCMC.interval，从 1000 提到 5000 甚至 10000。增大MCMC.burnin，让链充分预热。如果还不行，简化模型——去掉相关性高的项。可以用mcmc.diagnostics()看自相关图，确认自相关降到 0.1 以下再采样。

4.3 坑三：GOF 图显示模型拟合差但不知道加什么项

现象：gof()输出中degree或espartners的观测值在模拟分布之外，但不确定该加哪个统计量。

原因：模型设定遗漏了重要结构特征。比如度分布拟合差说明没有捕捉度异质性，espartners 拟合差说明聚类模式没捕捉到。

解决：按 GOF 维度对应加项。degree差加gwdegree。espartners差加gwesp或调衰减参数。distance差加gwdegree或gwdistance。nodematch差加对应的属性匹配项。每次加一项，重新跑 GOF，直到所有维度都拟合可接受。

4.4 坑四：网络对象构建时节点属性丢失

现象：nodematch("group")报错 “vertex attribute not found”，或者估计结果中 nodematch 系数为 NA。

原因：network()函数构建对象时，vertices参数的 data.frame 行数必须和节点数一致，且 id 列要和边列表中的节点 id 匹配。如果节点属性表有重复 id 或缺失 id，属性会丢失。

解决：构建 network 对象前，先检查node_attr$id是否唯一，是否覆盖了边列表中所有节点。用set.vertex.attribute(net, "group", node_attr$group)手动设置属性，确保对齐。构建后用list.vertex.attributes(net)确认属性存在。

4.5 坑五：Python 和 R 结果不一致

现象：同样的数据和模型设定，Python 的pyERGM和 R 的ergm给出的系数估计差异较大。

原因：MCMC 采样器实现不同，随机种子不同，统计量计算方式有细微差异（比如 gwesp 的衰减实现）。另外，R 的ergm默认用 “Stochastic Approximation” 优化，Python 可能用不同的优化器。

解决：不要期望完全一致。检查两边的统计量计算是否数学等价——用同一个网络对象，分别计算 edges、mutual、gwesp 的值，确认数值一致。如果统计量一致但估计不同，增大 samplesize 和 burnin，让两边都充分收敛。如果还差异大，以 R 的ergm为准，因为它的 MCMC 诊断工具更成熟。

5. 进阶技巧：用模拟网络验证模型与自定义统计量

5.1 从估计模型生成模拟网络：验证与预测

跑完 ERGM 后，一个常被忽略的步骤是用估计的模型生成模拟网络，看模拟网络是否重现观测网络的关键特征。这比 GOF 更灵活，因为你可以自定义比较维度。

# 从估计模型模拟 100 个网络 set.seed(123) sim_nets <- simulate(ergm_fit, nsim = 100, seed = 123) # 提取每个模拟网络的统计量 sim_stats <- t(sapply(sim_nets, function(net) { c( edges = network.edgecount(net), mutual = summary(net ~ mutual), gwesp = summary(net ~ gwesp(0.5, fixed = TRUE)), nodematch = summary(net ~ nodematch("group")) ) })) # 观测网络的统计量 obs_stats <- c( edges = network.edgecount(net), mutual = summary(net ~ mutual), gwesp = summary(net ~ gwesp(0.5, fixed = TRUE)), nodematch = summary(net ~ nodematch("group")) ) # 比较：观测值在模拟分布中的分位数 for (s in names(obs_stats)) { pct <- mean(sim_stats[, s] <= obs_stats[s]) cat(s, ": observed =", obs_stats[s], ", percentile =", round(pct * 100, 1), "%\n") }

这段代码的逻辑是：用simulate()从估计模型生成 100 个网络，计算每个网络的统计量，然后看观测值落在模拟分布的哪个分位。如果观测值在 5%-95% 之外，说明模型在该维度上拟合不好。

这个方法的优势是你可以自定义任何比较维度——比如网络中的三元组数量、特定节点的度、属性混合矩阵的某个单元格。GOF 只覆盖预定义的维度，模拟验证更灵活。

5.2 自定义统计量：当内置项不够用

R 的ergm允许用户自定义统计量项。假设你想检验“跨组边是否倾向于形成三角形”，可以定义一个组合统计量：

# 自定义统计量：跨组三角形数量 # 注册到 ergm 的统计量系统 InitErgmTerm.crosstriangle <- function(nw, arglist, ...) { a <- check.ErgmTerm(nw, arglist, varnames = c("attrname"), vartypes = c("character"), defaultvalues = list(NULL), required = c(TRUE)) attrname <- a$attrname node_attr <- get.vertex.attribute(nw, attrname) list( name = "crosstriangle", coef.names = paste0("crosstriangle.", attrname), inputs = c(length(unique(node_attr))), dependence = TRUE, minval = 0 ) } # C 语言层面的计算函数需要额外实现 # 这里只展示 R 层面的注册逻辑

自定义统计量需要同时实现 R 层面的注册和 C 层面的计算函数，门槛较高。更实用的替代方案是用ergm的userterms接口，或者用statnet的ergm.userterms包。如果只是探索性分析，建议先用内置项的线性组合近似，比如用nodematch+gwesp的交互来近似跨组三角形。

Python 的pyERGM自定义统计量更简单，因为统计量就是 Python 函数，直接传给ERGM类即可。但要注意计算效率——Python 循环比 C 慢很多，大网络可能跑不动。

5.3 模型比较：AIC/BIC 与交叉验证

多个候选模型怎么选？ergm提供 AIC 和 BIC：

# 拟合多个模型 fit1 <- ergm(net ~ edges + mutual) fit2 <- ergm(net ~ edges + mutual + gwesp(0.5, fixed = TRUE)) fit3 <- ergm(net ~ edges + mutual + gwesp(0.5, fixed = TRUE) + gwdegree(0.5, fixed = TRUE)) # 比较 AIC/BIC anova(fit1, fit2, fit3) # 或 AIC(fit1, fit2, fit3) BIC(fit1, fit2, fit3)

AIC 惩罚参数数量，BIC 惩罚更重。对于 ERGM，AIC 可能偏向复杂模型，因为 MCMC 估计的似然本身有噪声。BIC 更保守。我一般两个都看，如果结论一致就选，不一致就优先看 GOF 和实质解释性。

交叉验证在 ERGM 里比较少见，因为网络数据通常只有一个观测。但可以用边留一法：随机去掉 10% 的边，用剩余边拟合模型，预测去掉的边是否存在。这个方法计算量大，适合小网络。

5.4 一个实用习惯：先跑空模型再逐步加项

我自己的习惯是：拿到数据后先跑edges空模型，看密度参数是否合理。然后加mutual，看互惠性是否显著。再加gwesp，看聚类是否显著。每加一项，检查 MCMC 诊断和 GOF。如果某一步退化或 GOF 变差，就回退，换统计量或调参数。

这个流程看起来慢，但比一上来跑完整模型然后花几小时调试退化要快得多。ERGM 的玄学在于，模型设定和 MCMC 诊断是耦合的——参数估计不可信时，你无法判断是模型错了还是采样不够。逐步加项能帮你定位问题出在哪一步。

希望帮到你。

本文还有配套的精品资源，点击获取