0

I am working on a model with 8 ODEs, and want to fit some of the parameters (not all) based on observed data I have on just 3 of the 8 state variables. This is my code:

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")

After running this code I get the following warning message:

Warning message:
In nls.lm(par = Pars, fn = Fun, control = Contr, ...) :
  lmdif: info = 0. Improper input parameters.

And summary(fit)gives me this error:

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

I just want to fit these parameters: us, ul, ms, ml, g, mz and ub. I am quite confident with the rest of the parameters. Any help or hint on how to do this would be much appreciated.

Tomas Marina
  • 73
  • 10

1 Answers1

4

A bit late maybe but one never knows. Regarding your code, this is not recommended to try to fit so many parameters at the same time. You will see in my code below that I use the sensFun() function in order to select the parameters that have the biggest impact for that simulation. This enables me to select 5 parameters, instead of your whole list. I would also add a part on the Collin() function, that can help you to decide whether a given parameter is identifiable or not, how many parameters you can simultaneously estimate... With the code below, I manage to get a correct fit.

Model Fitted


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)
Dav
  • 122
  • 9