rm(list = ls())
library(dplyr)
library(data.table)
library(readxl)
library(tidyverse)
library(knitr)
library(kableExtra)
library(reshape2)
library(ggplot2)
library(plotly)
library(gt)
library(tiff)
library(grid)

1 Introduction

This report is based mostly on WADOH DAP participation data, which are used to construct parameters and targets for program participation dynamics and viral load outcomes. There are also some parameters pulled from Darcy’s analysis of (MMP data?) that inform treatment trajectories and dynamics (see the “Other parameters” section).

2 Parameters

#### Reading data from updated WApopdata repository
msm.pop.totals_2014 <- WApopdata::msm.pop.totals_2014
msm.pop.totals_2015 <- WApopdata::msm.pop.totals_2015
msm.pop.totals_2016 <- WApopdata::msm.pop.totals_2016
msm.pop.totals_2017 <- WApopdata::msm.pop.totals_2017
msm.pop.totals_2018 <- WApopdata::msm.pop.totals_2018
msm.pop.totals_2019 <- WApopdata::msm.pop.totals_2019
pop.totals <- WApopdata::pop.totals # old object

## Reading DAP data from Zoe, fixing the region names in the files first if necessary
#source(here::here("MakeData", "MM", "update_dap_data.R"))

adap_data <- readRDS(here::here("Data", "Clean", "DAPdata",
                                "adap_data.rds"))
pdap_data <- readRDS(here::here("Data", "Clean", "DAPdata",
                                "pdap_data.rds"))

#### Reading WHAMP survey data objects
whamp <- readRDS(here::here("MakeData", "MM", "Old",
                            "ZKWhampDemog.rds"))

2.1 Description of the data objects

These data objects are used for creating the final whamp_data object that is used to set parameters and targets in EpiModelHIV for WHAMP.

2.1.1 Population estimates

The population estimates msm.pop.totals_20XX come from the WApopdata repository, for XX=year (14 thru 19). The objects are documented, so the details are available via ?WApopdata::msm.pop.totals_2019.

These are list objects, which include the MSM population estimates by demographic attributes (age, race, region) and HIV status. Summary of the 2019 object below.

summary(msm.pop.totals_2019)
##                     Length Class      Mode
## pop.all             5      tbl_df     list
## pop.pos             8      tbl_df     list
## pop.neg             6      tbl_df     list
## pop.age.all         5      tbl_df     list
## pop.age.neg         4      tbl_df     list
## pop.age.pos         9      tbl_df     list
## pop.race.all        5      tbl_df     list
## pop.race.neg        4      tbl_df     list
## pop.race.pos        9      tbl_df     list
## pop.region.all      3      tbl_df     list
## pop.region.neg      4      tbl_df     list
## pop.region.pos      6      tbl_df     list
## pop.racexregion.all 4      grouped_df list
## pop.racexregion.neg 5      grouped_df list
## pop.racexregion.pos 7      grouped_df list
## MSM_PLWH_WA         5      tbl_df     list
## msm.wts.all.age5    5      tbl_df     list
## msm.wts.all.age10   5      tbl_df     list
## msm.wts.15_65.age10 5      tbl_df     list

2.1.2 DAP data files

These files are derived from the PDAP program data. Construction workflow includes DataManagement.R and DAPoutputs.R.

2.1.2.1 ADAP data

adap_data is a list including the following elements.

kable(adap_data$des_table, 
      caption= "ADAP data") %>% 
  kable_styling(full_width = F, position = "center", bootstrap_options = c("striped"))
ADAP data
Parameters Value Description Calculation Source group
surv_coef NA coefficients of the survivial model Using exponential and Weibull survival regression ADAP claims data ADAP-related
pred_disenroll NA predicted probability of disenrollment of ADAP Calculated based on surv_coef ADAP claims data ADAP-related
adap_demog NA demographics of ADAP clients by year, race, and region ADAP claims data ADAP-related
adap_ins NA insurance type of ADAP clients by year, race, and region ADAP claims data ADAP-related
adap_new_clients NA demographics of new clients of ADAP clients by year, race and region ADAP claims data ADAP-related
ADAPcost_pd NA ADAP daily cost by insurance type by insurance type and race ADAP claims data ADAP-related
ADAPcost_tot NA ADAP total cost by insurance type by insurance type and race ADAP claims data ADAP-related
vl_race NA proportion of population with viral suppression by race and ADAP status ADAP claims data and viral load data ADAP-related

