Potential indexing issue in ODE?

Question

Following @tpetzoldt suggestion i'm opening this as a question following the previous discussion (Parameter values as a function of another vector. deSolve).

What i'm trying to achieve is to be able to integrate the model at each timestep over a vector of DailyTemperature and then the corresponding parameter values for each day are a function of values from a dataframe of other temperature outputs.

library(deSolve)
set.seed(1)

deriv <- function(t, state, pars) {
  
  pars <- parameters[match(DailyTemperature[floor(t + 1)],parameters$TraitTemperature),2:5]
  #print(pars)
  
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y
    list(c(x = d_x, y = d_y), alpha=alpha, beta=beta, gamma=gamma, delta=delta)
  })
}


state <- c(x = 1000, y = 10)
times = seq(0, 50, by = 1)

# Parameter datasets 
parameters <- data.frame(
  TraitTemperature = seq(0.1,40,0.1),
  alpha = rtruncnorm(40,a=0,b=1,mean = 1,sd=2),
  beta =  rtruncnorm(40,a=0,b=1,mean = 1,sd=2),
  delta =  rtruncnorm(40,a=0,b=1,mean = 1,sd=2),
  gamma = seq(0.025,1,0.025)
)

# random daily temperature dataset 
DailyTemperature <- round(runif(51, 0, 40),1) # one more because start zero
DailyTemperature

out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)
out

Im actually beginning to think this is an issue with parameter values now rather than code. Regardless, i'd be interested to know if my indexing is right?

Thanks for opening a new question with an updated code example, even if the title "indexing issue" is somewhat misleading. It is not an "issue". The right topic is, how to organize indexing (or table lookup). In addition, there went indeed something wrong with the definition of parameters: what is `rtruncnorm`? I guess you mean `trunc(rnorm())`. And, what is `a=0, b=1`? — tpetzoldt, Dec 21 '21 at 20:42

score 0 · Answer 1 · answered Dec 21 '21 at 21:20

The code above has several issues:

The line with rtruncnormseems completely mangled. Here I assume that alpha, beta, gamma are linearly dependent on temperature plus a random component (that I don't really understand in this case), but this does not matter here if we concentrate on the programming approach.
States and parameter values are relatively big. As the terms are essentially exponential (best seen with alpha * x that is exponential growth), the population is likely to explode or to die out. This can be prevented by balancing parameters and state variables carefully.
match can have problems due to rounding errors, even if round and truncare used. Instead of testing for equality, it is generally better to check for minimum distance. This can be done for example with which.min

Putting this together, we could do it like follows:

library(deSolve)
set.seed(1)

deriv <- function(t, state, pars) {
  pars <- parameters[
    which.min(abs(DailyTemperature[floor(t + 1)] - parameters$TraitTemperature)), 
    1:5
  ]
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y

    list(c(x = d_x, y = d_y), 
         alpha = alpha, beta = beta, gamma = gamma, delta = delta, 
         ## some extra output for debugging
         DailyTemperature = DailyTemperature[floor(t + 1)],
         TraitTemperature = TraitTemperature
    )
  })
}


state <- c(x = 10, y = 5)
times = seq(0, 150, by = 1)
TraitTemperature = seq(0.1, 40, 0.1)

## here we assume a linear increase of the parameters with temperature
parameters <- data.frame(
  TraitTemperature = TraitTemperature,
  alpha = 1.0 + 0.01 * TraitTemperature + rnorm(40, mean = 0, sd = 0.01),
  beta =  0.2 + 0.01 * TraitTemperature + rnorm(40, mean = 0, sd = 0.01),
  delta = 0.2 + 0.01 * TraitTemperature + rnorm(40, mean = 0, sd = 0.01),
  gamma = seq(0.025, 1, 0.025)
)

## which.min will also work without rounding the temperature

#DailyTemperature <- round(runif(length(times), 0, 40), 1)
DailyTemperature <- runif(length(times), 0, 40)
DailyTemperature

out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)
head(out)

I want to add that there are also other (considerably faster) ways for table lookup.

Potential indexing issue in ODE?

1 Answers1