如何在R中同时拟合多时间序列的Logistic增长ODE模型

要实现两组初始密度数据的统一参数拟合，核心是修改代价函数，让它同时计算两组数据的模拟值与观测值的误差总和，再基于总误差进行参数优化。以下是具体实现步骤：

1. 准备数据与模型 #

先加载所需包并确认数据：

library(ggplot2)
library(FME)

# 观测数据
N0 <- c(5, 200)
time <- seq(0, 200, 10)
df <- data.frame(
  time = rep(time, 2),
  N0 = rep(N0, each = 21),
  N = c(5.0,12.6,27.9,51.9,74.7,88.7,94.6,98.1,100.3,100.2,99.3,97.8,99.3,97.9,98.7,99.6,99.3,99.1,99.9,99.7,98.1,200.0,123.7,107.5,101.1,99.5,100.7,98.4,99.7,99.4,101.1,99.2,99.2,100.8,99.0,100.0,100.2,99.7,99.3,100.7,99.6,99.7)
)

# 定义Logistic ODE模型
logistic.ode <- function(t, N, parms) {
  r <- parms["r"]
  K <- parms["K"]
  dN <- r * N * (1 - N/K)
  list(dN)
}

2. 修改代价函数，同时处理两组数据 #

在代价函数中分别模拟两个初始条件的时间序列，然后将两组的误差合并：

# 定义同时计算两组误差的代价函数
mod.cost_unified <- function(parms) {
  # 模拟N0=5的种群动态
  sim_N0_5 <- ode(y = c(N = N0[1]), func = logistic.ode, parms = parms, times = time)
  cost1 <- modCost(sim_N0_5, df[df$N0 == 5, c("time", "N")])
  
  # 模拟N0=200的种群动态
  sim_N0_200 <- ode(y = c(N = N0[2]), func = logistic.ode, parms = parms, times = time)
  cost2 <- modCost(sim_N0_200, df[df$N0 == 200, c("time", "N")])
  
  # 返回两组误差的总和
  return(cost1 + cost2)
}

3. 统一参数拟合 #

设置初始参数猜测，调用modFit进行拟合：

# 初始参数猜测（根据观测数据调整，比如K明显接近100）
parms.ini <- c(r = 0.1, K = 100)

# 拟合统一参数，设置参数上下限（符合生物学意义）
mod.fit_unified <- modFit(
  f = mod.cost_unified,
  p = parms.ini,
  lower = c(r = 0, K = 0),
  upper = c(r = 0.5, K = 150)
)

# 查看拟合结果
summary(mod.fit_unified)

4. 验证拟合效果 #

模拟拟合后的两组数据，和观测值对比绘图：

# 用拟合得到的参数模拟两组数据
sim_5 <- ode(y = c(N = 5), func = logistic.ode, parms = mod.fit_unified$par, times = time)
sim_200 <- ode(y = c(N = 200), func = logistic.ode, parms = mod.fit_unified$par, times = time)

# 转换为数据框方便绘图
sim_df <- rbind(
  data.frame(time = sim_5[,1], N = sim_5[,2], N0 = 5, type = "simulated"),
  data.frame(time = sim_200[,1], N = sim_200[,2], N0 = 200, type = "simulated")
)

# 合并观测数据和模拟数据
plot_df <- rbind(
  df |> mutate(type = "observed"),
  sim_df
)

# 绘图对比
ggplot(plot_df, aes(time, N, color = factor(N0), shape = type)) +
  geom_point(size = 2) +
  geom_line(data = sim_df) +
  theme_bw() +
  labs(color = "初始密度N0", shape = "数据类型")

关键说明 #

• 代价函数中通过cost1 + cost2合并两组误差，确保拟合过程同时考虑两组数据的信息，得到统一的r和K。
• 初始参数猜测建议结合观测数据调整，能提升拟合效率和稳定性。
• 参数上下限可依据生物学意义设置（比如r、K均为非负值）。

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