2.1.2.2 PrEP-DAP data

pdap_data is a list including the following elements.

kable(pdap_data$des_table, 
      caption= "PrEP-DAP data") %>% 
  kable_styling(full_width=F, position="center", bootstrap_options = c("striped"))
PrEP-DAP data
Parameters Value Description Calculation Source group
pdap_demog NA demographics of PDAP clients by year, race and region PDAP claims data PDAP-related
pdap_ins NA insurance type of PDAP clients by year, race and region PDAP claims data PDAP-related
pdap_new_clients NA demographics of new clients of PDAP clients by year, race, region PDAP claims data PDAP-related
PDAPcost_pd NA PDAP daily cost by insurance type by race, and insurance type PDAP claims data PDAP-related
init_pdap_dat NA initial enrollment days of PDAP clients by race, insurance, region PDAP claims data PDAP-related
drop_prob_dat NA probability of disenrollment by race, insurance, and years on PDAP PDAP claims data PDAP-related

2.1.3 WHAMP cost statistics

This is a list object, with elements derived from the WHAMP survey. See the make_WhampCostStats.Rmd file for details (peport also up online here)

whamp$des_table %>%
  gt::gt(rowname_col = "Outputs",
     groupname_col = "group")
Description Source
Parameters: Demog
pred_inc predicted income multinomial regression on age, race, region
income_levels income FPL levels FPL guidelines
pred_ins predicted insurance multinomial regression on age, race, region, hiv status
Parameters: PrEP
pred_prep_aware predicted PrEP awareness logistic regression on race, region
pred_prep_ever predicted PrEP ever (interest) logistic regression on age, region
Targets
target_prep_current target HIV- on PrEP obs props overall and by age, race and region
target_pdap target PrEP users in PDAP obs props overall and by age, race and region
target_adap target ART patients in ADAP obs props overall and by race and region

2.3 PrEP-DAP parameters

2.3.1 pdap_demog

  • Demographics of PrEP-DAP clients.
  • Calculated the number of PrEPDAP clients by year, race, and region
  • Calculated the number of HIV-negative individuals by race and region using Darcy’s estimation and the total population in the DSdata. The object created is called pop.racexregion.neg. Reference scripts could be found in this section.
  • Calculated the proportion of HIV-negative who are on PrEP-DAP by race and region.
### PDAP
### we again use num_ub, which for HIV- will be a the lower value

pop.racexregion.neg <- bind_rows(msm.pop.totals_2014$pop.racexregion.neg,
                                 msm.pop.totals_2015$pop.racexregion.neg,
                                 msm.pop.totals_2016$pop.racexregion.neg,
                                 msm.pop.totals_2017$pop.racexregion.neg,
                                 msm.pop.totals_2018$pop.racexregion.neg) %>%
  bind_cols(year, .) %>%
  select(year, region, race, popsize.neg = num_ub) %>%
  data.table()

pdap_data$pdap_demog <- merge(pdap_data$pdap_demog, 
                              pop.racexregion.neg, 
                              by = c("year", "race", "region"),
                              all.x = T)

pdap_data$pdap_demog[, `:=` (prop = N / popsize.neg, 
                             year = year - 2013)]
pdap_data$pdap_demog <- pdap_data$pdap_demog[year >= 4]
setkeyv(pdap_data$pdap_demog, c("year", "race", "region"))

2.3.1.1 Plot

p <- ggplot(pdap_data$pdap_demog, 
            aes(x = region, y = prop, text = N)) + 
  geom_bar(aes(fill = race), alpha = 0.7,
           stat = "identity", size = 1.2, 
           position = position_dodge2(width = 0.9)) + 
  labs(title = "Demographics of PrEP-DAP clients", 
       ylab = "proportion") + 
  facet_wrap(~year)

ggplotly(p, tooltip = "text")

2.3.1.2 Table

kable(pdap_data$pdap_demog, 
      digits = c(0, 0, 0, 0, 0, 4), 
      caption= "Demographics of PrEP-DAP clients by year, race, and region") %>% 
  kable_styling(full_width=F, position="center", bootstrap_options = c("striped"))
