How to plot several curves of the same variable for different values of a parameter using R?

Question

I am working with an SIR model using R. I need to plot multiple curves of I on the same figure for different values of beta, say for the values of 0.001, 0.002, 0.003, 0.004, and 0.005. Given below is the code that I have been working with so far. I know this might be a very simple problem, but I am new to R and couldn't find anything helpful yet.

library(deSolve)
sir_model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <-  beta * I * S - gamma * I
    dR <-  gamma * I
    return(list(c(dS, dI, dR)))
  })
}
parameters <- c(beta  = 0.001, gamma = 0.3)
initial_values <- c(S = 999, I = 1, R = 0)
time_series <- seq(0, 100)

sir_model_1 <- ode(
  y = initial_values,
  times = time_series,
  func = sir_model,
  parms = parameters
)

sir_model_1 <- as.data.frame(sir_model_1)


with(sir_model_1, { 
  plot(time, I, type = "l", col = "black",
       xlab = "time (days)", ylab = "Number of infections")     
})

I tried to use a for loop, but I think I am not doing it right.

Answer 1

Package deSolve contains a plotting function that supports to add multiple scenarios. It works directly with the output of ode that is a matrix of class deSolve. The first argument of the plot function needs to be such a deSolve object and the second can be a list of such objects.

This way, it can be run and plotted as follows:

library(deSolve)
sir_model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <-  beta * I * S - gamma * I
    dR <-  gamma * I
    return(list(c(dS, dI, dR)))
  })
}
parameters <- c(beta  = 0.001, gamma = 0.3)
initial_values <- c(S = 999, I = 1, R = 0)
time_series <- seq(0, 100)

run_sir <- function(beta) {
  parameters["beta"] <- beta
  ode(y = initial_values, times = time_series,func = sir_model, parms = parameters)
}

## define scenarios
beta <- 0.001 * 1:5

## run default scenario
run0 <- run_sir(beta = beta[1])
plot(run0)

## run other 4 scenarios
## [-1] means all except the first, that we already have
runs <- lapply(beta[-1], run_sir)

plot(run0, runs, las=1)
legend("bottomright", legend = paste("beta = ", beta), lty=1:5, col=1:5)

The plot.deSolvefunction is highly configurable with respect to layout, colors, line types, selection of variables, etc., see help page ?plot.deSolve.

To plot only selected state variables, one can use the which argument, e.g. which="I" to plot only the infected.

plot(run0, runs, which="I")

Some examples about plotting deSolve outputs can be found in the tutorial slides in section "Plotting".

Answer 2

We could wrap the code in a function and loop over the sequence of beta values, and plot. If we need to do this in a single plot window, modify the par

par(mfrow = c(5, 1))
lapply(seq(0.001, length.out = 5, by = 0.001), f1)

where

f1 <- function(beta) {

parameters <- c(beta  = beta, gamma = 0.3)
initial_values <- c(S = 999, I = 1, R = 0)
time_series <- seq(0, 100)

sir_model_1 <- ode(
  y = initial_values,
  times = time_series,
  func = sir_model,
  parms = parameters
)

sir_model_1 <- as.data.frame(sir_model_1)


with(sir_model_1, { 
  plot(time, I, type = "l", col = "black",
       xlab = "time (days)", main = beta, ylab = "Number of infections")     
})

}

-output

---

In case we want to loop over the column names, use the formula method

par(mfrow = c(3, 1))
lapply(names(sir_model_1)[-1], function(nm) 
   plot(reformulate("time", response = nm), data = sir_model_1, main = nm,
      type = "l", col = "black",
        xlab = "time (days)", ylab = "Number of infections"))

-output