名师讲堂|Stata:如何根据专利申请数据构造城市间合作申请数量矩阵?

《经济评论》2024年第3期论文「高铁网络如何促进城市间合作创新——基于高铁网络通达性与合作专利的实证分析」中介绍了这样的一个指标:

论文全文:https://mp.weixin.qq.com/s/jDl_v0-gCHWgoTUaI86BCA

今天我们一起来学习下如何在 Stata 中高效的计算这一指标,也就是城市间合作申请专利数量。

前不久我们分享的「1985~2024年共同申请专利的所有申请人经纬度及其所属的省市区县」数据正好可以用来计算这一指标:

1985~2024年共同申请专利的所有申请人经纬度及其所属的省市区县:https://rstata.duanshu.com/#/brief/course/6ce40daa5ebe4164af66b07f68374c49

该数据是对所有申请人为企业的进行地理编码,由于直接对企业名称地理编码会错误率很高,所以先匹配工商注册信息,使用工商注册地址解析,然后匹配不到工商注册信息的使用企业名称解析获得经纬度,有了经纬度之后就可以根据经纬度判断每个企业所处的省市区县了。对于实在无法解析的企业,尽可能使用其名称中包含的城市信息。

在附件中我选择了 2015~2017 年的数据作为示例:

2015~2017年共同申请专利的所有申请人经纬度及其所属的省市区县.dta

首先读取该数据,去除重复的专利:

cd "~/Desktop/Stata:如何根据专利申请数据构造城市间合作申请数量矩阵?"
*> /Users/ac/Desktop/Stata:如何根据专利申请数据构造城市间合作申请数量矩阵?
use "2015~2017年共同申请专利的所有申请人经纬度及其所属的省市区县.dta", clear
*> (数据处理:微信公众号 RStata)
*- 删除重复专利
replace 公开公告号 = subinstr(公开公告号, "A", "", .)
*> (444,049 real changes made)
replace 公开公告号 = subinstr(公开公告号, "B", "", .)
*> (245,323 real changes made)
replace 公开公告号 = subinstr(公开公告号, "U", "", .)
*> (296,944 real changes made)
replace 公开公告号 = subinstr(公开公告号, "S", "", .)
*> (36,974 real changes made)
duplicates drop 公开公告号, force
*> Duplicates in terms of 公开公告号
*> (674,896 observations deleted)
replace 专利类型 = "发明" if index(专利类型, "发明")
*> (202,118 real changes made)
keep newipzlid
save unique_patent, replace
*> file unique_patent.dta saved
list in 1/10
*> +-------------+
*> | newipzlid |
*> |-------------|
*> 1. | 20150000221 |
*> 2. | 20150000222 |
*> 3. | 20150000223 |
*> 4. | 20150000224 |
*> 5. | 20150000225 |
*> |-------------|
*> 6. | 20150001365 |
*> 7. | 20150001366 |
*> 8. | 20150001385 |
*> 9. | 20150001389 |
*> 10. | 20150001420 |
*> +-------------+

这是因为该专利数据中同时包含了专利的引用与授权公告,直接删除授权公告又担心有些授权专利的申请公告没有包含在数据里面:

统计的时候可以先去除公开公告号里面的 A、B、U、S。其中 A 代表发明专利的申请公开,B 代表发明专利的授权公告,U 代表实用新型专利的授权公告,S 代表外观设计专利的授权公告。

上面获取的是所有互不重复的专利编号,再和原始数据匹配,保留匹配成功的,就得到了所有互不重复的专利申请人信息:

use "2015~2017年共同申请专利的所有申请人经纬度及其所属的省市区县.dta", clear
*> (数据处理:微信公众号 RStata)
merge m:1 newipzlid using unique_patent
*> Result Number of obs
*> -----------------------------------------
*> Not matched 235,927
*> from master 235,927 (_merge==1)
*> from using 0 (_merge==2)
*> Matched 787,363 (_merge==3)
*> -----------------------------------------
keep if _m == 3
*> (235,927 observations deleted)
drop _m
replace 专利类型 = "发明" if index(专利类型, "发明")
*> (457,348 real changes made)

然后按照论文中介绍的方法进行一些筛选:

