Introduction to TB Models

A shiny app reproducing the models used in the Introduction to Tuberculosis modelling course practicals, run by TB MAC at the 2017 Union conference. See the TB MAC website for course materials and further resources. The models used in this course, and reproduced in this shiny app, were based on one published by Lin et al.

Installing the shiny app locally

Manual Install

To install and run the shiny app locally on your own computer you will need to first install R, it is also suggested that you install Rstudio. After downloading the source code from this repository click on the intro_to_tb_models.Rprof file, this will open an Rstudio window. Type the following code into the command line;

install.packages("shiny")
install.packages("shinydashboard")
install.packages("DT")
install.packages("tidyverse")
install.packages("rmarkdown")
install.packages("plotly")

To run the app open the ui.R file and press run, depending on your computer this may take some time.

Using Docker

Docker is a container software that seeks to eliminate “works on my machine” issues. For installation and set up instructions see here.

This docker container is based on the shiny docker image, see here for instructions on use. To run the docker image run the following in a bash shell:

docker pull seabbs/intro_to_tb_models
docker run --rm -p 3838:3838 seabbs/intro_to_tb_models

The shiny server can be found on port :3838 at your local machines ip (or localhost on windows), fcdashboard can be found at your-ip:3838/intro_to_tb_models.

UI

Download

## Load packages
library(shiny)
library(shinydashboard)
library(tidyverse)
library(rmarkdown)
library(plotly)
library(DT)


##Load code
source("TB_model.R")
source("graph_tb_model.R")
source("TB_model_diag.R")
source("graph_tb_model_diag.R")


sidebar <- dashboardSidebar(
  hr(),
  sidebarMenu(id = "menu",
              menuItem("TB model", tabName = "tb-model", icon = icon("line-chart")),
              menuItem("TB model with diagnostics", tabName = "tb-model-diag", icon = icon("line-chart")),
              menuItem("About", tabName = "readme", icon = icon("mortar-board"), selected = TRUE),
              menuItem("Code",  icon = icon("code"),
                       menuSubItem("Github", href = "https://github.com/seabbs/intro_to_tb_models", icon = icon("github")),
                       menuSubItem("ui.R", tabName = "ui", icon = icon("angle-right")),
                       menuSubItem("server.R", tabName = "server", icon = icon("angle-right")),
                       menuSubItem("TB_model.R", tabName = "tb_model", icon = icon("angle-right")),
                       menuSubItem("TB_model_diag.R", tabName = "tb_model_diag", icon = icon("angle-right")),
                       menuSubItem("graph_tb_model.R", tabName = "graph_tb_model", icon = icon("angle-right"))
                       
              )
  ),
  conditionalPanel(condition = 'input.menu == "tb-model"',
                   sliderInput("ecr",
                               "Effective Contacts (per year)",
                               min = 0, 
                               max = 20,
                               value = 15,
                               step = 0.5),
                   sliderInput("wks_inf_p", 
                               "Avg. No. of Weeks Infectious (Positive Sputum Smear)",
                               min = 0,
                               max = 102,
                               value = 51 ),
                   sliderInput("wks_inf_n", 
                               "Avg. No. of Weeks Infectious (Negative Sputum Smear)",
                               min = 0,
                               max = 190,
                               value = 95),
                   sliderInput("prot_init_reint",
                               "Protection from reinfection relative to infection",
                               min = 0, 
                               max = 1, 
                               value = 0.65),
                   selectInput("timestep",
                               "Set timestep",
                               choices = list(
                                 Week = 52,
                                 Month = 12,
                                 Year = 1),
                               selected = 1
                   ),
                   checkboxInput("burn_in",
                                 "Show burn in",
                                 value = FALSE)
  ),
  conditionalPanel(condition = 'input.menu == "tb-model-diag"',
                   sliderInput("ecr_diag",
                               "Effective Contacts (per year)",
                               min = 0, 
                               max = 20,
                               value = 15,
                               step = 0.5),
                   checkboxInput("intervention",
                                 "Intervention",
                                 value = TRUE),
                   selectInput("timestep_diag",
                               "Set timestep",
                               choices = list(
                                 Week = 52,
                                 Month = 12,
                                 Year = 1),
                               selected = 1
                   ),
                   checkboxInput("burn_in_diag",
                                 "Show burn in",
                                 value = FALSE)
  ),
  hr(),
  helpText("Developed by ", a("Sam Abbott", href = "http://samabbott.co.uk"), 
           style = "padding-left:1em; padding-right:1em;position:absolute; bottom:1em; ")
)

