名师讲堂｜使用 R 语言测算各城市虚拟集聚程度

一、虚拟集聚度指标介绍

1.1 指标背景

虚拟集聚度（Virtual Agglomeration）是衡量城市在数字经济领域集聚程度的重要指标。与传统的地理集聚不同，虚拟集聚强调通过信息技术实现的空间联系和资源共享，反映了城市在信息传输、软件和信息技术服务业领域的相对优势。

该指标的计算方法参考了《虚拟集聚与城市经济韧性》（宋林等）一文，通过结合区位熵和空间距离权重来综合测度城市的虚拟集聚水平。

1.2 计算公式

虚拟集聚度的计算公式为：

1.3 计算步骤概述

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

1. 数据准备：读取新增企业和注销企业数据
2. 存续企业计算：通过累计新增减去累计注销得到各城市各行业的存续企业数
3. IT 行业筛选：提取”信息传输、软件和信息技术服务业”企业数据
4. 区位熵计算：计算各城市 IT 行业的区位熵
5. 距离矩阵计算：使用 sf 包计算城市间的地理距离
6. 虚拟集聚度计算：结合区位熵和距离权重计算最终指标
7. 可视化：绘制空间分布地图和趋势折线图

---

二、详细计算代码

2.1 环境准备与数据读取

首先加载必要的 R 包并读取原始数据：

# 加载必要的包
library(tidyverse)  # 数据处理
library(sf)         # 空间数据处理

# 设置工作目录
setwd("/Users/ac/Desktop/使用 R 语言测算各城市虚拟集聚程度")

# 读取新增企业数据
new_firms <- read_csv("新增企业数量.csv")
message(sprintf("新增数据: %d 行", nrow(new_firms)))

# 读取注销企业数据
closed_firms <- read_csv("注销公司数量.csv")
message(sprintf("注销数据: %d 行", nrow(closed_firms)))

数据说明：

• 新增企业数量.csv：包含各城市各行业每年的新增注册企业数量
• 注销公司数量.csv：包含各城市各行业每年的注销企业数量
• 数据时间范围：2010-2023年
• 行业分类：按照国民经济行业分类标准

上面代码中读取的两个 csv 文件分别来源于平台上的工商注册信息和注销信息：

use "1949～2023年各省市区县、行业新增企业数量统计.dta", clear
export delimited using "新增企业数量.csv", replace

use "1970～2023年各年各省市区县、各行业注销公司数量面板数据.dta", clear
export delimited using "注销公司数量.csv", replace

2.2 数据预处理——计算累计存续企业

计算每个城市各行业的累计存续企业数，需要考虑企业的进入（新增）和退出（注销）：

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

关键函数说明：

• cumsum()：计算累计和
• arrange()：按指定变量排序
• group_by() + ungroup()：分组计算后解除分组

2.3 合并计算存续企业数

将新增和注销数据合并，计算最终的存续企业数量：

# 合并计算存续企业
firms <- new_summary %>%
  full_join(
    closed_summary %>% select(年份, 省份, 城市, 行业门类, 累计注销),
    by = c("年份", "省份", "城市", "行业门类")
  ) %>%
  mutate(
    累计新增 = coalesce(累计新增, 0),
    累计注销 = coalesce(累计注销, 0),
    存续数量 = pmax(累计新增 - 累计注销, 0)
  )

关键函数说明：

• coalesce()：返回第一个非缺失值
• pmax()：并行取最大值，确保存续数量不为负数

2.4 筛选 IT 行业并汇总

提取”信息传输、软件和信息技术服务业”数据，并按城市-年份汇总：

# 定义 IT 行业名称
it_sector <- "信息传输、软件和信息技术服务业"

# 按城市-年份汇总
city_data <- firms %>%
  group_by(年份, 城市) %>%
  summarise(
    总企业数 = sum(存续数量, na.rm = TRUE),
    IT企业数 = sum(ifelse(行业门类 == it_sector, 存续数量, 0), na.rm = TRUE),
    .groups = "drop"
  ) %>%
  filter(总企业数 > 0)

2.5 计算区位熵

区位熵（Location Quotient）衡量某城市 IT 行业的专业化程度相对于全国的水平：

# 计算区位熵
city_data <- city_data %>%
  group_by(年份) %>%
  mutate(
    全国IT = sum(IT企业数),
    全国总 = sum(总企业数)
  ) %>%
  ungroup() %>%
  mutate(
    全国IT集聚度 = 全国IT / 全国总,
    城市IT集聚度 = IT企业数 / 总企业数,
    区位熵 = 城市IT集聚度 / 全国IT集聚度,
    区位熵 = ifelse(is.finite(区位熵), 区位熵, 0)
  )

区位熵解读：

• 区位熵 > 1：该城市 IT 行业集聚度高于全国平均水平
• 区位熵 = 1：与全国平均水平相当
• 区位熵 < 1：低于全国平均水平

2.6 计算城市间距离矩阵

使用 sf 包读取行政区划数据并计算城市质心间的距离：

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

2.7 计算虚拟集聚度

