世界各国新发呼吸道传染病的风险评估系统

技术领域

本发明涉及传染病风险评估技术领域，具体涉及一种世界各国新发呼吸道传染病的风险评估系统。

背景技术

新发呼吸道传染病具有潜伏期长、传播迅速、致死数量多等显著特点，给全球卫生防疫系统带来极大的威胁与挑战，同时也对各国经济发展造成灾难性打击。呼吸道传染病亦是多发常见疾病，如甲型流感、中东呼吸综合征、SARS等，如果防控不及时，造成的危害是十分巨大的。因此，建立世界各国新发呼吸道传染病风险评估系统，全面动态的识别各国家自身疫情暴发风险和对我国的境外输入风险，具有十分重要的防疫价值。

由于不同国家医疗水平、人口学特征、交通特征等差异较大，其自身疫情的暴发趋势、防控水平会有所差异，对我国的输入风险亦不尽相同。识别疫情高风险国家、高风险输入来源国家至关重要。因此，亟需设计一种新的技术方案，以综合解决现有技术中存在的问题。

发明内容

本发明的目的是提供一种世界各国新发呼吸道传染病的风险评估系统，能全面动态的识别各国家自身疫情暴发风险和对我国的境外输入风险。

为解决上述技术问题，本发明采用了以下技术方案：

一种世界各国新发呼吸道传染病的风险评估系统，包括通讯连接的用于数据采集的数据采集模块、用于数据分析的数据处理模块以及用于实时显示分析数据的显示模块；其中数据采集模块采集各国的人口学数据、全国医疗水平数据、经济发展数据、疫情数据和全球健康安全指数GHS，数据处理模块接收数据采集模块采集的数据，并对数据进行以下处理：

(1)首先评估各国新发呼吸道传染病疫情爆发风险，评估公式如下：

式中：I

I

I

(2)之后采用加权校正法评估各国家与中国的交通流量，评估公式如下：

I

式中：I

I

I

(3)最后按照下式评估世界各国新发呼吸道传染病对中国的潜在境外输入风险：

Risk

式中：Risk

显示模块实时显示数据处理模块的分析处理结果。

其中，若人口学数据、全国医疗水平数据、经济发展数据、疫情数据和全球健康安全指数GHS的采集指标中包含多个数据信息，则使用多元统计方法主成分分析进行降维，提取第一主成分作为评价指标。

其中，人口学数据包括各国人口密度、人口数、华侨华人常住人口数、在外中国留学生数、海外务工人员数、来华留学生数以及航班旅客数。

其中，疫情数据包括每日累计确诊病例数、死亡数以及境外输入累计确诊病例数。

上述技术方案中提供的世界各国新发呼吸道传染病的风险评估系统，利用各国家人口密度、经济发展水平、医疗投入占比、医疗卫生评级、传染病全球病例数等信息基于传染病动力学模型预测各国家新发呼吸道传染病疫情暴发风险，基于上述风险，结合来华航班乘客数估计、全球健康安全指数、各国华人华侨、出国留学与来华留学生等大数据，实时动态地预测未来每个国家的境外输入风险，本发明能够直观各国家新发传染病疫情风险、来华乘客数动态估计、未来输入风险值预测等信息，为相关防疫部门提供重要的参考价值。

本发明的罹患率指数借助于全球健康安全指数，由于各国家经济卫生水平不尽相同，疫情暴发时计算出的罹患率仅能代表各国家的检出确诊人数比例，对于发展水平较为落后的国家可能存在大量漏报现象，对于全球健康安全指数较低的国家，漏诊数亦较多，因此本发明基于GHS指数对罹患率进行校正，罹患率校正后与GHS指数无关联，可排除混杂因素干扰。

具体实施方式

为了使本发明的目的及优点更加清楚明白，以下结合实施例对本发明进行具体说明。应当理解，以下文字仅仅用以描述本发明的一种或几种具体的实施方式，并不对本发明具体请求的保护范围进行严格限定。

本发明采取的技术方案如下，一种世界各国新发呼吸道传染病的风险评估系统，包括通讯连接的用于数据采集的数据采集模块、用于数据分析的数据处理模块以及用于实时显示分析数据的显示模块；其中数据采集模块采集各国的人口学数据、全国医疗水平数据、经济发展数据、疫情数据和全球健康安全指数GHS，数据处理模块接收数据采集模块采集的数据，并对数据进行以下处理：

