用 SIR 模型拟合疫情感染情况

用 SIR 模型拟合疫情感染情况

整理了一下大三统计计算的课程作业

目录

数据抓取

数据来源：

抓取工具：

抓取记录：

数据处理：

数据可视化：

拟合模型：

模型选择：

模型解释

模型求解：

---

数据抓取 数据来源：

由于提供疫情数据的网站，本文选用的是腾讯新闻的实施更新，网站链接如下：实时更新：新冠肺炎疫情最新动态

https://news.qq.com/zt2020/page/feiyan.htm#/global

抓取工具：

python和jupyter notebook

抓取记录：

由于该网页展现的只是当日的数据，所以我们要先通过查找网页源代码，找出包含中国及世界历史疫情数据的网页，以下有三个链接分别表示的是包含中国各省市当日实时数据、包含中国历史数据及每日新增数据、包含世界历史数据及每日新增数据的链接

https://view.inews.qq.com/g2/getOnsInfo?name=disease_h5

https://view.inews.qq.com/g2/getOnsInfo?name=disease_other

https://view.inews.qq.com/g2/getOnsInfo?name=disease_foreign

 

以下只提供了获取湖北与非湖北的历史数据的代码

# 获取湖北与非湖北历史数据
def get_data_1():
    with open(filename, "w+", encoding="utf_8_sig", newline="") as csv_file:
        writer = csv.writer(csv_file)
        header = ["date", "dead", "heal", "now/confirm/i", "deadRate", "healRate"]  # 定义表头
        writer.writerow(header)
        for i in range(len(hubei_notHhubei)):
            data_row = [hubei_notHhubei[i]["date"], hubei_notHhubei[i][w]["dead"], hubei_notHhubei[i][w]["heal"],
                        hubei_notHhubei[i][w]["now/confirm/i"], hubei_notHhubei[i][w]["deadRate"],
                        hubei_notHhubei[i][w]["healRate"]]
            writer.writerow(data_row)

数据处理：

由于数据中存在数据缺失的情况，在画图中会存在一些问题，因此利用R对数据进行清洗；同时我们用SIR模型进行拟合时，并不会用到全部的数据，因此对不同时期的数据进行了选择。

data_washing <- function(input_dt){
  require(lubridate)
  input_dt$date_new <- format(as.Date(input_dt$date,format='%Y-%m-%d'),format='%m-%d')
  input_dt$date <- as.Date(input_dt$date,format = '%Y-%m-%d')
  #更改数据格式
  date_start <- ymd(input_dt$date[1])
  date_end <- ymd(input_dt$date[length(input_dt$date)])
  complete_date <- date_start + days(0:as.numeric(date_end-date_start))
  #将数据合并
  df1 <- data.frame(date = input_dt$date,total_confirm = input_dt$total_/confirm/i)
  df2 <- data.frame(date = complete_date,c = c(1:length(complete_date)))
  #df2 <- data.frame(date = complete_date,c = rep(0,length(complete_date)))
  merged_data <- merge(df1,df2,by = "date",all  =TRUE)
  return(merged_data)
}

数据可视化：

图一展示的是1月28号到6月3号世界总确诊人数和1月20号到6月3号中国确诊人数的对比散点图（注意这个图中的确诊人数是总的确诊人数，没有计算治愈和死亡，因此散点图是递增的，单独把中国的散点图拉出来，如图二，可以看出中国总的确诊人数在四月中旬至六月初呈现收敛的趋势，即新增确诊人数已经变得很少）

 

 

图三展示的是1月1月20号到6月3号中国、湖北、非湖北确诊人数的对比散点图

 

图四展示的湖北和非湖北的确的治愈率、死亡率的变化曲线

 

图五展示的是我国治愈率和死亡率的变化曲线（可以看出，我国的治愈率和死亡率在5月和六月趋于稳定）

 

拟合模型： 模型选择：

在对传染病模型的研究上有很多模型比如：SI、SIS、SERS、SIR等。本文利用SIR模型对这次新型冠状病毒的发展情况进行研究。

模型解释