Demographics of PrEP-DAP clients by year, race, and region
year race region N popsize.neg prop
4 B EasternWA 0 252 0.0000
4 B King 9 4658 0.0019
4 B WesternWA 7 1980 0.0035
4 H EasternWA 3 2191 0.0014
4 H King 60 5478 0.0110
4 H WesternWA 12 3390 0.0035
4 O EasternWA 20 9676 0.0021
4 O King 429 56734 0.0076
4 O WesternWA 81 35822 0.0023
5 B EasternWA 0 256 0.0000
5 B King 43 4813 0.0089
5 B WesternWA 13 2061 0.0063
5 H EasternWA 18 2249 0.0080
5 H King 276 5620 0.0491
5 H WesternWA 66 3556 0.0186
5 O EasternWA 50 9765 0.0051
5 O King 1005 57656 0.0174
5 O WesternWA 195 36230 0.0054

2.3.2 pdap_ins

  • Distribution of insurance type among PrEP-DAP clients

  • We calculated total number of PrEP-DAP clients by race, region, and insurance type in different year.

  • The insurance type includes individual, employer, and none.

    • In the original data, there were only 17 Medicare records so we removed the Medicare insurance in the data.
    • There were no Medicaid clients in the PrEP-DAP claims data

  • We only kept the data in 2017 and 2018.

  • Calculated the distribution of insurance type by race and region in each year.

2.3.2.1 Table

pdap_data$pdap_ins[, sum_N := sum(n), by = list(year, race, region)]
pdap_data$pdap_ins[, prop := n / sum_N]
pdap_data$pdap_ins[, sum_N := NULL]
pdap_data$pdap_ins[, year := year - 2013]
pdap_data$pdap_ins <- pdap_data$pdap_ins[year >= 4]

pdap_ins1 <- data.table(expand.grid(year = c(4, 5), 
                                    race = c("B", "H", "O"), 
                                    region = c("EasternWA", "King", "WesternWA"), 
                                    insurance = c("employer", "individual", "none")))

pdap_data$pdap_ins <- merge(pdap_ins1, pdap_data$pdap_ins, 
                            by = c("year", "race", "region", "insurance"), 
                            all.x = T)
pdap_data$pdap_ins[is.na(n)]$n <- 0
pdap_data$pdap_ins[is.na(prop)]$prop <- 0

kable(pdap_data$pdap_ins, digits = c(0, 0, 0, 0, 0, 3), 
      caption= "Insurance type of PrEP-DAP clients") %>% 
  kable_styling(full_width=F, position="center", bootstrap_options = c("striped"))
Insurance type of PrEP-DAP clients
year race region insurance n prop
4 B EasternWA employer 0 0.000
4 B EasternWA individual 0 0.000
4 B EasternWA none 0 0.000
4 B King employer 7 0.778
4 B King individual 1 0.111
4 B King none 1 0.111
4 B WesternWA employer 6 0.857
4 B WesternWA individual 1 0.143
4 B WesternWA none 0 0.000
4 H EasternWA employer 3 1.000
4 H EasternWA individual 0 0.000
4 H EasternWA none 0 0.000
4 H King employer 36 0.600
4 H King individual 16 0.267
4 H King none 8 0.133
4 H WesternWA employer 7 0.583
4 H WesternWA individual 4 0.333
4 H WesternWA none 1 0.083
4 O EasternWA employer 8 0.400
4 O EasternWA individual 10 0.500
4 O EasternWA none 2 0.100
4 O King employer 261 0.608
4 O King individual 127 0.296
4 O King none 41 0.096
4 O WesternWA employer 42 0.519
4 O WesternWA individual 23 0.284
4 O WesternWA none 16 0.198
5 B EasternWA employer 0 0.000
5 B EasternWA individual 0 0.000
5 B EasternWA none 0 0.000
5 B King employer 23 0.535
5 B King individual 4 0.093
5 B King none 16 0.372
5 B WesternWA employer 8 0.615
5 B WesternWA individual 3 0.231
5 B WesternWA none 2 0.154
5 H EasternWA employer 10 0.556
5 H EasternWA individual 2 0.111
5 H EasternWA none 6 0.333
5 H King employer 141 0.511
5 H King individual 24 0.087
5 H King none 111 0.402
5 H WesternWA employer 28 0.424
5 H WesternWA individual 6 0.091
5 H WesternWA none 32 0.485
5 O EasternWA employer 30 0.600
5 O EasternWA individual 12 0.240
5 O EasternWA none 8 0.160
5 O King employer 662 0.659
5 O King individual 196 0.195
5 O King none 147 0.146
5 O WesternWA employer 113 0.579
5 O WesternWA individual 31 0.159
5 O WesternWA none 51 0.262