body <- dashboardBody(
  tabItems(
    tabItem(tabName = "tb-model",
            fluidRow(
              tabBox(width = 12, 
                     title = "Model Plots", 
                     side = "right",
                     tabPanel(title = "Rates",
                              plotlyOutput("plot_rates")
                              ),
                     tabPanel(title = "Annual Risk",
                              plotlyOutput("plot_annual_risk")
                     ),
                     tabPanel(title = "Proportions",
                              plotlyOutput("plot_props")
                     )
              ),
                     tabBox(width = 12, 
                            title = "Model Statistics", 
                            side = "right",
                            tabPanel(title = "Rates and Annual Risk",
                                     DT::dataTableOutput("rates_table")
                            ),
                            tabPanel(title = "Proportions", 
                                     DT::dataTableOutput("props_table")
                            ),
                            tabPanel(title = "Trajectory",
                                     DT::dataTableOutput("traj_table")
                            )
                            )
              )
            ),
    tabItem(tabName = "tb-model-diag",
            fluidRow(
              tabBox(width = 12, 
                     title = "Model Plots", 
                     side = "right",
                     tabPanel(title = "Rates",
                              plotlyOutput("plot_rates_diag")
                     ),
                     tabPanel(title = "Annual Risk",
                              plotlyOutput("plot_annual_risk_diag")
                     ),
                     tabPanel(title = "Cumulative Counts",
                              plotlyOutput("plot_cum")
                     ),
                     tabPanel(title = "Diagnoses",
                              plotlyOutput("plot_diag")
                     ),
                     tabPanel(title = "Proportions",
                              plotlyOutput("plot_props_diag")
                     )
              ),
              tabBox(width = 12, 
                     title = "Model Statistics", 
                     side = "right",
                     tabPanel(title = "Rates and Annual Risk",
                              DT::dataTableOutput("rates_table_diag")
                     ),
                     tabPanel(title = "Cumulative Counts",
                              DT::dataTableOutput("cum_table")
                     ),
                     tabPanel(title = "Diagnoses",
                              DT::dataTableOutput("diag_table")
                     ),
                     tabPanel(title = "Proportions", 
                              DT::dataTableOutput("props_table_diag")
                     ),
                     tabPanel(title = "Trajectory",
                              DT::dataTableOutput("traj_table_diag")
                     )
              )
            )
    ),
    tabItem(tabName = "readme",
            withMathJax(), 
            includeMarkdown("README.md")
    ),
    tabItem(tabName = "ui",
            box( width = NULL, status = "primary", solidHeader = FALSE, title = "UI",
                 downloadButton('downloadData2', 'Download'),
                 br(),br(),
                 pre(includeText("ui.R"))
            )
    ),
    tabItem(tabName = "server",
            box( width = NULL, status = "primary", solidHeader = FALSE, title = "Server",
                 downloadButton('downloadData3', 'Download'),
                 br(),br(),
                 pre(includeText("server.R"))
            )
            
    ),
    tabItem(tabName = "tb_model",
            box( width = NULL, status = "primary", solidHeader = FALSE, title = "TB Model",
                 downloadButton('downloadData4', 'Download'),
                 br(),br(),
                 pre(includeText("TB_model.R"))
            )
            
    ),
    tabItem(tabName = "graph_tb_model",
            box( width = NULL, status = "primary", solidHeader = FALSE, title = "TB Model Graphs",
                 downloadButton('downloadData5', 'Download'),
                 br(),br(),
                 pre(includeText("graph_tb_model.R"))
            )
            
    ),
    tabItem(tabName = "tb_model_diag",
            box( width = NULL, status = "primary", solidHeader = FALSE, title = "TB Model with Diagnostics",
                 downloadButton('downloadData6', 'Download'),
                 br(),br(),
                 pre(includeText("TB_model_diag.R"))
            )
            
    ),
    tabItem(tabName = "graph_tb_model_diag",
            box( width = NULL, status = "primary", solidHeader = FALSE, title = "TB Model with Diagnostics Graphs",
                 downloadButton('downloadData7', 'Download'),
                 br(),br(),
                 pre(includeText("graph_tb_model_diag.R"))
            )
            
    )
  )
)

dashboardPage(
  dashboardHeader(title = "Intro to TB Models"),
  sidebar,
  body,
  skin = "black"
)

Server

Download

#Load packages
library(shiny)
library(shinydashboard)
library(tidyverse)
library(plotly)
library(DT)
library(deSolve)


##Load code
source("TB_model.R")
source("graph_tb_model.R")
source("TB_model_diag.R")
source("graph_tb_model_diag.R")

## Stop spurious warnings
options(warn = -1)

