名师讲堂｜使用 R 语言测算城市数字经济发展的工具变量

3015 字 约 9 分钟 名师讲堂 区域经济R 语言

今天给大家分享如何使用 R 语言 + sf 包，基于上市公司数据和行政区划边界，测算城市数字经济发展的地理工具变量（IV）。工具变量定义参考杨本建、唐金汶《数字经济与区域产业布局》（2026 年第 3 期），其采用「各城市到杭州的球面距离 × 全国与杭州上市公司数字化转型程度之比」作为城市数字经济发展的地理工具变量，以缓解内生性。

基于早期杭州在数字经济发展方面的引领地位，计算各城市到杭州的球面距离。本文将这一距离乘以全国上市公司与位于杭州的上市公司数字化转型程度之比来反映杭州数字经济辐射作用随时间的变化，并以这一变量作为城市数字经济发展的工具变量（地理工具变量）。其中，上市公司数字化转型程度参考吴非等（2021）的方法计算而得。

处理流程概览：

步骤 脚本位置 任务
Step 0 第 0 节 sf::st_read() 读取 2021行政区划/市.shp，用 sf::st_centroid() 生成各市域质心坐标
Step 1 第 1 节 haven::read_dta() 读入上市公司数字化转型关键词总词频（公司-年度）
Step 2 第 2 节 haven::read_dta() 读入注册/办公地址文件，提取「注册地址所在市」「办公地址所在市」
Step 3 第 3 节 dplyr::group_by() + summarise() 计算年度全国均值、杭州（两口径）均值，并求「全国/杭州」之比
Step 4 第 4 节 sf::st_distance() 计算各城市到杭州的 WGS84 测地距离（km）
Step 5 第 5 节 dplyr::cross_join() 构造注册地址 / 办公地址两套城市-年度面板与工具变量
Step 6 + (B) 第 6 节及标签段 attr() 添加中文变量标签，haven::write_dta() 保存为 Stata 14 格式

数据来源说明：测算依赖三类数据——上市公司数字化转型关键词总词频（吴非 2021，公司-年度）、上市公司注册与办公地址（含所在地级市，公司-年度）、2021 年行政区划市域边界（shapefile，用于生成质心）。其中上市公司数据与地址数据来自 RStata 平台分享的「2001～2024 年各上市公司数字转型关键词总词频」与「2000～2024 年上市公司注册地址与办公地址」数据集；行政区划来自「2021 行政区划」数据。

上市公司数字化转型关键词总词频（吴非 2021，公司-年度）数据来自之前分享的：

2001～2024 年各上市公司数字转型关键词总词频及TF-IDF指标计算结果（吴非2021版本，基于管理层讨论分析文本）: https://rstata.duanshu.com/#/brief/course/91b932a04b774f2785e4159ed56a3cd4

上市公司数据与地址数据来自之前分享的：

2000～2024 年上市公司注册地址与办公地址（含经纬度及其所处的省市区县）：https://rstata.duanshu.com/#/brief/course/746d4f595eba41a2a116b5a6edc3ac0a

一、指标来源与计算原理

1.1 核心指标：地理工具变量

变量 含义
dist_hz_km 各城市到杭州的球面距离（km，sf::st_distance() 在 EPSG:4326 下返回 WGS84 测地距离）
ratio_natl_hz 全国 / 杭州 上市公司数字化转型程度之比（随年份变化、对所有城市相同）
iv_product 地理工具变量 = dist_hz_km × ratio_natl_hz

两种“公司所在地”口径：为便于稳健性检验，本结果同时按注册地址所在市（注册地址_市）与办公地址所在市（办公地址_市）各算一套。两者仅在「哪些城市纳入、以及杭州均值取注册地还是办公地在杭州的公司」上不同；全国均值与球面距离口径一致。

1.2 数字化转型程度的测算原理

“上市公司数字化转型程度”由 吴非等（2021） 的方法判定：从上市公司年报的「管理层讨论与分析（MD&A）」文本中，提取数字转型相关关键词并统计其总词频。本讲义取公司-年度「总词频」的年度均值作为该群体的数字化转型程度：

  • 全国：年份 t 全部上市公司的总词频均值（两口径相同）；
  • 杭州（注册口径）：年份 t 注册地位于杭州市的上市公司总词频均值；
  • 杭州（办公口径）：年份 t 办公地位于杭州市的上市公司总词频均值。

注意：原始数字化文件仅含总词频 ≥ 1 的公司-年度（无 0 值记录），且 99.8% 为 A 股；比值在全国与杭州间采用一致口径，时间维度上的变化有效。

