今天给大家分享使用 Python 计算专利颠覆性创新指数及筛选颠覆性专利的方法。

不过需要注意，直接对每个专利计算和之前年份（即使是 3 年内的）所有专利的相似度是非常非常耗时的，建议还是限定在某个领域内，例如只计算该专利与自己所在的部、大类、小类的。

使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。

重要说明（避免”已初始化”报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错：

ERROR: The requested version of Python cannot be used, as another version has already been initialized.

因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。

安装 reticulate（仅首次）

# 设置 CRAN 镜像（knit 时 R 处于非交互模式，不会自动选择镜像）
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

# 仅在尚未安装时才安装，避免每次 knit 都重装
if (!requireNamespace("reticulate", quietly = TRUE)) {
  install.packages("reticulate")
  message("reticulate 安装完成！")
} else {
  message("reticulate 已安装，版本：", packageVersion("reticulate"))
}

虚拟环境初始化原理（已在 setup chunk 中完成）

本文档的 setup chunk（隐藏运行）包含如下逻辑：

library(reticulate)

.venv_name   <- ".venv"
.venv_python <- virtualenv_python(.venv_name)

# 虚拟环境不存在时自动创建
if (!file.exists(.venv_python)) {
  virtualenv_create(.venv_name)
  .venv_python <- virtualenv_python(.venv_name)
}

# 通过环境变量抢先锁定 Python（优先级最高，早于任何 {python} chunk）
Sys.setenv(RETICULATE_PYTHON = .venv_python)
use_virtualenv(.venv_name, required = TRUE)

这样做的关键在于：knitr 在处理第一个 {python} chunk 时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv，不会再去碰 Anaconda。

在虚拟环境中安装 Python 包（仅首次）

# 检查关键包是否已安装，缺失的才安装
py_pkgs <- c("numpy", "pandas", "scipy", "jieba")
installed <- py_list_packages(".venv")$package
need_install <- setdiff(py_pkgs, installed)

if (length(need_install) > 0) {
  virtualenv_install(".venv", packages = need_install)
  message("已安装缺失的包：", paste(need_install, collapse = ", "))
} else {
  message("所有 Python 包已就绪，无需安装")
}

验证激活状态

# 验证当前绑定的 Python 路径（应指向 .venv 目录）
py_config()

1. 指标来源与计算方法

参照冉征、刘修岩、陈露（2025）发表于《中国工业经济》2025 年第 10 期的论文《技术集群结构与颠覆式创新——兼论关联性”陷阱”的突破路径》，本文基于 Kelly et al.（2021）设计的文本相似度法，通过专利文本的前后对比测度每项专利的颠覆性指标。

1.1 核心思想

颠覆式创新聚焦突破原有技术模式，形成全新的生产和价值创造模式，即一方面突破已有技术的限制，另一方面对后续的创新活动产生深远影响（Kerr, 2010）。基于这一内涵，一项专利如果：

与此前的专利关联性较小（新颖性高）
与此后的专利关联性较大（影响力大）

则说明该项专利突破了已有技术的框架，形成了新的技术模式，具有更强的颠覆性。

1.2 计算步骤

1.3 两种计算方案

1.4 颠覆性专利筛选

考虑到 radical 指数是连续指标且可以跨时期对比，按照论文方法设置前 5% 作为颠覆式创新的识别区间，即全样本中颠覆性指数位于前 5% 的专利定义为颠覆式创新专利。

2. 数据读取与预处理

2.1 导入 Python 包

import pandas as pd
import numpy as np
import re
import os
from collections import Counter
from scipy import sparse
import jieba

print("Python 包加载完成")

2.2 读取专利数据

首先读取专利数据样本文件夹中的所有 CSV 文件。专利数据包含如下变量：newipzlid（专利 ID，前 4 位为年份）、标题、摘要、公开公告号、申请号。

# 数据处理：微信公众号 RStata
data_dir = "专利数据样本"
csv_files = sorted([f for f in os.listdir(data_dir) if f.endswith(".csv")])

print("\n===== 读取专利数据 =====")
df_raw = pd.concat(
    [pd.read_csv(os.path.join(data_dir, f)) for f in csv_files],
    ignore_index=True
)

