【新冠肺炎】SIR模型预测与数据分析之代码篇

【新冠肺炎】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_equations <- function(time, variables, parameters) {with(as.list(c(variables, parameters)), {#使用百分比的fraction方程dS <- -beta * I * SdI <-  beta * I * S - gamma * IdR <-  gamma * Ireturn(list(c(dS, dI, dR)))})}#----------定义变量parameters_values <- c(beta  = beta, gamma = gamma)initial_values <- c(S = S0, I = I0, R = R0)#---------通过ode模拟疫情out <- ode(initial_values, times, sir_equations, parameters_values)#----------返回结果as.data.frame(out)
}

2.清洗数据以及数据探索
我们开始编写测试主体，这里以加拿大的新冠肺炎数据为例。
2.1数据概览
读入数据：

setwd("~ /SIR_Model/")
source("model_fit.R")#<---SIR模型的主体
cov_df<-read.csv("covid19_5_5.csv")

首先我们概览一下数据长什么样

最左边的pruid是每个省的编号，其中1代表加拿大整体。剩下的variable为
numconf:确诊数
numprob:疑似数
numdeath:死亡数
numtotal:确诊数+疑似数
numtested:测试数
这里我们感兴趣的是加拿大单位时间的总确诊数，也就是numtotal。
2.2数据清洗
打开数据可以发现前期加拿大政府并不是很重视疫情的统计，四月之前有一些日期是不连续，那么我们写一个方程来填补时间以及根据对应省id提取出该省的数据：

#date
#date
data_washing<-function(input_dt,prov_id){require(lubridate)proc_data<-cov_df[cov_df$pruid==prov_id,]#formationproc_data$date_new<-format(as.Date(proc_data$date,format='%Y-%m-%d'),format='%m-%d')proc_data$date<-as.Date(proc_data$date,format='%Y-%m-%d')#complete datedate_start<-ymd(proc_data$date[1])date_end<-ymd(proc_data$date[length(proc_data$date)])complete_date<-date_start+days(0:as.numeric(date_end-date_start))#mergedf1<-data.frame(date=proc_data$date,numtotal=proc_data$numtotal)df2<-data.frame(date=complete_date,c=c(1:length(complete_date)))merged_proc_data<-merge(df1,df2,by="date",all=TRUE)return(merged_proc_data)
}

2.3透视数据

#data exploring
#<------canada total confirmed
canada_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

t<-merged_plot(input_data=trt_total,model_duration=seq(1,200),#时间周期200天population = 14570000,#安省人口一千四百多万beta = 0.25,gamma = 0.11)

来一张局部放大图，我们设置时间周期为100天：

根据实际数据显示，截至5月5日，只有0.125%的安省人被感染，也就是18213人左右。而根据模型来看，已感染的人数在1.125%左右，也就是16万多，整整高出10倍。这里有3种可能：1)参数有误差；2)检测样本小，无法体现总确诊；3)人们通过各种手段延缓、抑制了疫情的发展。
可以通过一下方案解决
1）调参：修改beta为平均值0.328，这时的R0为3，可以看到疫情基本上要过去了：

给大家推荐一个网站，可以在线对SIR模型进行模拟，有兴趣的朋友可以试试不同的参数对疫情的情况有什么影响，这里就不多介绍了：
Hans Nesse - Global Health - SIR Model

2）样本容量，“真实”数据也不一定真实。加拿大对于有症状人群的检测方式：

所以这里我们只说模型的应用，不讨论数据。

3.2 预测魁北克省情况
魁省的数据3月以后才出现，而短短15天飙升到几万人，可见这里的beta是偏高的，所以我们取平均值偏低一点儿0.3。这里魁省的人口为260万人。

#qbc
merged_plot(input_data=qbc_total,model_duration=seq(1,200),population = 2600000,beta = 0.3,gamma = 0.11)

强调一点，模型的结果并不是一定可靠的，虽然前期很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

【新冠肺炎】SIR模型预测与数据分析之代码篇相关推荐

最新！兰州大学发布对上海市的新冠肺炎疫情预测！
这段时间,上海市疫情牵动着所有人的心.据数据显示,自 2022 年 3 月 1 日上海市报告新冠肺炎本土确诊病例和本土无症状感染者以来,截至 2022 年 4 月 10 日 24 时,上海市已累计报告 ...
Python小白的数学建模课-B5. 新冠疫情 SEIR模型
传染病的数学模型是数学建模中的典型问题,常见的传染病模型有 SI.SIR.SIRS.SEIR 模型. 考虑存在易感者.暴露者.患病者和康复者四类人群,适用于具有潜伏期.治愈后获得终身免疫的传染病. 本 ...
Python小白的数学建模课-B3. 新冠疫情 SIS模型
传染病的数学模型是数学建模中的典型问题,常见的传染病模型有 SI.SIR.SIRS.SEIR 模型. SIS 模型型将人群分为 S 类和 I 类,考虑患病者可以治愈而变成易感者,但不考虑免疫期. 本文 ...
Python小白的数学建模课-B2. 新冠疫情 SI模型
传染病的数学模型是数学建模中的典型问题,常见的传染病模型有 SI.SIR.SIRS.SEIR 模型. SI 模型是最简单的传染病模型,适用于只有易感者和患病者两类人群. 我们就从 SI 模型开始吧,从 ...
AAAI 2021：一种跨城市迁移的新冠肺炎高危社区发现框架
新冠肺炎已经在世界范围内广泛传播,严重影响着人们的日常生活.面对新冠肺炎,人为干预的空间隔离手段(如限制出行或集中隔离)已经被证明其有效性.但是,确诊病例的统计往往是滞后且粗粒度的,比如对于尚未确诊的 ...
【大数据平台】基于Spark的美国新冠肺炎疫情数据分析及预测
(本实验系中国地质大学(武汉)2022年秋期大数据平台及应用课程设计) 一.选题背景新型冠状病毒疫情是由严重急性呼吸系统综合征冠状病毒2(SARS-CoV-2)导致的2019冠状病毒病(COVID- ...
Matlab基于SEIRD模型，NSIR预测模型，AHP层次分析法新冠肺炎预测与评估分析
全文链接:http://tecdat.cn/?p=32175 分析师:Jiahui Zhao 新型冠状病毒肺炎COVID-19 给中国乃至全世界都带来了深重的灾难,对世界经济也造成了不可逆的影响(点击 ...
【cov-19】新冠肺炎的SIR模型补充与应用
[cov-19]SIR模型应用上一篇文章中我们讲了SIR的基本概念,不熟悉的朋友可以看这篇: https://blog.csdn.net/weixin_45265581/article/detail ...
Python数据分析高薪实战第十天 EDA实战-全球新冠肺炎确诊病例趋势分析
27 初识 EDA:全球新冠肺炎确诊病例趋势分析从本讲开始,我们会通过四个具体的案例来将我们之前学习的 Python 数据分析方面的知识全都串起来.一方面能够融会贯通,另一方面也能帮你掌握数据分析基 ...
使用SAP Analytics Cloud显示全球新冠肺炎确诊人数和发展趋势的预测
注:本文只是借用新冠肺炎全球确诊人数作为历史数据,来介绍SAP Analytics Cloud基于机器学习的Time Series Forecasting功能,并没有对现实世界中新冠肺炎的发展趋势做出 ...

【新冠肺炎】SIR模型预测与数据分析之代码篇

【新冠肺炎】SIR模型预测之代码篇

【新冠肺炎】SIR模型预测与数据分析之代码篇相关推荐

最新文章

热门文章