2.3.2.2 Plot

p <- ggplot(pdap_data$pdap_ins %>% mutate(year = year+2013), 
       aes(x = factor(year), y = prop, text = prop)) + 
  geom_line(aes(group = race, color = race)) + 
  facet_grid(region ~ insurance) + 
  labs(title = "Insurance type for PDAP clients") +
  xlab("year") + ylab("proportion") +
  theme(axis.text.x = element_text(size = 9, angle = 90, hjust = 1))


# This plot takes a lot of tweaking to display properly using
# ggplotly

gp <- ggplotly(p, tooltip = "text")

#str(gp[['x']][['layout']][['annotations']]) 

layout_ggplotly <- function(
  gg, x = -0.1, y = -0.05, 
  x_legend=1.05, y_legend=0.95, 
  mar=list(l=50, r=150)){
  # The 1 and 2 goes into the list that contains the options for the x and y axis labels respectively. last component is the legend
  gg[['x']][['layout']][['annotations']][[1]][['y']] <- x
  gg[['x']][['layout']][['annotations']][[2]][['x']] <- y
  gg[['x']][['layout']][['annotations']][[9]][['x']] <- x_legend
  gg[['x']][['layout']][['legend']][['y']] <- y_legend
  gg[['x']][['layout']][['legend']][['x']] <- x_legend
  gg %>% layout(margin = mar)
}

gp %>% layout_ggplotly

2.3.3 PDAPcost_pd

  • The calculation method is very similar to the ADAPcost_pd.
  • The primary difference is that we focus on years 2017 and 2018 and insurance types: none, employer, and individual.
  • Calculated the number of clients and the total number of days on PrEP-DAP among all the clients in each race and insurance category.
  • Caculated the total cost of a cost category among all PrEP-DAP active clients in each race and insurance category.
  • The daily cost for a cost category was its total cost divided by the number of PrEP-DAP active days among all clients in the race and insurance category.
  • The daily costs were calculated for 2017 and 2018 respectively.
tmpDF <- data.table(expand.grid(year = c(2017, 2018), 
                                race = c("B", "H", "O"), 
                                insurance = c("none", "employer", "individual")))
pdap_data$PDAPcost_pd <- pdap_data$PDAPcost_pd %>%
  select(year, race, insurance, pdtruvCost, Overallpdcost) %>%
  rename(cTruv = pdtruvCost, 
         cTotCost = Overallpdcost)

pdap_data$PDAPcost_pd <- merge(tmpDF, pdap_data$PDAPcost_pd, 
                               by = c("year", "race", "insurance"), 
                               all.x = T)
pdap_data$PDAPcost_pd[is.na(pdap_data$PDAPcost_pd)] <- 0

2.3.3.1 Plots

2.3.3.1.1 Truvada
p <- ggplot(pdap_data$PDAPcost_pd, 
            aes(x = insurance, y = cTruv, text = round(cTruv,2))) + 
  geom_bar(aes(fill = race), stat = "identity", 
           alpha = 0.7, size = 1.2, 
           position = position_dodge2(width = 0.9)) +
  facet_wrap(~ year, nrow = 2) + 
  labs(title = "Truvada cost per day", y = "cost ($)") 

ggplotly(p, tooltip = "text")
2.3.3.1.2 Total costs
p <- ggplot(pdap_data$PDAPcost_pd, 
            aes(x = insurance, y = cTotCost, text = round(cTotCost,2))) + 
  geom_bar(aes(fill = race), stat = "identity", 
           alpha = 0.7, size = 1.2, 
           position = position_dodge2(width = 0.9)) +
  facet_wrap(~ year, nrow = 2) + 
  labs(title = "Total cost per day", y = "cost ($)") 

ggplotly(p, tooltip = "text")

2.3.3.2 Table

kable(pdap_data$PDAPcost_pd, digits = c(0, 0, 0, 2, 2, 2, 2), 
      caption = "Daily PrEP-DAP cost by race and insurance") %>% 
  kable_styling(full_width=F, position="center", bootstrap_options = c("striped"))
