今天给大家分享使用 R 语言测算各城市数字产业集聚程度的方法。该方法参考屠西伟、史丹（2025）《数字产业集聚与企业能源效率改进》，通过区位熵来综合测度城市的数字产业集聚水平。

附件中提供了该参考文献的 PDF 文件，感兴趣的小伙伴可以阅读原文。

指标来源与计算原理

数字产业集聚度（Location Quotient）

区位熵的经济含义

DL > 1：该城市数字产业集聚度高于全国平均水平，具有相对专业化优势
DL = 1：与全国平均水平相当
DL < 1：低于全国平均水平

两种测算方法

本文介绍两种测算方式，主要区别在于分子分母的衡量单位不同：

方法	Xct	Sct	优点	局限
注册资本版（论文方法）	数字产业注册资本（万元）	全部企业注册资本（万元）	反映资本密度，与论文一致	大城市分母稀释效应明显
企业数量版（备选方法）	数字产业企业数量（家）	全部企业数量（家）	不受极值影响，城市间对比更直观	无法区分大企业与小企业的贡献

计算步骤概述

整个计算流程分为以下几个步骤：

读取行业分类：加载《数字经济及其核心产业统计分类（2021）》代码表
构建注销查找表：从注销企业 CSV 中提取 newgcid → exit_year 映射
单年聚合：逐年读取注册企业 CSV，先过滤已注销企业，再按城市聚合
面板累计：跨年累加，得到各城市各年的存量企业指标
计算区位熵：按公式计算 DL，并进行 Winsorize 极端值处理
输出结果：保存为 .dta 文件

数据说明

数据来源

工商注册信息：

1949～2023 年工商企业注册信息数据（含经纬度及其所属的省市区县）(版本2)：https://rstata.duanshu.com/#/brief/course/6d38a3f10cdb467492f3204d1ebdd313

注销企业信息：

1970～2023 年各年各省市区县、各行业注销公司工商信息及数量统计面板数据：https://rstata.duanshu.com/#/brief/course/bcdf21ad0e614645b8449e69342e0851

数字经济核心产业分类：数字经济及其核心产业统计分类.dta，提取自《数字经济及其核心产业统计分类（2021）》。

注销企业处理逻辑

本文采用个体层面过滤的方法处理注销企业：

第 t 年存量企业 = t 年及之前注册的且 t 年及之前未注销的企业

具体实现：预先从注销企业 CSV 提取 newgcid → exit_year 查找表，在每年聚合前关联该表，直接过滤掉 exit_year <= t 的企业，再对存活企业进行聚合。

环境准备

加载依赖包

library(data.table)   # 高性能数据处理
library(haven)        # 读写 Stata .dta 文件
library(future)       # 并行计划设置
library(furrr)        # 多线程 map

以上四个包各司其职：data.table 负责大数据量下的高速读写与聚合；haven 负责输出 Stata 可读的 .dta 文件；future + furrr 组合实现多线程并行，显著加快跨年文件的处理速度。

路径与参数配置

dir_reg  <- file.path(getwd(), "工商注册信息_精简")   # 注册企业 CSV 目录
dir_exit <- file.path(getwd(), "注销企业_精简")        # 注销企业 CSV 目录
dta_file <- file.path(getwd(), "数字经济及其核心产业统计分类.dta")

year_min <- 2010   # 计算起始年份
year_max <- 2012   # 计算截止年份

out_file <- file.path(getwd(), "数字产业集聚度_注册资本_各城市.dta")

# 多线程设置
n_cores   <- suppressWarnings(parallel::detectCores())
if (is.na(n_cores) || n_cores < 1) n_cores <- 4
n_workers <- min(6, max(1, n_cores - 1))

# 放宽 future 全局对象传输上限（默认约 500 MiB，全量注销查找表可能超过）
options(future.globals.maxSize = 4 * 1024^3)   # 4 GB

plan(multisession, workers = n_workers)

其中：

year_min / year_max：控制计算年份范围，便于缩小到示例数据快速验证
n_workers：并行线程数，上限 6，至少保留 1 个核心给主进程
plan(multisession)：启动多进程并行计划，在 Windows/macOS 上均适用

步骤一：读取数字经济产业行业代码

加载分类标准

digi_dta <- read_dta(dta_file) |>
  filter(!is.na(国民经济行业代码) & 国民经济行业代码 != "")

