Using R ode models as functions in nls

user_12370404 0
Viewed 74

I try to estimate parameters for an ode model in R using nls (and later nlme). My testing code gives me error messages.

library(deSolve)
seed=2423
dat2<-data.frame(days=(runif(20)+1)*10, X1=runif(20), X2=runif(20)) 
dat2$y<-0.4*exp(dat2$X1)+0.6*exp(dat2$X2)+rnorm(20, sd=0.3)
# example intentionally simple. I would usually solve it analysically
#***************************************************
#*Model definition
#***************************************************
decomp<-function(t, state, parameters){
  with(as.list(c(state, parameters)), {
    dX1<-a1*X1
    dX2<-a2*X2#+a1*X1
    list(c(dX1, dX2))
  } )
}
# Testing the code to demonstrate that it works
parameters<-c(a1=0.05,a2=0.05) 
state<-c(X1=0.5, X2=0.5)
times<-seq(0,100, by=1)
out<-ode(y=state, times=times, func=decomp, parms=parameters)
out[100,]
#*****************************
#* Wrapper function to be passed to nls or nlme
#************************************** 
calcdecom<-function(a1,a2,t,x1, x2)
{ 
  state<-c(X1=x1, X2=x2)
  times<-c(0,t)
  parameters<-c(a1=a1,a2=a2)
  out<-ode(y=state, times=times, func=decomp, parms=parameters)
  return(as.numeric(out[2,2]+out[2,3]))
}
# ******************** test
calcdecom(0.1,0.1,5,0.3,0.3)


test<-nls(y~calcdecom(a1, a2, days, X1, X2 ), 
          start=list(a1=0.02, a2=0.4), data=dat2)

My error messages for the nls function is:

Error in lsoda(y, times, func, parms, ...) : illegal input detected before taking any integration steps - see written message

nls ode r
1 Answers

Here a possible approach. One remaining question is, whether the data should really come from independent cases (simulations or measurements) and not from a time series. If the first is intended, then ode must be called for each case separately. This can be done in a for-loop or with an applyfunction. Note also the correction of the set.seed-call, the plot function and the reduced standard deviation in the data generating process. It seems to me that the ode model and the data generating process do not match yet, so that the fitted parameters differ. It would be nice if the OP could post a corrected version by editing the question.

library("deSolve")

## use set.seed to make example reproducible
set.seed(2423)

## simulated data
dat2 <- data.frame(
  days = (runif(20) + 1) + 10,
  X1 = runif(20),
  X2 = runif(20)
)

## reduced error for testing
dat2$y <- 0.4 * exp(dat2$X1) + 0.6 * exp(dat2$X2) + rnorm(20, sd = 0.1)

plot(dat2)

decomp <- function(t, state, parameters){
  with(as.list(c(state, parameters)), {
    dX1 <- a1 * X1
    dX2 <- a2 * X2 + a1 * X1
    list(c(dX1, dX2))
  } )
}

## test the ode moel to demonstrate that it works
parameters <- c(a1 = 0.05, a2 = 0.05)
state <- c(X1 = 0.5, X2 = 0.5)
times <- seq(0, 100, by = 1)
out <- ode(y = state, times = times, func = decomp, parms = parameters)
plot(out)
out[100,]

## Wrapper function to be passed to nls or nlme
calcdecom <- function(a1, a2, t, x1, x2) {
  ret <- numeric(length(t))
  parameters <- c(a1 = a1, a2 = a2)
  for (i in 1:length(t)) {
    times <- c(0, t[i])
    state <- c(X1 = x1[i], X2 = x2[i])
    out <- ode(y = state, times = times, func = decomp, parms = parameters)
    ret[i] <- out[2, 2] + out[2, 3]
  }
  return(ret)
}

## test wrapper function
calcdecom(0.1, 0.1, 5, 0.3, 0.3)

## with time and state as vectors
calcdecom(0.1, 0.1, 1:10, dat2$X1, dat2$X2)

test <- nls(y ~ calcdecom(a1, a2, days, X1, X2 ),
          start = list(a1 = 0.02, a2 = 0.4), data = dat2)

summary(test)

Instead of nls one may also consider modFit from package FME, that has some more flexibility for this kind of models. Details are found in the package vignettes and the following JSS paper: https://doi.org/10.18637/jss.v033.i03

user_36775760
Related