继续上次的内容，今天我们再来学习如何基于真实专利数据计算数实融合水平。
仿照示例数据，我们这里也需要准备三份数据：
- 专利的 IPC 分类号；
- IPC 与国民经济行业的对照；
- 数字经济与实体经济的行业分类。
准备专利数据
首先我们准备专利的 IPC 分类号数据。专利的原始数据可以从之前分享的这里下载：
1985～2024 年专利申请与授权数据（版本 3，含申请人所处的省市区县）：https://rstata.duanshu.com/#/brief/course/2397451274c546d3a36e156ffc865988
原始数据非常大，我们仅仅保留需要的，由于提供的专利数据是 dta 格式的，所以还是使用 Stata 读取处理更快速。
cap mkdir "patent_small"
forval y = 1985/2024 {
use "/Volumes/C16T/大数据/newIP专利数据分年/`y'", clear
replace 公开公告号 = ustrregexs(1) if ustrregexm(公开公告号, "(.*)[A-Z]$")
duplicates drop 年份 公开公告号, force
duplicates drop 年份 申请号, force
keep newipzlid 年份 IPC 省 省代码 市 市代码 县 县代码
export delimited using patent_small/`y'.csv, replace
}
|
不过由于部分年份的专利数据非常大，所以大家的电脑可能无法直接运行上面的循环。不过也可以通过每次读取部分行来解决：
Stata 读取数据的大小受到电脑内存的限制，因此对于超大文件可以通过下面的方式读取：
use 命令有种用法是：
use [varlist] [if] [in] using filename [, clear nolabel]
|
例如每次读 100 条：use in 1/100 using 文件名.dta, clear
不过读取最后一部分的时候需要知道文件的总行数。
获取文件的总行数可以使用 describe 命令，该命令有种用法是：
describe [varlist] using filename [, file_options]
|
例如获取某个名为 2012.dta 数据的总行数：
由于全部年份的数据非常大，所以附件中只提供了 2012 的：patent_small2/2012.csv
按照参考文献的方法，删除只有一个 IPC 号的专利，然后每个专利的分类号只使用前 10 个：
准备专利 IPC 数据
library(tidyverse)
library(readxl)
library(haven)
patent_data <- read_csv("patent_small/2012.csv") %>%
filter(str_detect(IPC, ";")) %>%
select(newipzlid, IPC) %>%
mutate(IPC = str_replace_all(IPC, " ", "")) %>%
mutate(
ipc_list = str_split(IPC, ";"),
ipc = map_chr(ipc_list, ~ paste(head(.x, 10), collapse = ";"))
) %>%
select(newipzlid, ipc)
patent_data
|
#> # A tibble: 466,947 × 2
#> newipzlid ipc
#> <dbl> <chr>
#> 1 20120565869 G06F19/00;G01R29/08
#> 2 20120196633 G01S7/40;G01R29/02
#> 3 20120195174 G06F21/62;G06F12/14
#> 4 20121092201 H04L29/06;G01R31/00
#> 5 20120565160 B23P9/04;C21D7/04
#> 6 20120565868 G06F17/00;G08G5/04;G01C21/00
#> 7 20121034880 G01H11/06;H04R1/44
#> 8 20120554752 B01J13/02;C08K5/10;C08J7/04;C23F11/00
#> 9 20121681282 G01M9/04;F25D17/02
#> 10 20120565870 G06F17/00;G01C21/00
#> # ℹ 466,937 more rows
|
IPC 与 国民经济行业小类对照
这里我们使用的是「国际专利分类与国民经济行业分类参照关系表（2018）」，该文件下载自国家知识产权局官网。首先我们需要借助外部工具把 pdf 文件转换成 xlsx 文件，我使用的是 Adobe Acrobat。
处理整齐：
uniq_ipc.dta 数据是历年所有专利涉及的 IPC 分类号。例如 2012 年的可以这么处理得到：
patent_data %>%
select(ipc) %>%
distinct() %>%
mutate(ipc_split = str_split(ipc, ";")) %>%
unnest(ipc_split) %>%
rename(uniq_ipc = ipc_split) %>%
select(uniq_ipc) %>%
distinct()
|
#> # A tibble: 45,547 × 1
#> uniq_ipc
#> <chr>
#> 1 G06F19/00
#> 2 G01R29/08
#> 3 G01S7/40
#> 4 G01R29/02
#> 5 G06F21/62
#> 6 G06F12/14
#> 7 H04L29/06
#> 8 G01R31/00
#> 9 B23P9/04
#> 10 C21D7/04
#> # ℹ 45,537 more rows
|
附件中我给大家提前生成好了历年所有专利涉及的 IPC 分类号：
haven::read_dta("uniq_ipc.dta") -> expanded_ipc
expanded_ipc
|
#> # A tibble: 85,187 × 1
#> uniq_ipc
#> <chr>
#> 1 A47G21/10
#> 2 A47G21/00
#> 3 B28B3/22
#> 4 B28B3/20
#> 5 B22D11/06
#> 6 B22D11/10
#> 7 B22D11/18
#> 8 A47L13/24
#> 9 A47L13/20
#> 10 G02B21/14
#> # ℹ 85,177 more rows
|
下面我们需要给每个 IPC 对应上行业小类：
industry_ipc_clean %>%
rename(IPC = 国际专利分类号) %>%
mutate(
IPC = str_replace_all(IPC, "[* ]", ""),
len = str_length(IPC)
) -> ipc_reference
lapply(3:11, function(i){
ipc_reference %>%
filter(len == i) %>%
select(-len) -> len_data
expanded_ipc %>%
mutate(IPC = str_sub(uniq_ipc, 1, i)) %>%
inner_join(len_data, by = "IPC", relationship = "many-to-many") %>%
select(uniq_ipc, 国民经济行业代码)
}) %>%
bind_rows() %>%
distinct(uniq_ipc, 国民经济行业代码) -> all_matched
write_dta(all_matched, "专利数据与行业小类代码简易对照表.dta")
|
可以看到这个数据的格式其实就类似上次课说的：IPC与产业对照表_示例2.dta
数字经济产业与实体产业
按照论文的介绍，实体产业包括下面四类：
industry_ipc_clean %>%
mutate(
行业门类 = str_sub(国民经济行业代码, 1, 1),
实体产业 = case_when(
行业门类 == "C" ~ "制造业",
行业门类 == "A" ~ "农业",
行业门类 %in% c("B", "D", "E") ~ "建筑业及其他工业",
行业门类 %in% c("I", "O") ~ "服务业",
TRUE ~ ""
)
) %>%
filter(实体产业 != "") %>%
select(国民经济行业代码, 实体产业) %>%
distinct() -> entity_industry
write_dta(entity_industry, "实体产业分类.dta")
|
这里其实论文里面说的有些含糊，例如农业到底是特指 A01 还是 A，也就是农林牧渔？
数字经济产业分类使用的是「数字经济及其核心产业统计分类（2021）」，该文件也是来源于国家知识产权局。
library(tidyverse)
library(docxtractr)
doc <- read_docx("数字经济及其核心产业统计分类（2021）.docx")
doc %>%
docx_extract_all_tbls() %>%
.[[1]] %>%
slice(-1, -2) %>%
set_names("大类", "中类", "小类", "名称", "说明", "国民经济行业代码及名称") %>%
type_convert() %>%
mutate(国民经济行业代码 = str_extract(国民经济行业代码及名称, "\\d{4}")) %>%
haven::write_dta("数字经济及其核心产业统计分类.dta")
|
然后再处理：
read_dta("数字经济及其核心产业统计分类.dta") %>%
filter(!is.na(小类), 小类 != "") %>%
filter(str_sub(小类, 1, 2) != "05") %>%
select(国民经济行业代码及名称, 小类) -> digital_economy
digital_economy
|
#> # A tibble: 114 × 2
#> 国民经济行业代码及名称 小类
#> <chr> <chr>
#> 1 3911 计算机整机制造 010101
#> 2 3912 计算机零部件制造 010102
#> 3 3913 计算机外围设备制造 010103
#> 4 3914 工业控制计算机及系统制造 010104
#> 5 3915 信息安全设备制造 010105
#> 6 3919 其他计算机制造 010106
#> 7 3921 通信系统设备制造 010201
#> 8 3922 通信终端设备制造 010202
#> 9 3940 雷达及配套设备制造 010203
#> 10 3931 广播电视节目制作及发射设备制造 010301
#> # ℹ 104 more rows
|
digital_economy %>%
mutate(
v1 = str_extract(国民经济行业代码及名称, "\\d{4}"),
国民经济行业代码及名称 = str_remove(国民经济行业代码及名称, v1),
v2 = str_extract(国民经济行业代码及名称, "\\d{4}"),
国民经济行业代码及名称 = str_remove(国民经济行业代码及名称, v2),
v3 = str_extract(国民经济行业代码及名称, "\\d{4}"),
国民经济行业代码及名称 = str_remove(国民经济行业代码及名称, v3),
v4 = str_extract(国民经济行业代码及名称, "\\d{4}"),
国民经济行业代码及名称 = str_remove(国民经济行业代码及名称, v4)
) %>%
mutate(
数字经济产业 = case_when(
str_sub(小类, 1, 2) == "01" ~ "数字产品制造业",
str_sub(小类, 1, 2) == "02" ~ "数字产品服务业",
str_sub(小类, 1, 2) == "03" ~ "数字技术应用业",
str_sub(小类, 1, 2) == "04" ~ "数字要素驱动业",
TRUE ~ ""
)
) %>%
select(v1, v2, v3, v4, 数字经济产业) %>%
pivot_longer(cols = c(v1, v2, v3, v4), names_to = "var", values_to = "industry2") %>%
filter(!is.na(industry2)) %>%
select(industry2, 数字经济产业) %>%
distinct() -> digital_economy_clean
write_dta(digital_economy_clean, "数字经济核心产业.dta")
|
然后就可以区分数字与实体了：
industry_ipc_clean %>%
select(国民经济行业代码) %>%
distinct() %>%
mutate(industry2 = str_sub(国民经济行业代码, 2, -1)) %>%
left_join(digital_economy_clean, by = "industry2") %>%
left_join(entity_industry, by = "国民经济行业代码") %>%
select(国民经济行业代码, 数字经济产业, 实体产业) %>%
mutate(class = as.numeric(!is.na(数字经济产业))) %>%
group_by(国民经济行业代码) %>%
summarize(class = sum(class)) %>%
mutate(
类别 = ifelse(class > 0, "数字经济产业", "实体产业")
) %>%
select(国民经济行业代码, 类别) -> industry_classification
industry_classification %>%
count(类别)
|
#> # A tibble: 2 × 2
#> 类别 n
#> <chr> <int>
#> 1 实体产业 732
#> 2 数字经济产业 89
|
write_dta(industry_classification, "产业数实分类.dta")
|
这里感觉论文里面说的也比较含糊，所以我采取的方案是如果该小类同时属于数字和实体，那么该小类就归类为数字经济行业。
计算 IPC 融合矩阵
#> # A tibble: 466,947 × 2
#> newipzlid ipc
#> <dbl> <chr>
#> 1 20120565869 G06F19/00;G01R29/08
#> 2 20120196633 G01S7/40;G01R29/02
#> 3 20120195174 G06F21/62;G06F12/14
#> 4 20121092201 H04L29/06;G01R31/00
#> 5 20120565160 B23P9/04;C21D7/04
#> 6 20120565868 G06F17/00;G08G5/04;G01C21/00
#> 7 20121034880 G01H11/06;H04R1/44
#> 8 20120554752 B01J13/02;C08K5/10;C08J7/04;C23F11/00
#> 9 20121681282 G01M9/04;F25D17/02
#> 10 20120565870 G06F17/00;G01C21/00
#> # ℹ 466,937 more rows
|
计算行业共现矩阵
分别使用 ipc 匹配各自的行业即可：
ipc_matrix %>%
rename(uniq_ipc = ipc1) %>%
left_join(all_matched, by = "uniq_ipc", relationship = "many-to-many") %>%
rename(industry1 = 国民经济行业代码) %>%
select(-uniq_ipc) %>%
rename(uniq_ipc = ipc2) %>%
left_join(all_matched, by = "uniq_ipc", relationship = "many-to-many") %>%
rename(industry2 = 国民经济行业代码) %>%
select(-uniq_ipc) %>%
filter(!is.na(industry1) & !is.na(industry2)) %>%
mutate(value = ifelse(industry1 == industry2, 0, value)) %>%
group_by(industry1, industry2) %>%
summarize(value = sum(value), .groups = "drop") %>%
filter(!is.na(industry1) & !is.na(industry2)) -> industry_matrix
industry_matrix
|
#> # A tibble: 299,131 × 3
#> industry1 industry2 value
#> <chr> <chr> <dbl>
#> 1 A0111 A0111 0
#> 2 A0111 A0112 3858
#> 3 A0111 A0113 3858
#> 4 A0111 A0119 3856
#> 5 A0111 A0121 3856
#> 6 A0111 A0122 3856
#> 7 A0111 A0123 3861
#> 8 A0111 A0131 3859
#> 9 A0111 A0132 3856
#> 10 A0111 A0133 3856
#> # ℹ 299,121 more rows
|
write_dta(industry_matrix, "industrymat.dta")
|
计算数实融合
然后就可以计算数实融合水平了：
由此我们便计算得到了 2012 年的数实融合水平，不过这个结果比参考文献的大挺多，但是仔细检查了很多遍也没发现计算过程的问题，也有可能是原作者的一些处理细节我我的不一样。
循环各年的就可以计算得到各年的结果了。下次课我们再讲解如何计算分城市和产业的～
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算数实融合水平（二）：基于真实专利数据
评论