SIRD模型拟合死亡时序数据报错排查与代码修正请求

问题概述 #

使用R实现SIRD模型拟合死亡时序数据，目标是最小化预测值与观测值的RMSE，但运行时出现错误：

Error in while (tol < End) { : missing value where TRUE/FALSE needed

以下是可复现的数据集、原代码，以及问题分析与修正方案：

可复现数据集 #

time <- 1:40
D_obs <- round(seq(1, 200, length.out = 40))
hmod <- data.frame(time, D_obs)
print(hmod)

原模型代码 #

hmod = (assign_df)  # 未定义变量错误

sir_model = function (t, y, parms)
{
  # label state variables from within the y vector of initial values
  S = y[1]        # susceptibles
  I = y[2]        # infectious
  R = y[3]        # recovered
  D = y[4]        # dead
  N = S + I + R + D
  
  # Pull parameter values from the parms vector
  beta = parms["beta"]
  gamma = parms["gamma"]
  alpha = parms["alpha"]
  
  # Define equations (dS etc. are simplified notation for dS/dt)
  dS = -(beta*S*I)/N
  dI = (beta*S*I)/N - gamma*I
  dR = alpha*gamma*I
  dD = (1-alpha)*gamma*I
  
  # return list of gradients   
  res = c(dS, dI, dR, dD)
  list(res)
}

par_init = c(beta=1, gamma=1, alpha =0.5)
init_sav = c(S=1079, I=1, R=0, D=0)
times_sav = hmod$time

n_obs <- length(hmod$D_obs)

opt_fun_sav = function(par){
  
  ## The following lines prevent 'optim_nm' from wandering into negative parameter space
  if (par[1]<0){
    par[1] = -par[1]
  } 
  if(par[2]<0){
    par[2] = -par[2]
  }
  ## Run the SIR model with the current parameter values 
  output_sav_try = ode(
    y=init_sav, 
    times = times_sav, 
    func = sir_model, 
parms= c(
      beta = par_init[1], 
      gamma = par_init[2],
      alpha = par_init[3]
    )
  )
  n_obs = nrow(hmod)
  
  output_savtry_df = as.data.frame(output_sav_try) # coerce to a data.frame object

  
  ## Create a data frame with observed and model prediction data for Savoonga
  sav_rmse_df = data.frame(
    time = hmod$time,
    obs = hmod$D_obs,
    # the following lines find the cumulative I prediction for each time in the mummps data
    pred = sapply(1:n_obs, function(i){output_savtry_df$D[
      output_savtry_df$time == hmod$time[i]
    ]})
  )
  # Calculate the RMSE for this prediction
  rmse_try = sqrt((1/n_obs)*(sum((sav_rmse_df$pred - sav_rmse_df$obs)^2)))  
  
  return(rmse_try)
  
}

optnm_result=optim_nm(opt_fun_sav,start=c(0.6,0.6),k=2, trace=TRUE)

问题分析 #

• 未定义变量错误：第一行hmod = (assign_df)中assign_df未定义，直接覆盖了之前生成的有效数据集，导致后续计算出现缺失值。
• 参数传递逻辑错误：优化函数opt_fun_sav调用ode时，始终使用固定的par_init而非输入的优化参数par，模型无法进行参数迭代，最终导致RMSE计算异常，触发optim_nm循环条件的缺失值错误。
• 优化参数数量不匹配：optim_nm起始参数仅传入2个，但模型有3个待优化参数（beta、gamma、alpha），参数维度不一致。
• 冗余时间匹配逻辑：output_sav_try的时间序列与hmod$time完全对齐，无需用sapply逐一匹配，直接提取output_savtry_df$D即可。

修正后的代码 #

# 加载必要依赖包
library(deSolve)
library(dfoptim)

# 生成数据集
time <- 1:40
D_obs <- round(seq(1, 200, length.out = 40))
hmod <- data.frame(time, D_obs)

sir_model <- function(t, y, parms) {
  S <- y[1]
  I <- y[2]
  R <- y[3]
  D <- y[4]
  N <- S + I + R + D
  
  beta <- parms["beta"]
  gamma <- parms["gamma"]
  alpha <- parms["alpha"]
  
  dS <- -(beta * S * I) / N
  dI <- (beta * S * I) / N - gamma * I
  dR <- alpha * gamma * I
  dD <- (1 - alpha) * gamma * I
  
  list(c(dS, dI, dR, dD))
}

# 初始参数与状态设置
par_init <- c(beta = 1, gamma = 1, alpha = 0.5)
init_sav <- c(S = 1079, I = 1, R = 0, D = 0)
times_sav <- hmod$time

opt_fun_sav <- function(par) {
  # 限制参数为正数，避免无效值
  par <- pmax(par, 1e-8)
  
  # 使用当前优化参数运行SIRD模型
  output_sav_try <- ode(
    y = init_sav,
    times = times_sav,
    func = sir_model,
    parms = c(beta = par[1], gamma = par[2], alpha = par[3])
  )
  
  output_savtry_df <- as.data.frame(output_sav_try)
  
  # 直接提取预测值，时间序列完全匹配
  rmse_try <- sqrt(mean((output_savtry_df$D - hmod$D_obs)^2))
  
  return(rmse_try)
}

# 传入3个参数的初始值，对应beta、gamma、alpha
optnm_result <- optim_nm(opt_fun_sav, start = c(0.6, 0.6, 0.5), k = 2, trace = TRUE)

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