上次课程的代码里面计算的时候忽略了“如果某个专利 j 没有引用专利 i 但是引用了其后项引用，则 cd = 0” 的情况，导致计算结果里面 1 和 -1 出现的太多。所以赶紧更正了下。
在《Nature》论文：「Papers and patents are becoming less disruptive over time」 中使用了 CD 指数数据，该指标的计算公式如下：
其中当专利 j 引用了专利 i 时 fit 取1，否则取 0；若专利 j 引用了专利 i 的后向引用专利时 bit 取 1，否则取 0；
附件中也提供了该论文的 pdf 文件。
论文中也提供了一个示意图帮助我们理解这个指数：
看起来很复杂，其实只需要关注三个：
- 若后面的专利既引用专利 i 又引用专利 i 的后向引用专利，那么专利 i 的该项引用的 CD1 指数 -1；
- 若后面的专利只引用专利 i 但未引用专利 i 的后向引用专利，那么专利i的该项引用的 CDI 指数为 1；
- 若后面的专利只引用专利 i 的后向引用专利，但是没有引用专利 i，那么专利 i 的该项引用的 CDI 指数为 0；
这样直接计算得到的是每个专利引用关系的 CD 指数，平均之后就得到了某个专利的 CD 指数。
今天的课程中我们将以上市公司专利为例进行讲解。
在计算之前我们需要准备两组数据：
- 一个是专利的引用关系，也就是每个专利引用了其他的哪些专利；
- 第二个是专利的属性信息，也就是每个公司每年申请了哪些专利。
这两个数据我们都分享过。全部专利的引用与被引用关系数据在这里：
1985~2024 年全部专利引用与被引用详细信息：https://rstata.duanshu.com/#/brief/course/225ad0b59a9945e1831d2e8b96ca1001；
由于全部的数据非常大，所以附件中仅仅提供了含有 newipzlid 和 newipzlid2 的数据作，这两个变量就足够用来计算专利的 CD 指数了。通过类似如下代码即可处理得到：
这部分代码就不用运行了，“全部专利引用与被引用信息分年.rds” 是一个超过 40 GB 的大数据，代码需要电脑性能足够好才能运行。
library(tidyverse)
read_rds("全部专利引用与被引用信息分年.rds") -> df1
df1 %>%
mutate(公开公告号_original = str_remove(公开公告号_original, "[A-Z]$"),
公开公告号 = str_remove(公开公告号, "[A-Z]$")) %>%
distinct(公开公告号_original, 公开公告号, .keep_all = T) -> df2
bind_rows(
df2 %>%
filter(引用或被引用 %in% c("他引信息", "自引信息")) %>%
set_names(c("newipzlid1", "申请日1", "公开公告号1",
"省1", "省代码1", "市1", "市代码1", "县1", "县代码1",
"引用或被引用",
"newipzlid2", "申请日2", "公开公告号2",
"省2", "省代码2", "市2", "市代码2", "县2", "县代码2")) %>%
select(-引用或被引用),
df2 %>%
filter(引用或被引用 %in% c("被他引信息", "被自引信息")) %>%
set_names(c("newipzlid2", "申请日2", "公开公告号2",
"省2", "省代码2", "市2", "市代码2", "县2", "县代码2",
"引用或被引用",
"newipzlid1", "申请日1", "公开公告号1",
"省1", "省代码1", "市1", "市代码1", "县1", "县代码1")) %>%
select(-引用或被引用)
) %>%
distinct() -> df0a
df0a
df0a %>%
write_rds("全部专利引用信息（去重后）.rds")
df0a %>%
select(newipzlid = newipzlid1, newipzlid2,
date = 申请日1, date2 = 申请日2) -> patent_citations
patent_citations %>%
write_rds("patent_citations.rds")
|
附件中提供的“patent_citations.rds”文件夹就是根据这个代码处理得到的。
考虑专利被引用情况通常需要限定窗口，例如这里限定 5 年的引用窗口，也就是只保留引用间隔 5 年内的：
library(tidyverse)
read_rds("patent_citations.rds") -> patent_citations
patent_citations %>%
filter(date - date2 <= dyears(5) & date - date2 >= 0) -> patent_citations2
patent_citations2 %>%
write_rds("patent_citations2.rds")
patent_citations2 %>%
slice_sample(n = 10)
|
#> # A tibble: 10 × 4
#> newipzlid newipzlid2 date date2
#> <dbl> <dbl> <date> <date>
#> 1 20141301234 20140334274 2014-09-15 2014-04-02
#> 2 20183713037 20172740680 2018-12-25 2017-11-10
#> 3 20213287175 20201944456 2021-09-15 2020-06-06
#> 4 20211043088 20171503410 2021-04-12 2017-07-21
#> 5 20211160131 20162258019 2021-04-23 2016-10-13
#> 6 20221814582 20182316620 2022-06-02 2018-08-30
#> 7 20173363949 20160236524 2017-12-28 2016-02-15
#> 8 20191239857 20170241169 2019-05-29 2017-02-28
#> 9 20191095681 20190968608 2019-04-24 2019-04-24
#> 10 20194387760 20190298138 2019-12-31 2019-01-15
|
附件中也提供了 “patent_citations2.rds”，因此可以直接读取这个。
由于需要考虑这种情况：“若后面的专利只引用专利 i 的后向引用专利，但是没有引用专利 i，那么专利 i 的该项引用的 CDI 指数为 0”， 所以我们就不能使用之前旧版本课程里面的方法，因为那个方法只能获取引用专利 i 的情况。也就是我们必须得循环所有的专利来计算了，效率很低，但是没特别好的办法了。
以某个专利为例：
id <- 20190809613
patent_citations2 %>%
filter(newipzlid2 == id) -> testpat
testpat
|
#> # A tibble: 1 × 4
#> newipzlid newipzlid2 date date2
#> <dbl> <dbl> <date> <date>
#> 1 20210978229 20190809613 2021-03-02 2019-04-22
|
patent_citations2 %>%
filter(newipzlid == id) %>%
select(newipzlid2) %>%
pull(newipzlid2) -> backrefs
backrefs
|
patent_citations2 %>%
filter(date >= testpat$date2[1]) -> part0
part0 %>%
filter(newipzlid2 == id) %>%
select(newipzlid) %>%
distinct(newipzlid) -> part1
part1
|
#> # A tibble: 1 × 1
#> newipzlid
#> <dbl>
#> 1 20210978229
|
part0 %>%
filter(newipzlid2 %in% backrefs) %>%
select(newipzlid) %>%
distinct(newipzlid) -> part2
part2
|
#> # A tibble: 0 × 1
#> # ℹ 1 variable: newipzlid <dbl>
|
# 由此就可以分出三种情况：
# 此处代码需要下载讲义材料查看～
|
这样就得到了一个专利的 CD 指数。循环所有的即可。
不过所以需要计算的专利总数超过了 4000 万，也就是需要循环四千次上面的代码，效率非常低。
我们可以进一步思考一种特殊的情况，也就是 如果某个专利 i，也就是 patent_citations2 里面的 newipzlid2 一列的专利，没有引用其他专利，那么就是说这个专利 i 没有后向引用，因此它就不存在 CD 指数为 0 或者为 -1 的情况，都是 1，平均起来也都是 1。所以这些专利就可以筛选出来不必计算了：
附件中也提供了“nobackrefs.rds”文件。
那么我们需要计算的其实只有这些：
read_rds("nobackrefs.rds") -> nobackrefs
nobackrefs
|
#> # A tibble: 8,277,416 × 2
#> newipzlid2 cd_index
#> <dbl> <dbl>
#> 1 19850027 1
#> 2 19857435 1
#> 3 19850937 1
#> 4 19854634 1
#> 5 19859532 1
#> 6 19850147 1
#> 7 19850143 1
#> 8 19851158 1
#> 9 19850131 1
#> 10 19851493 1
#> # ℹ 8,277,406 more rows
|
patent_citations2 %>%
distinct(newipzlid2) %>%
anti_join(nobackrefs) -> withbackrefs
withbackrefs
|
#> # A tibble: 7,028,753 × 1
#> newipzlid2
#> <dbl>
#> 1 198608410
#> 2 19850568
#> 3 19851877
#> 4 19854601
#> 5 19856362
#> 6 198601453
#> 7 198600518
#> 8 198605384
#> 9 198604263
#> 10 19857890
#> # ℹ 7,028,743 more rows
|
然后读取上市公司专利信息：
read_csv("2010年上市公司与专利数据匹配结果.csv") -> sspat
sspat %>%
distinct(newipzlid) %>%
filter(newipzlid %in% withbackrefs$newipzlid2) -> testmat
testmat
|
#> # A tibble: 19,385 × 1
#> newipzlid
#> <dbl>
#> 1 2010000002
#> 2 2010000004
#> 3 2010000005
#> 4 2010000006
#> 5 2010000091
#> 6 2010000147
#> 7 2010000281
#> 8 2010000287
#> 9 2010000346
#> 10 2010000353
#> # ℹ 19,375 more rows
|
testmat %>%
write_rds("待计算.rds")
|
不过这里我们就不完全计算了，而是选择 20 个：
#> Time difference of 2.751031 mins
|
通过这个结果我们就可以估算计算全部的需要多长时间了。
显然单线程计算这个是个巨耗时的工作，改用多线程：
这部分代码我单独存放的，可以方便直接全选运行。电脑性能不好的小伙伴还是建议老老实实单线程。
合并计算结果：
library(tidyverse)
fs::dir_ls("res1") %>%
lapply(readr::read_csv, col_names = F) %>%
bind_rows() %>%
set_names("newipzlid", "cd_index") -> cdindex
read_rds("nobackrefs.rds") %>%
set_names("newipzlid", "cd_index") -> nobackrefs
|
然后我们再次读取和处理上市公司专利，这次我们就需要进行去重了，不过我们需要尽可能保留出现在 patent_citations 里面的专利：
附件中也提供了这个文件：“具有引用与被引用信息的专利unique_newipzlid.rds”。
read_rds("/Volumes/ADT/常用大数据/具有引用与被引用信息的专利unique_newipzlid.rds") -> unique_newipzlid
readr::read_csv("2010年上市公司与专利数据匹配结果.csv") -> df1
df1 %>%
left_join(unique_newipzlid %>% mutate(value = 1)) -> df1b
df1b %>%
mutate(公开公告号 = str_remove(公开公告号, "[A-Z]$")) %>%
mutate(value = if_else(is.na(value), 0, value)) %>%
group_by(股票代码, 年份, 公开公告号) %>%
filter(value == max(value, na.rm = T)) %>%
ungroup() %>%
distinct(股票代码, 年份, 公开公告号, .keep_all = T) %>%
select(-公开公告号) %>%
group_by(股票代码, 年份, 申请号) %>%
filter(value == max(value, na.rm = T)) %>%
ungroup() %>%
distinct(股票代码, 年份, 申请号, .keep_all = T) %>%
select(-申请号, -value) %>%
select(股票代码, newipzlid, 年份) -> patent_info
patent_info
|
和上述计算结果链接：
patent_info %>%
left_join(
bind_rows(cdindex, nobackrefs)
) -> dfres
dfres %>%
mutate(cd_index = if_else(is.na(cd_index), 0, cd_index)) -> dfres
dfres
|
由于上面并没有计算所有的专利 CD 指数，这里这里直接用 0 代替缺失的其实是不对的，但是如果你计算了全部可计算的专利，这样做是可以的。没有被引用的专利 CD 指数是 0。
分公司汇总就可以得到公司层面的结果了：
dfres %>%
group_by(股票代码, 年份) %>%
summarise(cd_index = mean(cd_index, na.rm = T), .groups = "drop")
|
仔细学习的小伙伴也会发现，其实这个代码是可以计算全部专利 CD 指数的，不过由于运算过于耗时，还是只计算自己需要的就好。
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算专利的创新突破度：CD指数（更新）
评论