r - 具有随机效应和 lsoda 的非线性回归

标签 r ode differential-equations non-linear-regression nlme

我面临一个我无法解决的问题。我想使用 nlmenlmODE使用具有固定系数(阻尼振荡器)的二阶微分方程的解作为模型,执行具有随机效应的非线性回归。
我设法使用 nlme用简单的模型,但似乎使用deSolve生成微分方程的解会导致问题。下面是一个例子,以及我面临的问题。
数据和功能
这是使用 deSolve 生成微分方程解的函数:

library(deSolve)
ODE2_nls <- function(t, y, parms) {
  S1 <- y[1]
  dS1 <- y[2]
  dS2 <- dS1
  dS1 <- - parms["esp2omega"]*dS1  - parms["omega2"]*S1 + parms["omega2"]*parms["yeq"]
  res <- c(dS2,dS1)
  list(res)}

solution_analy_ODE2 = function(omega2,esp2omega,time,y0,v0,yeq){
  parms  <- c(esp2omega = esp2omega,
              omega2 = omega2,
              yeq = yeq)
  xstart = c(S1 =  y0, dS1 = v0)
  out <-  lsoda(xstart, time, ODE2_nls, parms)
  return(out[,2])
}
我可以为给定的周期和阻尼因子生成一个解决方案,例如这里的周期为 20,阻尼为 0.2:

# small example:
time <- 1:100
period <- 20 # period of oscillation
amort_factor <- 0.2
omega <- 2*pi/period # agular frequency
oscil <- solution_analy_ODE2(omega^2,amort_factor*2*omega,time,1,0,0)
plot(time,oscil)

enter image description here
现在我生成一个由 10 个人组成的面板,具有随机的起始阶段(即不同的起始位置和速度)。目标是使用 执行非线性回归。对起始值的随机影响
library(data.table)
# generate panel
Npoint <- 100 # number of time poitns
Nindiv <- 10 # number of individuals
period <- 20 # period of oscillation
amort_factor <- 0.2
omega <- 2*pi/period # agular frequency
# random phase
phase <- sample(seq(0,2*pi,0.01),Nindiv)
# simu data:
data_simu <- data.table(time = rep(1:Npoint,Nindiv), ID = rep(1:Nindiv,each = Npoint))

# signal generation
data_simu[,signal := solution_analy_ODE2(omega2 = omega^2,
                                         esp2omega = 2*0.2*omega,
                                         time = time,
                                         y0 = sin(phase[.GRP]),
                                         v0 = omega*cos(phase[.GRP]),
                                         yeq = 0)+ 
            rnorm(.N,0,0.02),by = ID]
如果我们看一下,我们有一个正确的数据集:
library(ggplot2)
ggplot(data_simu,aes(time,signal,color = ID))+
  geom_line()+
  facet_wrap(~ID)
enter image description here
问题
使用 nlme
使用 nlme使用类似的语法处理更简单的示例(不使用 deSolve 的非线性函数),我尝试了:
fit <- nlme(model = signal ~ solution_analy_ODE2(esp2omega,omega2,time,y0,v0,yeq), 
     data = data_simu,
     fixed = esp2omega + omega2 + y0 + v0 + yeq ~ 1,
     random = y0 ~ 1 ,
     groups = ~ ID, 
     start = c(esp2omega = 0.08, 
               omega2 = 0.04,
               yeq = 0,
               y0 = 1,
               v0 = 0))
我得到:

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


追溯:
12. stop(paste("The number of derivatives returned by func() (", length(tmp[[1]]), ") must equal the length of the initial conditions vector (", length(y), ")", sep = ""))
11. checkFunc(Func2, times, y, rho)
10. lsoda(xstart, time, ODE2_nls, parms)
9. solution_analy_ODE2(omega2, esp2omega, time, y0, v0, yeq)
.
.
我看起来像 nlme正在尝试将起始条件向量传递给 solution_analy_ODE2 ,并导致 checkFunc 中的错误来自 lasoda .
我尝试使用 nlsList :
test <- nlsList(model = signal ~ solution_analy_ODE2(omega2,esp2omega,time,y0,v0,yeq) | ID, 
        data = data_simu, 
        start = list(esp2omega = 0.08, omega2 = 0.04,yeq = 0,
                     y0 = 1,v0 = 0),
        control = list(maxiter=150, warnOnly=T,minFactor = 1e-10), 
        na.action = na.fail, pool = TRUE)
head(test)
Call:
  Model: signal ~ solution_analy_ODE2(omega2, esp2omega, time, y0, v0, yeq) | ID 
   Data: data_simu 

