如何以tidyverse风格为deSolve包组织系列模拟场景

问题背景 #

需要为deSolve包构建状态与参数可变的系列模拟场景，要求：

• 适配管道式编程风格
• 不使用嵌套数据框，采用清晰的表格结构
• 输出直接支持ggplot可视化
• 替代原有base-R风格（apply+列表）或带for循环/do.call的过渡方案

完整解决方案 #

以下以经典Logistic种群增长模型为例，演示纯Tidyverse风格的实现流程：

1. 加载依赖包

library(deSolve)
library(tidyverse)

2. 定义微分方程模型

编写标准的deSolve模型函数，输入时间、状态变量和参数，返回状态变量的变化率：

logistic_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N * (1 - N / K)  # Logistic增长核心方程
    return(list(dN))
  })
}

3. 创建多场景参数表格

用expand_grid生成所有参数组合的场景表格，每个场景对应一行，参数清晰可追溯：

scenarios <- expand_grid(
  scenario_id = 1:6,          # 场景唯一标识，方便后续区分
  r = c(0.1, 0.2, 0.3),       # 可变参数：种群增长率
  K = c(100, 200),            # 可变参数：环境容纳量
  initial_N = 10,             # 固定初始种群数量
  times = list(seq(0, 50, 1)) # 时间序列用list存储，便于传递给模拟函数
)

4. 编写单场景模拟函数

该函数接收单个场景的参数，执行模拟并返回tidy格式的结果（包含场景标识、参数和模拟数据）：

simulate_scenario <- function(scenario_id, r, K, initial_N, times) {
  # 定义初始状态与参数向量
  state <- c(N = initial_N)
  parameters <- c(r = r, K = K)
  
  # 执行ode模拟并转换为tidy数据框
  ode(
    y = state,
    times = times,
    func = logistic_model,
    parms = parameters
  ) %>%
    as_tibble() %>%
    # 添加场景标识与参数信息，方便后续分组可视化
    mutate(scenario_id = scenario_id, r = r, K = K)
}

5. 批量执行多场景模拟

用purrr::pmap替代for循环/do.call，按行传递场景参数并批量执行，最终合并为扁平的tidy数据框：

simulation_results <- scenarios %>%
  # pmap自动按行匹配参数与函数输入
  mutate(sim_data = pmap(., simulate_scenario)) %>%
  # 展开模拟结果列表，移除冗余的times列
  unnest(sim_data) %>%
  select(-times)

6. 直接适配ggplot可视化

由于结果是标准tidy数据框，可直接传入ggplot进行分组可视化：

ggplot(simulation_results, aes(x = time, y = N, color = factor(K), linetype = factor(r))) +
  geom_line(linewidth = 1) +
  labs(
    x = "时间",
    y = "种群数量",
    color = "环境容纳量(K)",
    linetype = "增长率(r)",
    title = "Logistic模型多场景模拟对比"
  ) +
  theme_minimal()

方案优势 #

• 管道式流程：从场景定义→模拟执行→可视化全程采用Tidyverse管道，逻辑连贯易读
• 清晰的表格结构：场景参数与模拟结果均为扁平数据框，教学示例中便于理解参数组合与结果对应关系
• 无特殊技巧：用pmap替代for循环和do.call，代码简洁且符合Tidyverse编程范式
• 无缝适配可视化：输出直接满足ggplot对数据格式的要求，无需额外转换

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