名师讲堂｜使用 R 语言测算专利的创新突破度：CD 指数

在《Nature》论文：「Papers and patents are becoming less disruptive over time」 中使用了 CD 指数数据，该指标的计算公式如下：

其中当专利 j 引用了专利 i 时 fit 取1，否则取 0；若专利 j 引用了专利 i 的后向引用专利时 bit 取 1，否则取 0；

附件中也提供了该论文的 pdf 文件。

论文中也提供了一个示意图帮助我们理解这个指数：

看起来很复杂，其实只需要关注两个：

1. 若后面的专利既引用专利 i 又引用专利 i 的后向引用专利，那么专利 i 的该项引用的 CD1 指数 -1；
2. 若后面的专利只引用专利 i 但未引用专利 i 的后向引用专利，那么专利i的该项引用的 CDI 指数为 1；
3. 其他的情况都是 0 ；

这样直接计算得到的是每个专利引用关系的 CD 指数，平均之后就得到了某个专利的 CD 指数。

今天的课程中我们将以上市公司专利为例进行讲解。

在计算之前我们需要准备两组数据：

1. 一个是专利的引用关系，也就是每个专利引用了其他的哪些专利；
2. 第二个是专利的属性信息，也就是每个公司每年申请了哪些专利。

这两个数据我们都分享过。全部专利的引用与被引用关系数据在这里：

1985~2024 年全部专利引用与被引用详细信息：https://rstata.duanshu.com/#/brief/course/225ad0b59a9945e1831d2e8b96ca1001；

由于全部的数据非常大，所以附件中仅仅提供了 2010 年前的数据作为样本，通过类似如下代码即可处理得到：

这部分代码就不用运行了，“全部专利引用与被引用信息分年.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")

# 选择 2010 年前的样本
read_rds("全部专利引用信息（去重后）.rds") %>%
  filter(year(申请日1) <= 2010) %>%
  write_csv("2010年及之前专利引用信息.csv")

专利的属性信息就是上市公司专利匹配结果：

1985～2024 年上市公司与专利数据匹配结果（版本3， 含申请、授权信息）：https://rstata.duanshu.com/#/brief/course/04100321f88b411f90429be934934bff

同样这个数据也是非常巨大，附件中提供了 “2010年上市公司与专利数据匹配结果.csv” 样本。

需要注意，如果想要计算 2010 年的专利 CD 指数，需要使用全部的 2010 年及之前的专利引用信息。因为这里需要专利引用的引用信息，2010 年的专利会引用 2010 年之前的专利，然后这些专利还会引用更早期的专利。

下面我们开始使用 R 语言处理，首先准备数据：

library(tidyverse)

# 在计算 CD 指数前，我们需要准备两个数据，一个是专利的引用关系，也就是每个专利引用了其他的哪些专利，第二个是专利的属性信息，也就是每个公司每年申请了哪些专利。

# 首先是 2010 年全部专利的引用关系
read_csv("2010年及之前专利引用信息.csv") -> df

df %>%
  select(newipzlid = newipzlid1, newipzlid2) -> patent_citations

patent_citations

#> # A tibble: 2,327,648 × 2
#>    newipzlid newipzlid2
#>        <dbl>      <dbl>
#>  1  19850046   19850027
#>  2  19850046   19857435
#>  3  19850568   19850937
#>  4  19850908   19854634
#>  5  19850908   19859532
#>  6  19851074   19850147
#>  7  19851074   19850143
#>  8  19851742   19851158
#>  9  19851759   19850131
#> 10  19851877   19851493
#> # ℹ 2,327,638 more rows

# patent_citations 里面包含的专利编号
unique(c(patent_citations$newipzlid, patent_citations$newipzlid2)) %>%
  as_tibble() %>%
  set_names("newipzlid") -> unique_newipzlid

unique_newipzlid

#> # A tibble: 1,732,706 × 1
#>    newipzlid
#>        <dbl>
#>  1  19850046
#>  2  19850568
#>  3  19850908
#>  4  19851074
#>  5  19851742
#>  6  19851759
#>  7  19851877
#>  8  19851978
#>  9  19852054
#> 10  19852055
#> # ℹ 1,732,696 more rows

# 上市公司专利数据
readr::read_csv("2010年上市公司与专利数据匹配结果.csv") -> df1
df1 %>%
  left_join(unique_newipzlid %>% mutate(value = 1)) -> df1b