Daily PrEP-DAP cost by race and insurance
year race insurance cTruv cTotCost
2017 B employer 3.00 3.10
2017 B individual 0.00 0.00
2017 B none 0.09 0.09
2017 H employer 5.47 5.48
2017 H individual 11.52 11.52
2017 H none 25.30 25.49
2017 O employer 5.30 5.34
2017 O individual 9.92 9.93
2017 O none 10.59 10.65
2018 B employer 1.15 1.63
2018 B individual 2.02 2.02
2018 B none 4.84 6.48
2018 H employer 1.15 1.45
2018 H individual 6.43 6.66
2018 H none 3.35 5.90
2018 O employer 1.82 2.19
2018 O individual 5.46 5.62
2018 O none 8.56 10.14

2.3.4 pred_disenroll_pdap

  • Disenrollment from ADAP

  • We fitted a logistic regression to predict the probability of PrEP-DAP disenrollment at the end of year on PrEP-DAP.

  • The model controls for the years on PrEP-DAP (year 1, 2 or 3), race, and insurance type (non, individual, and employer).

2.3.4.1 Model

summary(pdap_data$dropout_model)
## 
## Call:
## glm(formula = event ~ factor(years) + RaceCat + Insurance, family = binomial(link = "logit"), 
##     data = PDAPexpand %>% filter(cateYear >= 2016 & Insurance != 
##         "Public"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0268  -1.2224   0.7484   1.0368   1.4471  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.58545    0.23153   2.529   0.0115 *  
## factor(years)2       0.23476    0.08720   2.692   0.0071 ** 
## factor(years)3+      0.77688    0.13538   5.739 9.55e-09 ***
## RaceCatHispanic      0.01032    0.24848   0.042   0.9669    
## RaceCatOther        -0.48017    0.23284  -2.062   0.0392 *  
## InsuranceIndividual -0.72000    0.09900  -7.273 3.52e-13 ***
## InsuranceNone        0.54403    0.10480   5.191 2.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3901.7  on 2864  degrees of freedom
## Residual deviance: 3728.0  on 2858  degrees of freedom
## AIC: 3742
## 
## Number of Fisher Scoring iterations: 4

2.3.4.2 Table

tmpDF <- pdap_data$drop_prob_dat

kable(tmpDF, digits = c(0, 0, 0, 3), 
      caption = "Probability of PDAP disenrollment at year end") %>% 
  kable_styling(full_width=F, position="center", bootstrap_options = c("striped"))
Probability of PDAP disenrollment at year end
race insurance years Prob
B employer 1 0.642
H employer 1 0.645
O employer 1 0.526
B individual 1 0.466
H individual 1 0.469
O individual 1 0.351
B none 1 0.756
H none 1 0.758
O none 1 0.657
B employer 2 0.694
H employer 2 0.696
O employer 2 0.584
B individual 2 0.525
H individual 2 0.528
O individual 2 0.406
B none 2 0.796
H none 2 0.798
O none 2 0.708
B employer 3+ 0.796
H employer 3+ 0.798
O employer 3+ 0.707
B individual 3+ 0.655
H individual 3+ 0.658
O individual 3+ 0.540
B none 3+ 0.871
H none 3+ 0.872
O none 3+ 0.806

2.3.4.3 Plot

p <- ggplot(tmpDF, 
            aes(x = years, y = Prob, 
                text = scales::percent(Prob,1))) + 
  geom_line(aes(group=race, color = race)) +
  facet_wrap(~ insurance) + 
  labs(title = "Probability of disenrollment from PrEP-DAP",
       y = "probability") 

ggplotly(p, tooltip = "text")

2.4 PrEP cascade

Predicted awareness and interest come from the WHAMP survey. Interest is proxied by the “ever use” variable. The model for awareness includes race and region; the model for interest includes age_group and region.

  • Note we don’t have data over time on either awareness or ever use/interest.
# whamp$pred_prep_ever must be by region only

# We only have one year of the WHAMP predicted probs,
# those are scaled to 2017 and 2018 popsizes for the PDAP
# program

prep_cascade <- merge(pdap_data$pdap_demog, 
                      whamp$pred_prep_aware, 
                      by = c("race", "region"), 
                      all.x = T) %>%
  mutate(n_aware = round(prob * popsize.neg)) %>%
  select(-c(prob)) %>%
  merge(., whamp$pred_prep_ever, by = "region", all.x = T) %>%
  mutate(n_interest = round(n_aware * prob)) %>% 
  select(-prob) %>%
  select(year, region, race, popsize.neg, n_aware, n_interest, n_pdap = N) %>%
  mutate(prop_pdap = n_pdap/n_interest) %>%
  data.table()

2.4.1 Plots