# 提取 3 位前缀（匹配 CSV 行业中类代码后 3 位数字）
digi_prefixes <- unique(substr(digi_dta$国民经济行业代码, 1, 3))
# 构建完整 4 位代码集合（备用于小类匹配）
digi_codes_4d <- unique(digi_dta$国民经济行业代码)

《数字经济及其核心产业统计分类（2021）》使用 4 位行业代码，而工商注册信息 CSV 中的行业代码字段格式为 I641（中类）或 I6411（小类），均带有字母前缀。

代码的匹配策略如下：

行业字段	提取规则	与分类标准对比
`行业小类代码`（4位数字）	`substr(code, 2, 5)`	与 `digi_codes_4d` 完全匹配
`行业中类代码`（3位数字）	`substr(code, 2, 4)`	与 `digi_prefixes` 前缀匹配

优先使用小类代码，兜底使用中类代码，确保最大覆盖率。

步骤二：构建注销查找表

注销企业处理是整个计算中最关键的一步。

预构建查找表

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

代码要点：

select = c(“newgcid”, “退出日期”)：只读两列，大幅节省内存
unique(dt[, .(newgcid, exit_year)], by = “newgcid”)：每个企业只保留最早的注销记录（一次注销即永久退出）
setkey(exit_list, newgcid)：对查找表建索引，与注册数据 merge 时速度接近哈希表查找

步骤三：单年聚合函数

函数设计

这是计算流程的核心函数，逻辑为：读取 → 过滤注销 → 识别数字产业 → 企业缩尾 → 按城市聚合。

注册资本版与企业数量版的唯一区别在于：注册资本版在聚合前需做企业层面 Winsorize，企业数量版（计数）无需缩尾。以下以注册资本版为例说明。