shinyServer(function(input, output) {
  
  ## Run TB model
  model_traj <- reactive({
    ##Check model runs
    initial <- init_TB_model()
    times <- seq(0, 2024*as.numeric(input$timestep), 1)
    params <- params_TB_model(ecr_pyr = input$ecr, wks_infect_n = input$wks_inf_n,
                              wks_infect_p = input$wks_inf_p,  prot_reinf = input$prot_init_reint,
                              timestep = as.numeric(input$timestep))
    
    model_traj <- deSolve::lsoda(initial, times, TB_model, params)
    
  })
  
  ## summarise TB model rates
  sum_rates <- reactive({
    if (input$burn_in) {
      start_time <- 0
    }else{
      start_time <- 2010
    }
    
    sum_model_rates(model_traj(), start_time = start_time)
  })
  
  sum_props <- reactive({
    if (input$burn_in) {
      start_time <- 0
    }else{
      start_time <- 2010
    }
    
    sum_model_proportions(model_traj(), start_time = start_time)
  })
  
  ## Plot rates
  output$plot_rates <- renderPlotly({
    plot <- sum_rates() %>% 
      graph_rates
    
    ggplotly(plot)
  })
  
  ##Annual risk plot
  output$plot_annual_risk <- renderPlotly({
    plot <- sum_rates() %>% 
      graph_annual_risk()
    
    ggplotly(plot)
  })
  
  ## Rates as table
  output$rates_table <- DT::renderDataTable({
    sum_rates()},
    options = list(
      pageLength = 5,
      scrollX = TRUE,
      scrollY = TRUE,
      orderClasses = TRUE)
  )
  
  ## Plot props
  output$plot_props <- renderPlotly({
    plot <- sum_props() %>% 
      graph_proportions
    
    ggplotly(plot)
  })
  ##Prop as table
  output$props_table <- DT::renderDataTable({
    sum_props()
  },
    options = list(
                  pageLength = 5,
                  scrollX = TRUE,
                  scrollY = TRUE,
                  orderClasses = TRUE)
  )
  
  ##Model Trajectories as table
  output$traj_table <- DT::renderDataTable({
    model_traj()
  },
  options = list(
    pageLength = 5,
    scrollX = TRUE,
    scrollY = TRUE,
    orderClasses = TRUE)
  )
  
  
  ## TB model with diagnostics
  
  ## Run TB model
  diag_params <- reactive({
    if (input$intervention) {
      intervention <- 1
    } else{
      intervention <- 0
    }
    
    params <- params_TB_model_diag(ecr_pyr = input$ecr_diag,  wks_health_service = 52,
                                   prot_reinf = 0.65, intervention = intervention, 
                                   intervention_start = 2014,
                                   timestep = as.numeric(input$timestep_diag))
  })
  
  model_traj_diag <- reactive({

    ##Check model runs
    initial <- init_TB_model()
    times <- seq(0, 2024*as.numeric(input$timestep_diag), 1)
    model_traj <- deSolve::lsoda(initial, times, TB_model_diag, diag_params())
    
  })
  
  ## summarise TB model rates
  sum_rates_diag <- reactive({
    if (input$burn_in) {
      start_time <- 0
    }else{
      start_time <- 2010
    }
    
    sum_model_rates(model_traj_diag(), start_time = start_time)
  })
  
  sum_props_diag <- reactive({
    if (input$burn_in_diag) {
      start_time <- 0
    }else{
      start_time <- 2010
    }
    
    sum_model_proportions(model_traj_diag(), start_time = start_time)
  })
  
  ## Plot rates
  output$plot_rates_diag <- renderPlotly({
    plot <- sum_rates_diag() %>% 
      graph_rates
    
    ggplotly(plot)
  })
  
  ##Annual risk plot
  output$plot_annual_risk_diag <- renderPlotly({
    plot <- sum_rates_diag() %>% 
      graph_annual_risk()
    
    ggplotly(plot)
  })
  
  ## Rates as table
  output$rates_table_diag <- DT::renderDataTable({
    sum_rates_diag()},
    options = list(
      pageLength = 5,
      scrollX = TRUE,
      scrollY = TRUE,
      orderClasses = TRUE)
  )
  
  ## Plot props
  output$plot_props_diag <- renderPlotly({
    plot <- sum_props_diag() %>% 
      graph_proportions
    
    ggplotly(plot)
  })
  ##Prop as table
  output$props_table_diag <- DT::renderDataTable({
    sum_props_diag()
  },
  options = list(
    pageLength = 5,
    scrollX = TRUE,
    scrollY = TRUE,
    orderClasses = TRUE)
  )
  
  ##Model Trajectories as table
  output$traj_table_diag <- DT::renderDataTable({
    model_traj_diag()
  },
  options = list(
    pageLength = 5,
    scrollX = TRUE,
    scrollY = TRUE,
    orderClasses = TRUE)
  )

  ## Cumulative measures
  cum_traj <- reactive({

    cum_measures(model_traj_diag(), diag_params())
  })
  
  ##Cum as table
  output$cum_table <- DT::renderDataTable({
    cum_traj()
  },
  options = list(
    pageLength = 5,
    scrollX = TRUE,
    scrollY = TRUE,
    orderClasses = TRUE)
  )
  
  ## Plot cum
  output$plot_cum <- renderPlotly({
    plot <- cum_traj() %>% 
      graph_cum
    
    ggplotly(plot)
  })
  
  ## Diagnoses measures
  sum_diag_traj <- reactive({
    if (input$burn_in_diag) {
      start_time <- 0
    }else{
      start_time <- 2010
    }
    
    sum_diag(model_traj_diag(), start_time = start_time)
  })
  
  ##Cum as table
  output$diag_table <- DT::renderDataTable({
    sum_diag_traj()
  },
  options = list(
    pageLength = 5,
    scrollX = TRUE,
    scrollY = TRUE,
    orderClasses = TRUE)
  )
  
  ## Plot cum
  output$plot_diag <- renderPlotly({
    plot <- sum_diag_traj() %>% 
      graph_diag
    
    ggplotly(plot)
  })
  
  output$downloadData2 <- downloadHandler(filename = "ui.R",
                                          content = function(file) {
                                            file.copy("ui.R", file, overwrite = TRUE)
                                            }
                                          )
  output$downloadData3 <- downloadHandler(filename = "server.R",
                                          content = function(file) {
                                            file.copy("server.R", file, overwrite = TRUE)
                                            }
                                          )
  output$downloadData4 <- downloadHandler(filename = "TB_model.R",
                                          content = function(file) {
                                            file.copy("TB_model.R", file, overwrite = TRUE)
                                          }
  )
  output$downloadData5 <- downloadHandler(filename = "graph_tb_model.R",
                                          content = function(file) {
                                            file.copy("graph_tb_model.R", file, overwrite = TRUE)
                                          }
  )
  
  output$downloadData6 <- downloadHandler(filename = "tb_model_diag.R",
                                          content = function(file) {
                                            file.copy("tb_model_diag.R", file, overwrite = TRUE)
                                          }
  )
  
  output$downloadData7 <- downloadHandler(filename = "graph_tb_model_diag.R",
                                          content = function(file) {
                                            file.copy("graph_tb_model_diag.R", file, overwrite = TRUE)
                                          }
  )
  
})

