修改 SIR 模型以包含随机性

问题描述

我正在尝试通过将真实流行曲线与随机 SIR 模型的模拟进行比较来建立一种估计传染病参数的方法。为了构建随机 SIR 模型，我使用了 deSolve 包，而不是使用固定参数值，我想从以原始参数值为中心的泊松分布中的每个时间点绘制方程中使用的参数值。

以参数 beta 为例，beta 表示人均传播事件的平均次数，是平均接触次数与接触时发生传播的概率的乘积。实际上，一个人将拥有的接触次数存在差异，并且由于传播也是一种概率事件，因此也存在差异。 因此，即使平均传播率为 2.4（例如），一个人也可以继续以不同的概率感染 0、1、2 或 3 名等人。

我尝试使用 rpois 函数将其合并到我下面的代码中，并将方程中使用的参数重新分配给 rpois 的输出。

我已多次使用相同的初始值和参数运行我的代码，并且所有曲线都不同，表明某些“随机”正在发生，但我不确定代码是否在每个时间点使用 rpois 进行采样或一开始就只有一次。我最近才开始编码，所以没有太多经验。

如果有比我更有经验的人可以验证我的代码实际上在做什么以及它是否在每个时间点使用 rpois 进行采样，我将不胜感激。如果不是，我将不胜感激任何实现这一目标的建议。也许需要一个循环？

library('deSolve')
library('reshape2')
library('ggplot2')

#MODEL INPUTS
 initial_state_values <- c(S = 10000,I = 1,R = 0)

 #ParaMETERS
  parameters <- c(beta = 2.4,gamma = 0.1)


 #POISSON MODELLING OF ParaMETERS
 #BETA
  beta_p <- rpois(1,parameters[1])

 #GAMMA
  infectIoUs_period_p <- rpois(1,1/(parameters[2]))
 
  gamma_p <- 1/infectIoUs_period_p

 #timestepS
  times <- seq(from = 0,to = 50,by = 1)

 #SIR MODEL FUNCTION 
  sir_model <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
  N <- S + I + R   
  lambda <- beta_p * I/N 
  dS <- -lambda * S
  dI <- lambda*S - gamma_P*I
 dR <- gamma_P*I

return(list(c(dS,dI,dR)))
 })
 }

 output<- as.data.frame(ode(y= initial_state_values,times = times,func = sir_model,parms = parameters))

解决方法

问题中给出的代码随着时间的推移以恒定参数运行模型。这是一个参数随时间变化的示例。但是，此设置假设对于给定的时间步长，总体中的所有个体的参数都相同。如果您想获得个体变异性，可以对不同的亚群使用矩阵公式，也可以改用单个模型。

具有波动人口参数的模型：

library('deSolve')

initial_state_values <- c(S = 10000,I = 1,R = 0)

parameters <- c(beta = 2.4,gamma = 0.1)

times <- seq(from = 0,to = 50,by = 1) # note time step = 1!

# +1 to add one for time = zero
beta_p <- rpois(max(times) + 1,parameters[1])

infectious_period_p <- rpois(max(times) + 1,1/(parameters[2]))

gamma_p <- 1/infectious_period_p

sir_model <- function(time,state,parameters) {
  # cat(time,"\n") # show time steps for debugging
  with(as.list(c(state,parameters)),{
    
    # this overwrites the parms passed via parameters
    beta  <- beta_p[floor(time)  + 1]
    gamma <- gamma_p[floor(time) + 1]
    
    N <- S + I + R   
    lambda <- beta * I/N 
    
    dS <- -lambda * S
    dI <- lambda * S - gamma * I
    dR <- gamma * I
    
    list(c(dS,dI,dR))
  })
}

output <- ode(y = initial_state_values,times = times,func = sir_model,parms = parameters)

plot(output)

,

这是另一个稍微通用的版本。添加它作为第二个答案，以保持原始版本的紧凑和简单。新版本在以下方面有所不同：

• 泛化，以便它可以使用固定参数和随机强迫
• 将参数作为列表传递
• 运行基本的蒙特卡罗模拟

library('deSolve')

sir_model <- function(time,parameters) {
  with(as.list(c(state,{

    # this overwrites the parms passed via parameters
    if (time_dependent) {
      beta  <- beta_p[floor(time)  + 1]
      gamma <- gamma_p[floor(time) + 1]
    }
    N <- S + I + R
    lambda <- beta * I/N

    dS <- -lambda * S
    dI <- lambda * S - gamma * I
    dR <- gamma * I

    list(c(dS,dR))
  })
}

initial_state_values <- c(S = 10000,R = 0)
times <- seq(from = 0,by = 1) # note time step = 1!

## (1) standard simulation with constant parameters
parameters <- c(beta = 2.4,gamma = 0.1)
out0 <- ode(y= initial_state_values,func  = sir_model,parms = c(parameters,time_dependent = FALSE))
plot(out0)

## (2) single simulation with time varying parameters
beta_p <- rpois(max(times) + 1,parameters[1])
infectious_period_p <- rpois(times + 1,1/(parameters[2]))
gamma_p <- 1/infectious_period_p

## here we need pass the vectorized parameters globally
## for simplicity,it can also be done as list
out1 <- ode(y = initial_state_values,parms = c(time_dependent = TRUE))
plot(out0,out1)

## (3) a sample of simulations
monte_carlo <- function(i) {
  #parameters <- c(beta = 2.4,gamma = 0.1)
  beta_p <- rpois(max(times) + 1,parameters[1])
  infectious_period_p <- rpois(max(times) + 1,1 / (parameters[2]))
  gamma_p <- 1/infectious_period_p
  ode(y = initial_state_values,parms = list(beta_p  = beta_p,gamma_p = gamma_p,time_dependent = TRUE))
}

## run 10 simulations
out_mc <- lapply(1:10,monte_carlo)
plot(out0,out_mc,mfrow=c(1,3))