在附件文献「地区间技术互补与创新驱动型经济增长_郑江淮」中提到了地区间技术互补指数的指标：
今天的课程中我们将一起使用 R 语言复现这个指标。
使用的专利数据就是之前给大家分享的这个：
1985～2024 年专利申请与授权数据（版本 3，含申请人所处的省市区县）：https://rstata.duanshu.com/#/brief/course/2397451274c546d3a36e156ffc865988
在附件中我提供了 2001 年的专利数据作为示例。
首先处理专利数据，主要是去除重复专利、筛选发明专利、拆分 IPC 等：
library(tidyverse)
haven::read_dta("2001.dta") -> df
df
|
#> # A tibble: 146,970 × 49
#> newipzlid 年份 申请日 标题 摘要 申请人 公开公告号 公开公告日 申请号
#> <chr> <dbl> <date> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 2001101496 2001 2001-05-30 模拟实战射击… "" "" CN0110104… 2002-05-30 CN011…
#> 2 2001037217 2001 2001-03-12 一种碲镉汞晶… "" "" CN0110096… 2004-04-14 CN011…
#> 3 2001082750 2001 2001-10-15 小型无刷直流… "一种小… "上海元山… CN2504411Y 2002-08-07 CN012…
#> 4 2001012438 2001 2001-04-25 扑克(千叶红… "1.省… "上海宇琛… CN3210531D 2001-11-21 CN013…
#> 5 2001005665 2001 2001-03-14 扑克(回形针… "1.平… "上海宇琛… CN3205260D 2001-10-17 CN013…
#> 6 2001073966 2001 2001-12-25 箱式变电站(… "右视图… "上海瑞华… CN3254449D 2002-09-11 CN013…
#> 7 2001033740 2001 2001-05-09 内压缩式倍速… "本发明… "凌岗" CN1170063C 2004-10-06 CN011…
#> 8 2001043836 2001 2001-07-27 用位相反转的… "一种用… "中国科学… CN2554765Y 2003-06-04 CN012…
#> 9 2001034900 2001 2001-11-15 医药用氧化铁… "本发明… "上海氧化… CN1162477C 2004-08-18 CN011…
#> 10 2001042879 2001 2001-12-14 基因芯片分子… "本发明… "宋克" CN1424405A 2003-06-18 CN011…
#> # ℹ 146,960 more rows
#> # ℹ 40 more variables: 专利类型 <chr>, 公开国别 <chr>, 首项权利要求 <chr>,
#> # 独立权利要求 <chr>, 文献页数 <chr>, IPC主分类 <chr>, IPC <chr>,
#> # 洛迦诺分类号 <chr>, 当前权利人 <chr>, 申请人类型 <chr>,
#> # 申请人国家_地区 <chr>, 申请人地址 <chr>, 当前专利权人地址 <chr>,
#> # 工商注册地址 <chr>, 工商公司类型 <chr>, 工商成立日期 <chr>,
#> # 工商统一社会信用代码 <chr>, 工商注册号 <chr>, 工商上市代码 <chr>, …
|
由于专利数据中同时包含了专利的申请和授权公告，所以有重复的，在统计数量前需要进行去重。通常公开公告号的结尾字母用以区分公告类别，例如 A 代表发明专利的申请公开，B 代表发明专利的授权公告，U 代表实用新型专利的授权公告，S 代表外观设计专利的授权公告。因此可以通过去除公开公告号结尾的字母进行去重，另外申请号也存在重复的，也需要进行去重。
df %>%
mutate(公开公告号 = str_remove_all(公开公告号, "[A-Z]$")) %>%
distinct(公开公告号, .keep_all = T) %>%
distinct(申请号, .keep_all = T) -> df
df %>%
filter(str_detect(专利类型, "发明")) %>%
select(year = 年份, city = 市, ipc = IPC, newipzlid) %>%
filter(!is.na(city), city != "", !is.na(ipc), ipc != "") -> df
df
|
#> # A tibble: 27,425 × 4
#> year city ipc newipzlid
#> <dbl> <chr> <chr> <chr>
#> 1 2001 上海市 F04C18/14 2001033740
#> 2 2001 上海市 C09C1/24 2001034900
#> 3 2001 上海市 C12Q1/68; B01L3/00; C40B40/06 2001042879
#> 4 2001 上海市 G21K7/00; G03H5/00; G02B21/00 2001011181
#> 5 2001 上海市 H01S3/08; H01S3/11; H01S3/04 2001038650
#> 6 2001 上海市 G02B6/34 2001127808
#> 7 2001 上海市 H05G2/00; H01S4/00 2001080181
#> 8 2001 上海市 A61B3/12; A61B6/00 2001034395
#> 9 2001 上海市 G02B6/35 2001029535
#> 10 2001 上海市 G01B11/24; G01B9/021 2001001280
#> # ℹ 27,415 more rows
|
IPC 分类号使用分号分隔，可以使用 tidytext::unnest_tokens() 函数进行拆分和转置：
df %>%
tidytext::unnest_tokens("ipc", "ipc", token = stringr::str_split,
pattern = "; ", to_lower = F) %>%
mutate(ipc = str_remove_all(ipc, "[\\s.]"),
ipc = str_sub(ipc, 1, 4)) -> patent_data
patent_data
|
#> # A tibble: 88,407 × 4
#> year city ipc newipzlid
#> <dbl> <chr> <chr> <chr>
#> 1 2001 上海市 F04C 2001033740
#> 2 2001 上海市 C09C 2001034900
#> 3 2001 上海市 C12Q 2001042879
#> 4 2001 上海市 B01L 2001042879
#> 5 2001 上海市 C40B 2001042879
#> 6 2001 上海市 G21K 2001011181
#> 7 2001 上海市 G03H 2001011181
#> 8 2001 上海市 G02B 2001011181
#> 9 2001 上海市 H01S 2001038650
#> 10 2001 上海市 H01S 2001038650
#> # ℹ 88,397 more rows
|
下面就可以开始一步步的计算了。
步骤1：计算每个城市-年份-IPC的专利数量
patent_data %>%
count(year, city, ipc, name = "patent_count") -> patent_count
|
步骤2：计算每个城市每年所有 IPC 的专利总数
patent_count %>%
group_by(year, city) %>%
summarise(total_city = sum(patent_count), .groups = "drop") -> city_year_total
|
步骤3：计算每年每个IPC的全国专利总数
patent_count %>%
group_by(year, ipc) %>%
summarise(total_ipc = sum(patent_count), .groups = "drop") -> ipc_year_total
|
步骤4：计算 RPCA（显性专利比较优势指数）
patent_count %>%
left_join(city_year_total, by = c("year", "city")) %>%
left_join(ipc_year_total, by = c("year", "ipc")) %>%
mutate(
share_city = patent_count / total_city,
share_national = total_ipc / sum(patent_count),
RPCA = share_city / share_national
) %>%
mutate(x_id = if_else(RPCA > 1, 1, 0)) -> rpc_data
rpc_data
|
#> # A tibble: 11,158 × 10
#> year city ipc patent_count total_city total_ipc share_city share_national
#> <dbl> <chr> <chr> <int> <int> <int> <dbl> <dbl>
#> 1 2001 七台河市… A61K 5 14 9936 0.357 0.112
#> 2 2001 七台河市… A61P 4 14 9825 0.286 0.111
#> 3 2001 七台河市… C07C 4 14 1846 0.286 0.0209
#> 4 2001 七台河市… G09B 1 14 140 0.0714 0.00158
#> 5 2001 三明市…… A01N 2 35 879 0.0571 0.00994
#> 6 2001 三明市…… A23L 2 35 2819 0.0571 0.0319
#> 7 2001 三明市…… A61K 2 35 9936 0.0571 0.112
#> 8 2001 三明市…… A61P 6 35 9825 0.171 0.111
#> 9 2001 三明市…… B21C 1 35 73 0.0286 0.000826
#> 10 2001 三明市…… B27N 1 35 29 0.0286 0.000328
#> # ℹ 11,148 more rows
#> # ℹ 2 more variables: RPCA <dbl>, x_id <dbl>
|
步骤5：构建技术关联矩阵 Φ（基于所有年份的共现）
注意这里是需要把所有年份的发明专利数据放在一起处理，不过这里只提供了 2001 年的样本，所以还是以 2001 年的数据为例进行演示。
patent_data %>%
distinct(newipzlid, ipc) %>%
group_by(newipzlid) %>%
filter(n() > 1) %>%
summarise(
ipc_list = list(ipc),
.groups = "drop"
) %>%
mutate(ipc_pairs = map(ipc_list, ~expand.grid(.x, .x))) %>%
unnest(ipc_pairs) %>%
filter(Var1 != Var2) %>%
rename(ipc1 = Var1, ipc2 = Var2) %>%
count(ipc1, ipc2, name = "C_ij") -> ipc_pairs
patent_data %>%
distinct(newipzlid, ipc) %>%
count(ipc, name = "total_patents") -> ipc_total_patents
ipc_total_patents
|
#> # A tibble: 603 × 2
#> ipc total_patents
#> <chr> <int>
#> 1 A01B 16
#> 2 A01C 51
#> 3 A01D 10
#> 4 A01F 6
#> 5 A01G 246
#> 6 A01H 117
#> 7 A01K 90
#> 8 A01M 11
#> 9 A01N 481
#> 10 A01P 71
#> # ℹ 593 more rows
|
# 构建技术关联矩阵 Φ
# 此处代码需下载讲义材料查看～
|
步骤6：构建每个城市每年的技术优势向量 M_mat
rpc_data %>%
select(year, city, ipc, x_id) %>%
filter(ipc %in% unique(phi_mat$ipc1)) %>%
spread(ipc, x_id, fill = 0) -> M_mat
|
步骤7：计算两个城市 c 和 d 之间的技术互补指数 comp_cd
year_select <- 2001
city_c <- "北京市"
city_d <- "上海市"
M_mat %>%
filter(year == year_select, city == city_c) %>%
select(-year, -city) %>%
as.numeric() %>%
setNames(colnames(M_mat)[-c(1:2)]) -> M_c_vec
M_mat %>%
filter(year == year_select, city == city_d) %>%
select(-year, -city) %>%
as.numeric() %>%
setNames(colnames(M_mat)[-c(1:2)]) -> M_d_vec
A_cd <- pmin(M_c_vec, M_d_vec)
M_c_tilde <- M_c_vec - A_cd
M_d_tilde <- M_d_vec - A_cd
techs <- names(M_c_vec)
phi_sub <- phi_mat_wide[techs, techs]
numerator <- t(M_c_tilde) %*% phi_sub %*% M_d_tilde
denominator <- sqrt(sum(M_c_tilde) * sum(M_d_tilde))
comp_cd <- as.numeric(numerator / denominator)
print(comp_cd)
|
由此就可以编写一个计算两个城市间技术互补指数的函数了：
获取所有城市对并计算技术互补指数：
unique(M_mat$city) %>%
crossing(unique(M_mat$city)) %>%
set_names("city1", "city2") %>%
filter(city1 != city2) %>%
slice_sample(n = 100) %>%
mutate(comp_cd = map2_dbl(city1, city2, comp_cd_fun,
M_mat = M_mat,
phi_mat_wide = phi_mat_wide)) -> dfres
dfres
mean(dfres$comp_cd, na.rm = T)
|
不过完全计算所有的使用单线程会非常缓慢，可以考虑使用多线程：
这种多线程方法非常适合使用 map 族函数的代码。
如果想要获取所有年份的数据，循环年份即可，不过需要注意技术关联矩阵是需要基于所有年份的专利数据重新计算。
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言计算地区间技术互补指数
评论