# 由于专利数据里面包含重复的，所以我们要进行去重，不过去重的时候要尽可能保留 unique_newipzlid 里面的
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) %>%
  rename(firm_id = 股票代码, patent_id = newipzlid, year = 年份) -> patent_info

patent_info

#> # A tibble: 63,434 × 3
#>    firm_id  patent_id  year
#>    <chr>        <dbl> <dbl>
#>  1 002705  2010000002  2010
#>  2 000333  2010000004  2010
#>  3 000333  2010000005  2010
#>  4 000333  2010000006  2010
#>  5 002705  2010000016  2010
#>  6 603630  2010000091  2010
#>  7 600467  2010000147  2010
#>  8 002728  2010000188  2010
#>  9 301372  2010000262  2010
#> 10 000825  2010000264  2010
#> # ℹ 63,424 more rows

patent_citations

#> # A tibble: 2,327,648 × 2
#>    newipzlid newipzlid2
#>        <dbl>      <dbl>
#>  1  19850046   19850027
#>  2  19850046   19857435
#>  3  19850568   19850937
#>  4  19850908   19854634
#>  5  19850908   19859532
#>  6  19851074   19850147
#>  7  19851074   19850143
#>  8  19851742   19851158
#>  9  19851759   19850131
#> 10  19851877   19851493
#> # ℹ 2,327,638 more rows

准备后向引用关系，也就是每个专利引用了哪些专利：

patent_citations %>%
  rename(patent = newipzlid, backward_cited = newipzlid2) -> backward_refs

backward_refs

#> # A tibble: 2,327,648 × 2
#>      patent backward_cited
#>       <dbl>          <dbl>
#>  1 19850046       19850027
#>  2 19850046       19857435
#>  3 19850568       19850937
#>  4 19850908       19854634
#>  5 19850908       19859532
#>  6 19851074       19850147
#>  7 19851074       19850143
#>  8 19851742       19851158
#>  9 19851759       19850131
#> 10 19851877       19851493
#> # ℹ 2,327,638 more rows

计算每次引用的 CD1 值：

# 当专利 j 引用了专利 i 时 fit 取1，否则取 0；
# 若专利 j 引用了专利 i 的后向引用专利时 bit 取 1，否则取 0；
# 若后面的专利既引用专利 i 又引用专利 i 的后向引用专利，那么专利 i 的该项引用的 CD1 指数 -1
# 若后面的专利只引用专利 i 但未引用专利 i 的后向引用专利，那么专利i的该项引用的 CDI 指数为 1
patent_citations %>%
  rename(citing_patent = newipzlid, cited_patent = newipzlid2) %>%
  # 获取被引用专利的后向引用专利（即它引用的专利）
  left_join(backward_refs, by = c("cited_patent" = "patent"), relationship = "many-to-many") %>%
  # 检查引用专利是否也引用了这些后向引用专利
  left_join(
    patent_citations %>%
      rename(citing_patent = newipzlid,
             check_cited = newipzlid2),
    by = "citing_patent",
    relationship = "many-to-many"
  ) %>%
  # 标记是否同时引用了后向引用专利
  mutate(
    cites_backward = (backward_cited == check_cited) &
      (backward_cited != cited_patent) &
      !is.na(backward_cited)
  ) %>%
  # 这里运行下面这句就明白上面的代码在做什么准备了：
  # filter(cites_backward == 1)
  # 计算 CD1
  group_by(citing_patent, cited_patent) %>%
  summarise(
    cd1 = if_else(any(cites_backward), -1, 1),
    .groups = "drop"
  ) -> cd1_values

cd1_values

#> # A tibble: 2,327,648 × 3
#>    citing_patent cited_patent   cd1
#>            <dbl>        <dbl> <dbl>
#>  1      19850046     19850027     1
#>  2      19850046     19857435     1
#>  3      19850568     19850937     1
#>  4      19850908     19854634     1
#>  5      19850908     19859532     1
#>  6      19851074     19850142     1
#>  7      19851074     19850143     1
#>  8      19851074     19850146     1
#>  9      19851074     19850147     1
#> 10      19851742     19851158     1
#> # ℹ 2,327,638 more rows

cd1_values 是专利引用对的结果，分组汇总就得到了专利层面的结果：

cd2_values <- cd1_values %>%
  group_by(cited_patent) %>%
  summarise(
    cd2 = mean(cd1),
    n_citations = n(),
    .groups = "drop"
  )

再分企业平均就可以计算企业年度的 CD 指数：

