今天给大家分享使用 R 语言测算地区产业专业化指标的方法。该方法参考自杨本建、唐金汶（2022）《数字经济与区域产业布局》中的公式 (2)，其思想源于 Kalemli-Ozcan et al. (2003) 与 Du et al. (2022) 的 Krugman 式专业化指数。
附件中提供了该参考文献的 PDF 文件（数字经济与区域产业布局.pdf），感兴趣的小伙伴可以阅读原文。
本文使用的原始数据为工商企业注册信息（已在本项目内裁剪为测算所需的 9 个变量），以 2000–2005 年为例演示完整计算过程。
指标来源与计算过程
地区产业专业化指数（公式）
行业范围：制造业 C13–C43，剔除 C39
按国民经济行业分类（GB/T 4754）的制造业门类（行业门类 = “制造业”），取 2 位行业大类代码 C13–C43，并剔除 C39（计算机、通信和其他电子设备制造业）。剔除 C39 的依据是论文第 248 页附录 1 的说明（该行业受数字经济影响特殊，在相关研究中通常单独处理）。
产业规模的代理变量
论文以企业的注册资本作为行业规模的代理变量（主指标）；同时以企业数量作为稳健性口径。两者计算逻辑完全一致，仅在汇总时替换聚合字段。
存续企业的界定（进入与退出）
参考李磊等（2023）的做法，按”进入—退出”口径统计每年各城市各行业的存续企业：
- 进入：以企业的成立年份作为进入年份；
- 退出：当经营状态含”注销/吊销”时视为退出企业，退出年份取核准日期的年份；
- 存续：在年份 t 满足”成立年份 ≤ t，且（未退出 或 退出年份 > t）”的企业。
计算步骤概述
整体计算分为以下几个步骤：
- 读取与清洗：读取各年工商注册数据，仅保留制造业、有效城市代码、必需变量，并界定进入/退出年份；
- 构建存续面板：对每个目标年份 t，筛选存活企业，按”城市 × 行业”汇总规模，得到城市×行业×年份存续企业规模面板；
- 计算专业化指数：对每一年，按公式计算各城市的 spec 指数（注册资本口径 + 企业数量口径）；
- 输出与可视化：导出结果 CSV，并绘制趋势、分布与 Top 城市等图表。
详细讲解计算代码
下面按步骤完整展示 R 代码，每一段均可直接运行。
0. 路径与参数
首先加载 tidyverse，设定工程目录、输出目录、数据目录，以及制造业行业代码与所需变量。
library(tidyverse)
proj_dir <- dirname(knitr::current_input(dir = TRUE))
out_dir <- file.path(proj_dir, "输出")
data_dir <- file.path(proj_dir, "工商注册信息2025-sample")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
file_years <- 2000:2005
target_years <- 2000:2005
mfg_codes <- setdiff(sprintf("C%02d", 13:43), "C39")
mfg_codes
need_cols <- c("注册资本", "实缴资本", "行业门类", "行业大类代码",
"经营状态", "成立年份", "核准日期", "市", "市代码")
|
1. 读取并清洗单个年份文件
read_one_year() 负责把一年的 CSV 读入并清洗为”企业级”明细：
- 仅保留行业门类 == “制造业”且大类在 C13–C43 且非 C39 的记录；
- 仅保留市代码为 6 位数字行政区划码的有效城市；
- 把成立年份、核准日期解析为年份，按经营状态判定是否退出企业并得到退出年份；
- 剔除成立年份或注册资本缺失、注册资本非正的样本，以及”退出早于成立”的逻辑异常样本。
此处代码需要下载讲义材料查看～
将企业级明细缓存为 RDS，便于后续复算或单独调试：
# 缓存企业级明细，便于复算
write_rds(firms, file.path(out_dir, "firms_manufacturing_2000_2005.rds"), compress = "gz")
|
2. 构建城市×行业×年份存续企业规模面板
build_year_scale(t) 对目标年份 t 筛选”存续企业”（成立年份 ≤ t，且未退出或退出年份 > t），按城市 × 行业汇总该年的注册资本总额（output_reg）与企业数（n_firm）。逐年份构建后即得到城市×行业×年份规模面板。
message("步骤 2/4：构建城市×行业×年份存续企业规模面板 ……")
build_year_scale <- function(t) {
firms |>
filter(entry_year <= t, is.na(exit_year) | exit_year > t) |>
group_by(city_code, city, industry) |>
summarise(
year = t,
output_reg = sum(cap_reg, na.rm = TRUE),
n_firm = n(),
.groups = "drop"
)
}
city_ind_year <- map(target_years, build_year_scale) |> list_rbind() |>
relocate(year)
write_csv(city_ind_year,
file.path(out_dir, "城市_行业_年份_制造业规模.csv"))
city_ind_year
|
3. 计算 Krugman 式地区产业专业化指数
compute_spec_one_year(df, value_col) 是核心：对某一年的城市×行业规模数据，计算各城市的专业化指数。
关键点：
- 计算每个城市该口径的总规模 city_total，并求各行业份额 share；
- 通过 crossing() 补全”城市 × 行业”全网格，该城市该行业无企业时份额记为 0；
- 对每个行业求全部城市份额之和 sum_share，则”其他城市平均份额”为 (sum_share - share)/(J-1)；
- 计算差的平方 (share - other_avg)^2，跨行业求和即得该城市的 spec。
message("步骤 3/4：计算地区产业专业化指数 ……")
spec_all <- spec_reg |>
left_join(select(spec_cnt, year, city_code, spec_count),
by = c("year", "city_code")) |>
select(year, city_code, city, n_city, spec_reg, spec_count) |>
arrange(year, desc(spec_reg))
write_csv(spec_all,
file.path(out_dir, "地区产业专业化指数_2000_2005.csv"))
message("步骤 4/4：输出结果 ……")
message(sprintf(" 各年城市数(J)：\n%s",
paste(capture.output(print(count(spec_all, year, n_city = n_city))),
collapse = "\n")))
message("\n2005 年产业专业化指数最高的 10 个城市：")
spec_all |>
filter(year == 2005) |>
slice_max(spec_reg, n = 10) |>
select(city, spec_reg, spec_count) |>
print(n = 10)
message("\n完成！结果已写入：", out_dir)
|
结果预览
读取结果，展示 2005 年专业化指数最高的若干城市，便于核对：
spec_all <- read_csv(file.path(out_dir, "地区产业专业化指数_2000_2005.csv"),
show_col_types = FALSE)
spec_all |>
filter(year == 2005) |>
slice_max(spec_reg, n = 15) |>
select(year, city, n_city, spec_reg, spec_count) |>
knitr::kable(digits = 4, caption = "2005 年产业专业化指数最高的 15 个城市")
|
数据可视化
下面使用 ggplot2 绘制三张图表，直观展示专业化指数的整体趋势、年份分布与头部城市。
图1：各城市平均专业化指数随年份变化趋势
计算全部城市在各年份的平均值与中位数，并叠加”平均值 ± 1 倍标准差”阴影带，刻画整体趋势。
# 计算全部城市在各年份的平均值与中位数，刻画整体趋势
# 此处代码需要下载讲义材料查看～
|
图2：各年专业化指数分布（箱线图）
p_box <- spec_all |>
mutate(year = factor(year)) |>
ggplot(aes(year, spec_reg, fill = year)) +
geom_boxplot(alpha = 0.8, outlier.alpha = 0.3, show.legend = FALSE) +
scale_fill_manual(values = c("#0055aa", "#c40003", "#00c19b", "#eac862",
"#7fd2ff", "#007ed3")) +
labs(title = "各年份城市产业专业化指数分布 (2000-2005)",
subtitle = "箱线图展示全部城市的专业化指数分布及其演变",
x = NULL, y = "产业专业化指数 spec",
caption = "数据处理 & 绘图：微信公众号 RStata") +
theme(plot.title = element_text(face = "bold"))
ggsave(file.path(out_dir, "图2_专业化指数分布.png"), device = png,
p_box, width = 10, height = 6, dpi = 300, bg = "white")
|
图3：2005 年专业化程度最高的 20 城市
p_top <- spec_all |>
filter(year == 2005) |>
slice_max(spec_reg, n = 20) |>
mutate(city = fct_reorder(city, spec_reg)) |>
ggplot(aes(spec_reg, city)) +
geom_col(fill = "#7fd2ff", alpha = 0.85) +
labs(title = "2005 年制造业产业专业化程度最高的 20 个城市",
subtitle = "数据处理 & 绘图：微信公众号 RStata",
x = "产业专业化指数 spec", y = NULL) +
theme(plot.title = element_text(face = "bold"))
ggsave(file.path(out_dir, "图3_2005专业化top20.png"), device = png,
p_top, width = 10, height = 7, dpi = 300, bg = "white")
|
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算地区产业专业化指标
评论