*- 删除申请人数量超过 10 个的
bysort newipzlid: gen count = _N
drop if count < 2 | count > 10
*> (922 observations deleted)
*- 有些只剩一个的是因为另外一个没解析成功
tab count
*> count | Freq. Percent Cum.
*> ------------+-----------------------------------
*> 2 | 554,622 70.52 70.52
*> 3 | 162,960 20.72 91.24
*> 4 | 47,452 6.03 97.28
*> 5 | 14,035 1.78 99.06
*> 6 | 5,424 0.69 99.75
*> 7 | 1,498 0.19 99.94
*> 8 | 280 0.04 99.98
*> 9 | 90 0.01 99.99
*> 10 | 80 0.01 100.00
*> ------------+-----------------------------------
*> Total | 786,441 100.00
drop count
save tempdata, replace
*> file tempdata.dta saved
list in 1/10
*> +--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
*> | 年份 newipzlid 申请~号 企业名称 经度 纬度 省 省代码 市 市代码 县 县代码 公开公告号 专利类型 IPC 授权公告日 |
*> |--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
*> 1. | 2015 20150000221 申请人1 北京比动商贸有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280535S 外观设计 2015-07-15 |
*> 2. | 2015 20150000221 申请人2 北京比动广告有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280535S 外观设计 2015-07-15 |
*> 3. | 2015 20150000222 申请人1 北京比动商贸有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280536S 外观设计 2015-07-15 |
*> 4. | 2015 20150000222 申请人2 北京比动广告有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280536S 外观设计 2015-07-15 |
*> 5. | 2015 20150000223 申请人1 北京比动商贸有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280537S 外观设计 2015-07-15 |
*> |--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
*> 6. | 2015 20150000223 申请人2 北京比动广告有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280537S 外观设计 2015-07-15 |
*> 7. | 2015 20150000224 申请人1 北京比动商贸有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280538S 外观设计 2015-07-15 |
*> 8. | 2015 20150000224 申请人2 北京比动广告有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280538S 外观设计 2015-07-15 |
*> 9. | 2015 20150000225 申请人1 北京比动商贸有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280539S 外观设计 2015-07-15 |
*> 10. | 2015 20150000225 申请人2 北京比动广告有限公司 116.40105 39.914867 北京市 110000 北京市 110000 东城区 110101 CN303280539S 外观设计 2015-07-15 |
*> +--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

在循环计算多年的数据前,我们先选择一个年份测试。以 2015 年为例:

use tempdata, clear
*> (数据处理:微信公众号 RStata)
keep if 年份 == 2015
*> (558,651 observations deleted)
duplicates drop newipzlid 申请人编号, force
*> Duplicates in terms of newipzlid 申请人编号
*> (180 observations deleted)
keep newipzlid 申请人编号 市
list in 1/10
*> +--------------------------------+
*> | newipzlid 申请~号 市 |
*> |--------------------------------|
*> 1. | 20150000221 申请人1 北京市 |
*> 2. | 20150000221 申请人2 北京市 |
*> 3. | 20150000222 申请人1 北京市 |
*> 4. | 20150000222 申请人2 北京市 |
*> 5. | 20150000223 申请人1 北京市 |
*> |--------------------------------|
*> 6. | 20150000223 申请人2 北京市 |
*> 7. | 20150000224 申请人1 北京市 |
*> 8. | 20150000224 申请人2 北京市 |
*> 9. | 20150000225 申请人1 北京市 |
*> 10. | 20150000225 申请人2 北京市 |
*> +--------------------------------+

在之前的课程「名师讲堂|使用 Stata 测算数实融合水平」中我们也介绍了类似的操作,也就是计算 IPC 融合矩阵。因此我们这里也仿照那个课程里面的样子,把数据整理成下面这样:

*- 此处代码需下载讲义材料查看~

汇总下:

contract city

由于不需要考虑自引,所以我们可以去除重复城市:

*- 去除重复的城市
egen city2 = ipc_unique(city), parse(;)
drop city
ren city2 city
replace city = subinstr(city, ";", " ", .)
*> (6,597 real changes made)
replace city = strtrim(city)
*> (7 real changes made)
replace city = subinstr(city, " ", " ", .)
*> (0 real changes made)
replace city = subinstr(city, " ", ";", .)
*> (6,595 real changes made)