print(f"原始数据总行数: {len(df_raw)}")

2.3 提取年份

从 newipzlid 的前 4 位提取专利申请年份：

# 数据处理：微信公众号 RStata
df = df_raw.copy()
df["年份"] = df["newipzlid"].astype(str).str[:4].astype(int)

print(f"年份范围: {df['年份'].min()} - {df['年份'].max()}")

2.4 专利去重

专利数据中可能存在重复专利，需要按年份和申请号进行去重。先去重公开公告号（去掉末尾的类型标识字母），再去重申请号：

# 数据处理：微信公众号 RStata
print("\n===== 专利去重 =====")
print(f"去重前: {len(df)} 行")

# 按 年份 + 公开公告号 去重（去掉末尾字母）
df["公开公告号_clean"] = df["公开公告号"].str.replace(r"[A-Z]$", "", regex=True)
df = df.drop_duplicates(subset=["年份", "公开公告号_clean"], keep="first")
df = df.drop(columns=["公开公告号_clean"])

# 按 年份 + 申请号 去重（保留第一次出现的记录）
df = df.drop_duplicates(subset=["年份", "申请号"], keep="first")

print(f"去重后: {len(df)} 行")

关键函数对照：

R 函数	Python 函数	说明
`str_replace(x, "[A-Z]$", "")`	`str.replace(r"[A-Z]$", "", regex=True)`	去掉末尾大写字母
`distinct(年份, 列, .keep_all=TRUE)`	`drop_duplicates(subset=[...], keep="first")`	去重保留首条
`select(-列)`	`drop(columns=[...])`	删除列

2.5 文本清洗

合并专利的标题和摘要作为完整文本，去除标点、数字、英文及多余空格：

# 数据处理：微信公众号 RStata
# 此处代码需下载讲义材料查看～

2.6 中文分词

使用 jieba 进行中文分词，加载停用词表和用户词典，过滤长度小于 2 的词和非中文字符：

# 数据处理：微信公众号 RStata
# 此处代码需下载讲义材料查看～

分词逻辑对照：

R (jiebaR)	Python (jieba)	说明
`worker(stop_word="stopwords.txt", user="dictionary.txt")`	`jieba.load_userdict()` + 手动停用词过滤	加载词典和停用词
`segment(txt, seg)`	`jieba.cut(txt)`	中文分词
`nchar(tokens) >= 2`	`len(t) >= 2`	过滤短词
`!str_detect(tokens, "[^\\u4e00-\\u9fff]")`	`re.match(r"^[\u4e00-\u9fff]+$", t)`	仅保留纯中文

3. 辅助函数：TF-IDF 矩阵与颠覆性指数

3.1 构建时变 TF-IDF 矩阵

# 数据处理：微信公众号 RStata
# 此处代码需下载讲义材料查看～

TF-IDF 计算对照：

R (tidytext)	Python	说明
`bind_tf_idf(word, newipzlid, n)`	手动计算 tf, idf	TF-IDF 权重
`cast_sparse(newipzlid, word, tf_idf)`	`sparse.csr_matrix(...)`	构建稀疏 DTM
`dtm / row_norms`	`sparse.diags(1/norms) @ dtm`	L2 行归一化

3.2 计算颠覆性指数

对每一年，确定前后窗口内的专利索引，计算后向相似度（BPS）和前向相似度（FPS），并同时输出两种方案的 radical 指数。其中求和法是论文中使用的：

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

相似度计算对照：

R (Matrix)	Python (scipy.sparse)	说明
`dtm[target_idx, ] %*% t(dtm[backward_idx, ])`	`dtm[target_idx, :] @ dtm[backward_idx, :].T`	矩阵乘法求余弦相似度
`rowSums(sim_mat)`	`np.asarray(sim.sum(axis=1)).flatten()`	行求和
`ifelse(BPS > 0, FPS/BPS, NA)`	`np.where(BPS > 0, FPS/BPS, np.nan)`	条件计算

4. 逐年计算颠覆性指数

# 数据处理：微信公众号 RStata
target_years = list(range(2003, 2008))

