ODE状态依赖参数调整问题：带疫苗阈值的SIR模型报错求助

问题分析 #

1. 导数数量不匹配：状态向量包含8个变量（S、I、R、P、Sn、In、Rn、alpha），但原模型仅返回7个变量的导数，缺少dalpha的定义。
2. Event函数返回格式错误：event_func的返回值不符合deSolve的要求，未正确返回完整的修改后状态向量。
3. 未启用root与event机制：原模型运行代码未指定rootfun和events参数，导致状态触发的事件从未执行。
4. 冗余时间依赖判断：模型中alpha <- ifelse(t >=100, 0, alpha)与状态触发的需求冲突，需移除。

修正后的完整代码 #

1. 修正Root函数（保留你定义的接种人口阈值逻辑）

root_func <- function(t, states, parameters) {
  with(as.list(c(states, parameters)), {
    N <- S + I + R + P + Sn + In + Rn
    # 触发条件：接种人口（P+Rn）达到总人口的50%
    return(P + Rn - 0.5 * N)  
  })
}

2. 修正Event函数

Event函数需返回完整的修改后状态向量，格式符合deSolve要求：

event_func <- function(t, states, parameters) {
  states["alpha"] <- 0  # 触发事件后将接种率设为0
  return(list(y = states))
}

3. 修正模型函数（补充导数、移除时间依赖）

为状态变量alpha添加导数dalpha=0（仅通过event修改，微分方程中无变化），同时移除时间依赖的alpha判断：

model_1_take <- function(t, states, parameters) {
  with(
    as.list(c(states, parameters)),
    {
      N <- S + I + R + P + Sn + In + Rn
      
      dS <- (mu * N) - 
        (beta * S * (I + In) / N) - 
        (alpha * epsilon * S) - 
        (alpha * (1 - epsilon) * S) - 
        (nu * S)
      dI <- (beta * S * (I + In) / N) - 
        (gamma * I) - 
        (nu * I)
      dR <- (gamma * I) - 
        (alpha * R) - 
        (nu * R)
      dP <- (alpha * epsilon * S) - 
        (nu * P)
      dSn <- (alpha * (1 - epsilon) * S) - 
        (beta * Sn * (I + In) / N) - 
        (nu * Sn)
      dIn <- (beta * Sn * (I + In) / N) - 
        (gamma * In) - 
        (nu * In)
      dRn <- (alpha * R) +
        (gamma * In) -
        (nu * Rn)
      dalpha <- 0  # alpha仅通过event修改，微分导数为0
      
      return(list(c(dS, dI, dR, dP, dSn, dIn, dRn, dalpha)))
    }
  )
}

4. 修正模型运行代码（启用root与event）

在ode调用中添加rootfun和events参数，指定状态触发的事件：

# 定义时间序列（示例）
times <- seq(0, 365, by = 1)

out <- ode(
  y = states_1_take,
  times = times, 
  func = model_1_take,
  parms = parameters_1_take,
  method = "lsodes",
  rootfun = root_func,
  events = list(func = event_func, root = TRUE)  # 由root函数触发事件
)

关键说明 #

• 导数匹配：状态向量的每个变量必须对应一个导数，添加dalpha=0确保返回的导数数量与状态数量一致。
• Event函数规范：必须返回包含修改后状态的列表（list(y = states)），保证所有8个状态变量都被正确返回。
• 状态触发逻辑：通过events=list(func=event_func, root=TRUE)指定事件由根函数触发，当root_func返回0时自动执行停止接种的操作。
• 逻辑一致性：移除时间依赖的alpha判断，确保仅由接种人口阈值触发停止接种行为。

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