修改 SIR 模型以包含随机性

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

以参数β为例，β表示人均传播事件数，是平均接触次数与接触传播概率的乘积。实际上，一个人接触的次数存在差异，而且由于传播也是一个概率事件，因此也存在差异。 因此，即使平均传播率为 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)