R语言SEIR模型

#定义初始状态各人数
N = 50000 #初始总人数
E = 0     #暴露人数
I = 1     #感染人数
S = N-I   #易感人数
R = 0     #治愈人数

# 常量赋值
r = 150    #r：感染患者每天接触的易感者数目
B = 0.03  #β：传染系数;由疾病本身的传播能力,人群的防控能力决定
a = 0.1   #α：潜伏者的发病概率;一般为潜伏期的倒数
y = 0.1   #γ：恢复系数;一般为病程的倒数

#设置时间观察(eg:到第100天)
T = 150

#构建上述方程组
for (i in 1:(T-1)){
  S[i+1] = S[i] - r*B*S[i]*I[i]/N
  E[i+1] = E[i] + r*B*S[i]*I[i]/N - a*E[i]
  I[i+1] = I[i] + a*E[i] - y*I[i]
  R[i+1] = R[i] + y*I[i]
}

#生成表格
seir <- data.frame(S, E, I, R)
seir$date<-seq(1,150,1)
head(seir,10)
inflection_point<-subset(seir[,"date"],I==max(I))#拐点
inflection_num<-ceiling(max(I))#向上取整
#长宽数据转换
library(tidyr)#gather
seir_long<-gather(seir,condition,num,S:R)
head(seir_long,10)

#作图
library(ggplot2)
library(ggsci)#配色用的包
#设置Theme
windowsFonts(myFont = windowsFont("微软雅黑")) #设定桌面字体，以备图形使用
theme2<-theme_bw()+theme(panel.border=element_blank(),#图形边界
                         panel.grid.major=element_blank(), #网格线
                         panel.grid.minor=element_blank(), #次级网格线
                         legend.title=element_text(family="myFont",size=14,colour='black'), #图例标题
                         legend.text=element_text(family="myFont",size=14,colour='black'), #图例文字
                         legend.position="bottom",#图例位置
                         legend.background =element_blank(),#图例背景
                         axis.line= element_line(colour = "black",size=0.6),#坐标轴样式：颜色黑色，粗细0.6
                         axis.title=element_text(family="myFont",size=14,colour='black'),#坐标轴标题
                         axis.text=element_text(family="myFont",size=14,colour='black'), #坐标轴文字
                         plot.title=element_text(family="myFont",size=20,colour='black',hjust=0.5),#图片标题居中
                         plot.caption=element_text(family="myFont",size=14,colour='black'))#图片脚注
#作图
label<-paste("拐点出现日期 : 第",inflection_point,"天\n累计感染人数 : ",inflection_num,"人",sep=" ")#在图中添加拐点标签
ggplot(seir_long)+#绘出基础图
  geom_line(aes(x=date,y=num,group=condition,colour=condition),size=1.2)+#四条曲线颜色及粗细
  geom_vline(xintercept=inflection_point,linetype=2,colour="darkgrey",size=1)+#添加拐点竖线
  geom_hline(yintercept=inflection_num,linetype=2,colour="darkgrey",size=1)+#添加拐点横线
  geom_label(label=label,x=90,y=3400,size=4,family="myFont") +#标签位置及大小设置
  scale_color_npg()+#调用NPG配色方案
  labs(x="观察时间",y="人数",title="SEIR模型",caption=" ")+#坐标轴及标题名称设置
  theme2#使用上边设置好的主题
#无标签的
ggplot(seir_long)+#绘出基础图
  geom_line(aes(x=date,y=num,group=condition,colour=condition),size=1.2)+#四条曲线颜色及粗细
  geom_vline(xintercept=inflection_point,linetype=2,colour="darkgrey",size=1)+#添加拐点竖线
  geom_hline(yintercept=inflection_num,linetype=2,colour="darkgrey",size=1)+#添加拐点横线
  scale_color_npg()+#调用NPG配色方案
  labs(x="观察时间",y="人数",title="SEIR模型",caption=" ")+#坐标轴及标题名称设置
  theme2#使用上边设置好的主题