cd_index <- cd2_values %>%
  left_join(patent_info, by = c("cited_patent" = "patent_id")) %>%
  group_by(firm_id, year) %>%
  summarise(
    CD = mean(cd2),
    n_patents = n(),
    .groups = "drop"
  )

cd_index

#> # A tibble: 1,049 × 4
#>    firm_id  year    CD n_patents
#>    <chr>   <dbl> <dbl>     <int>
#>  1 000002   2010 1             3
#>  2 000012   2010 1             1
#>  3 000016   2010 1             7
#>  4 000021   2010 1             4
#>  5 000030   2010 1             1
#>  6 000032   2010 1             5
#>  7 000039   2010 0.6          10
#>  8 000050   2010 0.778         3
#>  9 000055   2010 1             2
#> 10 000063   2010 0.686       351
#> # ℹ 1,039 more rows

也可以编写一个函数快速计算：

calculate_cd_index_verbose <- function(patent_citations, patent_info) {
  # 1. 准备后向引用关系：每个专利引用了哪些专利
  backward_refs <- patent_citations %>%
    rename(patent = newipzlid, backward_cited = newipzlid2)

  # 2. 计算每次引用的CD1值
  cd1_values <- patent_citations %>%
    rename(citing_patent = newipzlid, cited_patent = newipzlid2) %>%
    # 获取被引用专利的后向引用专利（即它引用的专利）
    left_join(backward_refs, by = c("cited_patent" = "patent")) %>%
    # 检查引用专利是否也引用了这些后向引用专利
    left_join(
      patent_citations %>% rename(citing_patent = newipzlid,
                                  check_cited = newipzlid2),
      by = "citing_patent",
      relationship = "many-to-many"
    ) %>%
    # 标记是否同时引用了后向引用专利
    mutate(
      cites_backward = (backward_cited == check_cited) &
        (backward_cited != cited_patent) &
        !is.na(backward_cited)
    ) %>%
    # 计算CD1
    group_by(citing_patent, cited_patent) %>%
    summarise(
      cd1 = ifelse(any(cites_backward), -1, 1),
      .groups = "drop"
    )

  # 3. 计算专利级别的CD2
  cd2_values <- cd1_values %>%
    group_by(cited_patent) %>%
    summarise(
      cd2 = mean(cd1),
      n_citations = n(),
      .groups = "drop"
    )

  # 4. 计算企业年度的CD
  cd_index <- cd2_values %>%
    left_join(patent_info, by = c("cited_patent" = "patent_id")) %>%
    group_by(firm_id, year) %>%
    summarise(
      CD = mean(cd2),
      n_patents = n(),
      .groups = "drop"
    )

  return(cd_index)
}

创建一个测试数据：

test_cites <- tibble(
  newipzlid = c("P2", "P3", "P4", "P5", "P6", "P7", "P8"),
  newipzlid2 = c("P1", "P1", "P1", "P2", "P2", "P3", "P3")
)

test_info <- tibble(
  patent_id = paste0("P", 1:8),
  firm_id = rep(c("F1", "F2"), each = 4),
  year = rep(2020:2021, 4)
)

calculate_cd_index_verbose(patent_citations = test_cites, patent_info = test_info)

#> # A tibble: 2 × 4
#>   firm_id  year    CD n_patents
#>   <chr>   <int> <dbl>     <int>
#> 1 F1       2020     1         2
#> 2 F1       2021     1         1

这样如果准备好类似上面结果的数据，就可以直接计算每个公司的 CD 指数了：

result <- calculate_cd_index_verbose(patent_citations, patent_info)
print(result)

#> # A tibble: 1,049 × 4
#>    firm_id  year    CD n_patents
#>    <chr>   <dbl> <dbl>     <int>
#>  1 000002   2010 1             3
#>  2 000012   2010 1             1
#>  3 000016   2010 1             7
#>  4 000021   2010 1             4
#>  5 000030   2010 1             1
#>  6 000032   2010 1             5
#>  7 000039   2010 0.6          10
#>  8 000050   2010 0.778         3
#>  9 000055   2010 1             2
#> 10 000063   2010 0.686       351
#> # ℹ 1,039 more rows

# 保存
result %>%
  rename(股票代码 = firm_id, 年份 = year, CD指数 = CD) %>%
  select(-n_patents) %>%
  writexl::write_xlsx("2010年上市公司CD指数计算结果.xlsx")

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言测算专利的创新突破度：CD指数（更新）