基于R语言deSolve绘制多参数组合的疾病模型600天时序曲线

下面是针对你需求的完整实现代码，涵盖所有参数组合并生成600天的趋势曲线：

1. 加载依赖包 #

首先安装并加载deSolve（求解微分方程）和ggplot2（绘图）：

# 首次运行需安装包
install.packages(c("deSolve", "ggplot2"))

# 加载包
library(deSolve)
library(ggplot2)

2. 定义疾病模型 #

这里采用包含**易感(S)、有症状感染(I)、无症状感染(A)、康复(R)**的SIA-R模型（适配你提供的beta_s和gamma_a参数），微分方程逻辑如下：

• 易感人群因接触感染者减少
• 有症状感染者由易感人群转化，同时以固定速率康复
• 无症状感染者传播力设为有症状的70%（可按需调整），以gamma_a速率康复
• 康复人群由两类感染者转化而来

模型函数定义：

disease_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta_s * S * (I + A)
    dI <- beta_s * S * (I + A) - gamma_i * I
    dA <- 0.7 * beta_s * S * (I + A) - gamma_a * A
    dR <- gamma_i * I + gamma_a * A
    
    return(list(c(dS, dI, dA, dR)))
  })
}

3. 设置参数与初始条件 #

# 所有参数交叉组合
beta_s_values <- c(0.1, 0.3, 0.5)
gamma_a_values <- c(0.2, 0.3, 0.4)
params_grid <- expand.grid(beta_s = beta_s_values, gamma_a = gamma_a_values)

# 固定其他参数（可根据你的模型调整）
gamma_i <- 0.1  # 有症状感染者康复率

# 初始状态：总人口10000，1个有症状感染者，其余为易感
initial_state <- c(S = 9999, I = 1, A = 0, R = 0)

# 时间序列：0到600天
times <- seq(0, 600, by = 1)

4. 运行模型并收集结果 #

# 初始化结果容器
all_results <- list()

# 循环求解每个参数组合的模型
for (i in 1:nrow(params_grid)) {
  params <- c(params_grid[i, ], gamma_i = gamma_i)
  output <- ode(y = initial_state, times = times, func = disease_model, parms = params)
  # 转换为数据框并标记参数
  df <- as.data.frame(output)
  df$beta_s <- params_grid$beta_s[i]
  df$gamma_a <- params_grid$gamma_a[i]
  all_results[[i]] <- df
}

# 合并所有结果
combined_df <- do.call(rbind, all_results)

5. 绘制多组曲线 #

采用分面图区分不同gamma_a值，颜色区分beta_s值，同时展示有症状和无症状感染者的趋势：

ggplot(combined_df, aes(x = time, color = factor(beta_s))) +
  geom_line(aes(y = I), linewidth = 1) +
  geom_line(aes(y = A), linewidth = 1, linetype = "dashed") +
  facet_wrap(~gamma_a, labeller = label_both) +
  labs(
    x = "时间（天）",
    y = "人数",
    color = "β_s",
    title = "不同参数组合下的疾病感染趋势（600天）",
    subtitle = "实线：有症状感染者(I)；虚线：无症状感染者(A)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

自定义调整建议 #

• 若你的模型结构不同（如包含暴露期E），直接修改disease_model中的微分方程即可
• 如需展示易感/康复人群，添加geom_line(aes(y = S))或geom_line(aes(y = R))
• 可将facet_wrap替换为facet_grid(gamma_a ~ beta_s)，实现行列分别展示两个参数
• 颜色、线型、主题等可根据论文示例图的样式自由调整

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