Follow up: building a for-loop to run ODE simulation through different parameter sets in R

Question

I ended up trying to build a for-loop to achieve my end goal (running an ODE function modeling disease transmission through different scenarios/parameter sets in the most transparent way possible) but it's throwing the following warnings.

Note: Thank you to everyone for their initial help on my original post: How to loop differential equations through multiple parameter sets in R?. The suggestions I got worked well for my test code, but were either not transparent enough or were not generalizable to a larger set of parameters. I couldn't get either to work for my real data set, so I instead built a for-loop.

DLSODA-  At T (=R1), too much accuracy requested  
      for precision of machine..  See TOLSF (=R2) 
In above message, R1 = 1, R2 = nan
 
DLSODA-  At T (=R1), too much accuracy requested  
      for precision of machine..  See TOLSF (=R2) 
In above message, R1 = 1, R2 = nan
 
DLSODA-  At T (=R1), too much accuracy requested  
      for precision of machine..  See TOLSF (=R2) 
In above message, R1 = 1, R2 = nan

and

Warning messages:
1: In lsoda(y, times, func, parms, ...) :
  Excessive precision requested.  scale up `rtol' and `atol' e.g by the factor 10
2: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go
3: In lsoda(y, times, func, parms, ...) :
  Excessive precision requested.  scale up `rtol' and `atol' e.g by the factor 10
4: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go
5: In lsoda(y, times, func, parms, ...) :
  Excessive precision requested.  scale up `rtol' and `atol' e.g by the factor 10
6: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go

I have tried increasing rtol and atol by factors of 10, 100, and 1000, but that doesn't change the warning messages. Can someone help me get this for-loop running properly?

---

Here's the reproducible mock-up code I'm using to troubleshoot code mechanics before plugging the for-loop into my much larger model (note: the real model uses 13 parameters, so I am trying to make this for-loop able to run through parameter sets of any size).

library(deSolve)
library(ggplot2)
library(dplyr)
library(tidyverse)

parms = c(
"beta" = 0.00016,
"gamma" = 0.12
)

CoVode = function(t, x, parms) {
S =  x[1]  # Susceptible
I =  x[2]  # Infected
R =  x[3]  # Recovered

beta = parms["beta"]
gamma  = parms["gamma"]

dSdt <- -beta*S*I
dIdt <- beta*S*I-gamma*I
dRdt <- gamma*I

output = c(dSdt,dIdt,dRdt)
names(output) = c('S', 'I', 'R')

return(list(output))
}

# Initial conditions

init = numeric(3)
init[1] = 10000
init[2] = 1
names(init) = c('S','I','R')

# Run the model

odeSim = ode(y = init, 0:50, CoVode, parms)

# Plot results

Simdat <- data.frame(odeSim)

Simdat_long <-
Simdat %>%
pivot_longer(!time, names_to = "groups", values_to = "count")

ggplot(Simdat_long, aes(x= time, y = count, group = groups, colour = groups)) +
geom_line()

Now I need the code to run the model through scenarios A-I and save the model outputs in a data frame (so I can later use cowplot and ggplot to create a combined figure comparing the S-I-R dynamics for each scenario).

parmdf <- data.frame(scenario=c("A", "B", "C", "D", "E", "F", "G", "H", "I"),
beta=c(0.0001,0.0001, 0.0001, 0.001, 0.001, 0.001, 0.124, 0.124, 0.124),                       
gamma=c(0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3) )


# Attempted for-loop

for (i in parmdf) {
parms = as.numeric(unlist(parmdf[i,-1]))
names(parmdf) = names(parms)
odeSim5 = ode(y = init, 0:50, CoVode, parms)
}

print(odeSim5)

You're overwriting odeSim5 with each iteration of your loop. Surely you want to append to it? or, equivalently, use lapply to create a list of the individual scenario results and the bind them together? — Limey
– Limey, Commented Mar 21 at 15:21
for (i in parmdf) loops through columns of parmdf, not through rows. parmdf[i,-1] makes no sense, try to print(i) as the loop's 1st instruction. — Rui Barradas
– Rui Barradas, Commented Mar 21 at 19:32
@RuiBarradas I changed the scenarios from letters to numbers (i.e. "A" became 1), transposed the now-all-numerical dataframe (parmdf <- t(parmdf)) and updated the code such that for (i in parmdf) {print(i) parms = (parmdf[-1,i]) ode(y = init, 0:50, CoVode, parms) }. I'm still getting the error message DLSODA- At T (=R1), too much accuracy requested for precision of machine.. See TOLSF (=R2) In above message, R1 = 1, R2 = nan but 36 times! I can see from the using print(i) that the model is pulling the correct parameter values during each loop, so thank you. Any tips for errors? — Kit McLean
– Kit McLean, Commented Mar 21 at 21:24
@Limey I do need to figure that out! My main concern atm is getting the for-loop working correctly. It looks like /@RuiBarradas accomplished both in their answer! — Kit McLean
– Kit McLean, Commented Mar 21 at 21:35

Rui Barradas · Accepted Answer · 2025-03-21 19:51:44Z

Here is a way of running all scenarios and saving them in a named list.

First, create an empty list of length equal to the number of scenarios;
Name the list after column parmdf$scenario;
Then, for each row get the parameters in the named vector parms;
Fit the models.

parmdf <- data.frame(scenario=c("A", "B", "C", "D", "E", "F", "G", "H", "I"),
                     beta=c(0.0001,0.0001, 0.0001, 0.001, 0.001, 0.001, 0.124, 0.124, 0.124),                       
                     gamma=c(0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3) )

parms <- setNames(numeric(2L), c("beta", "gamma"))
models_list <- vector("list", nrow(parmdf)) %>% setNames(parmdf$scenario)
for(i in seq_len(nrow(parmdf))) {
  parms["beta"] <- parmdf$beta[i]
  parms["gamma"] <- parmdf$gamma[i]
  models_list[[i]] <- ode(y = init, 0:50, CoVode, parms)
}

plot_one_scenario <- function(odeSim, scenario) {
  lbl <- paste("Scenario:", scenario)
  odeSim %>%
    as.data.frame() %>%
    pivot_longer(!time, names_to = "groups", values_to = "count") %>%
    mutate(groups = factor(groups, levels = c('S','I','R'))) %>%
    ggplot(aes(x= time, y = count, group = groups, colour = groups)) +
    geom_line() +
    ggtitle(label = lbl)
}

map2(models_list, names(models_list), plot_one_scenario)

Collectives™ on Stack Overflow

Follow up: building a for-loop to run ODE simulation through different parameter sets in R

1 Answer 1

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related