这里的 ipc_unique() 函数来自之前的课程:「名师讲堂|使用 Stata 处理专利数据的分类号」

名师讲堂|使用 Stata 处理专利数据的分类号: https://rstata.duanshu.com/#/brief/course/1cfb2e8e8e5f4716959fc78bbaf6d446

对编写过程感兴趣的小伙伴可以回去再看看。

没有互相合作的也是我们不关心的,直接删除:

*- 删除没有互相合作的
drop if !index(city, ";")
*> (321 observations deleted)
order city
save tempdata2, replace
*> file tempdata2.dta saved

然后我们就可以继续借鉴「名师讲堂|使用 Stata 测算数实融合水平」中的代码统计城市共现数量了:

*- 统计各城市对的合作数量
use tempdata2, clear

mata:
city = st_sdata(., "city")
value = st_data(., "_freq")
fullcombinations = J(0, 3, "")
for (r = 1; r <= rows(city); r++) {
class_list = colshape(ustrsplit(city[r, 1], ";"), 1)
n_classes = rows(class_list)
if (n_classes == 2) {
fullcombinations = fullcombinations (class_list[1,1], class_list[2,1], strofreal(value[r,1]))
}
if (n_classes > 2) {
combinations = J(n_classes * n_classes, 3, "")
row_index = 1

for (i = 1; i <= n_classes; i++) {
for (j = 1; j <= n_classes; j++) {
combinations[row_index, 1] = class_list[i,1]
combinations[row_index, 2] = class_list[j,1]
combinations[row_index, 3] = strofreal(value[r,1])
row_index++
}
}
fullcombinations = fullcombinations combinations
}
}

st_matrix("obs", rows(fullcombinations))
stata("clear")
st_addvar("str100", "city1")
st_addvar("str100", "city2")
st_addvar("str100", "value")
stata("set obs `=obs[1,1]'")
st_sstore(., ("city1", "city2", "value"), fullcombinations)
end