TB Model

Download

#Packages
library(tidyverse)

## Model equations
TB_model <- function(t, x, params) {
  
  ## Specify model compartments
  S <- x[1]
  H <- x[2]
  L <- x[3]
  L_r <- x[4] 
  I_n <- x[5]
  I_p <- x[6]
  R <- x[7]
  
  with(as.list(params),{
    
    ## Specify total population
    N = S + H + L + L_r + I_n + I_p + R
    
    ## Force of infection
    foi = ecr * (rel_inf_n * I_n + I_p)/N

    ## Total infected
    total_infected <-  prim_dis_onset * H + sec_dis_onset * L + reinf_dis_onset * L_r + relapse * R
    
    ##Births
    births <- mu * N + mu_n * I_n + mu_p * I_p
    
    ## Derivative Expressions
    dS = births - foi * S - mu * S
    dH = foi * S - prim_dis_onset * H - rate_low_latent * H - mu * H
    dL = rate_low_latent * (H + L_r) - foi * L - sec_dis_onset * L - mu * L
    dL_r = foi * L + foi * R - reinf_dis_onset * L_r - rate_low_latent * L_r - mu * L_r
    dI_n = prop_n * total_infected - natural_cure * I_n - detect_recover_n * I_n - (mu + mu_n) * I_n
    dI_p = prop_p * total_infected - natural_cure * I_p - detect_recover_p * I_p - (mu + mu_p) * I_p
    dR = natural_cure * (I_n + I_p) + detect_recover_n * I_n + detect_recover_p * I_p - relapse * R - foi * R - mu * R
    
    ## derivatives
    derivatives <- c(dS, dH, dL, dL_r, dI_n, dI_p, dR)
    
    ##summary measures
    summary_measures <- c(
      year = t/timestep,
      ann_risk_inf = 100 * foi * timestep,
      total_deaths = births,
      total_inf = I_n + I_p,
      total_new_cases = prim_dis_onset * H + sec_dis_onset * L + reinf_dis_onset * L_r,
      total_new_deaths = mu_n * I_n + mu_p * I_p,
      total_new_infs = foi * S,
      total_new_prim_inf = prim_dis_onset * H,
      total_new_reinf_inf = reinf_dis_onset * L_r,
      total_new_react = sec_dis_onset * L)
    
      summary_measures <- c(summary_measures,
                            prop_dis_recent_inf = summary_measures["total_new_prim_inf.H"]/summary_measures["total_new_cases.H"],
                            prop_dis_recent_reinf = summary_measures["total_new_reinf_inf.L_r"]/summary_measures["total_new_cases.H"],
                            prop_dis_react = summary_measures["total_new_react.L"]/summary_measures["total_new_cases.H"],
                            new_cases_react_pHK = sec_dis_onset * L * 100000 * timestep / N,
                            new_cases_recent_ing_pHK = prim_dis_onset * H * 100000 * timestep / N,
                            new_cases_recent_reinf_pHK = reinf_dis_onset * L_r * 100000 * timestep / N,
                            ann_TB_incidence_pHK = summary_measures["total_new_cases.H"] * 100000 * timestep / N,
                            TB_mort_rate_pHK = summary_measures["total_new_deaths.I_n"] * 100000 * timestep / N,
                            TB_prevalence = (I_n + I_p)/N
      )
  
    ## output
    list(derivatives, summary_measures)
  })
}