1.3 球面距离的测算原理

sf::st_distance() 在地理坐标系（EPSG:4326，WGS84）下直接计算两点间测地距离（geodesic），单位米，比手写 haversine 更精确，且与 Stata 的 geodist 命令等价（geodist 使用 WGS84 椭球，sf 默认经 s2 球面实现，两者差异极小，可忽略）。

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

其中 centroid_lat / centroid_lon 为各城市市域质心经纬度，hz_pt 为杭州市质心点。

1.4 数据来源与口径

数据 角色
2001～2024年各上市公司数字转型关键词总词频（吴非2021，基于管理层讨论与分析部分文本）.dta 公司-年度数字化转型程度（总词频）
2000～2024年上市公司注册地址与办公地址（含经纬度、所处的省市区县及搬迁距离）.dta 公司-年度所在地级市（注册地址_市 / 办公地址_市）
2021行政区划/市.shp 各地级市行政边界，用于求市域质心坐标

二、数据准备工作

2.1 Step 0：生成城市质心坐标

关键函数说明：

  • sf::st_read()：读取 shapefile，得到 sf 对象（含 geometry 列）。
  • sf::st_centroid()：直接计算每个多边形的面积加权质心（自动处理含岛屿/飞地的多部件多边形），与 Stata shp2dta 的 gencentroid() 等价。返回的 sfc_POINT 用 st_coordinates() 提取经纬度（X=经度，Y=纬度）。
library(tidyverse)
library(sf)
library(haven)

# 路径设置（与脚本同目录）
data_dir <- "."
out_dir  <- "输出"
if (!dir.exists(out_dir)) dir.create(out_dir)

shp_file <- file.path(data_dir, "2021行政区划/市.shp")

# 读取 shapefile，保留省/市名称，并计算面积加权质心
city_sf <- st_read(shp_file, quiet = TRUE) %>%
  rename(city = 市, province = 省) %>%
  select(province, city)

cent <- st_centroid(st_geometry(city_sf))      # 面积加权质心
coords <- st_coordinates(cent)                 # X=经度, Y=纬度

centroid_df <- city_sf %>%
  mutate(
    centroid_lon = coords[, "X"],
    centroid_lat = coords[, "Y"]
  ) %>%
  st_drop_geometry() %>%
  select(province, city, centroid_lat, centroid_lon)

head(centroid_df)

2.2 Step 1–2：读取数字化转型文件与地址文件

主程序先用 haven::read_dta() 读入两张公司-年度表，统一把中文变量名改为英文（year / stkcd / digi / city_reg / city_off），规避后续合并与计算中的字符变量名问题。

# 数字化转型文件（公司-年度）
dt <- read_dta(file.path(data_dir, "2001～2024年各上市公司数字转型关键词总词频（吴非2021，基于管理层讨论与分析部分文本）.dta")) %>%
  rename(year = 年份, stkcd = 股票代码, digi = 总词频) %>%
  select(year, stkcd, digi)

# 地址文件（公司-年度）
addr <- read_dta(file.path(data_dir, "2000～2024年上市公司注册地址与办公地址（含经纬度、所处的省市区县及搬迁距离）.dta")) %>%
  rename(year = 年份, stkcd = 股票代码,
         city_reg = 注册地址_市, city_off = 办公地址_市) %>%
  select(year, stkcd, city_reg, city_off)

说明：select() 仅保留后续需要的变量，减小内存占用；两表通过 stkcd + year 唯一对应。

三、主程序：年度测算与面板构造

以下代码对应 计算地理工具变量.R 第 3–5 节。

3.1 年度数字化转型程度：全国均值与杭州均值（两口径）

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

说明：group_by(year) %>% summarise(mean(...)) 对应 Stata collapse (mean)；注册/办公两口径通过 filter(city_? == "杭州市") 分别筛选；full_join 把三组年度指标横向合并；filter(...) 自动剔除杭州双口径都无样本的年份（2006 年原始数据无样本），保证后续面板年份一致。

3.2 城市到杭州的距离（sf::st_distance）

# 构造 sf 点对象（经纬度顺序：先经度，后纬度）
pts <- centroid_df %>%
  st_as_sf(coords = c("centroid_lon", "centroid_lat"), crs = 4326)

# 杭州市质心点
hz_pt <- pts %>% filter(city == "杭州市")

# 计算各城市到杭州的测地距离（米），再换算为千米
dmat <- st_distance(pts, hz_pt)

cent_df <- centroid_df %>%
  mutate(dist_hz_km = as.numeric(dmat) / 1000)

