attr_names <- EpiModelWHAMPDX::attr_names
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(magrittr))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(ggrepel))
library(EpiModelWHAMPDX)
sim_dat <- readRDS("../../EpiModel/AE/sim_epimodel3/sim_at_2033_feb16b.rds")
options(dplyr.summarise.inform = FALSE)

print_many_plots <- function(testing_plots, sim_dat, num_hash = 2){
  sml_temp <- paste0(c(rep("#", num_hash), " %s ", "
", " 
 "), collapse = "")
  sub_plot_info <- testing_plots
  if (is.null(sub_plot_info$plt_type)) {
    plt_type <- "line" }else{ 
      plt_type <- sub_plot_info$plt_type }
  meas_data <- make_epi_plot_data(sim_dat, sub_plot_info)
  for (ct_idx in seq(attr_names)) {
    this_cat <- names(attr_names)[ct_idx]
    targs <- EpiModelWHAMPDX::WHAMP.targs %>% 
      filter(measure == sub_plot_info$name)
    if (nrow(targs) == 0) {
      targs <- NULL
    }else{
      sub_trg <- targs %>% filter(cat == this_cat)
      if (nrow(sub_trg) == 0) {
        catgs <- EpiModelWHAMPDX::attr_names[[this_cat]]
        sub_trg <- do.call(rbind, rep(list(targs), length(catgs)))
        sub_trg$sub_cat <- catgs
        sub_trg$cat <- this_cat
      }
      targs <- sub_trg
    }
    if (this_cat == "ovr") {
      cat_meas_name <- "overall" 
    }else if (this_cat == "age.grp") {
      cat_meas_name <- "age"
    }else{
      cat_meas_name <- this_cat
    }
    cat(sprintf(sml_temp, cat_meas_name))
    print(EpiModelWHAMPDX::plot_epi(meas_data, this_cat,
                                    plot_type = plt_type, 
                                    targets = targs,
                                    year_range = c(1990, 2030)))
    cat("\n\n")
  }
}

1 Introduction to measures

There are three different measures studied in this document.

Proportion of DX HIV+ on ART

  • This measure is the fraction:

\[ \frac{\text{Number of individuals who have HIV who are currently on ART}}{ \text{Number of individuals who have been diagnosed with HIV }} \]

  • The attribute (in EpiModel) that tracks if an individual has been diagnosed is diag.status == 1. ART status is determined by tx.status.

Proportion of those on ART Suppressed

  • This measure is the fraction:

\[ \frac{\text{Number of individuals currently on ART who are virally suppressed}}{ \text{Number of individuals currently on ART }} \]

  • Note that individuals in the simulation alternate on and off of ART treatment, and are only included in this calculation when they are on ART.
    • Before starting ART : not counted
    • Directly after starting ART : counted
    • After stopping ART : not counted (even though they may be virally suppressed for a short time after discontinuation)
    • After starting ART again : counted
  • The attribute (in EpiModel) that determines ART status is tx.status and the attribute that determines viral suppression in vl == 1.5.

Mean time on treatment

  • This mean is broken down across all individuals on ART, and all individuals on ART who are partial suppressors.

The construction of this mean is somewhat counterintuitive:

  • Before an individual is placed on treatment, their time on and off treatment are both zero. One exception to this is that during initialization HIV positive individuals are a time off treatment equal to their (randomly assigned) time infected. Additionally, all HIV positive individuals are assigned a time on treatment of zero.
  • After first starting treatment, at each time step, either time on treatment or time off treatment is increased by one depending on if the individual is on treatment at the current time step.

The mean time on and off treatment is just the mean value of these times among the specified subgroup.

2 Proportion of DX HIV+ on ART

2.1 Construction rule

Open the first code chunk below to see the construction rule

# prop.of.hiv.on.art <- tx.i.num / diag.i.num =
#   (status == 1 & tx.status == 1) / (status == 1 & diag.status == 1) = 
# (tx.status == 1) / (diag.status == 1)
# Note that tx.status == 1 requires that diag.status == 1 and both
# tx.status and diag.status cannot be 1 without status == 1
prop.art <- list(name = "prop.art",
                 plot_name = "Proportion of DX HIV+ on ART",
                 plt_type = "line",
                 vars = c("tx.i.num.", "diag.i.num."),
                 sum_fun = function(x, y) x / y)

print_many_plots(prop.art, sim_dat)

2.2 overall

2.3 race

2.4 region

2.5 age

3 Proportion of those on ART Suppressed

3.1 Construction rule

Open the first code chunk below to see the construction rule

# prop.of.hiv.supr <- supr.on.tx. / tx.i.num =
#   (vl == 1.5 & tx.status == 1)  / (tx.status == 1)

# Note that both tx.status cannot be 1 without status == 1.
prop.supr = list(name = "prop.supr",
                 plot_name = "Proportion of HIV+ Suppressed", 
                 plt_type = "line", 
                 vars = c("supr.on.tx.", "tx.i.num."),
                 sum_fun = function(x, y) x/y)
print_many_plots(prop.supr, sim_dat)

3.2 overall

3.3 race

3.4 region

3.5 age

4 Mean time on treatment

4.1 All on ART

4.1.1 On treatment

time.on.treat <- list(name = "mean.tx.on",
                      plot_name = "Mean Time on Treatment",
                      vars = "mean.tx.on.",
                      sum_fun = function(x) x / 4,
                      plt_type = "line")

print_many_plots(time.on.treat, sim_dat, num_hash = 4)

4.1.1.1 overall

4.1.1.2 race

4.1.1.3 region

4.1.1.4 age

4.1.2 Off treatment

time.off.treat <- list(name = "mean.tx.off",
                      plot_name = "Mean Time off Treatment",
                      vars = "mean.tx.off.",
                      sum_fun =  function(x) x / 4,
                      plt_type = "line")

print_many_plots(time.off.treat, sim_dat, num_hash = 4)

4.1.2.1 overall

4.1.2.2 race

4.1.2.3 region

4.1.2.4 age

4.2 Partial suppressors (tt.traj == 1)

4.2.1 On treatment

time.on.treat.part <- list(
  name = "time.on.part",
  plot_name = "Mean Time on Treatment in Months (Partial Supressors)",
  plt_type = "line",
  vars = "mean.tx.on.part.",
  sum_fun = function(x) x / 4
  )

print_many_plots(time.on.treat.part, sim_dat, num_hash = 4)

4.2.1.1 overall

4.2.1.2 race

4.2.1.3 region

4.2.1.4 age

4.2.2 Off treatment

time.off.treat.part <- list(
  name = "mean.tx.off.part",
  plot_name = "Mean Time off Treatment (Partial Suppressors)",
  vars = "mean.tx.off.part.",
  sum_fun = function(x) x / 4,
  plt_type = "line")

print_many_plots(time.off.treat.part, sim_dat, num_hash = 4)

4.2.2.1 overall

4.2.2.2 race

4.2.2.3 region

4.2.2.4 age

5 Proportion ever on treatment after N weeks

dx_df <- data.frame(
  status = sim_dat$attr[[1]]$cuml.time.on.tx > 0,
  time = sim_dat$control$nsteps - sim_dat$attr[[1]]$diag.time
)

dx_df %>% group_by(time) %>% summarise(txed = mean(status)) %>%
    ggplot(aes(time, txed)) + geom_point() + ggtitle("Proportion Treated X time steps (weeks) ago") + geom_vline(xintercept = 8) + 
  scale_x_log10()
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 1 rows containing missing values (geom_point).

# prp_late_tester <- list(name = "prop.late.tester",
#       plot_name = "Proportion who are late testers", 
#       plt_type = "line", 
#       vars = c("num.late.testers.", "num."),
#       sum_fun = function(x, y) x / y
#  )
# 
# prp_late_tester_arrive <- list(name = "prop.late.tester.arriv",
#       plot_name = "Proportion who enter who are late testers", 
#       plt_type = "line", 
#       vars = c("num.arrived.", "new.late.testers."),
#       sum_fun = function(x, y) {
#         data.table::frollsum(y[!is.na(y)], 200) /
#         data.table::frollsum(x[!is.na(x)], 200)}
#  )
#  
# late_test_rate <- list(name = "late.test.rate",
#       plot_name = "Late testing rate", 
#       plt_type = "line", 
#       vars = c("num.late.testers.", "num.get.test.late."),
#       sum_fun = function(x, y, z) data.table::frollmean(y / (x),
#                                                      200, na.rm = TRUE)
#  )
#  
# prp_late_test_dat <- make_plot_data(dat, prp_late_tester)
#  
#  make_cmpr_plot(prp_late_test_dat, year_range = c(1940, 2030), 
#                 targets = plt_targ)
#  late_test_rate_dat <- make_plot_data(dat, late_test_rate)
#  make_cmpr_plot(late_test_rate_dat, year_range = c(1940, 2030), 
#                 targets = ltr_targ)
# 
#  prp_late_arv_test_dat <- make_plot_data(dat, prp_late_tester_arrive)
#  
#  make_cmpr_plot(prp_late_arv_test_dat, year_range = c(1940, 2000))
#  
#  sub_plot_info <- testing_plots[[ms_idx]]
#      if (is.null(sub_plot_info$plt_type)) {
#        plt_type <- "line" }else{ 
#          plt_type <- sub_plot_info$plt_type }
#      meas_data <- make_plot_data(sim_dat, sub_plot_info)
#      this_meas_name <- sub_plot_info$plot_name
#      cat(sprintf(big_temp, this_meas_name))
#      for (ct_idx in seq(attr_names)) {
#        this_cat <- names(attr_names)[ct_idx]
#        cat(sprintf(sml_temp, this_cat))
#        print(make_cmpr_plot(meas_data, this_cat,
#                             plot_type = plt_type))
#        cat("\n\n")
# 
#   big_temp <- "# %s {.tabset .tabset-fade .tabset-pills}
# 
# "
#   sml_temp <- "## %s
# 
# "