在 R 中的 ODE 模型中拟合特定参数集
Fit specific set of parameters in ODEs model in R
我正在研究一个具有 8 个 ODE 的模型，并希望根据我在 8 个状态变量中的 3 个上的观察数据来拟合一些参数（不是全部）。
这是我的代码：
library("FME")
library("deSolve")
library("lattice")
# Model construction and definition of derivatives
model.sal <- function(time, y, param)
{
N <- y[1]
NH4 <- y[2]
Ps <- y[3]
Pl <- y[4]
Z <- y[5]
B <- y[6]
DON <- y[7]
D <- y[8]
with(as.list(param), {
dNdt <- nit*NH4*B - us*(N/(N+kns))*Ps - ul*(N/(N+knl))*Pl
dNH4dt <- fraz*Z + exb*B - us*(NH4/(NH4+kas))*Ps - ul*(NH4/(NH4+kal))*Pl - ub*(NH4/(NH4+kb))*B
dPsdt <- Ps*(us*((N/N+kns)*(NH4/NH4+kas)*(exp(-((S-Sop)^2)/ts^2))*(tanh(alfa*Im/Pm))) - exs - ms*(Ps/kms+Ps) - g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2)*Z)
dPldt <- Pl*(ul*((N/N+knl)*(NH4/NH4+kal)*(exp(-((S-Sop)^2)/tl^2))*(tanh(alfa*Im/Pm))) - exl - ml*(Pl/kml+Pl) - g*(pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2)*Z)
dZdt <- Z*(ge*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2)) - frdz
- fraz - mz*(Z/kmz+Z))
dBdt <- B*(ub*(NH4/(NH4+kb))*(DON/(DON+kb)) - exb - g*(pfb*B^2/B*pfb*kg+pfb*B^2)*Z)
dDONdt <- frdz*Z + exs*Ps + exl*Pl + bd*D - ub*(DON/(DON+kb))
dDdt <- (1-ge)*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2))
+ ms*(Ps/kms+Ps) + ml*(Pl/kml+Pl) + mz*(Z/kmz+Z) - bd*D
return(list(c(dNdt, dNH4dt, dPsdt, dPldt, dZdt, dBdt, dDONdt, dDdt)))
})
}
# Observed data on 3 of the 8 state variables
dat <- data.frame(
time = seq(0, 8, 1),
N = c(11.54, 16.6, 7.86, 6.73, 5.6, 5.2, 4.81, 4.18, 3.55),
Pl = c(3.85, 6.25, 3.41, 6.16, 8.92, 12.79, 16.26, 19.21, 22.36),
Ps = c(0.09, 0.33, 0.18, 0.06, 0.12, 0.4, 0.84, 0.7, 0.48))
# Parameters
param.gotm <- c(nit=0.1, us=0.7, kns=0.5, kas=0.5, exs=0.05, ms=0.05,
kms=0.2, ul=0.7, knl=0.5, kal=0.5, exl=0.02, ml=0.05,
kml=0.2, ge=0.625, g=0.35, kg=1, pfs=0.55, pfl=0.3, pfb=0.1,
pfd=0.05, frdz=0.1, fraz=0.7, mz=0.2, kmz=0.2, ub=0.24,
kb=0.05, exb=0.03, bd=0.33, alfa=0.1, Im=100, Pm=0.04,
Sop=34, S=34, ts=2, tl=1)
# Time options, initial values and ODE solution
times <- seq(0, 10, length=200)
y0 <- c(N=7, NH4=0.01, Ps=0.17, Pl=0.77, Z=0.012, B=0.001, DON= 0.001, D=0.01)
out1 <- ode(y0, times, model.sal, param.gotm)
plot(out1, obs = dat)
# Definition of the cost function
cost <- function(p)
{
out <- ode(y0, times, model.sal, p)
modCost(out, dat, weight = "none")
}
fit <- modFit(f = cost, p = param.gotm, method = "Marq")
在 运行 这段代码之后，我收到以下警告消息：
Warning message:
In nls.lm(par = Pars, fn = Fun, control = Contr, ...) :
lmdif: info = 0. Improper input parameters.
并且summary(fit)给我这个错误：
Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix
In addition: Warning message:
In summary.modFit(fit) : Cannot estimate covariance; system is singular
我只想拟合这些参数：us、ul、ms、ml、g、mz 和 ub。我对其余参数非常有信心。任何有关如何执行此操作的帮助或提示将不胜感激。
也许晚了一点，但谁也说不准。
关于您的代码，不建议尝试同时适应这么多参数。您将在下面的代码中看到，我使用 sensFun() 函数来 select 对模拟影响最大的参数。这使我能够 select 5 个参数，而不是你的整个列表。我还会在 Collin() 函数上添加一个部分，它可以帮助您确定给定参数是否可识别，您可以同时估计多少个参数......使用下面的代码，我设法得到一个正确合身。
library("FME")
library("deSolve")
library("lattice")
# # # # # # # # # # # # # # # # #
#
# 1) Preliminary functions
#
# # # # # # # # # # # # # # # # #
# Parameters
pars <- c(
nit=0.1, us=0.7, kns=0.5, kas=0.5, exs=0.05, ms=0.05,
kms=0.2, ul=0.7, knl=0.5, kal=0.5, exl=0.02, ml=0.05,
kml=0.2, ge=0.625, g=0.35, kg=1, pfs=0.55, pfl=0.3, pfb=0.1,
pfd=0.05, frdz=0.1, fraz=0.7, mz=0.2, kmz=0.2, ub=0.24,
kb=0.05, exb=0.03, bd=0.33, alfa=0.1, Im=100, Pm=0.04,
Sop=34, S=34, ts=2, tl=1
)
# Model construction and definition of derivatives
model.sal <- function(t, state, pars){
with(as.list(c(state, pars)), {
dNdt <- nit*NH4*B - us*(N/(N+kns))*Ps - ul*(N/(N+knl))*Pl
dNH4dt <- fraz*Z + exb*B - us*(NH4/(NH4+kas))*Ps - ul*(NH4/(NH4+kal))*Pl - ub*(NH4/(NH4+kb))*B
dPsdt <- Ps*(us*((N/N+kns)*(NH4/NH4+kas)*(exp(-((S-Sop)^2)/ts^2))*(tanh(alfa*Im/Pm))) - exs - ms*(Ps/kms+Ps) - g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2)*Z)
dPldt <- Pl*(ul*((N/N+knl)*(NH4/NH4+kal)*(exp(-((S-Sop)^2)/tl^2))*(tanh(alfa*Im/Pm))) - exl - ml*(Pl/kml+Pl) - g*(pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2)*Z)
dZdt <- Z*(ge*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2)) - frdz
- fraz - mz*(Z/kmz+Z))
dBdt <- B*(ub*(NH4/(NH4+kb))*(DON/(DON+kb)) - exb - g*(pfb*B^2/B*pfb*kg+pfb*B^2)*Z)
dDONdt <- frdz*Z + exs*Ps + exl*Pl + bd*D - ub*(DON/(DON+kb))
dDdt <- (1-ge)*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2))
+ ms*(Ps/kms+Ps) + ml*(Pl/kml+Pl) + mz*(Z/kmz+Z) - bd*D
return(list(c(dNdt, dNH4dt, dPsdt, dPldt, dZdt, dBdt, dDONdt, dDdt)))
})
}
# wrapper
solve_model <- function(pars, times = seq(0, 10, length=200)) {
# initial values
state <- c(N=7, NH4=0.01, Ps=0.17, Pl=0.77, Z=0.012, B=0.001, DON= 0.001, D=0.01)
out <- ode(y = state, times = times, func = model.sal, parms = pars)
return(out)
}
# Definition of the cost function
Objective <- function(x, parset = names(x)) {
pars[parset] <- x
tout <- seq(0, 10, length=200)
out <- solve_model(pars, tout)
modCost(out, dat, weight = "none")
}
# # # # # # # # # # # # # # # # #
#
# 2) Preliminary data
#
# # # # # # # # # # # # # # # # #
# Observed data on 3 of the 8 state variables
dat <- data.frame(
time = seq(0, 8, 1),
N = c(11.54, 16.6, 7.86, 6.73, 5.6, 5.2, 4.81, 4.18, 3.55),
Pl = c(3.85, 6.25, 3.41, 6.16, 8.92, 12.79, 16.26, 19.21, 22.36),
Ps = c(0.09, 0.33, 0.18, 0.06, 0.12, 0.4, 0.84, 0.7, 0.48))
# # # # # # # # # # # # # # # # # # # # #
#
# 3) Select the good parameters to fit
#
# # # # # # # # # # # # # # # # # # # # #
# Determine what are the best parameters to fit
Sfun <- sensFun(Objective, pars)
plot(summary(Sfun))
# from the mean plot, I see that ul, us, pfl, ge and knl have the most influence for the simulation
# I will do the optimization on them so.
# # # # # # # # # # # # #
#
# 4) Optimization
#
# # # # # # # # # # # # #
# set up the subset of parameters
parToFit <- c(ul = 0.7, us = 0.7, pfl = 0.3, ge = 0.625, knl = 0.5)
# run the beast
Fit <- modFit(
f = Objective,
p = parToFit,
lower = 0,
upper = Inf,
method = "Marq",
jac = NULL,
control = list(
#maxiter = 100,
ftol = 1e-06,
ptol = 1e-06,
gtol = 1e-06,
nprint = 1
),
hessian = TRUE
)
# # # # # # # # # # # # #
#
# 5) Rerun simulations and plot
#
# # # # # # # # # # # # #
# recover the optimized parameters and plot the results
# You could also plot the non optimized curves to compare
pars[names(parToFit)] <- Fit$par
optim <- solve_model(pars, times = seq(0, 10, length=200))
par(mfrow = c(2, 2))
plot(optim[, "time"], optim[, "N"], xlab = "Time (min)", ylab = "N", lwd = 2, type = "l", col = "red")
points(dat[, "time"], dat[, "N"], cex = 2, pch = 18)
plot(optim[, "time"], optim[, "Pl"], xlab = "Time (min)", ylab = "Pl", lwd = 2, type = "l", col = "red")
points(dat[, "time"], dat[, "Pl"], cex = 2, pch = 18)
plot(optim[, "time"], optim[, "Ps"], xlab = "Time (min)", ylab = "Ps", lwd = 2, type = "l", col = "red")
points(dat[, "time"], dat[, "Ps"], cex = 2, pch = 18)
我正在研究一个具有 8 个 ODE 的模型，并希望根据我在 8 个状态变量中的 3 个上的观察数据来拟合一些参数（不是全部）。 这是我的代码：
library("FME")
library("deSolve")
library("lattice")
# Model construction and definition of derivatives
model.sal <- function(time, y, param)
{
N <- y[1]
NH4 <- y[2]
Ps <- y[3]
Pl <- y[4]
Z <- y[5]
B <- y[6]
DON <- y[7]
D <- y[8]
with(as.list(param), {
dNdt <- nit*NH4*B - us*(N/(N+kns))*Ps - ul*(N/(N+knl))*Pl
dNH4dt <- fraz*Z + exb*B - us*(NH4/(NH4+kas))*Ps - ul*(NH4/(NH4+kal))*Pl - ub*(NH4/(NH4+kb))*B
dPsdt <- Ps*(us*((N/N+kns)*(NH4/NH4+kas)*(exp(-((S-Sop)^2)/ts^2))*(tanh(alfa*Im/Pm))) - exs - ms*(Ps/kms+Ps) - g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2)*Z)
dPldt <- Pl*(ul*((N/N+knl)*(NH4/NH4+kal)*(exp(-((S-Sop)^2)/tl^2))*(tanh(alfa*Im/Pm))) - exl - ml*(Pl/kml+Pl) - g*(pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2)*Z)
dZdt <- Z*(ge*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2)) - frdz
- fraz - mz*(Z/kmz+Z))
dBdt <- B*(ub*(NH4/(NH4+kb))*(DON/(DON+kb)) - exb - g*(pfb*B^2/B*pfb*kg+pfb*B^2)*Z)
dDONdt <- frdz*Z + exs*Ps + exl*Pl + bd*D - ub*(DON/(DON+kb))
dDdt <- (1-ge)*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2))
+ ms*(Ps/kms+Ps) + ml*(Pl/kml+Pl) + mz*(Z/kmz+Z) - bd*D
return(list(c(dNdt, dNH4dt, dPsdt, dPldt, dZdt, dBdt, dDONdt, dDdt)))
})
}
# Observed data on 3 of the 8 state variables
dat <- data.frame(
time = seq(0, 8, 1),
N = c(11.54, 16.6, 7.86, 6.73, 5.6, 5.2, 4.81, 4.18, 3.55),
Pl = c(3.85, 6.25, 3.41, 6.16, 8.92, 12.79, 16.26, 19.21, 22.36),
Ps = c(0.09, 0.33, 0.18, 0.06, 0.12, 0.4, 0.84, 0.7, 0.48))
# Parameters
param.gotm <- c(nit=0.1, us=0.7, kns=0.5, kas=0.5, exs=0.05, ms=0.05,
kms=0.2, ul=0.7, knl=0.5, kal=0.5, exl=0.02, ml=0.05,
kml=0.2, ge=0.625, g=0.35, kg=1, pfs=0.55, pfl=0.3, pfb=0.1,
pfd=0.05, frdz=0.1, fraz=0.7, mz=0.2, kmz=0.2, ub=0.24,
kb=0.05, exb=0.03, bd=0.33, alfa=0.1, Im=100, Pm=0.04,
Sop=34, S=34, ts=2, tl=1)
# Time options, initial values and ODE solution
times <- seq(0, 10, length=200)
y0 <- c(N=7, NH4=0.01, Ps=0.17, Pl=0.77, Z=0.012, B=0.001, DON= 0.001, D=0.01)
out1 <- ode(y0, times, model.sal, param.gotm)
plot(out1, obs = dat)
# Definition of the cost function
cost <- function(p)
{
out <- ode(y0, times, model.sal, p)
modCost(out, dat, weight = "none")
}
fit <- modFit(f = cost, p = param.gotm, method = "Marq")
在 运行 这段代码之后，我收到以下警告消息：
Warning message:
In nls.lm(par = Pars, fn = Fun, control = Contr, ...) :
lmdif: info = 0. Improper input parameters.
并且summary(fit)给我这个错误：
Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix
In addition: Warning message:
In summary.modFit(fit) : Cannot estimate covariance; system is singular
我只想拟合这些参数：us、ul、ms、ml、g、mz 和 ub。我对其余参数非常有信心。任何有关如何执行此操作的帮助或提示将不胜感激。
也许晚了一点，但谁也说不准。
关于您的代码，不建议尝试同时适应这么多参数。您将在下面的代码中看到，我使用 sensFun() 函数来 select 对模拟影响最大的参数。这使我能够 select 5 个参数，而不是你的整个列表。我还会在 Collin() 函数上添加一个部分，它可以帮助您确定给定参数是否可识别，您可以同时估计多少个参数......使用下面的代码，我设法得到一个正确合身。
library("FME")
library("deSolve")
library("lattice")
# # # # # # # # # # # # # # # # #
#
# 1) Preliminary functions
#
# # # # # # # # # # # # # # # # #
# Parameters
pars <- c(
nit=0.1, us=0.7, kns=0.5, kas=0.5, exs=0.05, ms=0.05,
kms=0.2, ul=0.7, knl=0.5, kal=0.5, exl=0.02, ml=0.05,
kml=0.2, ge=0.625, g=0.35, kg=1, pfs=0.55, pfl=0.3, pfb=0.1,
pfd=0.05, frdz=0.1, fraz=0.7, mz=0.2, kmz=0.2, ub=0.24,
kb=0.05, exb=0.03, bd=0.33, alfa=0.1, Im=100, Pm=0.04,
Sop=34, S=34, ts=2, tl=1
)
# Model construction and definition of derivatives
model.sal <- function(t, state, pars){
with(as.list(c(state, pars)), {
dNdt <- nit*NH4*B - us*(N/(N+kns))*Ps - ul*(N/(N+knl))*Pl
dNH4dt <- fraz*Z + exb*B - us*(NH4/(NH4+kas))*Ps - ul*(NH4/(NH4+kal))*Pl - ub*(NH4/(NH4+kb))*B
dPsdt <- Ps*(us*((N/N+kns)*(NH4/NH4+kas)*(exp(-((S-Sop)^2)/ts^2))*(tanh(alfa*Im/Pm))) - exs - ms*(Ps/kms+Ps) - g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2)*Z)
dPldt <- Pl*(ul*((N/N+knl)*(NH4/NH4+kal)*(exp(-((S-Sop)^2)/tl^2))*(tanh(alfa*Im/Pm))) - exl - ml*(Pl/kml+Pl) - g*(pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2)*Z)
dZdt <- Z*(ge*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2)) - frdz
- fraz - mz*(Z/kmz+Z))
dBdt <- B*(ub*(NH4/(NH4+kb))*(DON/(DON+kb)) - exb - g*(pfb*B^2/B*pfb*kg+pfb*B^2)*Z)
dDONdt <- frdz*Z + exs*Ps + exl*Pl + bd*D - ub*(DON/(DON+kb))
dDdt <- (1-ge)*(g*(pfs*Ps^2/Ps*pfs*kg+pfs*Ps^2) + (pfl*Pl^2/Pl*pfl*kg+pfl*Pl^2) + (pfb*B^2/B*pfb*kg+pfb*B^2))
+ ms*(Ps/kms+Ps) + ml*(Pl/kml+Pl) + mz*(Z/kmz+Z) - bd*D
return(list(c(dNdt, dNH4dt, dPsdt, dPldt, dZdt, dBdt, dDONdt, dDdt)))
})
}
# wrapper
solve_model <- function(pars, times = seq(0, 10, length=200)) {
# initial values
state <- c(N=7, NH4=0.01, Ps=0.17, Pl=0.77, Z=0.012, B=0.001, DON= 0.001, D=0.01)
out <- ode(y = state, times = times, func = model.sal, parms = pars)
return(out)
}
# Definition of the cost function
Objective <- function(x, parset = names(x)) {
pars[parset] <- x
tout <- seq(0, 10, length=200)
out <- solve_model(pars, tout)
modCost(out, dat, weight = "none")
}
# # # # # # # # # # # # # # # # #
#
# 2) Preliminary data
#
# # # # # # # # # # # # # # # # #
# Observed data on 3 of the 8 state variables
dat <- data.frame(
time = seq(0, 8, 1),
N = c(11.54, 16.6, 7.86, 6.73, 5.6, 5.2, 4.81, 4.18, 3.55),
Pl = c(3.85, 6.25, 3.41, 6.16, 8.92, 12.79, 16.26, 19.21, 22.36),
Ps = c(0.09, 0.33, 0.18, 0.06, 0.12, 0.4, 0.84, 0.7, 0.48))
# # # # # # # # # # # # # # # # # # # # #
#
# 3) Select the good parameters to fit
#
# # # # # # # # # # # # # # # # # # # # #
# Determine what are the best parameters to fit
Sfun <- sensFun(Objective, pars)
plot(summary(Sfun))
# from the mean plot, I see that ul, us, pfl, ge and knl have the most influence for the simulation
# I will do the optimization on them so.
# # # # # # # # # # # # #
#
# 4) Optimization
#
# # # # # # # # # # # # #
# set up the subset of parameters
parToFit <- c(ul = 0.7, us = 0.7, pfl = 0.3, ge = 0.625, knl = 0.5)
# run the beast
Fit <- modFit(
f = Objective,
p = parToFit,
lower = 0,
upper = Inf,
method = "Marq",
jac = NULL,
control = list(
#maxiter = 100,
ftol = 1e-06,
ptol = 1e-06,
gtol = 1e-06,
nprint = 1
),
hessian = TRUE
)
# # # # # # # # # # # # #
#
# 5) Rerun simulations and plot
#
# # # # # # # # # # # # #
# recover the optimized parameters and plot the results
# You could also plot the non optimized curves to compare
pars[names(parToFit)] <- Fit$par
optim <- solve_model(pars, times = seq(0, 10, length=200))
par(mfrow = c(2, 2))
plot(optim[, "time"], optim[, "N"], xlab = "Time (min)", ylab = "N", lwd = 2, type = "l", col = "red")
points(dat[, "time"], dat[, "N"], cex = 2, pch = 18)
plot(optim[, "time"], optim[, "Pl"], xlab = "Time (min)", ylab = "Pl", lwd = 2, type = "l", col = "red")
points(dat[, "time"], dat[, "Pl"], cex = 2, pch = 18)
plot(optim[, "time"], optim[, "Ps"], xlab = "Time (min)", ylab = "Ps", lwd = 2, type = "l", col = "red")
points(dat[, "time"], dat[, "Ps"], cex = 2, pch = 18)