如何在ODE中搭建蒙特卡洛模拟估计不确定性？适配R3.4.4环境

我想给我的常微分方程（ODE）做蒙特卡洛模拟，查看结果的方差，还要绘制时间序列解的可视化图。但搭建模拟时遇到了困难，找不到适配我问题的方案。

之前用过pksensi包，但当前R版本是3.4.4没法更新，用不了这个包。也看过FME包的modCRL函数示例，但搞不懂逻辑和输出。还查过简单模型的蒙特卡洛示例，还是没解决问题。

附上我写的报错代码，我知道代码有问题，但希望思路能被理解。除了用pksensi，我从没自己搭过蒙特卡洛模拟，求帮忙指导怎么搭建。

library(deSolve)

months <- seq(0, 144, 1)
all.df <- data.frame(months, c.weight = c(3.5,  4,  5,  5,  6, 7,  7,  7,  8,  8,  8,  9,  9,  9,  9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 37, 37, 37, 38, 38, 38, 39, 39, 39, 40, 40, 40, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 51, 51, 52, 54, 55))

model <- function(times, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    if (times <= 5) {
      volume <- (((0.2 * times) * c.weight[times+1]) * 30)
      transferred <- pm * volume
    } else if (times > 5 & times <= 14) {
      volume <- ((0.1 * c.weight[times+1]) * 30)
      transferred <- pm * volume
    } else {
      transferred <- 0
    }
    
    intake.c <- (dose * c.weight[times+1])
    elimination.c <- concentration * vd * c.weight[times+1] * log(2) / (halflife * 12)
    
    concentration <- (intake.c + transferred - elimination.c) / (vd * c.weight[times+1])
    
    list(c(concentration))
    
  })
}

#Following you can see my attempt to perform a Monte Carlo simulation with specific parameters using rnorm()
#The state can vary as well as the different parameters in params
state <- c(concentration = 0.5 * rnorm(1000, mean = 0.7, sd = 0.2)) #parameter to be tested
params <- list(vd = rnorm(1000, mean = 0.2, sd = 0.02), #parameter to be tested
               pm = rnorm(1000, mean = 0.05, sd = 0.03), #parameter to be tested
               halflife = rnorm(1000, mean = 4, sd = 2), #parameter to be tested
               dose = 0.001,
               c.weight = all.df$c.weight)
out <- as.data.frame(ode(y = state, times = months, func = model, parms = params))

原代码核心问题 #

1. ODE模型逻辑错误：deSolve的ode函数要求模型返回状态变量的导数（变化率），而非直接计算当前浓度值，你的代码直接赋值concentration，不符合要求。
2. 批量模拟方式错误：ode一次只能处理一组参数和初始状态，无法直接传入1000组参数，需循环或批量处理每组模拟。
3. 参数匹配错误：每组模拟需对应一组参数（vd、pm、halflife）和一个初始浓度，原代码的参数是长度1000的向量，无法被ode正确识别。

修正后的完整代码 #

library(deSolve)
library(ggplot2)

# 时间序列与体重数据
months <- seq(0, 144, 1)
all.df <- data.frame(
  months, 
  c.weight = c(3.5,  4,  5,  5,  6, 7,  7,  7,  8,  8,  8,  9,  9,  9,  9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 37, 37, 37, 38, 38, 38, 39, 39, 39, 40, 40, 40, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 51, 51, 52, 54, 55)
)

# 修正后的ODE模型：返回浓度的变化率dC/dt
model <- function(times, state, parameters) {
  with(as.list(c(state, parameters)), {
    # 根据当前时间获取对应体重
    t_idx <- which(all.df$months == times)
    current_weight <- all.df$c.weight[t_idx]
    
    # 计算transferred
    if (times <= 5) {
      volume <- (0.2 * times) * current_weight * 30
      transferred <- pm * volume
    } else if (times > 5 & times <= 14) {
      volume <- 0.1 * current_weight * 30
      transferred <- pm * volume
    } else {
      transferred <- 0
    }
    
    # 计算摄入、消除的量
    intake_c <- dose * current_weight
    elimination_c <- concentration * vd * current_weight * log(2) / (halflife * 12)
    
    # 关键：返回浓度的变化率dC/dt
    d_concentration <- (intake_c + transferred - elimination_c) / (vd * current_weight)
    
    list(c(d_concentration))
  })
}

# 蒙特卡洛模拟设置
n_sim <- 1000  # 模拟次数

# 生成每组模拟的参数：每行一组参数
mc_params <- data.frame(
  vd = pmax(rnorm(n_sim, mean = 0.2, sd = 0.02), 0.01), # 避免负值
  pm = pmax(rnorm(n_sim, mean = 0.05, sd = 0.03), 0.001),
  halflife = pmax(rnorm(n_sim, mean = 4, sd = 2), 0.5),
  dose = 0.001  # 固定参数
)

# 生成每组模拟的初始浓度
initial_conc <- 0.5 * rnorm(n_sim, mean = 0.7, sd = 0.2)

# 批量运行蒙特卡洛模拟
mc_results <- lapply(1:n_sim, function(i) {
  # 获取第i组参数和初始状态
  params_i <- as.list(mc_params[i, ])
  state_i <- c(concentration = initial_conc[i])
  
  # 运行ODE
  out <- ode(y = state_i, times = months, func = model, parms = params_i)
  
  # 转成数据框并添加模拟编号
  df <- as.data.frame(out)
  df$sim_id <- i
  return(df)
})

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

# 计算每个时间点的均值和方差
summary_stats <- aggregate(concentration ~ months, data = all_results, 
                          FUN = function(x) c(mean = mean(x), var = var(x)))
summary_stats <- do.call(data.frame, summary_stats)
colnames(summary_stats) <- c("months", "mean_conc", "var_conc")

# 可视化1：所有模拟曲线+均值曲线
ggplot(all_results, aes(x = months, y = concentration, group = sim_id)) +
  geom_line(alpha = 0.1, color = "blue") +
  geom_line(data = summary_stats, aes(y = mean_conc), color = "red", size = 1) +
  labs(x = "月份", y = "浓度", title = "ODE蒙特卡洛模拟结果") +
  theme_bw()

# 可视化2：各时间点浓度方差变化
ggplot(summary_stats, aes(x = months, y = var_conc)) +
  geom_line(color = "darkgreen", size = 1) +
  labs(x = "月份", y = "浓度方差", title = "浓度方差随时间变化") +
  theme_bw()

代码说明 #

1. 修正ODE模型：将直接计算浓度改为计算浓度的变化率d_concentration，这是deSolve的核心要求。
2. 批量模拟处理：用lapply循环处理每组参数和初始值，每次运行一次ode，并给结果添加模拟编号方便后续合并。
3. 结果统计：用aggregate计算每个时间点的浓度均值和方差，实现你查看方差的需求。
4. 可视化：用ggplot2绘制所有模拟曲线（调低透明度避免重叠）和均值曲线，同时绘制方差随时间的变化趋势。

注意事项 #

• 参数生成时用pmax限制最小值，避免出现不符合实际意义的负值（如半衰期不能为负）。
• 若R3.4.4安装ggplot2有问题，可替换为基础绘图函数（如plot、lines）实现可视化。

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