使用R包fitode拟合简单ODE时出现optim错误的原因排查

问题 #

使用R包fitode拟合肿瘤生长微分方程模型时，遇到如下错误：

Error in optim(par = c(log.r1 = 0, log.X0 = 0), fn = function (p)  : 
  initial value in 'vmmin' is not finite

即使简化到最简复现案例，错误依然存在：

• 调用包并模拟指数增长的细胞计数数据：

library(fitode)
days <- c(1, 3, 6, 8, 9, 11, 13, 14)
data_test <- data.frame(day = days, size = rpois(length(days), exp(days)))

• 定义指数增长模型：

model_test <- odemodel(
  name="test",
  model=list(
    X1 ~ r1 * X1
  ),
  observation=list(
    size ~ dpois(lambda=X1)
  ),
  initial=list(
    X1 ~ X0
  ),
  par=c("r1", "X0")
)

• 已知真实参数r1=1、X0=1，模拟数据正常，但用该初始值拟合时仍报错：

start_test <- c(r1=1, X0=1)
fit_test <- fitode(
  model=model_test,
  data=data_test,
  start=start_test
)

解决方案 #

问题核心是调用fitode时未指定时间变量，程序无法正确匹配观测数据的时间点与ODE求解时间，进而引发数值计算错误。

修改拟合代码，添加time参数指定数据中的时间列（即day）即可解决：

fit_test <- fitode(
  model=model_test,
  data=data_test,
  start=start_test,
  time = "day"  # 指定时间变量
)

也可以在定义odemodel时提前指定时间变量，让模型结构更清晰：

model_test <- odemodel(
  name="test",
  model=list(
    X1 ~ r1 * X1
  ),
  observation=list(
    size ~ dpois(lambda=X1)
  ),
  initial=list(
    X1 ~ X0
  ),
  par=c("r1", "X0"),
  time="day"  # 提前定义时间变量
)

修改后即可正常拟合，得到符合预期的参数估计结果。

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