名师讲堂|使用 R 语言计算上市公司供应链地理加权距离

指标来源

供应链地理加权距离指标来自邹颖、石福安、祁亚发表在《世界经济》2026 年第 1 期的论文《以数促联:公共数据开放与企业供应链地理布局》。

该论文采用堆叠双重差分模型考察公共数据开放对企业供应链地理布局的影响,其中核心被解释变量即为供应链地理加权距离,包括两个维度:

  • 供应商加权距离(Disws):衡量上市企业与主要供应商之间的地理距离(加权);
  • 客户加权距离(Diswc):衡量上市企业与主要客户之间的地理距离(加权)。

论文发现,公共数据开放能打破地理距离约束,拓宽企业供应链分布范围。

指标定义与计算公式

数据来源

该指标的测算需要三类数据:

  1. 上市公司注册地址与办公地址数据:包含 2000~2024 年所有沪深 A 股上市公司的注册地址和办公地址经纬度信息,以及所处省市区县。
  2. 上市公司前 5 大供应商数据:包含 2001~2024 年上市公司前 5 大供应商的工商注册信息匹配结果,含供应商经纬度、采购额及采购额占比。
  3. 上市公司前 5 大客户数据:包含 2001~2024 年上市公司前 5 大客户的工商注册信息匹配结果,含客户经纬度、销售额及销售额占比。

计算公式

地理距离计算方法

地理距离采用 sf 包的 st_distance() 函数计算大圆距离(Haversine 公式)。具体而言:

  1. 将上市公司办公地址的经纬度和供应商/客户地址的经纬度分别转换为 sf 的 POINT 几何对象,使用 WGS84 坐标系(CRS 4326);
  2. st_distance() 在 WGS84 坐标系下自动计算大圆距离(单位为米),然后转换为千米(km);
  3. 上市公司地址使用的是办公地址(而非注册地址),这与论文中基于实际经营地址的选择一致。

加权方式说明

权重采用的是供应商采购额占前五大供应商采购总额的比例(而非供应商采购额占企业总采购额的比例)。这是因为上市公司年报中通常只披露前五大供应商/客户的名称和交易金额,无法获取完整的采购/销售总额数据。

计算过程

第 1 步:参数设置与数据读取

首先加载所需的 R 包并设置数据路径:

library(tidyverse)
library(sf)
library(haven)

dir <- "~/Desktop/使用 R 语言计算上市公司供应链地理加权距离/"

gys_file <- paste0(dir, "2001~2024年上市公司前5大供应商工商注册信息匹配结果(含经纬度及所处的省市区县).dta")
kh_file <- paste0(dir, "2001~2024年上市公司前5大客户工商注册信息匹配结果(含经纬度及所处的省市区县).dta")
addr_file <- paste0(dir, "2000~2024年上市公司注册地址与办公地址(含经纬度、所处的省市区县及搬迁距离).dta")

然后读取三类数据:

cat("读取供应商数据...\n")
gys <- read_dta(gys_file)
cat(sprintf(" 供应商记录数: %d, 唯一公司数: %d\n", nrow(gys), gys %>% distinct(股票代码) %>% nrow()))

cat("读取客户数据...\n")
kh <- read_dta(kh_file)
cat(sprintf(" 客户记录数: %d, 唯一公司数: %d\n", nrow(kh), kh %>% distinct(股票代码) %>% nrow()))

cat("读取上市公司地址数据...\n")
addr <- read_dta(addr_file)
cat(sprintf(" 地址记录数: %d, 唯一公司数: %d\n", nrow(addr), addr %>% distinct(股票代码) %>% nrow()))

第 2 步:准备上市公司办公地址

从地址数据中提取上市公司办公地址的经纬度,并过滤缺失值:

addr_sub <- addr %>%
select(股票代码, 统计年份 = 年份, 办公地址_经度, 办公地址_纬度) %>%
filter(!is.na(统计年份), !is.na(办公地址_经度), !is.na(办公地址_纬度))

第 3 步:定义地理距离计算函数

核心函数使用 sf 包将经纬度转换为 POINT 几何对象,然后用 st_distance() 计算大圆距离:

