修改 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))

修改 SIR 模型以包含随机性

Modifying SIR model to include stochasticity

parameters

r

poisson

stochastic