## initial model
init_TB_model <- function(total_pop = 100000) {
  H <- 10
  L <- 10
  L_r <- 10
  I_n <- 50
  I_p <- 50
  R <- 10
  S <- total_pop - H - L - L_r - I_n - I_p - R
  
  return(c(S = S, H = H, L = L, L_r = L_r, I_n = I_n, I_p = I_p, R = R))
}


## Set parameters for model
params_TB_model <- function(ecr_pyr = 15, wks_infect_n = 95,
                            wks_infect_p = 51, prot_reinf = 0.65,
                            timestep = 52) {
  
  ## yearly rate of disease onset
  onset_yr <- 0.03

  ## proportion that are smear negative
  prop_p <- 0.7
  
  ## Mortality rates
  mu <- 0.021
  mu_n <- 0.19
  mu_p <- 0.22
  
  ## Cure rate
  natural_cure <- 0.2
  
  ## parameters
  params <- list(ecr = ecr_pyr / timestep,
                 rate_low_latent = 1 / (5 * timestep),
                 rel_inf_n = 0.22,
                 prim_dis_onset = onset_yr / timestep, 
                 sec_dis_onset  = 0.0003/timestep,
                 reinf_dis_onset = (onset_yr * (1 - prot_reinf))/timestep,
                 prop_n = 1 - prop_p, 
                 prop_p = prop_p,
                 mu = mu/timestep,
                 mu_n = mu_n/timestep,
                 mu_p = mu_p/timestep,
                 detect_recover_n = 52/(wks_infect_n * timestep) - (natural_cure + mu + mu_n)/timestep,
                 detect_recover_p = 52/(wks_infect_p * timestep)  - (natural_cure + mu + mu_p)/timestep,
                 natural_cure = natural_cure/timestep,
                 relapse = 0.001/timestep,
                 timestep = timestep
                 )
  return(params)
}

                      
##Summarise model rates
sum_model_rates <- function(model_traj, start_time = 2010, round = NULL) {
  model_traj <- as.data.frame(model_traj) %>% as_tibble
  ## Filter based on start time
  model_traj <- model_traj[unlist(model_traj[["year"]]) >= start_time, ]
  
  with(model_traj,{
    df <- data_frame(Year = year,
                     `Annual TB incidence rate` = ann_TB_incidence_pHK.total_new_cases.H,
                     `TB mortality rate` = TB_mort_rate_pHK.total_new_deaths.I_n,
                     `Annual risk of infection` = ann_risk_inf.I_n)
  })
  
}


##Sumarise model proportions
sum_model_proportions <- function(model_traj, start_time = 2010) {
  model_traj <- as.data.frame(model_traj) %>% as_tibble
  
  ## Filter based on start time
  model_traj <- model_traj[unlist(model_traj[["year"]]) >= start_time, ]
  
  with(model_traj,{
    df <- data_frame(Year = year,
                     `Recent infection` = prop_dis_recent_inf.total_new_prim_inf.H,
                     `Recent reinfection` = prop_dis_recent_reinf.total_new_reinf_inf.L_r,
                     `Reactivation` = prop_dis_react.total_new_react.L)
  })
  
}

                                                                                                                                                 
                                                                                                                                                 
                                                                                                                                                 
                                                                                                                                                 
                                                                                                                                                 
                                                                                                                                                 
                                                                                                                                                 

TB Model Graphs

Download

library(tidyverse)

graph_rates <- function(df){
  plot <- df %>%
    gather(key = "Summary measure", value = "Incidence rate (per 100,000)",
           `Annual TB incidence rate`, `TB mortality rate`) %>% 
    ggplot(aes(x = Year, y = `Incidence rate (per 100,000)`,
                      colour = `Summary measure`, group = `Summary measure`)) +
    geom_line() +
    theme_minimal() +
    theme(legend.position = "bottom") +
    expand_limits(y = 0)
  
  return(plot)  
}

graph_annual_risk <- function(df) {
  plot <-  df %>% 
    ggplot(aes(x = Year, y = `Annual risk of infection`)) +
    geom_line() +
    theme_minimal() +
    theme(legend.position = "none") +
    expand_limits(y = 0)
  
  return(plot)
}