说明：先构造 sf 点对象并赋 EPSG:4326；st_distance() 在地理坐标系下返回测地距离（米），除以 1000 得到千米。as.numeric() 会剥离 units 单位并返回米值。

3.3 构造两种口径的城市-年度面板

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

说明：cross_join() 把「城市 × 距离」与「年份 × 比值」做笛卡尔积，得到城市-年度面板；mutate(iv_product = ...) 即地理工具变量；两口径仅比值变量不同（注册用 ratio_natl_hz_reg，办公用 ratio_natl_hz_off），其余完全一致。随后第 6 节保存年度比值与质心文件。

四、变量标签与结果可视化

4.1 添加中文变量标签并保存

tidyverse + sf 的结果是 R 数据框，中文变量标签可通过 attr(df$col, "label") <- "中文含义" 设置。haven::write_dta() 在保存时会保留这些标签，供 Stata 用户 describe 查看。

# 小工具：批量添加变量标签
add_var_labels <- function(.data, ...) {
  lbls <- list(...)
  for (nm in names(lbls)) {
    if (nm %in% names(.data)) attr(.data[[nm]], "label") <- lbls[[nm]]
  }
  .data
}

# 构造输出数据框
yearly_out <- yearly %>%
  select(year, digi_natl, digi_hz_reg, digi_hz_off,
         ratio_natl_hz_reg, ratio_natl_hz_off)

cent_out <- cent_df %>%
  select(province, city, centroid_lat, centroid_lon, dist_hz_km)

# 给各表添加变量标签
# 此处代码需下载讲义材料查看～

# 保存为 Stata 14 格式（UTF-8 中文标签/中文城市名）
out_dir_dta <- "城市数字经济发展地理工具变量_结果"
if (!dir.exists(out_dir_dta)) dir.create(out_dir_dta)

write_dta(panel_reg, file.path(out_dir_dta, "城市数字经济发展地理工具变量_注册地址_2001-2024.dta"), version = 14)
write_dta(panel_off, file.path(out_dir_dta, "城市数字经济发展地理工具变量_办公地址_2001-2024.dta"), version = 14)
write_dta(yearly_out, file.path(out_dir_dta, "各年数字化转型程度与比值.dta"), version = 14)
write_dta(cent_out, file.path(out_dir_dta, "城市质心坐标_2021行政区划.dta"), version = 14)

说明：attr() 设置变量标签后，write_dta() 保存时会把这些标签写入 Stata 14 的 .dta 文件；变量名仍为英文（province/city/year/dist_hz_km/ratio_natl_hz/iv_product），便于跨软件使用，在 Stata 中 describe 即可看到中文标签。

4.2 结果可视化

以「各年全国/杭州数字化转型程度之比（两口径）趋势」「2024 年城市地理工具变量空间分布」以及「2024 年代表性城市 IV 柱状图」为例，使用 ggplot2 绘制并保存到 输出/ 目录。

# 图1：各年“全国/杭州”数字化转型程度之比趋势（两口径）
# 此处代码需下载讲义材料查看～
# 图2：2024 年代表性城市地理工具变量（IV）柱状图
reps <- c("杭州市", "上海市", "深圳市", "北京市", "广州市",
          "成都市", "武汉市", "西安市", "重庆市", "乌鲁木齐市")

sub2024 <- panel_reg %>%
  filter(year == 2024, city %in% reps) %>%
  mutate(city = fct_reorder(city, iv_product))

p2 <- ggplot(sub2024, aes(x = city, y = iv_product, fill = city)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = format(round(iv_product), big.mark = ",")),
                vjust = -0.3, family = "sc") +
  scale_fill_manual(values = c("#fed439", "#709ae1", "#8a9197", "#d2af81",
                               "#fd7446", "#d5e4a2", "#197ec0", "#f05c3b",
                               "#46732e", "#71d0f5")) +
  labs(title = "2024 年代表性城市地理工具变量(IV = 距离×比值)",
       x = "城市", y = "IV") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

ggsave(file.path(out_dir, "2024年代表性城市IV.png"), p2,
       width = 10, height = 5.5, dpi = 400, device = png)

说明：图 1 展示两口径比值均在 0.6–0.8 区间波动（杭州长期领先，故比值 < 1，随全国追赶比值上行）；图 2 用柱状图对比代表性城市的 IV，其排序与到杭州距离排序一致（同一年份比值对所有城市相同）。

4.3 结果展示

下图由上述脚本生成：

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算城市数字经济发展的工具变量

相关文章

评论