(1)首先评估各国新发呼吸道传染病疫情爆发风险，评估公式如下：

式中：I

I

I

(2)之后采用加权校正法评估各国家与中国的交通流量，评估公式如下：

I

式中：I

I

I

(3)最后按照下式评估世界各国新发呼吸道传染病对中国的潜在境外输入风险：

Risk

式中：Risk

显示模块实时显示数据处理模块的分析处理结果。

其中，若人口学数据、全国医疗水平数据、经济发展数据、疫情数据和全球健康安全指数GHS的采集指标中包含多个数据信息，则使用多元统计方法主成分分析进行降维，提取第一主成分作为评价指标。

本实施例中人口学数据中的各国人口数(万)采集自快易数据“世界各国人口总数统计”，世界各国华侨华人常住人口数收集自《华人经济年鉴(2009～2010)》；在外中国留学生数收集自《中国留学发展报告(2016)》；海外务工人员数收集自《中国统计年鉴2014》：“按国别(地区)分对外经济合作”，务工人员估计值为承包工程年末在外人数与劳务合作年末在外人数之和；各国家来华留学生数据收集自《来华留学生简明统计2014》。

全球医疗水平根据世界各国的医疗可及性和质量指数(Healthcare Access andQuality Index，简称HAQ指数)。各国指数分布在0到100之间，分值越高，说明该国个人获得的医疗服务机会和质量越高，数据收集自权威医学杂志《柳叶刀》在2018年发布了《2016年全球疾病、伤害和风险因素负担研究报告》(Global Burden of Diseases,Injuries,andRisk Factors Study 2016,GBD)，对195个国家和地区在1990到2016年间个人医疗服务的获取情况和质量进行了评价。