print("\n===== 计算颠覆性创新指数 =====")
print(f"目标年份: {', '.join(map(str, target_years))}")
print(f"计算窗口: tau = 3 年")

all_results = []

for yr in target_years:
    print(f"\n--- 正在计算 {yr} 年 ---")
    result = compute_radical(df, yr, tau=3)
    print(f"  专利数: {len(result)}")
    print(f"  有效 radical_sum 指数: {result['radical_sum'].notna().sum()}")
    print(f"  有效 radical_avg 指数: {result['radical_avg'].notna().sum()}")
    all_results.append(result)

#>
#> --- 正在计算 2003 年 ---
#>   专利数: 998
#>   有效 radical_sum 指数: 980
#>   有效 radical_avg 指数: 980
#>
#> --- 正在计算 2004 年 ---
#>   专利数: 995
#>   有效 radical_sum 指数: 974
#>   有效 radical_avg 指数: 974
#>
#> --- 正在计算 2005 年 ---
#>   专利数: 997
#>   有效 radical_sum 指数: 972
#>   有效 radical_avg 指数: 972
#>
#> --- 正在计算 2006 年 ---
#>   专利数: 997
#>   有效 radical_sum 指数: 979
#>   有效 radical_avg 指数: 979
#>
#> --- 正在计算 2007 年 ---
#>   专利数: 998
#>   有效 radical_sum 指数: 984
#>   有效 radical_avg 指数: 984

5. 汇总结果与筛选颠覆性专利

5.1 标记颠覆性专利

在全样本期内，将颠覆性指数位于前 5% 的专利定义为颠覆式创新，分别对两种方案计算阈值并标记：

# 数据处理：微信公众号 RStata
# 此处代码需下载讲义材料查看～

5.2 合并原始专利信息

# 数据处理：微信公众号 RStata
# 注意：统一 newipzlid 类型为字符型，避免 merge 时类型不匹配
df["newipzlid"] = df["newipzlid"].astype(str)
df_radical["newipzlid"] = df_radical["newipzlid"].astype(str)

df_final = df[["newipzlid", "年份", "标题", "摘要", "公开公告号", "申请号"]].merge(
    df_radical, on=["newipzlid", "年份"], how="left"
)

6. 结果展示

6.1 总体描述

# 数据处理：微信公众号 RStata
print("\n===== 结果摘要 =====")
print(f"目标年份: {min(target_years)} - {max(target_years)}")
print(f"总专利数（去重后，含窗口期）: {len(df)}")
print(f"计算了 radical 指数的专利数: {df_final['radical_avg'].notna().sum()}")

6.2 各年颠覆性专利分布（方案一：求和法）

# 数据处理：微信公众号 RStata
# 此处代码需下载讲义材料查看～

6.3 各年颠覆性专利分布（方案二：均值法）

# 数据处理：微信公众号 RStata
print("\n--- 各年颠覆性专利分布（方案二：均值法）---")
summary_avg = df_final[df_final["年份"].isin(target_years)].groupby("年份").agg(
    总专利数=("newipzlid", "count"),
    有效radical=("radical_avg", lambda x: x.notna().sum()),
    颠覆性专利数=("颠覆性_avg", "sum"),
    radical均值=("radical_avg", "mean"),
    radical中位数=("radical_avg", "median"),
).reset_index()
summary_avg["颠覆性占比"] = (summary_avg["颠覆性专利数"] / summary_avg["总专利数"] * 100).round(2)
summary_avg["radical均值"] = summary_avg["radical均值"].round(4)
summary_avg["radical中位数"] = summary_avg["radical中位数"].round(4)
print(summary_avg.to_string(index=False))

6.4 两种方案一致性比较

# 数据处理：微信公众号 RStata
print("\n--- 两种方案一致性比较 ---")
df_target = df_final[df_final["年份"].isin(target_years)]
both = ((df_target["颠覆性_sum"] == 1) & (df_target["颠覆性_avg"] == 1)).sum()
only_sum = ((df_target["颠覆性_sum"] == 1) & (df_target["颠覆性_avg"] == 0)).sum()
only_avg = ((df_target["颠覆性_sum"] == 0) & (df_target["颠覆性_avg"] == 1)).sum()
print(f"  两种方案均标记为颠覆性: {both}")
print(f"  仅求和法标记: {only_sum}")
print(f"  仅均值法标记: {only_avg}")
print(f"  一致率: {both / (both + only_sum + only_avg) * 100:.1f}%")

