SIR model in Rstudio shiny

Question

I´m trying to build the basic SIR model in Rstudio shiny. The model takes 2 parameters (beta = infection rate/day, gamma = recovery date/day), 3 initial values (S = numbers of susceptibles, I = infectious, R = recovered) and last variable is time (in days).

Here is the code of it just in R markdown:

library(deSolve)

sir_equations <- 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_values <- c(
  beta  = 0.05, # infectious rate/day
  gamma = 0.5    # recovery rate/day
)

initial_values <- c(
  S = 1000,  # susceptibles
  I =   1,  # infectious
  R =   0   # recovered (immune)
)

time_values <- seq(0, 10) #number of days (0-10)

sir_values_1 <- ode(
  y = initial_values,
  times = time_values,
  func = sir_equations,
  parms = parameters_values 
)

sir_values_1 <- as.data.frame(sir_values_1) # convert to data frame

with(sir_values_1, {
  plot(time, S, type = "l", col = "blue",
       xlab = "period (days)", ylab = "number of people")
  lines(time, I, col = "red")
  lines(time, R, col = "green")
})

legend("right", c("susceptibles", "infectious", "recovered"),
       col = c("blue", "red", "green"), lty = 1, bty = "n")

Now I want to add this into R shiny, where the user can input the beta, gamma and days value (sliderbar, or just input), then it will plot the result. I´m pretty new to R and tried some variations here, like putting the user input into ,,UI,, the calculating into ,,server,, then combine it in like this shinyApp(ui = ui, server = server). This code below I tried, but its not working. Can you guys help me, what I´m doing wrong, and what to follow to be able to put the code into R shiny?

library(deSolve)
library(shiny)


ui <- fluidPage(




  sliderInput(inputId = "time_values", label = "Dny", value = 10, min = 1, max = 100),
  sliderInput(inputId = "beta", label ="Míra nákazy", value = 0.05, min = 0.00, max = 1, step = 0.01),
  sliderInput(inputId = "gamma", label ="Míra uzdravení", value = 0.5, min = 0.00, max = 1, step = 0.1),

  plotOutput("plot")
)



server <- function(input, output) {
  sir_equations <- 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)))
  })
  }

  initial_values <- c(S = 1000, I = 1, R = 0)

  sir_values_1 <- ode(
  y = initial_values,
  times = time_values,
  func = sir_equations,
  parms = parameters_values 
)

  output$plot <- renderPlot({
    plot(rnorm(input$time_values))
    plot(rnorm(input$beta))
    plot(rnorm(input$gamma))
  })


}

shinyApp(ui = ui, server = server)

Thanks Michal

Answer 1

I guess it is something like this you want?

library(deSolve)
library(shiny)

ui <- fluidPage(
  sliderInput(inputId = "time_values", label = "Dny", value = 10, min = 1, max = 100),
  sliderInput(inputId = "beta", label ="Míra nákazy", value = 0.05, min = 0, max = 1, step = 0.01),
  sliderInput(inputId = "gamma", label ="Míra uzdravení", value = 0.5, min = 0, max = 1, step = 0.1),

  plotOutput("plot")
)

server <- function(input, output) {
  sir_equations <- 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)))
    })
  }

  initial_values <-  c(S = 1000, I = 1, R = 0)

  sir_values_1 <- reactiveValues(val = data.frame())

  observe({
    sir_values_1$val <- as.data.frame(ode(
      y = initial_values,
      times = seq(0, input$time_values),
      func = sir_equations,
      parms = c(beta=input$beta, gamma=input$gamma) 
    ))
  })

  output$plot <- renderPlot({
    with(sir_values_1$val, {
    plot(sir_values_1$val$time, sir_values_1$val$S, type = "l", col = "blue",
         xlab = "period (days)", ylab = "number of people")
    lines(sir_values_1$val$time, sir_values_1$val$I, col = "red")
    lines(sir_values_1$val$time, sir_values_1$val$R, col = "green")
    legend("right", c("susceptibles", "infectious", "recovered"),
           col = c("blue", "red", "green"), lty = 1, bty = "n")
    })
  })
}

shinyApp(ui = ui, server = server)

Answer 2

Here another solution without the need of an observer function. More about deSolve and shiny at: https://tpetzoldt.github.io/deSolve-shiny/deSolve-shiny.html

library("deSolve")

sir_equations <- 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)))
  })
}

ui <- fluidPage(
  sliderInput(inputId = "time_values", label = "Dny", value = 10, min = 1, max = 100),
  sliderInput(inputId = "beta", label ="Míra nákazy", value = 0.05, min = 0.00, max = 1, step = 0.01),
  sliderInput(inputId = "gamma", label ="Míra uzdravení", value = 0.5, min = 0.00, max = 1, step = 0.1),

  plotOutput("plot")
)

server <- function(input, output) {
  output$plot <- renderPlot({
    initial_values <- c(S = 1000, I = 1, R = 0)
    sir_values <- ode(
      y = initial_values,
      times = seq(0, input$time_values, length.out=1000),
      func = sir_equations,
      parms = c(beta=input$beta, gamma=input$gamma)
    )

    ## easiest is to use the deSolve plot function
    #plot(sir_values, mfrow=c(1,3))
    ## but you can also do it with own plot functions, e.g.:
    matplot(sir_values[,1], sir_values[,-1], type="l", xlab="time", ylab="S, I, R")
    legend("topright", col=1:3, lty=1:3, legend=c("S", "I", "R"))
  })
}

shinyApp(ui = ui, server = server)

Answer 3

Just look at the error:

Warning: Error in ode: objet 'time_values' introuvable

In ode(), you should replace time_values by input$time_values and put the full ode() function in a reactive environment since you use some inputs:

  sir_values_1 <- reactive({
    ode(
    y = initial_values,
    times = input$time_values,
    func = sir_equations,
    parms = parameters_values 
  )
  })

Then you have some errors in your plot but setting xlim and ylim should make it work. However, if you want to display multiple plots, you must define several plotOutput and renderPlot. Putting three plot in one renderPlot will not display the three of them but only the last one.