结合区位熵和距离权重，计算每个城市的虚拟集聚度：

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

2.8 保存结果

# 合并结果并保存
final_result <- city_data %>%
  left_join(vag_results, by = c("年份", "城市"))

final_result %>%
  write_csv("城市虚拟集聚度结果.csv")

---

三、数据可视化

这部分可以学习平台上的地图绘制相关课程：

使用 R 语言绘制地图课程汇总索引: https://mp.weixin.qq.com/s/lYyVFUzKWkzyBFRAInoZFg

3.1 空间分布地图绘制

library(ggspatial)
library(scico)

# 定义Albers投影坐标系
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 读取数据
vag_data <- read_csv("城市虚拟集聚度结果.csv", show_col_types = FALSE)
vag_2023 <- vag_data %>%
  filter(年份 == 2023) %>%
  select(城市, 虚拟集聚度, 区位熵, IT企业数, 总企业数)

# 设置绘图区域
st_bbox(c(xmin = -2725586, xmax = 2982768, ymax = 6000000, ymin = 1800655),
        crs = st_crs(mycrs)) %>% st_as_sfc() -> plotbbox

# 读取地图数据
citymap <- read_sf("mapdata/chinacity2021mini/chinacity2021mini.shp") %>%
  filter(!is.na(省代码))

citylinemap <- read_sf("mapdata/chinacity2021mini/chinacity2021mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class)

china_neighboring <- read_sf("mapdata/china_neighboring/china_neighboring.shp")

数据分组处理

# 合并数据
citymap %>%
  left_join(vag_2023, by = c("市" = "城市")) %>%
  mutate(虚拟集聚度 = if_else(is.na(虚拟集聚度), 0, 虚拟集聚度)) -> citymap_data

# 将连续变量分组为离散变量
vag_values <- citymap_data$虚拟集聚度[citymap_data$虚拟集聚度 > 0]
breaks <- quantile(vag_values, probs = c(0, 0.2, 0.4, 0.6, 0.8, 1), na.rm = TRUE)

citymap_data <- citymap_data %>%
  mutate(vag_group = case_when(
    虚拟集聚度 == 0 ~ "无数据",
    虚拟集聚度 <= breaks[2] ~ sprintf("%.2f-%.2f", breaks[1], breaks[2]),
    虚拟集聚度 <= breaks[3] ~ sprintf("%.2f-%.2f", breaks[2], breaks[3]),
    虚拟集聚度 <= breaks[4] ~ sprintf("%.2f-%.2f", breaks[3], breaks[4]),
    虚拟集聚度 <= breaks[5] ~ sprintf("%.2f-%.2f", breaks[4], breaks[5]),
    TRUE ~ sprintf(">%.2f", breaks[5])
  ))

citymap_data$vag_group <- factor(citymap_data$vag_group,
                                  levels = c("无数据",
                                            sprintf("%.2f-%.2f", breaks[1], breaks[2]),
                                            sprintf("%.2f-%.2f", breaks[2], breaks[3]),
                                            sprintf("%.2f-%.2f", breaks[3], breaks[4]),
                                            sprintf("%.2f-%.2f", breaks[4], breaks[5]),
                                            sprintf(">%.2f", breaks[5])))

citymap_proj <- citymap_data %>% st_transform(mycrs)
china_neighboring_proj <- china_neighboring %>% st_transform(mycrs)

绘制地图

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

3.2 历年趋势折线图

library(gghighlight)

# 筛选数据
此处代码需下载讲义材料查看～

# 绘制折线图
p_line <- ggplot(vag_trend, aes(x = 年份, y = 虚拟集聚度, group = 城市, color = 城市)) +
  geom_line(linewidth = 0.8, alpha = 0.7) +
  gghighlight(城市 %in% municipalities, label_key = 城市,
              label_params = list(family = cnfont, size = 4, nudge_x = 0.5, nudge_y = 0.1,
                                  segment.size = 0.3, segment.color = "gray50"),
              unhighlighted_params = list(color = "gray80", alpha = 0.3, linewidth = 0.4),
              use_direct_label = TRUE) +
  scale_color_manual(values = municipality_colors) +
  scale_x_continuous(breaks = c(seq(1990, 2020, 5), 2023), limits = c(1990, 2023)) +
  scale_y_continuous(breaks = seq(0, 4, by = 0.5), limits = c(0, 4)) +
  labs(title = "1990-2023年中国城市虚拟集聚度变化趋势",
       subtitle = "突出显示四个直辖市 | 其他城市以灰色显示",
       x = NULL, y = "虚拟集聚度",
       caption = "数据来源：全国工商企业注册数据 | 绘制：RStata") +
  hrbrthemes::theme_ipsum(base_family = cnfont, grid = "Y") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5, size = 11, color = "gray40"),
        plot.caption = element_text(hjust = 1, size = 9, color = "gray50"),
        legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank()) -> p2

ggsave("charts/09_城市虚拟集聚度历年趋势_直辖市高亮.png", plot = p_line, width = 12, height = 7, dpi = 300, device = png)

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