Coefficients:
   esp2omega     omega2           yeq         y0          v0
1  0.1190764 0.09696076  0.0007577956 -0.1049423  0.30234654
2  0.1238936 0.09827158 -0.0003463023  0.9837386  0.04773775
3  0.1280399 0.09853310 -0.0004908579  0.6051663  0.25216134
4  0.1254053 0.09917855  0.0001922963 -0.5484005 -0.25972829
5  0.1249473 0.09884761  0.0017730823  0.7041049  0.22066652
6  0.1275408 0.09966155 -0.0017522320  0.8349450  0.17596648
我们可以看到,非线性拟合在单个信号上效果很好。现在,如果我想对具有随机效应的数据集进行回归,语法应该是:
fit <- nlme(test, 
     random = y0 ~ 1 ,
     groups = ~ ID, 
     start = c(esp2omega = 0.08, 
               omega2 = 0.04,
               yeq = 0,
               y0 = 1,
               v0 = 0))
但我得到了完全相同的错误信息。
然后我尝试使用 nlmODE ,根据 Bne Bolker 对 a similar question I asked some years ago 的评论
使用 nlmMODE
library(nlmeODE)
datas_grouped <- groupedData( signal ~ time | ID, data = data_simu, 
                              labels = list (x = "time", y = "signal"), 
                              units = list(x ="arbitrary", y = "arbitrary"))

modelODE <- list( DiffEq = list(dS2dt = ~ S1,
                                dS1dt = ~ -esp2omega*S1  - omega2*S2 + omega2*yeq),
                  ObsEq = list(yc = ~ S2),
                  States = c("S1","S2"),
                  Parms = c("esp2omega","omega2","yeq","ID"), 
                  Init = c(y0 = 0,v0 = 0))

resnlmeode = nlmeODE(modelODE, datas_grouped) 
assign("resnlmeode", resnlmeode, envir = .GlobalEnv)
#Fitting with nlme the resulting function
model <- nlme(signal ~ resnlmeode(esp2omega,omega2,yeq,time,ID), 
              data = datas_grouped, 
              fixed = esp2omega + omega2 + yeq + y0 + v0  ~ 1, 
              random = y0 + v0 ~1,
              start = c(esp2omega = 0.08, 
                        omega2 = 0.04,
                        yeq = 0,
                        y0 = 0,
                        v0 = 0)) # 

我得到错误:

Error in resnlmeode(esp2omega, omega2, yeq, time, ID) : object 'yhat' not found


