In antononcube/ECMMon-R: Epidemiologic Compartmental Modeling Monad

library(flexdashboard)
library(deSolve)
library(dplyr)
library(tidyr)
library(magrittr)
library(ggplot2)
library(ECMMon)

# load data in 'global' chunk so it can be shared by all users of the dashboard
modelSIR <- SEI2HRModel()

Column {.sidebar data-width=300}

Rates

sliderInput( inputId = "contactRateISSP", label = "ISSP contact rate:", min = 0, max = 6, value = 0.56, step = 0.01 )

sliderInput( inputId = "contactRateINSP", label = "INSP contact rate:", min = 0, max = 6, value = 0.56, step = 0.01 )

sliderInput( inputId = "deathRateISSP", label = "ISSP death rate:", min = 0, max = 1, value = 0.035, step = 0.001 )

sliderInput( inputId = "deathRateINSP", label = "INSP death rate:", min = 0, max = 1, value = 0.015, step = 0.001 )

sliderInput( inputId = "sspf", label = "Severely symptomatic population fraction:", min = 0, max = 1, value = 0.2, step = 0.01 )

sliderInput( inputId = "aip", label = "Average infectious period:", min = 1, max = 90, value = 28, step = 1 )

sliderInput( inputId = "aincp", label = "Average incubation period:", min = 1, max = 90, value = 6, step = 1 )

Quarantine parameters

sliderInput( inputId = "qsd", label = "Quarantine start date:", min = 0, max = 365, value = 60, step = 1 )

sliderInput( inputId = "ql", label = "Quarantine length:", min = 0, max = 365, value = 56, step = 1 )

sliderInput( inputId = "qcrf", label = "Quarantine contact rate factor:", min = 0, max = 2, value = 0.25, step = 0.01 )

Time horizon

sliderInput( inputId = "maxTime", label = "Max time:", min = 1, max = 365, value = 365, step = 1 )

Stocks to plot

selectInput( inputId = "stocksToPlot", 
             label = "Stocks to plot:", 
             choices = setNames( GetStocks( modelSIR, "^.*Population" ), modelSIR$Stocks[ GetStocks( modelSIR, "^.*Population" ) ]), 
             selected = setdiff( GetStocks( modelSIR, "^.*Population" ), c("TPt") ), 
             multiple = TRUE )

Appearance

checkboxInput( inputId = "moneyPlotLog10", label = "Log 10 money plots?", value = TRUE )

Column {data-width=650}

Populations plot

sol <- 
  reactive({
    times <- seq(0, input$maxTime, 1) 

    localParams <- as.list( modelSIR[["RateRules"]] )

    localParams[ c( "contactRateISSP", "contactRateINSP", "deathRateISSP", "deathRateINSP", "sspf") ] <- c( input$contactRateISSP, input$contactRateINSP, input$deathRateISSP / input$aip, input$deathRateINSP / input$aip, input$sspf)

    localParams[["contactRateINSP"]] <- function(t) { input$contactRateINSP * ifelse( input$qsd <= t && t <= input$qsd + input$ql, input$qcrf, 1 ) }
    localParams[["contactRateISSP"]] <- function(t) { input$contactRateISSP * ifelse( input$qsd <= t && t <= input$qsd + input$ql, input$qcrf, 1 ) }

    deSolve::ode(y = modelSIR[["InitialConditions"]], times = times, func = modelSIR[["RHSFunction"]], parms = localParams, method = "rk4" ) 
  })

dfSol <- 
  reactive({

    res <- as.data.frame(sol()) 
    colnames(res) <- gsub( "time", "Time", colnames(res) )

    res
  })

renderPlot( expr = {

  p <-
    dfSol() %>%
    tidyr::pivot_longer( cols = colnames(dfSol())[-1], names_to = "Stock", values_to = "Value" )%>% 
    dplyr::filter( Stock %in% input$stocksToPlot ) %>%
    dplyr::mutate( Stock = paste0( Stock, ", ", modelSIR$Stocks[ Stock ] ) ) %>% 
    ggplot2::ggplot( ) +
    ggplot2::geom_line( ggplot2::aes( x = Time, y = Value, color = Stock), lwd = 1.2 ) +
    ggplot2::geom_vline( xintercept = c( input$qsd, input$qsd + input$ql), color = "gray", show.legend = TRUE )

  print(p)
})

Money plots

renderPlot( expr = {

  p <-
    dfSol() %>%
    tidyr::pivot_longer( cols = colnames(dfSol())[-1], names_to = "Stock", values_to = "Value" )%>% 
    dplyr::filter( Stock %in% GetStocks( modelSIR, "^Money" ),  ) %>% 
    dplyr::mutate( Stock = paste0( Stock, ", ", modelSIR$Stocks[ Stock ] ))

  p <-
    if( input$moneyPlotLog10 ) {
      ggplot2::ggplot(p) +
        ggplot2::geom_line( aes( x = Time, y = log10(Value), color = Stock), lwd = 1.2 )
    } else {
      ggplot2::ggplot(p) +
        ggplot2::geom_line( aes( x = Time, y = Value, color = Stock), lwd = 1.2 )
    }  

  p <- p  + ggplot2::geom_vline( xintercept = c( input$qsd, input$qsd + input$ql), color = "gray", show.legend = TRUE )
  print(p)
})

Column {data-width=350}

Stock histograms

renderPlot( expr = {

  p <-
    dfSol() %>%
    tidyr::pivot_longer( cols = colnames(dfSol())[-1], names_to = "Stock", values_to = "Value" )%>% 
    ggplot2::ggplot( ) +
    ggplot2::geom_histogram( ggplot2::aes( x = Value), bins = 30 ) +
    ggplot2::facet_wrap( ~Stock, scales = "free", ncol = 2 )

  print(p)
})