SIR模型：SIR模型是是一种传播模型，是信息传播过程的抽象描述，是传染病模型中最经典的模型，其中S表示易感者，I表示感染者，R表示移出者。

基本假设

(1) 假设给定数据真实，不考虑人口的出生、死亡、流动等，该城市人口总数N是常数，即每年流入人口等于流出人口

(2) 假设人口当日所处的健康状态只分为健康人群、感染新型冠状病毒确诊病人和恢复者三类;

(3) 假设人口不可流动，即无病人流入和流出该城市;

(4) 假设严格隔离的新型冠状病毒感染病人都不再传染他人;

符号说明

N——所研究地区的人口总数；

S——易感类，该类成员没有感染新型冠状病毒，也没有免疫能力，是有可能被传染上新型冠状病毒的，设易感类人群t时刻在总人口中所占比例为s（t），S（t）=N * s（t）；

I——感染者，该类成员已经感染上新型冠状病毒，属于确诊患者，并且能够把病毒传染给S类（即易感类）成员，设该类人群t时刻在总人口中的比例为i（t），I（t）=N * i（t）；

R——免疫类，该类成员为新型冠状病毒康复者或因患新冠肺炎死亡，已经具有免疫力，不再对其他成员产生任何影响，设该类人群t时刻在总人口中的比例为r（t），R（t）=N * r（t）；

beta——被感染系数，在t时刻单位时间内被所有病人传染的人数为beta*N*s(t)*i(t)；

gamma——恢复系数，在t时刻单位时间内恢复者数量为gamma*N*i(t)；

t——模型中的时间单位。

感染机制：

感染者总人数为原时刻即j时刻的感染人数加上i时刻由易感者被感染成感染者的人数

新的免疫者（移除者）是由i时刻感染者变为免疫者的人数

S（i）+ I（j） I（i）+I（j）

I（i） R（i）

为拟合疫情爆发的动力学模型，需要下列三个微分方程:

注：该微分方程是用来表示3类人中人群向其他类转换速度的3个比率

模型求解：

#建立sir模型
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 * S
      dI <-  beta * I * S - gamma * I
      dR <-  gamma * I
      return(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)
}

1.求解该微分方程需要知道S、I、R的初始值，以及参数beta和gamma的值，其中在病毒爆发的初期我们可以假设每个人都有可能被感染，即 S0 = N，感染人数可以由数据得出，此时R0就得到了，由于该病潜伏期最大14天，设定γ=1/14，因此求解模型参数的问题旧转化成了求解beta参数的问题了。

2.根据公式R0 = beta/gamma，R0即基本传染数，指在流行病学上，没有外力介入，同时所有人都没有免疫力的情况下，一个感染到某种传染病的人，会把疾病传染给其他多少个人的平均数。R0对于疾病的传播和防控具有重要意义。 根据之前的一些研究，新冠肺炎的R0值总体上在2到3左右。但在美国CDC期刊《新兴传染病》杂志发表的一项新的研究认为，新冠肺炎的R0值高达5.7。

若微分方程根据方法2进行求解，根据R0 = 5.7，算出beta = 0.407，代入微分方程中进行计算，可以画出如下图像：

由图像可知我们把S的初始值设为1400000000时，确诊人数的拟合与真实情况是明显不符合的，因此我尝试缩小S的值，可以得到下图，可以看出，图像后半段与真实情况是比较符合的，因此我再次尝试改变R0的值和S的初始值

---

由图像可以看出，当用R0是5.7画出来的图是十分不准确的，因此我们可以通过改变R0的值，来和真实数据逼近。

 

第二种方法：求beta可以采用最小二乘方法，通过计算拟合的beta与真实beta曲线之间的最小二乘残差来判断拟合的好坏，选取使最小二乘残差即SSE值最小的beta值带入到SIR模型中进行拟合。我们用作循坏比较的beta范围是0.142~0.407之间且步长为0.001的数，从下面的运行结果可知，最合适的beta值为0.356，对应的R0值为4.984。

由下图可知，确诊人数真实值和拟合值的最小二乘残差在不同beta下的变化

 图九是SIR模型的拟合散点图以及确诊人数真实值的散点