【新冠肺炎】SIR模型预测的编码关键词：新冠肺炎、SIR模型预测、数据分析
本篇使用r语言
以前介绍过SIR模型的基本理论及其微分方程，不熟悉的请看这里：
SIR SIRE 传染病预测模型与代码应用之概念篇
根据SIR模型的微分方程，了解了病毒感染的各种特征，包括何时终止、最大规模、疫情控制手段、疫苗和群体免疫问题等，不熟悉的朋友不妨回顾一下这篇文章：
【cov-19】新冠肺炎的SIR模型补充与应用
掌握这些基础知识后，可以基于微分方程在r上建立SIR模型，并根据模型预测病毒的整体流向。 其实，对于预测问题来说，建立模型非常简单，我们知道提高精度的关键主要是调参和概念上的修正。 此处简要介绍了模型的构建方法和比较预测。
首先强调，模型预测并不可靠，只能表现出一些数据上的趋势，真实情况还是需要根据事实自我判断。
1.建立模型
思路：先编辑SIR模型，再用ode完成仿真：
注意：这里使用fraction方程式，即百分比的数量。
sir _ model-function (beta，gamma，S0，I0，R0，times ) require(desolve ) #-------sir微分模型sir _ equation parameters )、{#百分比fraction方程ds-- bet a* I * SDI-bet a* I * s-gamm a * IDR-gamm a* ireturn (list (c (ds，di，di )， (dR ) ) )、伽玛) initial_values-c ) beta=beta，定义伽玛(initial_values-c ) s=S0 )
我们开始写测试主体了。 这里以加拿大新冠引起的肺炎数据为例。
2.清洗数据以及数据探索
读取数据：
setwd(((/sir_model/) ) source ) model_fit.r ) #---SIR模型的主体cov_df-read.CSV ) covid19_5_5.5
最左边的pruid是每个省的号码，其中1表示整个加拿大。 剩下的变量是
numconf:确诊数
numprob:伪数
numdeath:死亡人数
numtotal:确诊数疑似数
numtested:测试数
这里我们感兴趣的是加拿大单位时间的总确诊数，也就是numtotal。
2.1数据概览
打开数据，可以看到前期加拿大政府不太重视疫情统计，到4月份有几个日期是连续的。 那么，写出填补时间的方程式，根据对应的省id提取省的数据。
# date # date data _ washing-function (input _ dt，prov_id ) require ) lubridate ) proc _ data-cov _ df [ cov _ ] format='%Y-%m-%d ' )，format=format=' % y-% m-% d ' # completedate _ start-ymd (proc _ data $ date ) ) date_end-ymd ) proc_data ) date(1) date_end-ymd ) proc complete _ date-date _ start days (0: as.numeremert ) mergedf1-data.frame ) date=num total=proc _ data $ num total (df2-data.frame ) date=complete_date，c ) 110
date))) merged_proc_data<-merge(df1,df2,by="date",all=TRUE) return(merged_proc_data)}2.3透视数据
#data exploring#<------canada total confirmedcanada_total<-data_washing(cov_df,"1")trt_total<-data_washing(cov_df,"35")qbc_total<-data_washing(cov_df,"24")ggplot(data=canada_total,aes(x=canada_total$date,y=canada_total$numtotal))+ geom_point(aes(x=canada_total$date,y=canada_total$numtotal,colour="CA_Total"))+ geom_point(data=trt_total,aes(x=trt_total$date,y=trt_total$numtotal,colour="trt_total"))+ geom_point(data=qbc_total,aes(x=qbc_total$date,y=qbc_total$numtotal,colour="qbc_total"))+ theme(axis.text.x = element_text(angle=90,hjust=1))2.3.1加拿大总体的确诊走向已经进入疫情爆发的前期阶段。
2.3.2 Ontario（安省）政府一直对疫情比较重视，比较早的颁布了社交距离禁令，所以总体确诊要低于魁北克。
2.3.3 Quebec（魁北克）政府并不是很重视疫情监控，可以看到头铁的魁北克3月之前都没有记录数据，而在开始记录之后的短短15天直接指数起飞。
3.SIR模型数据预测
3.1 先看一下安省，人口为1457万，由于安省的社交禁令比较早，所以S（可能感染者）之间的接触系数也就会比较小。各媒体官方给出的R0为1.4-6.49之间，结合官方数据再加上一顿乱敲，我们取R0=1.9, gamma=0.13,beta=0.24
来一张局部放大图，我们设置时间周期为100天：
根据实际数据显示，截至5月5日，只有0.125%的安省人被感染，也就是18213人左右。而根据模型来看，已感染的人数在1.125%左右，也就是16万多，整整高出10倍。这里有3种可能：1)参数有误差；2)检测样本小，无法体现总确诊；3)人们通过各种手段延缓、抑制了疫情的发展。
可以通过一下方案解决
1）调参：修改beta为平均值0.328，这时的R0为3，可以看到疫情基本上要过去了：
给大家推荐一个网站，可以在线对SIR模型进行模拟，有兴趣的朋友可以试试不同的参数对疫情的情况有什么影响，这里就不多介绍了：
cjdqc Nesse - Global Health - SIR Model
2）样本容量，“真实”数据也不一定真实。加拿大对于有症状人群的检测方式：
所以这里我们只说模型的应用，不讨论数据。
3.2 预测魁北克省情况
魁省的数据3月以后才出现，而短短15天飙升到几万人，可见这里的beta是偏高的，所以我们取平均值偏低一点儿0.3。这里魁省的人口为260万人。
强调一点，模型的结果并不是一定可靠的，虽然前期很fit，但在实际的情况中，有很多人为的因素是模型无法应变的，比如人们恐慌后开始加强防护措施，保持社交距离等，都会使曲线变得更加平缓，甚至在大爆发前抑制住病毒。
第二，对于得到的数据我们所做的模型不应该很完美的fit，毕竟每天检测的sample size有限。打个比方，A市每天检测100名报名的人，而每天确诊的患者都会增加5人。那我们能说A市每天确诊增加了5人吗？真实数字当然远远大于5人，所以我们的模型结果所表示的更应该是一种整体的趋势，而不是用来与特定的样本来比较。且在数据中，我们更应该关心样本中确诊人数的占比。
Ontario比例：
Quebec 比例：
从确诊比例图中我们可以看出，Ontario似乎已经经过了一波玩命的云朵，目前已经趋于一个下降的过程。而Quebec还在持续的攀升中。
再次强调，模型预测并不可靠，只能表现出一些数据上的趋势，真实情况还是需要根据事实自我判断。
SIR模型到这里就结束了，需要完整代码可以评论。
Reference:
1. L.Eggertson cmaj News ‘COVID-19: Recent updates on the coronavirus pandemic’
url: https://cmajnews.com/2020/05/06/coronavirus-1095847/
2. Population of Cities in Canada (2020)
url: https://worldpopulationreview.com/countries/canada-population/cities/
3. some r code example
url: https://rpubs.com/choisy/sir