疫情数据包括(各国家(地区)每日累计确诊病例数、死亡数等信息。每日累计确诊病例数收集自霍普金斯大学“全球病例数报告”(https://coronavirus.jhu.edu/map.html)；境外输入病例累计确诊病例数收集自各省份卫生健康委员会官方网站。所有疫情相关数据收集每天实时更新。

全球健康安全指数(GHS)(2019)是由约翰斯·霍普金斯大学联合核威胁倡议和经济学人智库发布的全球各国家安全卫生水平指数，基于疾病预防、传染病检测和报告、快速响应等六个维度计算出的综合评分，介于0-100分之间，是目前能够较为全面准确衡量各国家公共卫生水平的指数。由于各国家经济卫生水平不尽相同，疫情暴发时计算出的罹患率仅能代表各国家的检出确诊人数比例，发展水平较为落后的国家可能存在大量漏报现象。因此，其很难用来衡量真实罹患率水平。本发明借助于全球健康安全指数(GHS)，对于GHS较低的国家，漏诊数亦较多。因此，基于GHS指数对罹患率进行校正：罹患率指数＝e

其中，如果有的指标包含多个数据的信息，如经济发展实力包括了GDP总量、人均GDP、恩格尔系数等，则使用多元统计方法主成分分析(PCA)进行降维，提取第一主成分作为其评价指标。

世界各国与中国交通流量估计主要基于两方面的指标：1、国际航班表，2、航班机型信息，使用网络爬虫技术，获取数据。

航班信息估计数据收集自每日计划的国际航班信息，数据收集自携程旅行“国际航班查询时刻表”(https://flights.ctrip.com/international/search/schedule/)，航班起飞或取消状态由携程旅行航班实时动态数据库进一步确认(https://flights.ctrip.com/actualtime/)。航班信息包括了境外各国家飞至全国内地各机场的起飞时间、落地时间、起飞机场、落地机场、飞机机型、起飞城市代码、终点城市代码等信息。所有航班相关数据收集每天实时更新。

此外，由于国际航班转机情况较为常见，亦需对航班入境人数进行校正。本发明提出一种加权校正法，基于输入来源国家的各机场与全球各机场的通航流量作为转机概率的权重，并把转机流量的总人数控制在30％左右，最终得到校正的航班入境流量。

为了估计不同航班的入境乘客数，本专利使用机型在常见舱位规格下搭载乘客数的近似值。对于常见机型的载客量估计值分别为：空客319(123)，空客320(152)，空客321(197)，空客330(334)，空客340(440)，空客350(323)，空客380(555)，波音737(147)，波音738(176)，波音747(407)，波音767(245)，波音777(361)，波音787(289)。

本专利中相关性分析使用Spearman秩相关进行分析，多重比较采用Bonferroni法进行校正。线性回归基于最小二乘法进行估计。对于非正态分布的数据(罹患率等)，模型中进行对数变换。本专利所有数据收集、数据分析均采用R软件(Version 3.6.1)完成。

部分参考计算机源程序如下：

setwd("D:\\OneDrive\\2020\\0831")

options(stringsAsFactors＝F)

#source("sourcecode.R")

require(readxl)

temp＝read.csv("time_series_confirmed_global.csv")

ic_latest＝aggregate(temp[,ncol(temp)],by＝list(temp$Country.Region),FUN＝sum)#最新确诊病例

ic_10dayago＝aggregate(temp[,(ncol(temp)-9)],by＝list(temp$Country.Region),FUN＝sum) #10日前确诊病例

colnames(ic_latest)＝c("country_eng","ic_latest")

ic_latest$ic_10dayago＝ic_10dayago$x

dict＝read_excel("Dict.xlsx")#中英文匹配

country_full＝intersect(dict$eng1,ic_latest$country_eng)

country_full_chinese＝dict$chinese[match(country_full,dict$eng1)]

ic_latest＝ic_latest[match(country_full,ic_latest$country_eng),]

risk_assess＝data.frame(country_full＝country_full_chinese,country_full_eng＝country_full,

ic_now＝ic_latest$ic_latest,ic_10dayago＝ ic_latest$ic_10dayago)

risk_assess$increase＝ (risk_assess$ic_now-risk_assess$ic_10dayago)/risk_assess$ic_10dayago#计算10日增长人数

#########整理数据

popul＝read_excel("世界人口.xlsx")

popul$eng＝dict$eng[match(popul$名称,dict$chinese)]

risk_assess$popul＝popul$`总人口(万)`[match(risk_assess$country_full,popul$名称)]

risk_assess$Prevalence＝risk_assess$ic_now/(risk_assess$popul*100)

dens＝read.csv("世界人口密度.csv")

risk_assess$dens＝dens$人口.平方公里[match(risk_assess$country_full,dens$国家.地区)]

###huaren

huaren＝read_excel("世界华人华侨统计.xlsx",col_types＝c("text","numeric"))

huaren$chinese＝unlist(lapply(huaren$Full_name,function(x)strsplit(x,"(",fixed＝T)[[1]][1]))

huaren$eng＝unlist(lapply(huaren$Full_name,function(x)strsplit(x,"(",fixed＝T)[[1]][2]))

huaren$eng＝substr(huaren$eng,1,nchar(huaren$eng)-1)

risk_assess$huaren＝huaren$Num[match(risk_assess$country_full,huaren$chinese)]

risk_assess$huaren[is.na(risk_assess$huaren)]＝0

###student

student＝read_excel("留学生数.xlsx")

risk_assess$student＝student$Num[match(risk_assess$country_full,student$Country)]*10000

risk_assess$student[is.na(risk_assess$student)]＝0

###worker

worker＝read_excel("中国统计年鉴2014按国别(地区)分对外经济合作(2013年).xls",sheet＝"clean")

risk_assess$worker＝worker$Num[match(risk_assess$country_full,worker$Country)]

risk_assess$worker[is.na(risk_assess$worker)]＝0

###outstudent

outstudent＝read_excel("来华留学生简明统计2014.xlsx")

risk_assess$outstudent＝ outstudent$Num[match(risk_assess$country_full,outstudent$Country)]

risk_assess$outstudent[is.na(risk_assess$outstudent)]＝0

###GHS

ghs<-read.csv("GHS_index.csv")

risk_assess$ghs＝ghs$ghs[match(dict$eng,ghs$country)]

risk_assess$realpre＝exp(log(risk_assess$Prevalence)*(risk_assess$ghs/100))##校正罹患率指数

###region大洲信息

region＝read_excel("Region.xlsx")

risk_assess$region_chinese＝ region$region_chinese[match(risk_assess$country_full_eng,region$country)]

###passenger旅客数估计

citylist＝read_excel("城市缩写.xlsx")

plane＝read_excel("plane.xlsx")

pas＝read.csv("9.02-9.03航班数据(国内21市+国外119市).csv")

pas$Depart_country＝citylist$Country[match(pas$Depart_city,citylist$City)]

pas$Destin_country＝citylist$Country[match(pas$Destin_city,citylist$City)]

pas1＝subset(pas,Destin_country％in％c("中国","中国香港")&！Depart_country％in％c(" 中国","中国香港"))

pas$passenger＝plane$passenger[match(pas$plane,plane$plane)]

pas$passenger[is.na(pas$passenger)]＝mean(pas$passenger,na.rm＝T)

passenger＝data.frame(country＝names(table(pas$Depart_country)))；passenger$num＝NA

for(i in 1:nrow(passenger))passenger$num[i]＝ sum(pas$passenger[pas$Depart_country＝＝passenger$country[i]])

risk_assess$passenger＝passenger$num[match(risk_assess$country_full,passenger$country)]

risk_assess$passenger[is.na(risk_assess$passenger)]＝0

risk_assess$passenger1＝risk_assess$passenger+sum(risk_assess$passenger)*0.05

#risk_assess$risk＝risk_assess$passenger*risk_assess$realpre

model＝lm(ic_now～Prevalence+dens+huaren+student+outstudent+worker,data＝ risk_assess)；summary(model)

####风险评分

risk_pca＝ risk_assess[,c("ic_now","increase","realpre","huaren","student","worker","outstudent")]

#colnames(risk_pca)＝c("境外输入","罹患率","华侨华人","留学生数")

rownames(risk_pca)＝risk_assess$country_full

risk_pca＝risk_pca[complete.cases(risk_pca),]

p＝princomp(log(risk_pca+1))；summary(p)

m＝p$scores

m；summary(p)

m＝as.data.frame(m)

m$region＝risk_assess$region_chinese[match(rownames(m),risk_assess$country_full)]

require(ggplot2)

ggplot()+

geom_point(aes(x＝Comp.1,y＝Comp.2,pch＝factor(region,levels＝c("西太平洋","美洲 ","欧洲","非洲","东地中海","东南亚")),

col＝factor(region,levels＝c("西太平洋","美洲","欧洲","非洲","东地中海","东南亚"))),data＝m,size＝3,position＝"jitter")+

geom_text(aes(x＝Comp.1,y＝Comp.2+0.5,label＝rownames(risk_pca)),data＝m,size＝ 3)+

theme(legend.title＝element_blank())+

xlab("PC1")+

ylab("PC2")+

theme(axis.text.y＝element_text(size＝15,color＝"black",angle＝0))+

theme(axis.text.x＝element_text(size＝15,color＝"black",angle＝45,hjust＝0.9))+

theme(axis.title.x＝element_text(size＝15))+

theme(axis.title.y＝element_text(size＝15))+

theme(panel.grid.major＝element_blank(),panel.grid.minor＝element_blank(),

panel.background＝element_blank(),axis.line＝element_line(colour＝"black"))

ggsave(paste0("./plot/PCA.jpg"),height＝8,width＝8,dpi＝300)

model_score＝lm(risk_assess$ic_now[match(rownames(m),risk_assess$country_full)]～m$Comp.1+m$Comp.2)；p rint(summary(model_score))

m$score＝coef(model_score)[2]*m$Comp.1

m$score＝m$score-min(m$score)

m$score＝m$score*100/max(m$score)

risk_assess$pca_score＝NA

risk_assess$pca_score[match(rownames(m),risk_assess$country_full)]＝m$score

risk_assess$risk＝log(risk_assess$passenger1+1)*risk_assess$pca_score

risk_assess$risk_norm＝risk_assess$risk

risk_assess$risk_norm[is.na(risk_assess$risk_norm)]＝0

risk_assess$risk_norm＝100*risk_assess$risk_norm/max(risk_assess$risk_norm)

######risk_assess$risk_norm即为标化后的风险评分

print(risk_assess$risk_norm)

上面结合实施例对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，对于本技术领域的普通技术人员来说，在获知本发明中记载内容后，在不脱离本发明原理的前提下，还可以对其作出若干同等变换和替代，这些同等变换和替代也应视为属于本发明的保护范围。