问题描述
我面临一个我无法解决的问题。我想使用 nlme 或 nlmODE 执行具有随机效应的非线性回归，使用具有固定系数的二阶微分方程（阻尼振荡器）的解作为模型。
我设法将 nlme 用于简单模型，但似乎使用 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,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,1,0)
plot(time,oscil)
现在我生成了一个由 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.02),by = ID]
如果我们看一看，我们有一个合适的数据集：
library(ggplot2)
ggplot(data_simu,aes(time,signal,color = ID))+
geom_line()+
facet_wrap(~ID)
问题
使用 nlme
使用具有类似语法的 nlme 处理更简单的示例（非线性函数不使用 deSolve），我尝试过：
fit <- nlme(model = signal ~ solution_analy_ODE2(esp2omega,omega2,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))
我得到：
错误 in checkFunc(Func2,times,rho) : func() 返回的导数 (2) 必须等于初始条件向量的长度 (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,rho)
10. lsoda(xstart,parms)
9. solution_analy_ODE2(omega2,yeq)
.
.
我看起来像 nlme 试图将起始条件向量传递给 solution_analy_ODE2，并导致 checkFunc 中的 lasoda 出错。
我尝试使用 nlsList：
test <- nlsList(model = signal ~ solution_analy_ODE2(omega2,yeq) | ID,start = list(esp2omega = 0.08,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,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,v0 = 0))
但我得到了完全相同的错误信息。
然后我尝试使用 nlmODE，遵循 Bne Bolker 对 a similar question I asked some years ago 的评论
使用 nlmMODE
library(nlmeODE)
datas_grouped <- groupedData( signal ~ time | ID,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,yeq,ID),data = datas_grouped,fixed = esp2omega + omega2 + yeq + y0 + v0 ~ 1,random = y0 + v0 ~1,y0 = 0,v0 = 0)) #
我收到错误：
resnlmeode(esp2omega,ID) 中的错误：找不到对象“yhat”
问题
- 你能重现这个问题吗？
- 有没有人想用
nlme或nlmODE解决这个问题？ - 如果没有，是否有使用其他软件包的解决方案？我看到了
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,out)
cat(apply(tmp,function(x) paste(paste(x,collapse = " "),"\n")),"\n")
return(out)
}
time <- 0:10*10
ramdom_y0 <- sample(seq(0,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.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,y0),fixed = b + A + y0 ~ 1,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 绑定 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,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,ODE2_nls,parms,atol=1e-4,rtol=1e-4,method="ode45")
return(out[,2])
}
然后设置（或计算）方程的数量，例如N <- 1 分别N <-1000 在通话之前。
模型通过这种方式运行，然后运行在数值问题中，但那是另一回事了......
然后您可以尝试使用另一个 ode 求解器（例如 vode），将 atol 和 rtol 设置为较低的值，调整 nmle 的优化参数，使用框约束...等等，就像在非线性优化中一样。
我发现了一个解决 nlme 行为的解决方案：如我的编辑所示，问题来自于 nlme 将 NindividualxNpoints 的向量传递给非线性函数的事实，假设函数为每个时间点一个值。但是 lsoda 不这样做，因为它沿时间积分一个方程（即它需要所有时间直到给定的时间点产生一个值）。
我的解决方案包括分解 nlme 传递给我的函数的参数，进行计算，并重新创建一个向量：
detect_id <- function(vec){
tmp <- c(0,diff(vec))
out <- tmp
out <- NA
out[tmp < 0] <- 1:sum(tmp < 0)
out <- na.locf(out,na.rm = F)
rleid(out)
}
detect_id 将时间向量分解为单个时间向量标识符：
detect_id(rep(1:10,3))
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3
然后，函数对每个个体进行数值积分循环，并将结果向量绑定在一起：
solution_analy_ODE2_modif = function(omega2,yeq){
tmp <- detect_id(time)
out <- lapply(unique(tmp),function(i){
idxs <- which(tmp == i)
parms <- c(esp2omega = esp2omega[idxs][1],omega2 = omega2[idxs][1],yeq = yeq[idxs][1])
xstart = c(S1 = y0[idxs][1],dS1 = v0[idxs][1])
out_tmp <- lsoda(xstart,time[idxs],parms)
out_tmp[,2]
}) %>% unlist()
return(out)
}
我做了一个测试，在那里我传递了一个类似于 nlme 传递给函数的向量：
omega2vec <- rep(0.1,30)
eps2omegavec <- rep(0.1,30)
timevec <- rep(1:10,3)
y0vec <- rep(1,30)
v0vec <- rep(0,30)
yeqvec = rep(0,30)
solution_analy_ODE2_modif(omega2 = omega2vec,esp2omega = eps2omegavec,time = timevec,y0 = y0vec,v0 = v0vec,yeq = yeqvec)
[1] 1.0000000 0.9520263 0.8187691 0.6209244 0.3833110 0.1321355 -0.1076071 -0.3143798
[9] -0.4718058 -0.5697255 1.0000000 0.9520263 0.8187691 0.6209244 0.3833110 0.1321355
[17] -0.1076071 -0.3143798 -0.4718058 -0.5697255 1.0000000 0.9520263 0.8187691 0.6209244
[25] 0.3833110 0.1321355 -0.1076071 -0.3143798 -0.4718058 -0.5697255
它有效。它不适用于@tpetzoldt 方法，因为时间向量从 10 传递到 0，这会导致积分问题。在这里，我真的需要破解 nlnme 的工作方式。
现在：
fit <- nlme(model = signal ~ solution_analy_ODE2_modif (esp2omega,omega2,yeq),data = data_simu,fixed = esp2omega + omega2 + y0 + v0 + yeq ~ 1,random = y0 ~ 1,groups = ~ ID,start = c(esp2omega = 0.5,omega2 = 0.5,yeq = 0,y0 = 1,v0 = 1))
就像一个魅力
summary(fit)
Nonlinear mixed-effects model fit by maximum likelihood
Model: signal ~ solution_analy_ODE2_modif(omega2,yeq)
Data: data_simu
AIC BIC logLik
-597.4215 -567.7366 307.7107
Random effects:
Formula: list(y0 ~ 1,v0 ~ 1)
Level: ID
Structure: General positive-definite,Log-Cholesky parametrization
StdDev Corr
y0 0.61713329 y0
v0 0.67815548 -0.269
Residual 0.03859165
Fixed effects: esp2omega + omega2 + y0 + v0 + yeq ~ 1
Value Std.Error DF t-value p-value
esp2omega 0.4113068 0.00866821 186 47.45002 0.0000
omega2 1.0916444 0.00923958 186 118.14876 0.0000
y0 0.3848382 0.19788896 186 1.94472 0.0533
v0 0.1892775 0.21762610 186 0.86974 0.3856
yeq 0.0000146 0.00283328 186 0.00515 0.9959
Correlation:
esp2mg omega2 y0 v0
omega2 0.224
y0 0.011 -0.008
v0 0.005 0.030 -0.269
yeq -0.091 -0.046 0.009 -0.009
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.2692477 -0.6122453 0.1149902 0.6460419 3.2890201
Number of Observations: 200
Number of Groups: 10