在这里我不明白错误来自哪里,也不知道如何解决它。
问题
  • 你能重现这个问题吗?
  • 有没有人有办法解决这个问题,使用nlmenlmODE ?
  • 如果没有,是否有使用其他软件包的解决方案?我看到了nlmixr ( https://cran.r-project.org/web/packages/nlmixr/index.html ),但我不知道,安装很复杂,最近从 CRAN 中删除了

  • 编辑
    @tpetzoldt 提出了一个很好的调试方法 nlme行为,这让我很惊讶。这是一个具有非线性函数的工作示例,其中我生成了一组 5 个个体,其中随机参数在个体之间变化:
    reg_fun = function(time,b,A,y0){
      cat("time : ",length(time)," b :",length(b)," A : ",length(A)," y0: ",length(y0),"\n")
      out <- A*exp(-b*time)+(y0-1)
      cat("out : ",length(out),"\n")
      tmp <- cbind(b,A,y0,time,out)
      cat(apply(tmp,1,function(x) paste(paste(x,collapse = " "),"\n")),"\n")
      return(out)
    }
    
    time <- 0:10*10
    ramdom_y0 <- sample(seq(0,1,0.01),10)
    Nid <- 5
    data_simu <- 
    data.table(time = rep(time,Nid),
               ID = rep(LETTERS[1:Nid],each = length(time)) )[,signal := reg_fun(time,0.02,2,ramdom_y0[.GRP]) + rnorm(.N,0,0.1),by = ID]
    
    函数中的猫在这里给出:
    time :  11  b : 1  A :  1  y0:  1 
    out :  11 
    0.02 2 0.64 0 1.64 
     0.02 2 0.64 10 1.27746150615596 
     0.02 2 0.64 20 0.980640092071279 
     0.02 2 0.64 30 0.737623272188053 
     0.02 2 0.64 40 0.538657928234443 
     0.02 2 0.64 50 0.375758882342885 
     0.02 2 0.64 60 0.242388423824404 
     0.02 2 0.64 70 0.133193927883213 
     0.02 2 0.64 80 0.0437930359893108 
     0.02 2 0.64 90 -0.0294022235568269 
     0.02 2 0.64 100 -0.0893294335267746
    .
    .
    .
    
    现在我用 nlme :
    nlme(model = signal ~ reg_fun(time,b,A,y0), 
         data = data_simu,
         fixed = b + A + y0 ~ 1,
         random = y0 ~ 1 ,
         groups = ~ ID, 
         start = c(b = 0.03, A = 1,y0 = 0))
    
    我得到:
    time :  55  b : 55  A :  55  y0:  55 
    out :  55 
    0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     
    time :  55  b : 55  A :  55  y0:  55 
    out :  55 
    0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
     0.03 1 0 0 0 
     0.03 1 0 10 -0.259181779318282 
     0.03 1 0 20 -0.451188363905974 
     0.03 1 0 30 -0.593430340259401 
     0.03 1 0 40 -0.698805788087798 
     0.03 1 0 50 -0.77686983985157 
     0.03 1 0 60 -0.834701111778413 
     0.03 1 0 70 -0.877543571747018 
     0.03 1 0 80 -0.909282046710588 
     0.03 1 0 90 -0.93279448726025 
     0.03 1 0 100 -0.950212931632136 
    ...
    
    所以nlme绑定(bind)5次(个体的数量)时间向量并将其传递给函数,参数重复相同的时间。这当然与 lsoda 的方式不兼容我的功能有效。

    最佳答案

    似乎使用错误的参数调用了 ode 模型,因此它得到了一个包含 2000 个状态变量而不是 2 个的向量。尝试以下操作来查看问题:

    ODE2_nls <- function(t, y, parms) {
      cat(length(y),"\n") # <----
      S1 <- y[1]
      dS1 <- y[2]
      dS2 <- dS1
      dS1 <- - parms["esp2omega"]*dS1  - parms["omega2"]*S1 + parms["omega2"]*parms["yeq"]
      res <- c(dS2,dS1)
      list(res)
    }
    
    编辑 :我认为分析函数有效,因为它是矢量化的,因此您可以尝试通过迭代 ode 模型或(更好地)在内部使用向量作为状态变量来对 ode 函数进行矢量化。如ode在求解具有几个 100k 方程的系统时速度很快,2000 应该是可行的。
    我猜想 nlme 中的状态和参数作为向量传递。 ode 模型的状态变量是一个“长”向量,参数可以实现为一个列表。
    这是一个示例(已编辑,现在将参数作为列表):
    ODE2_nls <- function(t, y, parms) {
      #cat(length(y),"\n")
      #cat(length(parms$omega2))
      ndx <- seq(1, 2*N-1, 2)
      S1  <- y[ndx]
      dS1 <- y[ndx + 1]
      dS2 <- dS1
      dS1 <- - parms$esp2omega * dS1  - parms$omega2 * S1 + parms$omega2 * parms$yeq
      res <- c(dS2, dS1)
      list(res)
    }
    
    solution_analy_ODE2 = function(omega2, esp2omega, time, y0, v0, yeq){
      parms  <- list(esp2omega = esp2omega, omega2 = omega2, yeq = yeq)
      xstart = c(S1 =  y0, dS1 = v0)
      out <-  ode(xstart, time, ODE2_nls, parms, atol=1e-4, rtol=1e-4, method="ode45")
      return(out[,2])
    }
    
    然后设置(或计算)方程的数量,例如N <- 1分别N <-1000在通话之前。
    该模型通过这种方式运行,然后运行在数值问题中,但这是另一个故事......
    然后您可以尝试使用另一个 ode 求解器(例如 vode ),设置 atolrtol要降低值,请调整 nmle的优化参数,使用框约束……等等,和非线性优化中的一样。

    关于r - 具有随机效应和 lsoda 的非线性回归,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/65778251/

    相关文章:

    r - 在对称数据帧中删除行和列满足条件

    python - 具有离散值的 ODE 积分

    python - 如何使用Python求解具有边界条件的微分方程,其中一个是不等式

    matlab - 不带计数器的自定义周期函数

    matlab - 积分微分方程

    r - 在ggplot中设置日期范围

    r - write.csv()不等大小的data.frames列表

    android - 减去 R 中数据框中的列。 A-B=已校正 A。这些是传感器中从数据帧 B 寻求校正的抽象值

    C++ - 系统的odeint解决方案差异很大

    执行 3d 图时 python 卡住