2.4.1.1 PrEP aware

tmpDF <- copy(prep_cascade)
tmpDF[, prop_aware := n_aware / popsize.neg]

ggplot(tmpDF[year == 5], aes(x = race, y = prop_aware)) + 
  geom_bar(aes(color = race), stat = "identity", 
           fill = "white", size = 1.2, 
           position = position_dodge2(width = 0.9)) +
  scale_color_brewer(palette = "Set1") + 
  facet_wrap(~ region) + 
  ylab("proportion") + 
  ggtitle("PrEP awareness (WHAMP Survey)") + 
  theme_bw() +
  theme(plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(size = 14), 
        axis.text.x = element_text(size = 10), 
        legend.position = "bottom")

2.4.1.2 PrEP interest

For PrEP interest (ever PrEP), ZK only considered the variation of region.

tmpDF[, prop_interest := n_interest / n_aware]

ggplot(tmpDF[year == 5 & race == "O"], aes(x = region, y = prop_interest)) + 
  geom_bar(color = "royalblue", stat = "identity", 
           fill = "white", size = 1.2, position = position_dodge2(width = 0.9)) +
  ylab("proportion") + 
  ggtitle("PrEP interest (ever_prep, WHAMP Survey)") + 
  theme_bw() +
  theme(plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(size = 14), 
        axis.text.x = element_text(size = 10), 
        legend.position = "bottom")

2.4.1.3 On PrEP-DAP

ggplot(tmpDF, aes(x = race, y = prop_pdap)) + 
  geom_bar(aes(x = race, y = prop_pdap, color = race, fill = factor(year)), 
           stat = "identity", size = 1.2, alpha = 0.6, 
           position = position_dodge2(width = 0.9)) +
  scale_color_brewer(palette = "Set1") + 
  scale_fill_brewer(palette = "Greys") + 
  facet_wrap(~ region) + 
  ylab("proportion") + 
  ggtitle("PDAP fraction of PrEP interested (ever_PrEP, WHAMP Survey)") + 
  theme_bw() +
  theme(plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(size = 14), 
        axis.text.x = element_text(size = 10), 
        legend.position = "bottom")

2.4.2 Table

kable(prep_cascade, digits = c(0, 0, 0, 0, 0, 0, 0, 3), 
      caption = "PrEP cascade") %>% 
  kable_styling(full_width=F, position="center", bootstrap_options = c("striped"))
PrEP cascade
year region race popsize.neg n_aware n_interest n_pdap prop_pdap
4 EasternWA B 252 195 41 0 0.000
5 EasternWA B 256 198 41 0 0.000
4 EasternWA H 2191 1241 259 3 0.012
5 EasternWA H 2249 1274 266 18 0.068
4 EasternWA O 9676 8353 1742 20 0.011
5 EasternWA O 9765 8430 1758 50 0.028
4 King B 4658 4491 1877 9 0.005
5 King B 4813 4641 1939 43 0.022
4 King H 5478 4991 2086 60 0.029
5 King H 5620 5121 2140 276 0.129
4 King O 56734 55613 23241 429 0.018
5 King O 57656 56517 23619 1005 0.043
4 WesternWA B 1980 1745 546 7 0.013
5 WesternWA B 2061 1816 569 13 0.023
4 WesternWA H 3390 2500 783 12 0.015
5 WesternWA H 3556 2623 821 66 0.080
4 WesternWA O 35822 33366 10447 81 0.008
5 WesternWA O 36230 33746 10566 195 0.018

2.5 WHAMP parameters

Please see the calculation of parameters here. (Note, this has been updated by MM) The parameters include the prediction of income distribution, insurance type, PrEP awareness, PrEP interest, PrEP current, PrEP-DAP, and ADAP active.

2.6 Other parameters

2.6.1 Treatment trajectory

  • The trajectory values were taken from Darcy’s report: statnet-WHAMP/DWR_dissertation/Parameter%20estimation/HIV_care_cascade.html
  • Because Sam’s PrevComb model considers 3 trajectory instead of 4, I assume that 50% of the full trajectory in Darcy’s estimates is durable suppression and the other 50% is really full suppression.
## Taken from "make_epi_data.R" file 
tt_traj_dt <- expand.grid(race = c("B", "H", "O"), 
                          region = c("EasternWA", "King", "WesternWA"), 
                          traj = c(1, 2, 3)) # 1: tt.part.supp; 2: tt.full.supp; 3: tt.dur.supp
