1

So I have a system of ode's and some data I am using the R packages deSolve and FME to fit the parameters of the ode system to data. I am getting a singular matrix result when I fit the full parameter set to the data. So I went back and looked at the collinearity of the parameters using a collinearity index cut-off of 20 as suggested in all the FME package documentation I then picked a few models with subsets of parameters to fit. Then when I run modFit I get this error:

Error in approx(xMod, yMod, xout = xDat) : need at least two non-NA values to interpolate

Can anyone enlighten me as to a fix for this. Everything else is working fine. So this is not a coding problem.

Here is a minimal working example (removing r=2 in modFit creates the error which I can fix in the minimal working example but not in my actual problem so I doubt a minimal working example helps here):

`## =======================================================================
## Now suppose we do not know K and r and they are to be fitted...
## The "observations" are the analytical solution
## =======================================================================

# You need these packages
library('deSolve')
library('FME')

## logistic growth model
TT <- seq(1, 100, 2.5)
N0 <- 0.1
r <- 0.5
K <- 100

## analytical solution

Ana <- cbind(time = TT, N = K/(1 + (K/N0 - 1) * exp(-r*TT)))

time <- 0:100
parms <- c(r = r, K = K)
x <- c(N = N0)

logist <- function(t, x, parms) {
  with(as.list(parms), {
    dx <- r * x[1] * (1 - x[1]/K)
    list(dx)
  })
}

## Run the model with initial guess: K = 10, r = 2

parms["K"] <- 10
parms["r"] <- 2
init <- ode(x, time, logist, parms)

## FITTING algorithm uses modFit
## First define the objective function (model cost) to be minimised
## more general: using modFit

Cost <- function(P) {
  parms["K"] <- P[1]
  parms["r"] <- P[2]
  out <- ode(x, time, logist, parms)
  return(modCost(out, Ana))
} 
(Fit<-modFit(p = c(K = 10,r=2), f = Cost))
summary(Fit)` 

1 Answers1

1

I think the problem is in your Cost function. If you don't provide both K and r, then the cost function will override the start value of r to NA. You can test this:

Cost <- function(P) {
  parms["K"] <- P[1]
  parms["r"] <- P[2]
  print(parms)
  #out <- ode(x, time, logist, parms)
  #return(modCost(out, Ana))
} 
Cost(c(K=10, r = 2))
Cost(c(K=10))

This function works:

Cost <- function(P) {
  parms[names(P)] <- P
  out <- ode(x, time, logist, parms)
  return(modCost(out, Ana))
} 

The vignette FMEDyna is very helpful: https://cran.r-project.org/web/packages/FME/vignettes/FMEdyna.pdf See page 14 on how to specify the Objective (Cost) function.

jtay
  • 53
  • 7
  • Thanks hadn't looked here in a while I had resolved this and forgot to report back yes, it's more or less as you say and the choice of solver makes a big difference needed to try a few more but I did get it working. Thanks. – Rodney Beard Jan 21 '16 at 01:03