graph_proportions <- function(df) {
  plot <- df %>% 
    gather(key = "Summary measures", value = "Proportions",
           `Recent infection`, `Recent reinfection`, Reactivation) %>% 
    ggplot(aes(x = Year, y = Proportions, colour = `Summary measures`, group = `Summary measures`)) +
    geom_line() + 
    theme_minimal() +
    theme(legend.position = "bottom") +
    expand_limits(y = 0)
  
  return(plot)
}

TB Model with Diagnostics

Download

#Packages
library(tidyverse)

## Model equations
TB_model_diag <- function(t, x, params) {
  
  ## Specify model compartments
  S <- x[1]
  H <- x[2]
  L <- x[3]
  L_r <- x[4] 
  I_n <- x[5]
  I_p <- x[6]
  R <- x[7]
  
  with(as.list(params),{
    
    ## Specify total population
    N = S + H + L + L_r + I_n + I_p + R
    
    ## Force of infection
    foi = ecr * (rel_inf_n * I_n + I_p)/N

    ## Total infected
    total_infected <-  prim_dis_onset * H + sec_dis_onset * L + reinf_dis_onset * L_r + relapse * R
    
    ##Births
    births <- mu * N + mu_n * I_n + mu_p * I_p
    
    ## Evaluate if intervention is in place
    if (intervention == 1 && t/timestep >= intervention_start) { 
      detect_recover_n <- interv_detect_recover_n 
      detect_only_n <- interv_detect_only_n
    } else{
        detect_recover_n <- preinterv_detect_recover_n 
        detect_only_n <- preinterv_detect_only_n 
    }
    
    ## Derivative Expressions
    dS = births - foi * S - mu * S
    dH = foi * S - prim_dis_onset * H - rate_low_latent * H - mu * H
    dL = rate_low_latent * (H + L_r) - foi * L - sec_dis_onset * L - mu * L
    dL_r = foi * L + foi * R - reinf_dis_onset * L_r - rate_low_latent * L_r - mu * L_r
    dI_n = prop_n * total_infected - natural_cure * I_n - detect_recover_n * I_n - (mu + mu_n) * I_n
    dI_p = prop_p * total_infected - natural_cure * I_p - detect_recover_p * I_p - (mu + mu_p) * I_p
    dR = natural_cure * (I_n + I_p) + detect_recover_n * I_n + detect_recover_p * I_p - relapse * R - foi * R - mu * R
    
    ## derivatives
    derivatives <- c(dS, dH, dL, dL_r, dI_n, dI_p, dR)
    
    ##summary measures
    summary_measures <- c(
      year = t/timestep,
      ann_risk_inf = 100 * foi * timestep,
      total_deaths = births,
      total_inf = I_n + I_p,
      total_new_cases = prim_dis_onset * H + sec_dis_onset * L + reinf_dis_onset * L_r,
      total_new_deaths = mu_n * I_n + mu_p * I_p,
      total_new_infs = foi * S,
      total_new_prim_inf = prim_dis_onset * H,
      total_new_reinf_inf = reinf_dis_onset * L_r,
      total_new_react = sec_dis_onset * L)
    
      summary_measures <- c(summary_measures,
                            prop_dis_recent_inf = summary_measures["total_new_prim_inf.H"]/summary_measures["total_new_cases.H"],
                            prop_dis_recent_reinf = summary_measures["total_new_reinf_inf.L_r"]/summary_measures["total_new_cases.H"],
                            prop_dis_react = summary_measures["total_new_react.L"]/summary_measures["total_new_cases.H"],
                            new_cases_react_pHK = sec_dis_onset * L * 100000 * timestep / N,
                            new_cases_recent_ing_pHK = prim_dis_onset * H * 100000 * timestep / N,
                            new_cases_recent_reinf_pHK = reinf_dis_onset * L_r * 100000 * timestep / N,
                            ann_TB_incidence_pHK = summary_measures["total_new_cases.H"] * 100000 * timestep / N,
                            TB_mort_rate_pHK = summary_measures["total_new_deaths.I_n"] * 100000 * timestep / N,
                            TB_prevalence = (I_n + I_p)/N
      )
      
      ## Intervention measures
      intervention_measures <- c(
        total_new_diag_n  = I_n * detect_only_n,
        total_new_diag_p = I_p * detect_only_p,
        smear_exams_p = I_p * test_TBcases,
        smear_exams_n = I_n * test_TBcases,
        smear_exams_no_TBcases = (1 - summary_measures["TB_prevalence.I_n"]) * test_non_TBcases * N
      )
      
      ## Time varying intervention measures
      if (intervention == 1 && t/timestep >= intervention_start) { 
        intervention_measures <- c(intervention_measures,
                                   Xray_exams_n = prop_access_Xray * intervention_measures["smear_exams_n.I_n"],
                                   Xray_exams_non_TBcases = prop_access_Xray * intervention_measures["smear_exams_non_TBcases.TB_prevalence.I_n"],
                                   new_test_n.smear_exams_n.I_n = 0,
                                   new_test_non_TBcases.smear_exams_non_TBcases.TB_prevalence.I_n = 0
                                   )
      } else{
        intervention_measures <- c(intervention_measures,
                                   Xray_exams_n.smear_exams_n.I_n = 0,
                                   Xray_exams_non_TBcases.smear_exams_non_TBcases.TB_prevalence.I_n = 0,
                                   new_test_n = intervention_measures["smear_exams_n.I_n"],
                                   new_test_non_TBcases = intervention_measures["smear_exams_non_TBcases.TB_prevalence.I_n"]
        ) 
      }
      

      ## Complex intervention measurs
      intervention_measures <- c(intervention_measures,
        total_new_diag = intervention_measures["total_new_diag_n.I_n"] + intervention_measures["total_new_diag_p.I_p"],
        total_smear_exams = intervention_measures["smear_exams_n.I_n"] + intervention_measures["smear_exams_p.I_p"] + intervention_measures["smear_exams_non_TBcases.TB_prevalence.I_n"],
        total_Xray_exams = intervention_measures["Xray_exams_n.smear_exams_n.I_n"] + intervention_measures["Xray_exams_non_TBcases.smear_exams_non_TBcases.TB_prevalence.I_n"],
        total_new_diag_tests = intervention_measures["new_test_n.smear_exams_n.I_n"] + intervention_measures["new_test_non_TBcases.smear_exams_non_TBcases.TB_prevalence.I_n"]
      )
      
    ## output
    list(derivatives, c(summary_measures, intervention_measures))
  })
}