aggregate_reg_year <- function(year_val, dir_reg, digi_prefixes, digi_codes_4d, exit_list) {
  library(data.table)   # 并行子进程中需重新加载包

  fpath <- file.path(dir_reg, paste0(year_val, ".csv"))
  if (!file.exists(fpath)) return(NULL)

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

五个步骤解析：

步骤	核心代码	说明
① 读取	`fread(..., select = ...)`	只读必要列，单年文件可达 200MB+，列选择可节省 60% 内存
② 过滤注销	`dt[is.na(exit_year) \| exit_year > year_val]`	`exit_year` 为 NA 表示未注销；`> year_val` 表示当年末尚未退出
③ 行业识别	`is_digi := ... %in% ...`	逻辑向量，TRUE/FALSE，`sum(is_digi)` 即数字企业计数
④ 企业缩尾	`pmax(pmin(注册资本, p95), p05)`	先于聚合，在企业个体层面截断极端注册资本，防止天价壳公司（如注册资本 1000 亿元的空壳）拉偏整个城市
⑤ 聚合	`sum(注册资本 * is_digi)`	利用逻辑值自动转 0/1，省去 `ifelse`

并行执行

reg_files <- list.files(dir_reg, pattern = "\\.csv$")
reg_years <- sort(as.integer(gsub("\\.csv$", "", reg_files)))
reg_years <- reg_years[reg_years >= year_min & reg_years <= year_max]

reg_results <- future_map(
  reg_years,
  ~ aggregate_reg_year(.x, dir_reg, digi_prefixes, digi_codes_4d, exit_list),
  .options = furrr_options(seed = TRUE),
  .progress = TRUE
)

reg_panel <- rbindlist(reg_results, use.names = TRUE, fill = TRUE)
plan(sequential)   # 并行任务结束后关闭多进程池

future_map 等价于 purrr::map，但自动将任务分发到多个子进程并行执行。.progress = TRUE 会显示进度条。完成后调用 plan(sequential) 释放资源。

步骤四：构建面板与计算累计存量

补全面板并计算累计值

注册数据是流量（某年新注册的企业），而区位熵需要存量（截至该年仍存活的企业总量）。由于已在聚合前过滤了注销企业，这里只需按城市做累计加总即可。

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

关键函数对照：

R 代码	含义
`CJ(a, b)`	data.table 的笛卡尔积，生成所有组合的完整面板
`set(dt, i, j, value)`	就地修改，避免整列复制，内存效率高于 `dt[i, j := value]`
`cumsum(...), by = city_code`	按城市分组计算逐年累积，等价于 Stata 的 `by city: gen cumsum = sum(x)`

步骤五：计算区位熵与极端值处理

计算 DL

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

fifelse 是 data.table 的向量化条件函数，等价于 dplyr::if_else，但速度更快。这里对零分母做安全保护，避免产生 Inf。

Winsorize 处理极端值

本项目采用两层 Winsorize 策略：

企业层面（步骤三已完成）：对每个年份内的企业注册资本做 5%/95% 截断，防止单个天价壳公司拉偏城市存量
城市层面：对最终 DL 指标做 5%/95% 截断，控制区位熵测量结果的极端异常

区位熵存在极端值（如某小城市有大量超高注册资本的数字产业企业），通常使用 5%/95% 分位数截断：

dl_valid <- panel_full[!is.na(DL) & is.finite(DL), DL]
p05 <- quantile(dl_valid, 0.05, na.rm = TRUE)
p95 <- quantile(dl_valid, 0.95, na.rm = TRUE)

panel_full[, DL_raw := DL]   # 保留原始值供对比
panel_full[!is.na(DL), DL := pmax(pmin(DL, p95), p05)]

pmin(DL, p95) 将超过上限的值截断为 p95，pmax(..., p05) 再将低于下限的值截断为 p05，两步合并即为 Winsorize。

步骤六：合并城市名称与输出

从注册文件反查城市名称

city_map <- data.table(city_code = character(), city_name = character())
for (y in sort(reg_years, decreasing = TRUE)) {
  fpath <- file.path(dir_reg, paste0(y, ".csv"))
  if (!file.exists(fpath)) next
  dt <- fread(fpath, select = c("市", "市代码"),
              encoding = "UTF-8", na.strings = "",
              colClasses = list(character = "市代码"))
  setnames(dt, c("市", "市代码"), c("city_name", "city_code"))
  dt <- unique(dt[city_code != "" & city_name != "" & city_name != " "],
               by = "city_code")
  dt <- dt[!city_code %in% city_map$city_code]   # 只添加尚未映射的城市
  if (nrow(dt) > 0) city_map <- rbind(city_map, dt)
}
output <- merge(output, city_map, by = "city_code", all.x = TRUE)

保存为 DTA 文件

haven::write_dta(
  as.data.frame(output),
  out_file,
  label = "数据处理：微信公众号 RStata"
)

haven::write_dta 的 label 参数设置数据集标签，在 Stata 中打开时可通过 notes 或 describe 看到该标注，方便数据溯源。

结果展示（2010–2012 示例）

以下直接读取已生成的结果文件进行展示，两种方法均以 2010–2012 年数据为例。

读取计算结果

library(data.table)
library(ggplot2)
library(hrbrthemes)
library(haven)

# 读取已生成的结果文件
dt_cap <- data.table(haven::read_dta(
  file.path(getwd(), "数字产业集聚度_注册资本_各城市.dta")
))
dt_cnt <- data.table(haven::read_dta(
  file.path(getwd(), "数字产业集聚度_企业数量_各城市.dta")
))

cat(sprintf("注册资本版：%d 行，%d 个城市，年份 %d–%d\n",
    nrow(dt_cap), length(unique(dt_cap$city_code)),
    min(dt_cap$year), max(dt_cap$year)))
cat(sprintf("企业数量版：%d 行，%d 个城市，年份 %d–%d\n",
    nrow(dt_cnt), length(unique(dt_cnt$city_code)),
    min(dt_cnt$year), max(dt_cnt$year)))

描述性统计

cat("========== 注册资本版 描述性统计 ==========\n")
print(dt_cap[, .(
  DL_mean   = round(mean(数字产业集聚度, na.rm = TRUE), 4),
  DL_median = round(median(数字产业集聚度, na.rm = TRUE), 4),
  DL_sd     = round(sd(数字产业集聚度, na.rm = TRUE), 4),
  DL_min    = round(min(数字产业集聚度, na.rm = TRUE), 4),
  DL_max    = round(max(数字产业集聚度, na.rm = TRUE), 4),
  N         = .N
)])

cat("\n========== 企业数量版 描述性统计 ==========\n")
print(dt_cnt[, .(
  DL_mean   = round(mean(数字产业集聚度, na.rm = TRUE), 4),
  DL_median = round(median(数字产业集聚度, na.rm = TRUE), 4),
  DL_sd     = round(sd(数字产业集聚度, na.rm = TRUE), 4),
  DL_min    = round(min(数字产业集聚度, na.rm = TRUE), 4),
  DL_max    = round(max(数字产业集聚度, na.rm = TRUE), 4),
  N         = .N
)])

Top-10 城市（2012 年）

cat("========== 注册资本版 Top-10（2012 年）==========\n")
print(dt_cap[year == 2012][order(-数字产业集聚度)][1:10, .(
  city_name,
  DL         = round(数字产业集聚度, 3),
  注册资本_亿元 = round(存量数字企业注册资本_万元 / 1e8, 2)
)], row.names = FALSE)

cat("\n========== 企业数量版 Top-10（2012 年）==========\n")
print(dt_cnt[year == 2012][order(-数字产业集聚度)][1:10, .(
  city_name,
  DL     = round(数字产业集聚度, 3),
  企业数 = 存量企业数_家
)], row.names = FALSE)

两种方法的 Top-10 城市排名存在差异，主要原因如下：

注册资本版更受大体量企业影响——若某城市少数几家高注册资本的数字企业占全国注册资本比例高，区位熵会被抬高
企业数量版更反映数字产业在城市中的普及度，对城市规模中性

区域对比

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

与论文结果的对比

参考论文（屠西伟、史丹 2025）的结论，数字产业集聚度应呈现 东部 > 中部 > 西部 > 东北 的梯度格局。以 2010–2012 示例数据计算的结果如下：

排名	论文预期	注册资本版	企业数量版
1	东部	东部（0.744）✅	东部（0.919）✅
2	中部	中部（0.653）✅	中部（0.845）✅
3	西部	东北（0.630）⚠️	东北（0.802）⚠️
4	东北	西部（0.562）⚠️	西部（0.749）⚠️

一致之处：东部排第一、中部排第二，与论文完全一致。两种方法均呈现”东高西低”的大格局。

差异之处：东北与西部的排序与论文相反（东北 > 西部，而非论文的西部 > 东北）。可能的原因：

数据年份不同：论文覆盖 2011–2021 年，本文仅以 2010–2012 示例；随着时间推移，西部数字产业可能加速发展、东北相对衰退，长面板下西部可能超过东北
统计口径差异：论文使用 GDP 加权均值（大城市的 DL 和 GDP 权重更大），本文使用简单均值；西部有成都、重庆等高 DL 大城市，加权后西部均值会被显著拉高
企业层面缩尾：本文在企业层面增加了 5%/95% 缩尾处理，论文是否做此处理未知，这可能影响注册资本版的排序
使用的数据不完全一样。

可视化：区域对比

library(ggplot2)

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

可视化：Top-20 城市趋势（企业数量版）

# 取 2012 年 Top-20 城市
top_cities <- dt_cnt[year == 2012][order(-数字产业集聚度)][1:20, city_code]
trend_dt   <- dt_cnt[city_code %in% top_cities]

ggplot(trend_dt, aes(x = year, y = 数字产业集聚度, color = city_name)) +
  geom_line(linewidth = 1.1, alpha = 0.85) +
  geom_point(size = 1.8, alpha = 0.85) +
  scale_x_continuous(breaks = 2010:2012) +
  scale_color_discrete(name = "城市") +
  labs(
    title    = "数字产业集聚度 Top-20 城市趋势（企业数量版，2010–2012）",
    subtitle = "数据来源：微信公众号 RStata",
    x        = NULL,
    y        = NULL
  ) +
  theme(
    legend.position  = "right",
    panel.grid.minor = element_blank()
  )

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算各城市数字产业集聚程度

指标来源与计算原理

数字产业集聚度（Location Quotient）

区位熵的经济含义

两种测算方法

计算步骤概述

数据说明

数据来源

注销企业处理逻辑

环境准备

加载依赖包

路径与参数配置

步骤一：读取数字经济产业行业代码

加载分类标准

步骤二：构建注销查找表

预构建查找表

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

步骤三：单年聚合函数

函数设计

并行执行

步骤四：构建面板与计算累计存量

补全面板并计算累计值

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

步骤五：计算区位熵与极端值处理

计算 DL

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

Winsorize 处理极端值

步骤六：合并城市名称与输出

从注册文件反查城市名称

保存为 DTA 文件

结果展示（2010–2012 示例）

读取计算结果

描述性统计

Top-10 城市（2012 年）

区域对比

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

与论文结果的对比

可视化：区域对比

可视化：Top-20 城市趋势（企业数量版）

相关文章

评论