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

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)

Solution

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)