## Set parameters for model
params_TB_model_diag <- function(ecr_pyr = 15, wks_health_service = 52,
                            prot_reinf = 0.65,
                            intervention = 1, intervention_start = 2014,
                            timestep = 52) {
  
  ## yearly rate of disease onset
  onset_yr <- 0.03

  ## proportion that are smear negative
  prop_p <- 0.7
  
  ## Mortality rates
  mu <- 0.021
  mu_n <- 0.19
  mu_p <- 0.22
  
  ## Cure rate
  natural_cure <- 0.2
  
  ## TB model parameters
  params <- list(ecr = ecr_pyr / timestep,
                 rate_low_latent = 1 / (5 * timestep),
                 rel_inf_n = 0.22,
                 prim_dis_onset = onset_yr / timestep, 
                 sec_dis_onset  = 0.0003/timestep,
                 reinf_dis_onset = (onset_yr * (1 - prot_reinf))/timestep,
                 prop_n = 1 - prop_p, 
                 prop_p = prop_p,
                 mu = mu/timestep,
                 mu_n = mu_n/timestep,
                 mu_p = mu_p/timestep,
                 natural_cure = natural_cure/timestep,
                 relapse = 0.001/timestep,
                 timestep = timestep
  )
  
  ## Treatment pathway parameters - weeks to treatmetns
  wks_to_health_services_n <- wks_health_service
  wks_to_health_services_p <- wks_health_service
  wks_HS_visit_to_treatment <- 5.4
  
  wks_to_treat_if_Xrayed_n <- wks_to_health_services_n + wks_HS_visit_to_treatment
  wks_to_treat_p <- wks_to_health_services_p + wks_HS_visit_to_treatment
  
  ## Xray related parameters
  prop_access_Xray <- 0.3
  X_ray_sensitivity <- 0.8
  prop_default_treatment <- 0.15
  prop_treated_successfully <- 0.75
  
  ## New Intervention adjusted weeks
  new_diag_sensitivity_n <- 0.7
  new_diag_wks_to_treat_n <- wks_to_treat_if_Xrayed_n
  
  ## Rates of detection and recovery for both postive and negative sputum smear
  intervention_rates <- list(
    prop_access_Xray = prop_access_Xray,
    preinterv_detect_recover_n = (52 * prop_access_Xray*X_ray_sensitivity*(1 - prop_default_treatment)*prop_treated_successfully) / (wks_to_treat_if_Xrayed_n * timestep),
    detect_recover_p = 52 * (1 - prop_default_treatment)*prop_treated_successfully / (wks_to_treat_p * timestep),
    preinterv_detect_only_n =  52 * prop_access_Xray*X_ray_sensitivity / (wks_to_treat_if_Xrayed_n * timestep),
    detect_only_p = 52 / (wks_to_treat_p * timestep),
    interv_detect_recover_n = 52 * new_diag_sensitivity_n * (1 - prop_default_treatment) * prop_treated_successfully / (new_diag_wks_to_treat_n * timestep),
    interv_detect_only_n = 52 * new_diag_sensitivity_n / (new_diag_wks_to_treat_n * timestep)
    )

  ## Intervention control parameters
  yrs_evaluate <- 10
  
  intervention_params <- list(
    intervention = intervention,
    intervention_start = intervention_start,
    yrs_evaulate = yrs_evaluate,
    yr_eval_impact = intervention_start + yrs_evaluate  
  )
  
  ## Cost parameters
  
  rate_ratio_attendance_non_TB_vs_TB <-  0.01
  
  costs_params <- list(
    rate_ratio_attendance_non_TB_vs_TB = rate_ratio_attendance_non_TB_vs_TB,
    test_TBcases = 52 / (wks_to_treat_p * timestep),
    test_non_TBcases = 52 * rate_ratio_attendance_non_TB_vs_TB / (timestep * wks_to_treat_p),
    new_diagnostic_cost = 20,                                                                                                                          
    smear_cost = 5,                                                                                                                                       
    X_ray_cost = 10 
  )
  return(c(params, intervention_rates, intervention_params, costs_params))
}
                                                                                                                                              
