
简介针对2018年全球主要贸易关系面向国际贸易网络研究者与经济学课程学习者提供一套基于R语言的网络分析脚本完整覆盖网络密度、平均路径长度、传递性、互易性、分类性、自旋玻璃社区检测、结构等效性及指数随机图模型等核心指标与算法。压缩包内含1个R脚本仅7KB轻量易用适合有一定R和网络分析基础的读者阅读、调试或二次开发也可作为经济地理或国际经济学课程的实验代码。目前已有233人学习脚本按指标分段组织函数与注释较完整能够帮助使用者理解各网络指标的计算逻辑与结果含义并据此复现2018年全球贸易格局的实证分析。通过研读这份脚本还可掌握ERGM等复杂网络模型的R实现要点为后续扩展研究提供参考。 做贸易研究的人早晚会遇到一个问题整个世界那么多国家彼此之间买进卖出这种千丝万缕的关系能不能像社交网络一样被分析答案是肯定的。TradeNetworkDistributionsAnalysis 这个项目做的就是这件事——用 R 语言网络分析把 2018 年主要国家之间的双边贸易关系整理成一张有向网络然后从密度、平均路径长度、传递性、互易性、分类性这些基础指标一直算到自旋玻璃社区检测算法、结构等效性和指数随机图模型ERGM。这篇文章会把每个指标的计算逻辑、R 代码和结果解读都展开讲适合想用网络视角研究贸易数据或者正在做类似 R 数据项目的读者参照。在我自己抠这个项目的过程中印象最深的是R 里做网络分析函数调用只是最简单的一层真正花时间的是数据清洗、阈值设定和结果解释。下面我会按项目实操的顺序来说尽量把容易踩的坑提前标出来。1. 数据准备与网络构建1.1 贸易数据长什么样别急着写代码先搞清楚数据格式。项目里用的数据源可以是 UN Comtrade 或者 CEPII-BACI都是公开的国际贸易数据库。通常下载下来是一个长表每一行代表一条双边贸易记录关键字段包括出口国、进口国、年份、贸易额可能还有产品类别。而且这类数据有几个共同毛病国家代码不统一有的用 ISO3有的用英文全名同一对国家的贸易额在不同表格里可能重复部分年份缺失严重。所以需要先做三步聚合和清洗。第一步如果数据是 HS 编码或者 SITC 编码的产品级数据要按 exporter-importer-year 聚合成国家层面的总贸易额第二步删除国家自己对自己的记录也就是排除掉 exporter importer 的行第三步把明显是缺失值或零值的记录统一处理成 0但要注意有些 0 是真的没有贸易有些是未报告。我自己的原则是单独保留一份“未清洗原始数据”所有清洗步骤都在新表里操作这样可以随时回溯哪一步出了问题。清洗完之后我们会得到一个 n×n 的贸易矩阵行表示出口国列表示进口国。矩阵里第 i 行第 j 列的数值就是 i 国出口到 j 国的贸易额。这个矩阵就是后续所有网络指标的基础。要注意矩阵的行列顺序必须保持一致否则后面算分类性、结构等效性时结果会出现莫名其妙的错位。1.2 构建有向加权网络R 里做网络分析我推荐优先用 igraph。原因是它经过大量工程优化像密度、互易性、社区检测这些计算都是底层 C 实现处理几十个国家的数据非常快图形输出也能跟 ggraph 无缝衔接。当然也可以只装 tidygraph 走数据管道风格但底层还是 igraph所以先熟悉 igraph 是更稳妥的路线。核心代码大概是这样的library(igraph) library(dplyr) library(tidyr) # 假设 trade_df 是清洗后的长表包含 exporter, importer, value g - graph_from_data_frame(trade_df, directed TRUE) # 删除自环和重复边 g - simplify(g, remove.loops TRUE, remove.multiple TRUE) # graph_from_data_frame 会自动把 value 作为边属性 E(g)$weight - E(g)$value如果你手里的数据是矩阵也可以用graph_from_adjacency_matrix(as.matrix(trade_mat), mode directed, weighted TRUE, diag FALSE)。这里diag FALSE可以直接忽略对角线省得之后还要清自环。有一个细节我踩过坑如果贸易额的单位差异很大比如有的数据源报美元有的报千美元那后面算权重相关指标时会非常离谱。最好在构建网络前统一单位并用log1p(value)取对数缓解长尾分布直接log(value)会在 0 值上报-Inf。2. 核心网络指标的计算与解读2.1 网络密度与平均路径长度网络密度是最直观的指标用来衡量这张网络有多“稠密”。在有向网络中n 个节点之间的最大可能边数是 n*(n-1)不考虑自环密度就是实际存在边数除以这个数。在 2018 年全球贸易场景下即使只取主要经济体密度一般也不会太高因为不是任意两个国家都有直接的大额贸易往来。R 代码很简单edge_density(g) # 整体密度 edge_density(g, loops FALSE) # 忽略自环效果同上 mean_distance(g, directed TRUE, unconnected TRUE)mean_distance算的是所有节点对之间的最短路径长度平均值。在贸易网络里如果平均路径长度在 2 左右说明任意两个国家平均只要通过一个中间国家就能建立贸易联系。这个“小世界”特征非常重要它意味着某个国家的供给冲击或需求波动很容易通过网络传导到其他国家而不是只影响直接贸易伙伴。顺带一提unconnected TRUE这个参数不能省。当网络里有孤立节点或者有向图里存在不可达的节点对时默认会计算距离为无穷大不设置的话要么报错要么会得到 Inf后面就全乱了。我在第一次跑的时候就是因为有一批小岛国没进入任何主要贸易关系导致结果直接变废。2.2 传递性与互易性传递性是三角关系的闭合程度。换成大白话A 国出口到 B 国B 国出口到 C 国那么 A 和 C 之间是否也会发生贸易如果这种“朋友的朋友也是朋友”的模式远高于随机水平说明网络中存在明显的三角循环或传递结构。transitivity(g, type global)这里有个细节transitivity默认处理的是无向图。如果你传的是有向图它会把有向边按无向边处理这其实是可以接受的因为它关心的只是“连接是否存在”。但如果想严格区分方向需要自己对邻接矩阵做逻辑运算比如计算 A→B、B→C、A→C 同时发生的概率。大多数贸易网络分析里无向的全局传递性已经够用。互易性则专门针对有向网络衡量双向贸易关系的比例我出口给你你是不是也出口给我reciprocity(g, ignore.loops TRUE)reciprocity()返回值在 0 到 1 之间0.7 说明有 70% 的贸易关系是双向对等的。现实世界中全球贸易网络的互易性通常比较高因为大多数国家之间既有出口也有进口只是金额不一定对等。这也是后面 ERGM 里mutual项往往显著的原因。要注意reciprocity如果不设置ignore.loops TRUE网络里的自环会被当成互易边结果虚高。我一开始就用默认值跑了半天才发现数字不太对。3. 分类性分析与结构等效性3.1 分类性强节点是否更倾向连接强节点分类性assortativity在网络分析里常用来看节点之间的“门当户对”程度。贸易网络里我们最关心两个版本一个是基于度的分类性另一个是基于属性的名义分类性。基于度的分类性看的是贸易规模大的国家是否更倾向于彼此连接。计算方法在 igraph 里是assortativity_degree(g, directed TRUE)返回值如果是正的说明“强强联手”明显如果是负的说明大型经济体更多是向小型经济体延伸连接比如典型的枢纽-辐条结构。现实中全球贸易网络往往呈现正的度分类性因为核心国家之间的贸易流量实在太大了。我跑 2018 年数据时得到正数和预期一致。如果想看区域属性可以用assortativity_nominal()例如按联合国区域分组V(g)$region - countries$region assortativity_nominal(V(g)$region, directed TRUE)正值表示相同区域的国家之间更容易形成贸易连接这就从网络层面验证了“区域化倾向”。这一步看起来很轻量但能为后面社区检测的结果提供交叉验证。如果分类性为正那么社区检测大概率会得到几个明显聚在一起的区域集团两者可以相互印证。3.2 结构等效性找网络中的“角色替补”结构等效性structural equivalence是一个和社区检测完全不同的分析。它不关心节点是否紧密连接而是关心两个节点是否拥有相似的“关系画像”如果我连接到的一堆国家和另外某个国家连接的几乎一样那这两个国家在网络中的角色就可能是同类的。实现思路是这样的用邻接矩阵的每一行代表一个节点的“出向关系画像”每一列代表“入向关系画像”然后计算节点两两之间的相似度或距离。常见做法是把行列拼成一个向量计算 Pearson 相关或欧氏距离。R 里可以这样跑adj - as_adjacency_matrix(g, sparse FALSE) profiles - cbind(adj, t(adj)) # 同时考虑出边和入边 sim_matrix - cor(t(profiles)) # 转换成距离并做层次聚类 dist_mat - as.dist(1 - sim_matrix) hc - hclust(dist_mat, method ward.D2) plot(hc, hang -1, cex 0.8)聚类出来之后可以按类别个数cutree(hc, k 4)给每个国家打上“角色标签”。这些标签和社区标签不一样一个小组里的国家可能分布在全球不同区域但因为它们都和同样的中间国家贸易所以被视为结构等效。这一步在实际项目中很适合做国家分类面板比如判断哪些国家是“转口贸易枢纽”哪些是“最终需求市场”。不过要注意如果有太多孤立节点会把相似度矩阵拉低建议先剔除完全没有连边的节点再跑。4. 自旋玻璃社区检测实战4.1 从统计物理到贸易集团自旋玻璃spin glass社区检测算法是我在这个项目里比较兴奋的一部分。它源自统计物理里的 Ising 模型和 Potts 模型基本思想是把每个节点看成一个“自旋”自旋的状态就是它所属的社区编号节点之间的边会对自旋状态产生约束如果一条边的两个端点被分到同一个社区系统能量就比较低如果分到不同社区能量就比较高。算法在随机扰动和局部更新中不断寻找使总能量最低的划分方式。用这个算法来识别贸易网络好处是它不要求社区内部特别稠密而是可以捕捉到一些由权重和链接模式共同决定的隐含集团。通俗点说它能帮我们看到 2018 年的贸易网络自然地形成了哪几个“派系”而且每个派系内部往往带有区域贸易协定的影子。相比 walktrap 或 label propagation自旋玻璃在面对权重分布不均的贸易数据时社区边界更干净。但这里必须提醒igraph 里的cluster_spinglass()在标准版本里主要支持无向图。有向加权网络要先决定是转成无向还是抽取其中的互惠边来跑。如果你硬传有向图有些版本会忽略方向或直接报错。我处理的方法是把有向边权取两条方向的平均值生成一个无向加权图这样信息损失最小。4.2 参数调节与结果可视化实操时我通常先用walktrap.community()或者fastgreedy.community()快速看一下社区大概有多少个再把数量作为spins的参考值传给cluster_spinglass()。set.seed(2024) sg - cluster_spinglass(g, weights E(g)$weight, spins 8, gamma 0.8) sg_membership - membership(sg) sg_modularity - modularity(sg)spins是社区数上限不是精确值。gamma控制社区粒度gamma 越大越容易分出更多小社区。根据我的经验gamma 在 0.5 到 1.5 之间变化比较敏感建议多试几组记录每个 gamma 下的模块度最后选模块度最高的那组参数。出来结果之后常规做法是画社区图plot(sg, g, vertex.label V(g)$name, vertex.size 8, edge.arrow.size 0.2, col membership(sg))不过直接用plot只是一张静态图。更好的方式是导出membership字段合并到原始数据框里用来做后续的统计或地图可视化。我还会顺便计算每个社区的总出口额占比看看哪些集团主导了当年的贸易流。这个结果往往比单纯看模块度更有业务价值。一个常见的坑自旋玻璃算法不是确定性的不同随机种子会得到不同社区划分。所以不论是用脚本还是跑分析都建议固定set.seed()并且在报告中注明种子值否则别人复现时对不上。5. 指数随机图模型ERGM建模5.1 ERGM 到底在做什么前面算的密度、互易性、分类性本质上都是“描述性指标”。它们告诉我们网络长什么样但不能告诉我们这些结构是不是显著更不能直接告诉我们是什么机制生成了这样的网络。指数随机图模型ERGM就是来解决这个问题的。ERGM 的建模思路是把当前观察到的网络视为众多可能网络中的一个实现用一批网络统计量当作解释变量比如边数、互惠边数、三角形数量、节点属性是否相似然后用最大似然估计或者 MCMC 估计出这些机制的回归系数。用贸易网络的话说就是在控制其他条件的情况下某种连接模式是否显著地比随机网络更常见。因为 ERGM 一般处理二值网络所以需要先把加权贸易网络二值化。我这里的做法是保留两条国家之间贸易额大于中位数的连接定义成 1否则定义成 0。这个阈值直接影响结果后面我会再强调。5.2 模型拟合与系数解读R 里做 ERGM 的主力军是statnet套件的ergm包。先要把 igraph 对象转成network对象然后加入节点属性。library(statnet) library(ergm) net - asNetwork(g) set.vertex.attribute(net, region, V(g)$region) set.vertex.attribute(net, gdp, V(g)$gdp) # 如果有 GDP 数据然后拟合模型fit - ergm(net ~ edges mutual triangles nodematch(region), control control.ergm(seed 1, MCMC.burnin 1000, MCMC.interval 200)) summary(fit)输出里最关键的是 Estimate 和 Pr(|z|)。比如edges系数为负说明在给定其他项的条件下形成边的概率比随机机会低所以网络整体比较稀疏mutual系数显著为正说明双向连接比单向连接更容易出现triangles显著为正说明网络中存在大量三角闭合结构nodematch(region)显著为正说明同一区域内的国家之间更容易建立贸易连接。还要看 AIC/BIC比较不同模型时越低越好。拟合之后最好再跑一下gof(fit)用拟合优度检验看模型能不能捕捉网络的关键特征比如边数分布、度数分布、最短路径分布。如果拟合样本和观察网络差太远说明模型漏掉了重要机制需要加项。这里我必须提醒ERGM 的参数估计对网络规模非常敏感。节点超过 100 个、边超过 1000 条之后MCMC 收敛会变得很慢尤其在带triangles的模型里。你可以用control.ergm(MCMC.samplesize 2000)或parallel参数来加速。如果还不行就要考虑用一个子网络或者先合并一些节点。6. 常见问题与避坑经验6.1 数据预处理中的隐性坑整理这个项目的过程中最常遇到的问题有三个一是原始数据的国家代码不统一比如有的数据源用 ISO3有的用中文名合并时会莫名丢失很多行二是贸易额为 0 的情况没有提前判断导致取对数时出现-Inf三是某些年份某条边只有单向记录另一方显示 0在构建有向网络时这个 0 和缺失值在后续 ERGM 中含义完全不同。我的处理原则是单独保存一份加权原始矩阵保留真实 0二值化矩阵时0 和缺失都处理成 0但要在文档里写明“未报告视为无贸易”。虽然这在统计上有点粗但至少可复现。更好一点的做法是如果某条边在数据库里完全没有记录可以先用 0 填充然后在 ERGM 里通过设置edge.missings来单独处理缺失但这会增加不少复杂性。6.2 指标计算不一致的说明igraph 和 statnet 很多函数名容易混淆比如 igraph 的reciprocity()和 statnet 的mutual都处理互惠但前者是比率后者是模型中一个参数igraph 的transitivity()和 statnet 的transitiveties也是两个概念。如果你同时加载两个包注意函数名前缀冲突建议用igraph::reciprocity()这种写法明确调用。另外edge_density()在有向网络里自动除以n*(n-1)不会把自环算进去。有人误以为密度应该除以n*n这是错的。写报告时最好标注清楚公式。平均路径长度也要注明是 “有向平均距离” 还是 “无向平均距离”两者结果差异很大尤其是在方向性明显的贸易网络中。6.3 常见问题速查现象可能原因解决方式密度异常低或高数据未聚合、自环未删按国家-年份聚合simplify 删除自环mean_distance 返回 Inf有孤立节点/不可达节点对设置 unconnected TRUE或剔除孤立点互易性接近 1未设置 ignore.loops设置 remove.loops / ignore.loopsspin glass 每次结果不同随机初始化set.seed 固定种子并注明种子值ERGM 一直不收敛网络太大或模型太复杂增大 MCMC、简化模型、抽样或合并节点ERGM 结果符号和预期相反二值化阈值不合适尝试不同阈值做敏感性分析做完 TradeNetworkDistributionsAnalysis 这个项目我的体会是R 网络分析的价值不在于一句函数跑出漂亮数字而在于你能从不同角度反复确认同一个故事。密度和路径长度告诉你网络骨架互易性和传递性告诉你关系模式自旋玻璃社区告诉你集团边界ERGM 则把这些现象还原成可检验的机制。如果你想在自己的数据上复现建议按这个顺序来不要跳过结构等效性它可以帮社区检测结果做解释。最后再分享一个小技巧二值化网络做 ERGM 前先把阈值画成曲线看不同阈值下网络密度和互惠比例的变化再选一个相对稳定的平台区间。这样模型结果不会因为某个随意定的阈值而变得不可信。本文还有配套的精品资源点击获取