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