使用nlme和lsoda拟合一阶方程

Question

使用nlme和lsoda拟合一阶方程

6

我正在尝试使用nlme和lsoda来拟合一阶常微分方程模型。基本思路是：首先定义函数以生成微分方程的解决方案：

library(deSolve)

ODE1 <- function(time, x, parms) {with(as.list(c(parms, x)), {
  import <- excfunc(time)
  dS <- import*k/tau - (S-yo)/tau 
  res <- c(dS)
  list(res)})}


solution_ODE1 = function(tau1,k1,yo1,excitation,time){
  excfunc <- approxfun(time, excitation, rule = 2)
  parms  <- c(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
  xstart = c(S = yo1)
  out <-  lsoda(xstart, time, ODE1, parms)
  return(out[,2])
}

我会生成两个ID的数据，遵循以下公式:

time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(signal = c(solution_ODE1(3,2,0.1,excitation,time)+rnorm(length(time),0,0.1),
                                   solution_ODE1(3.2,1.5,0.3,excitation,time)+rnorm(length(time),0,0.1)),
                        time = rep(time,2),
                        excitation = rep(excitation,2),
                        ID = rep(c("A","B"),each = length(time)))

这是它的外观：

library(ggplot2)
ggplot(simu_data)+
  geom_point(aes(time,signal,color = "signal"),size = 2)+
  geom_line(aes(time,excitation,color = "excitation"))+
  facet_wrap(~ID)

我正在尝试使用nlme进行拟合：

fit1 <- nlme(signal ~ solution_ODE1(damping,gain,eq,excitation,time),
             data = simu_data,
             fixed = damping + gain + eq ~1,
             random =  damping   ~ 1 ,
             groups = ~ ID,
             start = c(damping = 5, gain = 1,eq = 0))

我得到了这个错误，而我没有得到：

在评估（substitute（expr），data，enclos = parent.frame（））时出错：对象'k'未找到

traceback显示该错误来自ODE1模型，在生成值时有效。

16.    eval(substitute(expr), data, enclos = parent.frame()) 
15.    eval(substitute(expr), data, enclos = parent.frame()) 
14.    with.default(as.list(c(parms, x)), {
    import <- excfunc(time)
    dS <- import * k/tau - (S - yo)/tau
    res <- c(dS) ... 
13.    with(as.list(c(parms, x)), {
    import <- excfunc(time)
    dS <- import * k/tau - (S - yo)/tau
    res <- c(dS) ... 
12.    func(time, state, parms, ...) 
11.    Func2(times[1], y) 
10.    eval(Func2(times[1], y), rho) 
9.    checkFunc(Func2, times, y, rho) 
8.    lsoda(xstart, time, ODE1, parms) 
7.    solution_ODE1(damping, gain, eq, excitation, time) 
6.    eval(model, data.frame(data, pars)) 
5.    eval(model, data.frame(data, pars)) 
4.    eval(modelExpression[[2]], envir = nlEnv) 
3.    eval(modelExpression[[2]], envir = nlEnv) 
2.    nlme.formula(signal ~ solution_ODE1(damping, gain, eq, excitation, 
    time), data = simu_data, fixed = damping + gain + eq ~ 1, 
    random = damping ~ 1, groups = ~ID, start = c(damping = 5, 
        gain = 1, eq = 0)) 
1.    nlme(signal ~ solution_ODE1(damping, gain, eq, excitation, time), 
    data = simu_data, fixed = damping + gain + eq ~ 1, random = damping ~ 
        1, groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))

有没有人有想法，我应该如何继续？

编辑

我尝试根据mikeck的建议进行修改：

ODE1 <- function(time, x, parms) {
  import <- parms$excfunc(time)
  dS <- import*parms$k/parms$tau - (x["S"]-parms$yo)/parms$tau 
  res <- c(dS)
  list(res)}

生成数据没有问题。但使用 nlme 现在出现以下错误：

Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (0) must equal the length of the initial conditions vector (100)

并伴随着以下回溯信息：

> traceback()
10: stop(paste("The number of derivatives returned by func() (", 
        length(tmp[[1]]), ") must equal the length of the initial conditions vector (", 
        length(y), ")", sep = ""))
9: checkFunc(Func2, times, y, rho)
8: lsoda(xstart, time, ODE1, parms) at #5
7: solution_ODE1(damping, gain, eq, excitation, time)
6: eval(model, data.frame(data, pars))
5: eval(model, data.frame(data, pars))
4: eval(modelExpression[[2]], envir = nlEnv)
3: eval(modelExpression[[2]], envir = nlEnv)
2: nlme.formula(signal ~ solution_ODE1(damping, gain, eq, excitation, 
       time), data = simu_data, fixed = damping + gain + eq ~ 1, 
       random = damping ~ 1, groups = ~ID, start = c(damping = 5, 
           gain = 1, eq = 0))
1: nlme(signal ~ solution_ODE1(damping, gain, eq, excitation, time), 
       data = simu_data, fixed = damping + gain + eq ~ 1, random = damping ~ 
           1, groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))

- denis

你尝试过 nlmeODE 包吗？ - Ben Bolker

我确实在尝试。我遇到了一些困难，但也许这会奏效。我仍然很高兴能够得到关于这种奇怪行为的解决方案/说明。 - denis

我做了一些调整 - 参数应该使用list()而不是c()，并且我已经将xstart <- yo1 (然后在ODE1中直接引用x)，但我仍然收到“非法输入”消息... - Ben Bolker

1

你尝试过重新定义 ODE1()，不使用 with()，而是使用 parms$k 等吗？错误信息看起来可能是某种作用域问题引起的。 - mikeck

@mikeck 我试过了，它改变了错误信息。我编辑了我的问题。我不知道 nlme 在内部做什么，但它似乎给函数提供了初始条件的向量，从而导致错误。 - denis

@denis 这个链接可能会有帮助：https://stackoverflow.com/questions/40025139/solving-ode-with-desolve-in-r-number-of-derivatives-error - mikeck

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

在你的例子中，你的times向量并不单调递增。我认为这会影响到lsoda。在这里，时间的工作方式是什么上下文/意义？用两个组拟合随机效应模型没有太多意义。你是想将相同的曲线拟合到两个独立的时间序列中吗？

以下是一个简化的示例，进行了一些调整（不能将所有内容折叠成数字向量，否则会丢失必要的结构）：

library(deSolve)
ODE1 <- function(time, x, parms) {
    with(as.list(parms), {
        import <- excfunc(time)
        dS <- import*k/tau - (x-yo)/tau 
        res <- c(dS)
        list(res)
    })
}
solution_ODE1 = function(tau1,k1,yo1,excitation,time){
    excfunc <- approxfun(time, excitation, rule = 2)
    parms  <- list(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
    xstart = yo1
    out <-  lsoda(xstart, time, ODE1, parms)
    return(out[,2])
}
time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(time = rep(time,2),
                        excitation = rep(excitation,2))
svec <- c(damping = 3, gain = 1.75, eq = 0.2)

这是可行的：

with(c(simu_data, as.list(svec)),
     solution_ODE1(damping,gain,eq,excitation[1:50],time[1:50]))

但如果我们再加上一步（使时间重置为0），它将失败：

with(c(simu_data, as.list(svec)),
     solution_ODE1(damping,gain,eq,excitation[1:51],time[1:51]))

在执行任何积分步骤之前检测到非法输入 - 请参见书面消息。