calc_distance <- function(df, firm_lon_col, firm_lat_col, partner_lon_col, partner_lat_col) {
# 构建公司 POINT
firm_pts <- df %>%
select(firm_lon = all_of(firm_lon_col), firm_lat = all_of(firm_lat_col)) %>%
st_as_sf(coords = c("firm_lon", "firm_lat"), crs = 4326)

# 构建合作伙伴 POINT
partner_pts <- df %>%
select(part_lon = all_of(partner_lon_col), part_lat = all_of(partner_lat_col)) %>%
st_as_sf(coords = c("part_lon", "part_lat"), crs = 4326)

# st_distance 返回米为单位的矩阵,转 km
dist_m <- st_distance(firm_pts, partner_pts, by_element = TRUE)
as.numeric(dist_m) / 1000
}

注意:by_element = TRUE 表示逐元素计算距离(即第 i 个公司到第 i 个供应商),而非生成距离矩阵。st_as_sf() 中 crs = 4326 指定 WGS84 坐标系,st_distance() 会据此自动采用大圆距离公式。

第 4 步:计算供应商加权距离

筛选有效记录(年报、有经纬度、有采购额),匹配上市公司办公地址:

此处代码需下载讲义材料查看~

关键逻辑说明:首先在同一公司-年度组内计算前五大供应商的采购总额 raw_total,然后用每个供应商的采购额除以总额得到权重 ratio_s。最后用 log(1 + sum(dist_km * ratio_s)) 得到加权距离的对数值,这既压缩了极端值,又保证了取值非负。

第 5 步:计算客户加权距离

客户加权距离的计算逻辑与供应商完全对称,只是将”供应商采购额”替换为”客户销售额”:

此处代码需下载讲义材料查看~

第 6 步:合并结果并输出

将供应商加权距离和客户加权距离合并到上市公司信息表中:

info <- addr %>%
select(股票代码, 统计年份 = 年份, 股票简称, 行业代码C,
办公地址_经度, 办公地址_纬度, 办公地址_省, 办公地址_市)

result <- info %>%
left_join(gys_weighted, by = c("股票代码", "统计年份")) %>%
left_join(kh_weighted, by = c("股票代码", "统计年份")) %>%
arrange(股票代码, 统计年份)

# 去除 2000 年数据(供应商数据从 2001 年开始)
result <- result %>% filter(统计年份 >= 2001)

outfile <- paste0(dir, "2001~2024年上市公司供应链地理加权距离.dta")
write_dta(result, outfile, label = "数据处理:微信公众号 RStata")

result

描述性统计

论文中报告的描述性统计结果如下:

变量 样本量 均值 1/4 分位数 中位数 3/4 分位数 标准差
Disws 6,507 5.940 5.385 6.283 6.862 1.297
Diswc 9,518 5.967 5.528 6.412 6.915 1.416

可以看到:

  • 供应商加权距离(Disws)和客户加权距离(Diswc)的中位数分别为 6.283 和 6.412,均值分别为 5.940 和 5.967;
  • 由于取了对数变换,中位数大于均值说明存在明显的左偏分布;
  • 客户加权距离的标准差(1.416)略大于供应商加权距离(1.297),说明不同企业与客户的地理距离差异更大。

R 语言脚本

完整的 R 语言脚本如下(已整理为可直接运行的格式):

# ============================================================
# 计算上市公司与供应商/客户的加权地理距离
# 参考论文:以数促联:公共数据开放与企业供应链地理布局(邹颖等, 2026)
#
# 公式:
# Disw_s_it = ln(1 + Σ Dis_s_ipt × Ratio_s_ipt)
# Disw_c_it = ln(1 + Σ Dis_c_iqt × Ratio_c_iqt)
#
# 距离计算:sf::st_distance(大圆距离 / Haversine)
# 地址选择:上市公司使用办公地址
# 数据计算:微信公众号 RStata
# ============================================================

library(tidyverse)
library(sf)
library(haven)

# ---- 1. 参数设置 ----
dir <- "~/Desktop/使用 R 语言计算上市公司供应链地理加权距离/"

gys_file <- paste0(dir, "2001~2024年上市公司前5大供应商工商注册信息匹配结果(含经纬度及所处的省市区县).dta")
kh_file <- paste0(dir, "2001~2024年上市公司前5大客户工商注册信息匹配结果(含经纬度及所处的省市区县).dta")
addr_file <- paste0(dir, "2000~2024年上市公司注册地址与办公地址(含经纬度、所处的省市区县及搬迁距离).dta")

# ---- 2. 读取数据 ----
cat("读取供应商数据...\n")
gys <- read_dta(gys_file)
cat(sprintf(" 供应商记录数: %d, 唯一公司数: %d\n", nrow(gys), gys %>% distinct(股票代码) %>% nrow()))