## Cumulative measures                                                                                                                                                
cum_measures <- function(model_traj, params) {
  model_traj <- as.data.frame(model_traj) %>% as_tibble
  
  ## Filter based on start time
  model_traj <- model_traj[unlist(model_traj[["year"]]) >= params["intervention_start"] & 
                             unlist(model_traj[["year"]]) <= params["yr_eval_impact"], ]
  
  with(model_traj,{
    df <- data_frame(Year = year,
                     Cases = ifelse(params["intervention"] == 1, total_new_cases.H, 0),
                     Infections = ifelse(params["intervention"] == 1, total_new_infs.I_n, 0),
                     Deaths = ifelse(params["intervention"] == 1, total_new_deaths.I_n, 0),
                     Diagnoses = ifelse(params["intervention"] == 1, total_new_diag.total_new_diag_n.I_n, 0)
    ) %>% 
      mutate(cum_cases = cumsum(Cases),
             cum_inf = cumsum(Infections),
             cum_deaths = cumsum(Deaths),
             cum_diag = cumsum(Diagnoses)) %>% 
      mutate(Cases = round(cum_cases, digits = 0),
             Infections = round(cum_inf, digits = 0),
             Deaths = round(cum_deaths, digits = 0),
             Diagnoses = round(cum_diag, digits = 0)) %>% 
      select(-cum_cases, -cum_inf, -cum_deaths, -cum_diag)
    
    df
  })

}                                                                                                                                              
 
## Diagnosis summary
sum_diag <- function(model_traj, start_time = 2010) {
  model_traj <- as.data.frame(model_traj) %>% as_tibble
  
  ## Filter based on start time
  model_traj <- model_traj[unlist(model_traj[["year"]]) >= start_time, ]
  
  with(model_traj,{
    df <- data_frame(Year = year,
                     `Sputum negative` = total_new_diag_n.I_n,
                     `Sputum positive` = total_new_diag_p.I_p,
                     `TB cases` = total_new_diag.total_new_diag_n.I_n) %>% 
      mutate(`Sputum negative` = round(`Sputum negative`, digits = 0),
             `Sputum positive` = round(`Sputum positive`, digits = 0),
             `TB cases` = round(`TB cases`, digits = 0))
    
    df
    
  })
}
## Test                                                                                                                                              

TB Model with Diagnostics Graphs

Download

## Graphs Implemented
## Graph 1: Incidence in pops
## Graphs to make
## Graph 2: cum cases since the start of the intervention
## Graph 3: Total new diagnosis
## Graph 4: Number of smear exams and xrays in sputum negative populations
## Graph 5: total cumulative costs 
## Graph 6: cost per_diagnosed case

## Thoughts about graphs
## Graph 2, 5, and 6 only need data from the start of the intervention
## 

graph_cum <- function(df){
  plot <- df %>%
    gather(key = "Measure", value = "Cumulative Count",
           Cases, Infections, Deaths, Diagnoses) %>% 
    ggplot(aes(x = Year, y = `Cumulative Count`,
               colour = Measure, group = Measure)) +
    geom_line() +
    theme_minimal() +
    theme(legend.position = "bottom") +
    expand_limits(y = 0)
  
  return(plot)  
}

graph_diag <- function(df) {
  plot <- df %>%
    gather(key = "Disease", value = "New Diagnoses",
           `Sputum negative`, `Sputum positive`, `TB cases`) %>% 
    ggplot(aes(x = Year, y = `New Diagnoses`,
               colour = Disease, group = Disease)) +
    geom_line() +
    theme_minimal() +
    theme(legend.position = "bottom") +
    expand_limits(y = 0)
  
  return(plot)  
}