求解耦合ODE系统的R代码绘图异常问题排查

问题描述 #

作为R语言新手，需求解并绘制如下耦合ODE系统的S、I、R、D随时间变化曲线：

$$\frac{dS}{dT} = -\beta(1-\mu)S(t)\frac{I(t)}{1-D(t)}$$
$$\frac{dI}{dT} = \beta(1-\mu)S(t)\frac{I(t)}{1-D(t)}-\delta I(t)-\phi I(t)$$
$$\frac{dR}{dT}=\delta I(t)$$
$$\frac{dD}{dT}=\phi I(t)$$

参考deSolve手册编写代码后，生成绘图与预期不符，调整参数无明显改善。代码如下：

install.packages("deSolve")
library("deSolve")

sird_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
dS <- -beta * (1- mu) * S * (I/(1-D))
dI <- -beta * (1- mu) * S * (I/(1-D))-sigma * I - theta * I
dR <- sigma * I 
dD <- theta * I 
return(list(c(dS, dI, dR, dD)))
})
} 

parameter_values <- c( 
beta = 0.95, #infectious contact rate (/person/day)
sigma = 0.5, #recovery rate (/day)
theta = 0.043, #mortality rate (/day)
mu = 0.1 #adjustment for efficacy of countermeasures
)

inital_values <- c(
S = 126000000, #susceptible population (adjust based on your data)
I = 1, #infected on day 1 of data 
R = 0, #recovered on day 1
D = 0 #dead on day 1
)


time_values <- seq(0, 10) #days(adjust based on your data)
ls()

sird_values_1 <- ode(
y = inital_values,
times = time_values,
func = sird_equations,
parms = parameter_values 
)

sird_values_1 <- as.data.frame(sird_values_1)
with(sird_values_1, {
plot(time, S, type = "l", col = "blue",
xlab = "time (days)", ylab = "# of people")
lines(time, I, col = "black")
lines(time, R, col = "red")
lines(time, D, col = "orange")
})

legend("right", c("susceptible", "infectious", "recovered", "dead"),
col = c("blue", "black", "red", "orange"), lty = 1, bty = "n")

错误排查与修正 #

1. 核心符号错误（dI方程） #

原ODE中，$\frac{dI}{dT}$的第一项为正号——易感者被感染转化为感染者，但代码里写成了负号，直接导致感染者数量从初始值开始持续下降，完全违背疫情传播规律。

修正后的dI计算：

dI <- beta * (1- mu) * S * (I/(1-D)) - sigma * I - theta * I

2. 时间范围过短 #

代码中仅模拟10天，此时感染者还未达到传播峰值就结束模拟，无法呈现预期的曲线形态。建议延长时间范围：

time_values <- seq(0, 50) # 模拟50天

3. 绘图可视化问题 #

易感者S的初始值远大于I/R/D，直接绘制会导致后三者的曲线被压缩到X轴附近几乎不可见。推荐两种优化方式：

• 方式一：双Y轴绘图

with(sird_values_1, {
  # 绘制S的左侧Y轴
  plot(time, S, type = "l", col = "blue",
       xlab = "time (days)", ylab = "# of susceptible people",
       ylim = c(0, max(S)))
  # 添加右侧Y轴用于I/R/D
  par(new = TRUE)
  plot(time, I, type = "l", col = "black",
       xlab = "", ylab = "", axes = FALSE,
       ylim = c(0, max(I)*1.1))
  lines(time, R, col = "red")
  lines(time, D, col = "orange")
  axis(side = 4)
  mtext("# of infected/recovered/dead people", side = 4, line = 3)
})
legend("topright", c("susceptible", "infectious", "recovered", "dead"),
       col = c("blue", "black", "red", "orange"), lty = 1, bty = "n")

• 方式二：使用相对比例
  计算各群体占总人口的比例（总人口为S+I+R+D），使所有曲线都在0-1范围内，便于观察趋势。

完整修正后的代码 #

install.packages("deSolve")
library("deSolve")

sird_equations <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * (1 - mu) * S * (I / (1 - D))
    # 修正dI的符号
    dI <- beta * (1 - mu) * S * (I / (1 - D)) - sigma * I - theta * I
    dR <- sigma * I 
    dD <- theta * I 
    return(list(c(dS, dI, dR, dD)))
  })
} 

parameter_values <- c( 
  beta = 0.95,    # 传染接触率
  sigma = 0.5,    # 康复率
  theta = 0.043,  # 死亡率
  mu = 0.1        # 防控措施调整系数
)

inital_values <- c(
  S = 126000000,
  I = 1,
  R = 0,
  D = 0
)

# 延长模拟时间
time_values <- seq(0, 50)

sird_values_1 <- ode(
  y = inital_values,
  times = time_values,
  func = sird_equations,
  parms = parameter_values 
)

sird_values_1 <- as.data.frame(sird_values_1)

# 双Y轴绘图优化可视化
with(sird_values_1, {
  plot(time, S, type = "l", col = "blue",
       xlab = "time (days)", ylab = "# of susceptible people",
       lwd = 2, ylim = c(0, max(S)))
  par(new = TRUE)
  plot(time, I, type = "l", col = "black",
       xlab = "", ylab = "", axes = FALSE,
       lwd = 2, ylim = c(0, max(I)*1.2))
  lines(time, R, col = "red", lwd = 2)
  lines(time, D, col = "orange", lwd = 2)
  axis(side = 4, col = "gray")
  mtext("# of infected/recovered/dead people", side = 4, line = 3, col = "gray")
})
legend("topright", 
       c("susceptible", "infectious", "recovered", "dead"),
       col = c("blue", "black", "red", "orange"), 
       lty = 1, lwd = 2, bty = "n")

内容的提问来源于stack exchange，提问作者Tom FitzGerald-Jones