修改 SIR 模型以包含随机性

Question

我正在尝试通过将真实流行曲线与随机 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))

Answer 1

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

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

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)

Answer 2

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

• 泛化，使其可以与固定参数和随机强制一起工作
• 以列表形式传递参数
• 运行基本的蒙特卡罗模拟

library('deSolve')

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

    # 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, dI, dR))
  })
}

initial_state_values <- c(S = 10000, I = 1, R = 0)
times <- seq(from = 0, to = 50, 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,
            times = times,
            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, times = times,
          func = sir_model, 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, times = times,
      func = sir_model, 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))