cat("读取客户数据...\n")
kh <- read_dta(kh_file)
cat(sprintf(" 客户记录数: %d, 唯一公司数: %d\n", nrow(kh), kh %>% distinct(股票代码) %>% nrow()))

cat("读取上市公司地址数据...\n")
addr <- read_dta(addr_file)
cat(sprintf(" 地址记录数: %d, 唯一公司数: %d\n", nrow(addr), addr %>% distinct(股票代码) %>% nrow()))

# ---- 3. 准备上市公司办公地址 ----
cat("准备上市公司办公地址...\n")
addr_sub <- addr %>%
select(股票代码, 统计年份 = 年份, 办公地址_经度, 办公地址_纬度) %>%
filter(!is.na(统计年份), !is.na(办公地址_经度), !is.na(办公地址_纬度))
cat(sprintf(" 有效地址记录: %d\n", nrow(addr_sub)))

# ---- 4. 计算地理距离的核心函数 ----
# 使用 sf 的 st_distance 计算大圆距离(Haversine),结果单位为米
# 这里将经纬度转为 sf POINT(CRS 4326),然后用 st_distance
calc_distance <- function(df, firm_lon_col, firm_lat_col, partner_lon_col, partner_lat_col) {
# 构建公司 POINT
firm_pts <- df %>%
select(firm_lon = all_of(firm_lon_col), firm_lat = all_of(firm_lat_col)) %>%
st_as_sf(coords = c("firm_lon", "firm_lat"), crs = 4326)

# 构建合作伙伴 POINT
partner_pts <- df %>%
select(part_lon = all_of(partner_lon_col), part_lat = all_of(partner_lat_col)) %>%
st_as_sf(coords = c("part_lon", "part_lat"), crs = 4326)

# st_distance 返回米为单位的矩阵,转 km
dist_m <- st_distance(firm_pts, partner_pts, by_element = TRUE)
as.numeric(dist_m) / 1000
}

# ---- 5. 计算供应商加权距离 ----
# 此处代码需下载讲义材料查看~

# ---- 6. 计算客户加权距离 ----
# 此处代码需下载讲义材料查看~

# ---- 7. 合并结果 ----
cat("\n========== 合并结果 ==========\n")

info <- addr %>%
select(股票代码, 统计年份 = 年份, 股票简称, 行业代码C,
办公地址_经度, 办公地址_纬度, 办公地址_省, 办公地址_市)

result <- info %>%
left_join(gys_weighted, by = c("股票代码", "统计年份")) %>%
left_join(kh_weighted, by = c("股票代码", "统计年份")) %>%
arrange(股票代码, 统计年份)

result

cat(sprintf(" 最终样本: %d 个公司-年度观测值\n", nrow(result)))
cat(sprintf(" 唯一公司数: %d\n", result %>% distinct(股票代码) %>% nrow()))
cat(sprintf(" 年份范围: %d - %d\n",
min(result$统计年份, na.rm = TRUE),
max(result$统计年份, na.rm = TRUE)))

# ---- 8. 描述性统计 ----
cat("\n========== 描述性统计 ==========\n")
cat("--- Disw_s(供应商加权距离)---\n")
ds <- result$Disw_s %>% na.omit()
cat(sprintf(" N = %d\n", length(ds)))
cat(sprintf(" 均值 = %.3f, 中位数 = %.3f, 标准差 = %.3f\n",
mean(ds), median(ds), sd(ds)))
cat(sprintf(" 最小值 = %.3f, 最大值 = %.3f\n", min(ds), max(ds)))

cat("--- Disw_c(客户加权距离)---\n")
dc <- result$Disw_c %>% na.omit()
cat(sprintf(" N = %d\n", length(dc)))
cat(sprintf(" 均值 = %.3f, 中位数 = %.3f, 标准差 = %.3f\n",
mean(dc), median(dc), sd(dc)))
cat(sprintf(" 最小值 = %.3f, 最大值 = %.3f\n", min(dc), max(dc)))

# ---- 9. 去除 2000 年数据 ----
result <- result %>% filter(统计年份 >= 2001)

# ---- 10. 保存结果 ----
outfile <- paste0(dir, "2001~2024年上市公司供应链地理加权距离.dta")
write_dta(result, outfile)
cat(sprintf("\n结果已保存至: %s\n", outfile))
cat("完成!\n")

参考文献

邹颖、石福安、祁亚,2026:《以数促联:公共数据开放与企业供应链地理布局》,《世界经济》第 1 期。

点击这里跳转到 RStata 短书平台获取附件:名师讲堂|使用 R 语言计算上市公司供应链地理加权距离

评论