附件中的文献《生产性服务业集聚何以赋能技术扩散》中提到了城市的生产性服务业专业化集聚指标和多样化集聚指标，今天的课程中我们将会讲解如何使用 R 语言根据城市统计年鉴数据来测算这两个指标。
城市统计年鉴数据在这里：
1999～2025 年中国城市统计年鉴面板数据整理结果: https://rstata.duanshu.com/#/brief/course/c6b0aaf3bcba494faadbacacb012dcef
这里选择的是 2017～2019 年的样本作为演示使用。也就是附件中的 城市统计年鉴样本数据.dta。
一、指标来源与计算公式
1.1 专业化集聚指数 SP（公式4）
1.2 多样化集聚指数 DV（公式5）
二、数据与变量说明
2.1 数据来源
城市统计年鉴数据在这里：
1999～2025 年中国城市统计年鉴面板数据整理结果: https://rstata.duanshu.com/#/brief/course/c6b0aaf3bcba494faadbacacb012dcef
这里选择的是 2017～2019 年的样本作为演示使用。也就是附件中的 城市统计年鉴样本数据.dta。
由于计算所需的变量主要集中在 2003～2019 年，所以完整数据仅可计算得到 2003～2019 年的。
2.2 生产性服务业行业界定
参考 Ke et al.（2014）和韩峰和阳立高（2020），生产性服务业包括 7个细分行业：
| 代码 |
行业名称 |
说明 |
| s1 |
交通运输、仓储和邮政业 |
全市口径 |
| s2 |
金融业 |
全市口径 |
| s3 |
科学研究和技术服务业 |
2003-2016含地质勘探，2017-2019不含 |
| s4 |
租赁和商务服务业 |
全市口径 |
| s5 |
信息传输、计算机服务和软件业 |
全市口径 |
| s6 |
批发和零售业 |
全市口径 |
| s7 |
水利、环境和公共设施管理业 |
全市口径 |
注意：科学研究和技术服务业在2003-2016年统计口径含”地质勘探业”，2017-2019年不含。处理时用 coalesce() 合并两套口径。
2.3 就业数据说明
- 数据来自《中国城市统计年鉴》城镇单位从业人员（全市口径）
- 单位：人
- 总就业 = 城镇单位从业人员期末人数（全市）
三、R 语言代码实现
3.1 环境准备与数据读取
library(tidyverse)
library(haven)
df_raw <- read_dta("城市统计年鉴样本数据.dta")
|
3.2 变量选取与数据清洗
# ============================================================
# 二、选取变量并清洗
# 生产性服务业包括 7 个细分行业
# 就业数据来自《中国城市统计年鉴》城镇单位从业人员（全市口径）
# 数据处理：微信公众号 RStata
# ============================================================
# 此处代码需下载讲义材料查看～
|
3.3 异常值处理
jyg_fix <- df %>%
filter(city == "嘉峪关市", year %in% c(2012, 2014)) %>%
summarise(emp_total_fix = mean(emp_total, na.rm = TRUE)) %>%
pull(emp_total_fix)
cat("嘉峪关市 2013 年 emp_total 修复：\n")
cat(" 原始值:", df %>% filter(city == "嘉峪关市", year == 2013) %>% pull(emp_total), "\n")
cat(" 替换值（2012/2014均值）:", round(jyg_fix), "\n")
df <- df %>%
mutate(emp_total = if_else(city == "嘉峪关市" & year == 2013,
jyg_fix, emp_total))
cat("\n剔除鹤壁市 2013 年（科研技术异常值）:\n")
n_before <- nrow(df)
df <- df %>% filter(!(city == "鹤壁市" & year == 2013))
cat(" 剔除前:", n_before, "条 → 剔除后:", nrow(df), "条\n")
cat("\n异常值处理后样本量:", nrow(df), "\n")
|
3.4 计算全国层面基准（按年份）
# ============================================================
# 三、计算全国层面基准（按年份）
# 数据处理：微信公众号 RStata
# ============================================================
# 此处代码需下载讲义材料查看～
cat("全国层面基准计算完成，年份数:", nrow(nat_by_year), "\n")
|
3.5 测算专业化集聚指数 SP（公式4）
# ============================================================
# 四、测算专业化集聚指数 SP（公式4）
# SP_i = ∑_s(E_is/E_i) / ∑_s(E_s/E)
# 数据处理：微信公众号 RStata
# ============================================================
# 此处代码需下载讲义材料查看～
|
3.6 测算多样化集聚指数 DV（公式5）
# ============================================================
# 五、测算多样化集聚指数 DV（公式5）
# 关键：公式中的份额是"其余行业在剩余就业中的份额"
# DV_i = ∑_s (E_{i,s}/E_i) × [∑_{s'≠s}[E_{s'}/(E-E_s)]^2 / ∑_{s'≠s}[E_{i,s'}/(E_i-E_{i,s})]^2]
# 数据处理：微信公众号 RStata
# ============================================================
# 此处代码需下载讲义材料查看～
# 城市与全国就业宽表合并
# 此处代码需下载讲义材料查看～
|
3.7 缩尾处理与结果保存
cat(" 处理后范围:", round(min(df_result$`生产性服务业多样化集聚DV`), 4),
"～", round(max(df_result$`生产性服务业多样化集聚DV`), 4), "\n")
cat("\n=== 最终结果数据概览 ===\n")
cat("总行数:", nrow(df_result), "\n")
cat("年份范围:", min(df_result$year), "～", max(df_result$year), "\n")
cat("\nSP（专业化集聚）描述性统计:\n")
print(summary(df_result$SP))
cat("\nDV（多样化集聚，经1% Winsorize）描述性统计:\n")
print(summary(df_result$`生产性服务业多样化集聚DV`))
cat("\n样本展示（前10行）:\n")
print(
df_result %>%
select(year, prov, city, SP, `生产性服务业多样化集聚DV`) %>%
head(10),
digits = 4
)
|
year_min <- min(df_result$year)
year_max <- max(df_result$year)
output_filename <- paste0(year_min, "～", year_max, "年各城市生产性服务业集聚水平.dta")
df_out <- df_result %>%
rename(
`年份` = year,
`省` = prov,
`省代码` = prov_code,
`市` = city,
`市代码` = cityid,
`生产性服务业专业化集聚SP` = SP
)
write_dta(
df_out,
path = output_filename,
label = "数据处理：微信公众号 RStata"
)
cat("\n已保存:", output_filename, "\n")
|
四、结果验证：与论文对比
4.1 SP 指标验证
论文报告 SP 均值约为 0.83。
sp_mean <- mean(df_result$SP, na.rm = TRUE)
cat("SP 均值:", round(sp_mean, 4), "\n")
cat("论文报告 SP 均值: 0.83\n")
cat("差距:", round(sp_mean - 0.83, 4), "\n")
|
结论：SP 均值 r round(sp_mean, 4) 与论文报告的 0.83 高度吻合 ✅
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算各城市生产性服务业集聚水平
评论