如何让SEIR模型参数r随I舱室逐时间步发病率变化（ode实现）

你代码里的核心问题是动态参数r的定义方式错误：

• 定义parms时用的diff(I)调用的是初始状态的I=2000，而非随时间变化的I值，所以最终r被设为0，根本没实现“每个时间步依赖I的差值”的逻辑。
• deSolve的ode求解器中，parms是固定参数集合，无法直接引用随时间变化的状态变量（比如I），必须把动态参数的计算放到微分方程函数内部，或者改用离散时间模型。

下面给出两种针对性的解决方案：

方案一：离散时间SEIR模型（完全匹配“每个时间步用I差值”需求） #

如果你的需求是严格基于**离散时间步（比如每天）**计算当前I与前一步I的差值来更新r，直接用循环逐步计算更直观，不需要依赖deSolve的连续时间求解：

# 初始状态
start <- c(S=5000, E=3000, I=2000, R=0, CInc=0)
# 时间步：0到10000，步长1
times <- seq(0, 10000, 1)
# 固定参数
fixed_parms <- c(lambda=0.3, gamma_h=1/20, rho=1/(1*365))

# 初始化结果矩阵
result <- matrix(nrow=length(times), ncol=length(start))
colnames(result) <- names(start)
result[1,] <- start

# 逐时间步计算
for (t in 2:length(times)) {
  # 提取上一步的状态
  prev <- result[t-1,]
  
  # 计算I的差值：当前步（t）的I还没算，所以用前两步的I差值
  delta_I <- if (t == 2) 0 else (prev["I"] - result[t-2, "I"])
  # 计算r，限制范围避免不合理值（比如超过1或负数）
  r <- max(min(delta_I / 10000, 1), 0)
  
  # 更新各舱室状态
  result[t, "S"] <- prev["S"] + (-fixed_parms["lambda"] * prev["S"] + fixed_parms["rho"] * prev["R"])
  result[t, "E"] <- prev["E"] + (fixed_parms["lambda"] * prev["S"] - fixed_parms["gamma_h"] * prev["E"])
  result[t, "I"] <- prev["I"] + (fixed_parms["gamma_h"] * prev["E"] - r * prev["I"])
  result[t, "R"] <- prev["R"] + (r * prev["I"] - fixed_parms["rho"] * prev["R"])
  result[t, "CInc"] <- prev["CInc"] + fixed_parms["lambda"] * prev["S"]
}

# 转成数据框方便后续分析/绘图
result_df <- cbind(time=times, as.data.frame(result))

方案二：连续时间模型中动态计算r（基于I的变化率） #

如果你只是想让r依赖I的变化趋势（连续时间中的导数），而非严格的离散差值，可以在微分方程函数内部动态计算r：

library(deSolve)

seirs <- function(t, x, parms)  {
  with(as.list(c(parms, x)), {
    # 先计算I的变化率（近似离散差值的连续版本）
    dI_rate <- gamma_h * E
    # 计算r，同样限制范围
    r <- max(min(dI_rate / 10000, 1), 0)
    
    # 正式计算各舱室的导数
    dS <- -lambda*S + rho*R
    dE <- lambda*S - gamma_h*E
    dI <- gamma_h*E - r*I
    dR <- r*I - rho*R
    dCInc <- lambda*S
    
    return(list(c(dS, dE, dI, dR, dCInc)))
  })
}

start <- c(S=5000, E=3000, I=2000,R=0, CInc = 0 )
parms <- c(lambda = 0.3, gamma_h = 1/20, rho = 1/(1*365))
times <- seq(0, 10000, 1)
result <- ode(times=times, y=start, func=seirs, parms=parms)

关键说明 #

1. 离散模型的优势是完全贴合“时间步差值”的需求，但需要手动处理循环和状态更新；
2. 连续模型更符合SEIR的经典假设，但这里的r是基于I的瞬时变化率，而非离散步的差值；
3. 无论哪种方案，都要注意限制r的取值范围（比如max(min(r,1),0)），避免出现不合理的参数值导致模型崩溃。

内容的提问来源于stack exchange，提问作者ccc