compress
foreach i of varlist _all {
cap format `i' %10s
}
destring value, replace
drop if city1 == city2
collapse (sum) value, by(city1 city2)
save tempdata3, replace

list in 1/10

然而稍微仔细检查就会发现,这个数据里面 A-B 的数量和 B-A 的数量不完全一样,也就是说这个数据要是展开成矩阵不是对称的。

这里之所以有不对称的问题,其实主要是因为,如果申请人有两个,上面的 mata 程序里面只保留了 A-B-Value 的结果,没有同时生成 B-A-Value 的结果。

可以通过下面的方式把数据转换成对称的:

*- 处理成对称矩阵
use tempdata3, clear
gen id = _n
drop value
gather city1 city2
*> id
keep value
duplicates drop value, force
*> Duplicates in terms of value
*> (10,148 observations deleted)
save temp1b, replace
*> file temp1b.dta saved
ren value value2
cross using temp1b
ren value city1
ren value2 city2
merge 1:1 city1 city2 using tempdata3
*> Result Number of obs
*> -----------------------------------------
*> Not matched 117,251
*> from master 117,251 (_merge==1)
*> from using 0 (_merge==2)
*> Matched 5,249 (_merge==3)
*> -----------------------------------------
drop _m
replace value = 0 if mi(value)
*> (117,251 real changes made)
save temp1c, replace
*> file temp1c.dta saved
ren city1 temp
ren city2 city1
ren temp city2
ren value value2
merge 1:1 city1 city2 using tempdata3
*> Result Number of obs
*> -----------------------------------------
*> Not matched 117,251
*> from master 117,251 (_merge==1)
*> from using 0 (_merge==2)
*> Matched 5,249 (_merge==3)
*> -----------------------------------------
*- 因为对于每个专利,ab ba 的情况都统计了,所以这里选择最大的(统计最全面的)作为共现结果
egen value3 = rowmax(value value2)
replace value = value3
*> (118,258 real changes made)
drop value2 _m value3
drop if city1 == city2
*> (350 observations deleted)
gsort city1 city2
drop if value == 0
*> (115,256 observations deleted)
ren value 合作申请专利数量
save 2015年各城市对合作申请专利数量计算结果, replace
*> file 2015年各城市对合作申请专利数量计算结果.dta saved

这个也是那个课程里面讲解的方法。不过其实有更好的办法,也就是直接纠正 Mata 程序:

use tempdata2, clear
mata:
city = st_sdata(., "city")
value = st_data(., "_freq")
fullcombinations = J(0, 3, "")
for (r = 1; r <= rows(city); r++) {
class_list = colshape(ustrsplit(city[r, 1], ";"), 1)
n_classes = rows(class_list)
if (n_classes == 2) {
fullcombinations = fullcombinations (class_list[1,1], class_list[2,1], strofreal(value[r,1]))
fullcombinations = fullcombinations (class_list[2,1], class_list[1,1], strofreal(value[r,1]))
}
if (n_classes > 2) {
combinations = J(n_classes * n_classes, 3, "")
row_index = 1

for (i = 1; i <= n_classes; i++) {
for (j = 1; j <= n_classes; j++) {
combinations[row_index, 1] = class_list[i,1]
combinations[row_index, 2] = class_list[j,1]
combinations[row_index, 3] = strofreal(value[r,1])
row_index++
}
}
fullcombinations = fullcombinations combinations
}
}

st_matrix("obs", rows(fullcombinations))
stata("clear")
st_addvar("str100", "city1")
st_addvar("str100", "city2")
st_addvar("str100", "value")
stata("set obs `=obs[1,1]'")
st_sstore(., ("city1", "city2", "value"), fullcombinations)
end

compress
foreach i of varlist _all {
cap format `i' %10s
}
destring value, replace
drop if city1 == city2
collapse (sum) value, by(city1 city2)
ren value 合作申请专利数量
save 2015年各城市对合作申请专利数量计算结果2, replace

然后我们就可以循环所有的年份了:

*- 此处代码需下载讲义材料查看~

需要注意由于 Mata 程序以 end 结尾,然而 end 也会结束循环,所以我们不能直接把 Mata 程序放到循环里面,这里我是新建了一个 matacode.do,然后在循环里面使用 do matacode.do 调用。

这里我给每年保存了一个文件,使用我自编的 appendall 命令可以合并某个文件夹下面的所有 dta 文件:

*- 合并所有的
appendall res
*> (6,894 observations deleted)
foreach i of varlist _all {
*> 2. cap format `i' %10s
*> 3. }
order 年份
gsort 年份
label data "数据计算:微信公众号 RStata"
save data1, replace
*> (file data1.dta not found)
*> file data1.dta saved

类似的代码统计分类型的专利:

*- 分类型
foreach v in "发明" "外观设计" "实用新型" {
cap mkdir "res`v'"
forval y = 2015/2017 {
di "`v': `y'"
qui {
use tempdata, clear
keep if 专利类型 == "`v'"
keep if 年份 == `y'
if `=_N' > 0 {
duplicates drop newipzlid 申请人编号, force
keep newipzlid 申请人编号 市
spread 申请人编号 市
unite 申请人*, gen(city) sep(" ")
drop 申请人*
replace city = strtrim(city)
replace city = subinstr(city, " ", " ", .)
replace city = subinstr(city, " ", ";", .)
contract city

*- 去除重复的城市
egen city2 = ipc_unique(city), parse(;)
drop city
ren city2 city
replace city = subinstr(city, ";", " ", .)
replace city = strtrim(city)
replace city = subinstr(city, " ", " ", .)
replace city = subinstr(city, " ", ";", .)

*- 删除没有互相合作的
drop if !index(city, ";")
order city

do matacode.do
destring value, replace
drop if city1 == city2
collapse (sum) value, by(city1 city2)
ren value 合作申请`v'专利数量
gen 年份 = `y'
save res`v'/`y', replace

*- 删除文件
cap erase tempdata3
cap erase tempdata2
cap erase temp1b
cap erase temp1c
}
}
}
}

最后合并这四个文件即可:

*- 此处代码需下载讲义材料查看~

这样我们就解决了这个问题,区县和省份的计算代码也非常类似。

后面我们还会讲解如何使用 Stata 绘制合作网络图,类似这样:

点击这里跳转到 RStata 短书平台获取附件:名师讲堂|Stata:如何根据专利申请数据构造城市间合作申请数量矩阵?

评论