如何用R绘制SIR模型中不同beta值下的感染人数I曲线？

你可以通过遍历指定的beta值，循环计算每个beta对应的SIR模型结果，再将每条感染曲线叠加到同一张图上。修改后的完整代码如下：

library(deSolve)

# 定义SIR模型函数
sir_model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <-  beta * I * S - gamma * I
    dR <-  gamma * I
    return(list(c(dS, dI, dR)))
  })
}

# 固定参数与初始值
initial_values <- c(S = 999, I = 1, R = 0)
time_series <- seq(0, 100)
gamma <- 0.3
beta_values <- c(0.001, 0.002, 0.003, 0.004, 0.005)
# 定义对应beta的颜色，方便区分
line_colors <- c("black", "red", "blue", "green", "purple")

# 初始化空图，设置坐标轴范围和标签（type="n"表示不绘制任何点线）
plot(
  x = time_series,
  y = numeric(length(time_series)), # 空y值
  type = "n",
  xlab = "time (days)",
  ylab = "Number of infections",
  ylim = c(0, max(sapply(beta_values, function(b) {
    # 提前计算所有beta对应的最大感染数，用来设置y轴范围
    params <- c(beta = b, gamma = gamma)
    res <- ode(y = initial_values, times = time_series, func = sir_model, parms = params)
    max(res[, "I"])
  })))
)

# 循环遍历每个beta值，计算并绘制曲线
for (i in seq_along(beta_values)) {
  current_beta <- beta_values[i]
  parameters <- c(beta = current_beta, gamma = gamma)
  
  # 运行SIR模型
  sir_result <- ode(
    y = initial_values,
    times = time_series,
    func = sir_model,
    parms = parameters
  )
  sir_result_df <- as.data.frame(sir_result)
  
  # 将当前beta对应的感染曲线添加到图中
  lines(
    x = sir_result_df$time,
    y = sir_result_df$I,
    col = line_colors[i],
    lwd = 2 # 设置线宽，更清晰
  )
}

# 添加图例，说明每条曲线对应的beta值
legend(
  "topright",
  legend = paste("β =", beta_values),
  col = line_colors,
  lwd = 2
)

关键修改说明 #

• 先定义所有要测试的beta_values向量，以及对应的颜色向量line_colors，方便区分不同曲线
• 用plot(type="n")初始化空图，通过sapply提前计算所有beta对应的最大感染数，设置合适的y轴范围，避免部分曲线超出显示范围
• 在for循环中，每次替换parameters里的beta值，重新运行SIR模型，再用lines()函数将曲线叠加到已有图上（而不是重新开新图）
• 最后添加图例，让每条曲线对应的beta值一目了然

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