tt_traj_dt <- data.table(tt_traj_dt, key = c("race", "region"))
tt_traj_dt[, prop := c((1 - 0.759), 0.759 * 0.5, 0.759 * 0.5, 
                       (1 - 0.839), 0.839 * 0.5, 0.839 * 0.5, 
                       (1 - 0.811), 0.811 * 0.5, 0.811 * 0.5, 
                       (1 - 0.847), 0.847 * 0.5, 0.847 * 0.5, 
                       (1 - 0.902), 0.902 * 0.5, 0.902 * 0.5, 
                       (1 - 0.885), 0.885 * 0.5, 0.885 * 0.5, 
                       (1 - 0.874), 0.874 * 0.5, 0.874 * 0.5, 
                       (1 - 0.920), 0.920 * 0.5, 0.920 * 0.5, 
                       (1 - 0.905), 0.905 * 0.5, 0.905 * 0.5)]
tt_traj_dt <- reshape(tt_traj_dt, idvar = c("race", "region"), 
                      timevar = "traj", direction = "wide")
setkeyv(tt_traj_dt, c("race", "region"))

kable(tt_traj_dt, digits = c(0, 0, 3, 3, 3), 
      caption= "Distribution of treatment trajectories (partial/full/durable)") %>% 
  kable_styling(full_width=F, position="center")
Distribution of treatment trajectories (partial/full/durable)
race region prop.1 prop.2 prop.3
B EasternWA 0.241 0.380 0.380
B King 0.161 0.420 0.420
B WesternWA 0.189 0.406 0.406
H EasternWA 0.153 0.424 0.424
H King 0.098 0.451 0.451
H WesternWA 0.115 0.442 0.442
O EasternWA 0.126 0.437 0.437
O King 0.080 0.460 0.460
O WesternWA 0.095 0.452 0.452

2.6.2 Time to ART treatment initiation

## Taken from "make_epi_data.R" file 
tx_init_dt <- data.table(expand.grid(race = c("B", "H", "O"), 
                                region = c("King", "WesternWA", "EasternWA")))
tx_init_dt[, days := c(46, 43, 47, 56, 53, 57, 53, 50, 53)]
setkeyv(tx_init_dt, c("race", "region"))

kable(tx_init_dt, digits = c(0, 0, 0), caption= "Median days to ART initiation") %>% kable_styling(full_width=F, position="center")
Median days to ART initiation
race region days
B King 46
B WesternWA 56
B EasternWA 53
H King 43
H WesternWA 53
H EasternWA 50
O King 47
O WesternWA 57
O EasternWA 53

2.6.3 ART reinitiation

This is also based on Darcy’s estimation here

tx_reinit_full_dt <- data.table(expand.grid(race = c("B", "H", "O"), 
                                     region = c("King", "WesternWA", "EasternWA")))
tx_reinit_full_dt[, prop := c(0.0244, 0.0239, 0.0275, 0.0200, 0.0195, 0.0223, 0.0219, 0.0209, 0.0237)]
setkeyv(tx_reinit_full_dt, c("race", "region"))

kable(tx_reinit_full_dt, digits = c(0, 0, 3), caption= "ART reinitiation") %>% kable_styling(full_width=F, position="center")
ART reinitiation
race region prop
B King 0.024
B WesternWA 0.020
B EasternWA 0.022
H King 0.024
H WesternWA 0.020
H EasternWA 0.021
O King 0.028
O WesternWA 0.022
O EasternWA 0.024

3 Targets

3.1 ADAP program size

adap.size <- adap_data$adap_demog %>%
  group_by(year) %>%
  summarize(adapN = sum(N),
            popsize = sum(popsize.pos)) %>%
  mutate(prop = adapN/popsize) %>%
  data.table()
  
setkeyv(adap.size, c("year")) # check if we need this

p <- ggplot(data = adap.size %>% mutate(year = year+2013), 
            aes(x = factor(year), y = prop, text = adapN)) +
  geom_bar(fill = "blue", alpha = 0.5,
           stat="identity") + 
  labs(title = "ADAP program size: Proportion of HIV+ population",
       x = "year", y = "proportion")

ggplotly(p, tooltip = "text")

3.2 PrEP-DAP program size

pdap.size <- pdap_data$pdap_demog %>%
  group_by(year) %>%
  summarize(pdapN = sum(N),
            popsize = sum(popsize.neg)) %>%
  mutate(prop = pdapN/popsize) %>%
  data.table()
  
