用 nlme 和 lsoda 拟合一阶方程 [英] fitting first order equation with nlme and lsoda
问题描述
我尝试使用 nlme 和 lsoda 拟合一阶微分模型.这是基本思想:我首先定义允许生成微分方程解的函数:
library(deSolve)ODE1 <- function(time, x, parms) {with(as.list(c(parms, x)), {导入 <- excfunc(time)dS <- 进口*k/tau - (S-yo)/taures＜-c(dS)列表(资源)})}solution_ODE1 = function(tau1,k1,yo1,excitation,time){excfunc <- approxfun(时间，激发，规则= 2)parms <- c(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)xstart = c(S = yo1)out <- lsoda(xstart, time, ODE1, parms)返回(出[，2])}然后我按照两个 ID 的等式生成数据:
time <- 0:49激发 <- 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)+norm(length(time),0,0.1)),时间 = 代表(时间，2)，激发 = rep(excitation,2),ID = rep(c("A","B"),each = length(time)))这是它的样子:
库(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),数据 = simu_data，固定 = 阻尼 + 增益 + eq ~1，随机 = 阻尼 ~ 1 ,组 = ~ ID,开始 = c(阻尼 = 5, 增益 = 1,eq = 0))我收到此错误，但我没有收到:
<块引用>eval(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)), {导入 <- excfunc(time)dS <- 导入 * k/tau - (S - yo)/taures <- c(dS) ...13. with(as.list(c(parms, x)), {导入 <- excfunc(time)dS <- 导入 * k/tau - (S - yo)/taures <- 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(阻尼、增益、均衡、激励、时间)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,激励,时间)，数据 = simu_data，固定 = 阻尼 + 增益 + eq ~ 1，随机 = 阻尼 ~ 1，组 = ~ID，开始 = c(阻尼 = 5，增益 = 1，当量 = 0))1. nlme(信号~解_ODE1(阻尼、增益、均衡、激励、时间),数据 = simu_data，固定 = 阻尼 + 增益 + eq ~ 1，随机 = 阻尼 ~1、groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))有人知道我应该如何进行吗?
<小时>编辑
我尝试按照 mikeck 的建议进行修改:
ODE1 <- function(time, x, parms) {导入 <- parms$excfunc(time)dS <- import*parms$k/parms$tau - (x["S"]-parms$yo)/parms$taures＜-c(dS)列表(资源)}生成数据没有问题.但是使用 nlme 现在给出:
checkFunc(Func2, times, y, rho) 中的错误:func() 返回的导数个数(0) 必须等于初始条件向量的长度(100)
具有以下回溯:
<代码>>追溯()10: stop(paste("func() 返回的导数个数(",length(tmp[[1]]), ") 必须等于初始条件向量的长度 (",长度(y)，)"，sep ="))9: checkFunc(Func2, times, y, rho)8: lsoda(xstart, time, ODE1, parms) 在 #57:solution_ODE1(阻尼、增益、均衡、激励、时间)6: eval(model, data.frame(data, pars))5: eval(model, data.frame(data, pars))4: eval(modelExpression[[2]],envir = nlEnv)3:评估(模型表达[[2]]，环境= nlEnv)2: nlme.formula(signal ~ solution_ODE1(damping, gain, eq,激励,时间)，数据 = simu_data，固定 = 阻尼 + 增益 + eq ~ 1，随机 = 阻尼 ~ 1，组 = ~ID，开始 = c(阻尼 = 5，增益 = 1，当量 = 0))1:nlme(信号~解_ODE1(阻尼，增益，均衡，激励，时间)，数据 = simu_data，固定 = 阻尼 + 增益 + eq ~ 1，随机 = 阻尼 ~1、groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))在您的示例中，您的 times 向量不是单调运行的.我认为这与 lsoda 混乱.时间在这里运作的方式的上下文/含义是什么?将随机效应模型与两组拟合并没有真正意义.您是否试图将同一条曲线拟合到两个独立的时间序列?
这是一个精简的示例，进行了一些调整(并非所有内容都可以在不丢失必要结构的情况下折叠为数字向量):
library(deSolve)ODE1 <- 函数(时间，x，参数){with(as.list(parms), {导入 <- excfunc(time)dS <- 导入*k/tau - (x-yo)/taures＜-c(dS)列表(资源)})}solution_ODE1 = function(tau1,k1,yo1,excitation,time){excfunc <- approxfun(时间，激发，规则= 2)参数 <- 列表(tau = tau1，k = k1，yo = yo1，excfunc = excfunc)xstart = yo1out <- lsoda(xstart, time, ODE1, parms)返回(出[，2])}时间 <- 0:49激发 <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))simu_data <- data.frame(time = rep(time,2),激发 = rep(excitation,2))svec <- c(阻尼 = 3，增益 = 1.75，eq = 0.2)这有效:
with(c(simu_data, as.list(svec)),解决方案_ODE1(阻尼，增益，均衡，激励[1:50]，时间[1:50]))但是如果我们再包含一个步骤(这样时间重置为 0)，它就会失败:
with(c(simu_data, as.list(svec)),解决方案_ODE1(阻尼，增益，均衡，激励[1:51]，时间[1:51]))<块引用>
lsoda(xstart, time, ODE1, parms) 中的错误:在采取任何集成步骤之前检测到非法输入 - 请参阅书面消息
I a trying to fit a first order differential model using nlme and lsoda.
Here is the basic idea: I first define the function allowing to generate the solution of the differential equation:
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])
}
I then generate data following the equation for two IDs:
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)))
Here is what it looks like :
library(ggplot2)
ggplot(simu_data)+
geom_point(aes(time,signal,color = "signal"),size = 2)+
geom_line(aes(time,excitation,color = "excitation"))+
facet_wrap(~ID)
I am then trying to fit using 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))
I am getting this eror, that I don't get:
Error in eval(substitute(expr), data, enclos = parent.frame()) : object 'k' not found
The traceback shows that the error comes from the ODE1 model, which works when generating values.
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))
Does anyone have an idea How I should proceed ?
Edit
I tried to modify following the advise of 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)}
Generating the data works without problems. But use of nlme gives now:
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)
with the following traceback:
> 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))
In your example, your times vector doesn't run monotonically. I think that messes with lsoda. What is the context/meaning of the way that time works here? It doesn't really make sense to fit a random-effects model with two groups. Are you trying to fit the same curve to two independent time series?
Here's a stripped-down example, with some adjustments (not everything can be collapsed to a numeric vector without losing necessary structure):
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)
This works:
with(c(simu_data, as.list(svec)),
solution_ODE1(damping,gain,eq,excitation[1:50],time[1:50]))
But if we include one more step (so that time resets to 0), it fails:
with(c(simu_data, as.list(svec)),
solution_ODE1(damping,gain,eq,excitation[1:51],time[1:51]))
Error in lsoda(xstart, time, ODE1, parms) : illegal input detected before taking any integration steps - see written message
这篇关于用 nlme 和 lsoda 拟合一阶方程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
