R语言SIR模型

##加载所需要的包
library(deSolve) 
library(ggplot2)
library(ggsci)

##根据微分方程建模型
sir <- function(time, state, pars) {
  with(as.list(c(state, pars)), {
    dS <- -beta * S * I/N
    dI <- beta * S * I/N - gamma * I
    dR <- gamma * I
    return(list(c(dS, dI, dR)))
  })
}

##设置参数
N <- 10000 # 总人口
I0 <- 1 # 初始感染者数量
RM0 <- 0 # 初始移除人员数量
S0 <- N - I0 - RM0 # 初始易感人群数量
init <- c(S = S0, I = I0, R = RM0) # 初始值
# 以下参数在模型假定下是常量
pars <- c(
  beta = 0.55, # 有效接触率
  gamma = 0.2, # 移除率
  N = N # 人口
)

# 迭代次数，以天计
times <- seq(0, 75, by = 1)

#计算
res <- as.data.frame(ode(y = init, times = times, func = sir, parms = pars))

#作图
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),#坐标轴样式:黑色,size为线粗细
                         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=12,colour='black'))#图片脚注
ggplot(res) +
  geom_line(aes(x = time, y = S, col = '易感人群'),size=1.2)+
  geom_line(aes(x = time, y = I, col = '感染人群'),size=1.2)+
  geom_line(aes(x = time, y = R, col = '康复人群'),size=1.2)+
  scale_colour_manual("",values=c("易感人群"="cornflowerblue","感染人群" = "darkred", "康复人群" = "forestgreen") ) +#设置图例颜色
  labs(x="观察时间",y="人数",title="SIR模型 ",caption=" ")+theme2
#方法二
library(EpiModel)
param <- param.dcm(inf.prob = 0.2,act.rate = 1, rec.rate = 1/20,
                   a.rate = 1/95, ds.rate =1/100, di.rate = 1/80, dr.rate = 1/100)
init <- init.dcm(s.num = 1000, i.num =1, r.num = 0)
control <- control.dcm(type ="SIR", nsteps = 500, dt = 0.5)
mod <- dcm(param, init, control)
plot(mod, alpha = 0.5,
     lwd = 4, main = "Compartment Sizes")
#每个时刻具体数据
par(mfrow = c(1, 1))
comp_plot(mod, at = 50, digits = 1)