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