setkeyv(pdap.size, c("year")) # check if we need this

p <- ggplot(data = pdap.size %>% mutate(year = year+2013), 
            aes(x = factor(year), y = prop, text = pdapN)) +
  geom_bar(fill = "blue", alpha = 0.5,
           stat="identity") + 
  labs(title = "PDAP program size: Proportion of HIV- population",
       x = "year", y = "proportion")

ggplotly(p, tooltip = "text")

3.3 Viral suppression

  • We calculated the proportion of viral suppression by ADAP status and race
  • Individuals might get checked multiple times a year. We used the viral load of the last visit in a year for each individual. In addition, if an individual did not visit clinic for viral load in a year, we assumed that the individual is not virally suppressed.
  • The definition of viral suppression is viral load < 200 copies of HIV per milliliter of blood.
  • We calculated the proportion of individuals who were virally suppressed by their ADAP status and race in the last two years.

3.3.1 Plot

  • This is a plot ZK originally made on the TS so I can’t easily pull the N for the tooltip.
p <- ggplot(adap_data$vl_race %>%
              mutate(ADAP = ifelse(ADAP==0, "no", "yes")), 
            aes(x = race, y = avg_vl_supp, text = round(avg_vl_supp, 2),
                group = ADAP, fill = ADAP)) +
  geom_bar(stat = "identity", alpha = 0.7, position = "dodge") +
  labs(title = "Viral suppression by race and ADAP status",
       x = "Race", y = "proportion")

ggplotly(p, tooltip = "text")
# grid.raster(readTIFF(here::here("MakeData", "Zoe", "EconModelBook", "data", 
#                                 "vl_race.tiff")))

3.3.2 Table

kable(adap_data$vl_race %>%
        select(ADAP, race, avg_vl_supp), digits = c(0, 0, 2), 
      caption= "proportion of viral suppression") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped"))
proportion of viral suppression
ADAP race avg_vl_supp
0 B 0.72
1 B 0.91
0 H 0.74
1 H 0.94
0 O 0.81
1 O 0.96

4 Generate whamp_data object for EpiModelHIV simulation

### Creating whamp_data used in EpiModelHIV-p

#names(msm.pop.totals_2019) <- gsub("[.]", "_", names(msm.pop.totals_2019))

whamp_data <- list()
#whamp_data$wa_pop_data <- msm.pop.totals_2019
#whamp_data$pred_inc <- whamp$pred_inc
#whamp_data$pred_ins <- whamp$pred_ins
whamp_data$pred_disenroll <- adap_data$pred_disenroll
whamp_data$adap_cost <- adap_data$ADAPcost_pd
whamp_data$adap.size <- adap.size
whamp_data$adap_cli <- adap_data$adap_demog[year == 1] # Initial distribution of adap client 

whamp_data$adap_cliIns <- adap_data$adap_ins

# whamp_data$prep_aware <- whamp$prep_aware
# whamp_data$prep_interest <- whamp$prep_interest
# whamp_data$prep_current <- whamp$prep_current

whamp_data$pdap_disenroll <- pdap_data$drop_prob_dat[ , years := as.numeric(years)]
whamp_data$pdap_cost <- pdap_data$PDAPcost_pd
whamp_data$pdap.size <- pdap.size
whamp_data$pdap_cli <- pdap_data$pdap_demog
whamp_data$pdap_cliIns <- pdap_data$pdap_ins
whamp_data$init_pdap_dist <- pdap_data$init_pdap_dat

# Last thing to check
init_pdap_n <- pdap_data$init_pdap_dat[, .N, by = c("race", "region")]
init_pdap_n[, n_samp := N]
init_pdap_n <- merge(init_pdap_n, pop.racexregion.neg, 
                     by = c("race", "region"))
init_pdap_n[, N := NULL]
init_pdap_n[, Npop := Freq]
init_pdap_n[, Freq := NULL]
whamp_data$init_pdap_n <- init_pdap_n

whamp_data$tar_vl_suppress <- adap_data$vl_race 
whamp_data$tt_data$tt_traj_dt <- tt_traj_dt
whamp_data$tt_data$tx_init_dt <- tx_init_dt
whamp_data$tt_data$tx_reinit_full_dt <- tx_reinit_full_dt

### data used for EpiModelHIV
saveRDS(whamp_data, "whamp_data.rds")