#> 一致率: 99.2%

7. 保存结果

7.1 保存为 CSV 文件

# 数据处理：微信公众号 RStata
df_output = df_final[[
    "newipzlid", "年份", "标题", "摘要", "公开公告号", "申请号",
    "BPS_sum", "FPS_sum", "radical_sum", "颠覆性_sum",
    "BPS_avg", "FPS_avg", "radical_avg", "颠覆性_avg",
]].copy()

csv_file = "patent_disruptive_index.csv"
df_output.to_csv(csv_file, index=False, encoding="utf-8-sig")
print(f"结果已保存到: {csv_file}")

7.2 保存为 Stata .dta 文件

由于 Stata 变量名不支持中文字符，需将中文列名映射为英文变量名，并设置变量标签。最终通过 R 的 haven 包写入 .dta 文件，以保证与 R 版本的完全兼容：

# 数据处理：微信公众号 RStata
# Stata 变量名不支持中文，需映射为英文变量名
rename_map = {
    "newipzlid": "newipzlid",
    "年份": "year",
    "标题": "title",
    "摘要": "abstract",
    "公开公告号": "pub_id",
    "申请号": "app_id",
    "BPS_sum": "bps_sum",
    "FPS_sum": "fps_sum",
    "radical_sum": "radical_s",
    "颠覆性_sum": "disrup_s",
    "BPS_avg": "bps_avg",
    "FPS_avg": "fps_avg",
    "radical_avg": "radical_a",
    "颠覆性_avg": "disrup_a",
}

df_stata = df_output.rename(columns=rename_map).copy()

# 确保数值列类型正确
for col in ["bps_sum", "fps_sum", "radical_s", "bps_avg", "fps_avg", "radical_a"]:
    df_stata[col] = pd.to_numeric(df_stata[col], errors="coerce")
# 颠覆性标记：NaN 填充为 0 后转为整数
for col in ["disrup_s", "disrup_a"]:
    df_stata[col] = df_stata[col].fillna(0).astype(int)

# 保存为 pickle 供 R 读取
df_stata.to_pickle("df_stata_tmp.pkl")
print("数据准备完成，即将通过 R 写入 .dta 文件")

# 数据处理：微信公众号 RStata
library(haven)

# 从 Python 环境读取准备好的数据
df_stata <- py$df_stata

# 变量标签（中文）
variable_labels <- list(
  newipzlid = "专利ID",
  year      = "年份",
  title     = "标题",
  abstract  = "摘要",
  pub_id    = "公开公告号",
  app_id    = "申请号",
  bps_sum   = "后向相似度总和",
  fps_sum   = "前向相似度总和",
  radical_s = "颠覆性指数(求和法)",
  disrup_s  = "颠覆性标记(求和法)",
  bps_avg   = "后向平均相似度",
  fps_avg   = "前向平均相似度",
  radical_a = "颠覆性指数(均值法)",
  disrup_a  = "颠覆性标记(均值法)"
)

# 设置变量标签
for (v in names(variable_labels)) {
  attr(df_stata[[v]], "label") <- variable_labels[[v]]
}

# 写入全量 .dta 文件
write_dta(df_stata, "patent_disruptive_index.dta",
          label = "数据处理:微信公众号 RStata")
cat("Stata 文件已保存至: patent_disruptive_index.dta\n")

# 写入 2003-2007 年子集
df_target <- df_stata[df_stata$year %in% 2003:2007, ]

write_dta(df_target, "patent_disruptive_index_2003_2007.dta",
          label = "数据处理:微信公众号 RStata")
cat("2003-2007 年子集已保存至: patent_disruptive_index_2003_2007.dta\n")

# 清理临时文件
file.remove("df_stata_tmp.pkl")

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 计算专利颠覆性创新指数及筛选颠覆性专利