WHAMP DAP Descriptives

rm(list = ls())

library(dplyr)
library(ggplot2)
library(ggpubr)
library(plotly)
library(kableExtra)

runmode = "script" # to control how target summaries are built

# This report produces the WhampDapParams "adap.frac" and "pdap.frac"
# and it sources the makeDAPdesc.R file to construct a file with the
# DAP participation rates stratified by attrs.

# Subgroup sizes are too small when stratified to use as targets
# So stratified DAP participation rate targets are instead estimated
# from the DOH DAP data.


# Either make or read in data

dataset <- "readin"
#dataset <- "makeit"

if(dataset == "makeit") {
  source(here::here("MakeData", "MM", "Scripts",
                    "makeWideData.R"))
} else {
  wDF <- readRDS(here::here("MakeData", "MM", "Data",
                               "wideDF.RDS"))
}

Back to Outline

Introduction

This report presents the descriptive statistics for participation in the two DAPs (PDAP and ADAP) as reported in the WHAMP survey. All of the data used here come from the WHAMP Survey; ARTnet WA cases do not have information on DAP participation. The WHAMP survey was conducted over a three month period using social media (primarily FaceBook, with a small number of Grindr and Growlr respondents) from Sep 11 - Dec 11 2019.

---

The number of observations in the WHAMP sample to inform estimates of DAP participation is quite small, for both programs.

• For PDAP, the limiting factor is the sample number of MSM on PrEP: 161.

• For ADAP the limiting factor is the sample number of MSM who are Dx HIV+: 92.

So we only use them to estimate the overall DAP participation rate.

---

These are used as input parameters, but will be calibrated to DOH DAP program size estimates).

There are separate estimates of 2018 DAP participation stratified by attributes based on the WADOH program data – these are used for most other DAP parameter inputs. See make_DOHDAPdynamics.Rmd. The overall enrollment estimates from the two data sources are compared in that report.

---

Missing data

Missing values for ADAP and PDAP come primarily from survey attrition. For ADAP, 56 respondents had dropped out by that section (6.0%), and for PDAP the number was 84 (9.1%).

Missing values in the demographic breakdown variables are minimal (4 cases), so we use the complete case subset. There are a few missing insurance cases (6 cases) and we add these to the uninsured in the analyses, so we don’t lose them.

---

ADAP

• Enrollment fraction in the WA State ART-DAP (ADAP) program.

• Sample restricted to HIV+ currently on ART: 87.

• Too few cases to model or stratify, so we use only the overall fraction of ART users in ADAP as a target.

ART use

• ART use is nearly universal for HIV+ in the WHAMP survey

hivposDF <- cDF %>%
  filter(diag.status == 1) %>%
  mutate(artCurrent = factor(art.current,
                             levels = c(0, 1), 
                             labels = c("not on ART", "on ART")),
         ADAP = factor(ADAP,
                       levels = c(0, 1), 
                       labels = c("no ADAP", "ADAP")))

tempDF <- hivposDF %>%
  group_by(artCurrent) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

tempDF %>% kable(caption= "ART use among HIV+") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste("Sample restricted to HIV+, n =",
                 sum(tempDF$nobs)))

ART use among HIV+

artCurrent	nobs	n.wtd	prop.wtd
not on ART	1	0.9	0.01
on ART	87	89.7	0.94
NA	4	4.6	0.05
Note:
Sample restricted to HIV+, n = 92

Enrolled in ADAP

• This is restricted to HIV+ on ART, resulting in a very small sample size: 87.

• Because virtually all HIV+ respondents in WHAMP are on ART, there is no real difference in the proportion of ADAP participants among HIV+ or among ART users.

We show breakdowns by attributes, but the small cell sizes make these unreliable.

Overall

adapDF <- hivposDF %>% filter(art.current == 1)

adap_prop_art <- adapDF %>%
  group_by(ADAP) %>%
  summarise(nobs = n(),
            n.wtd = sum(ego.wawt)) %>%
  mutate(prop.wtd = n.wtd / sum(n.wtd))

adap_prop_art %>% kable(caption= "ADAP enrollment",
                    digits = c(0,0,1,3)) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste("Sample restricted to HIV+ on ART, n =",
                 sum(adap_prop_art$nobs)))

ADAP enrollment

ADAP	nobs	n.wtd	prop.wtd
no ADAP	58	62.9	0.701
ADAP	29	26.8	0.299
Note:
Sample restricted to HIV+ on ART, n = 87

# ADAP enrollment fraction starting parameter value
adap.frac.art <- adap_prop_art$prop.wtd[2]

# Calculate the adap fraction of HIV+ (will be almost the same)
adap_prop_pos <- hivposDF %>%
  group_by(ADAP) %>%
  summarise(nobs = n(),
            n.wtd = sum(ego.wawt)) %>%
  mutate(prop.wtd = n.wtd / sum(n.wtd))

# ADAP enrollment fraction starting parameter value
adap.frac.pos <- adap_prop_pos$prop.wtd[2]

Age

Note: No 15-24yo on ART

adap.prop.age <- adapDF %>%
  filter(art.current == 1) %>%
  group_by(ADAP, age_group) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(age_group) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

adap.prop.age %>% 
  tidyr::pivot_wider(names_from = ADAP,
              values_from = c(nobs,  n.wtd, prop.wtd)) %>%
  kable(col.names = c("Age", rep(c("no ADAP", "ADAP"), 3)),
        caption= paste("ADAP enrollment by age")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current ART users; n = ",
                 sum(adap.prop.age$nobs)))

ADAP enrollment by age

	Obs N		Wtd N		Wtd Prop
Age	no ADAP	ADAP	no ADAP	ADAP	no ADAP	ADAP
25-34	7	1	6.3	1.0	0.86	0.14
35-44	11	7	14.1	6.4	0.69	0.31
45-54	12	6	14.0	5.6	0.71	0.29
55-65	28	15	28.5	13.9	0.67	0.33
Note:
Restricted to current ART users; n = 87

Race

adap.prop.race <- adapDF %>%
  group_by(ADAP, race) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(race) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

adap.prop.race %>% 
  tidyr::pivot_wider(names_from = ADAP,
              values_from = c(nobs,  n.wtd, prop.wtd)) %>%
  kable(col.names = c("Race", rep(c("no ADAP", "ADAP"), 3)),
        caption= paste("ADAP enrollment by race")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current ART users; n = ",
                 sum(adap.prop.race$nobs)))

ADAP enrollment by race

	Obs N		Wtd N		Wtd Prop
Race	no ADAP	ADAP	no ADAP	ADAP	no ADAP	ADAP
B	1	2	2.4	3.0	0.44	0.56
H	5	2	9.4	1.9	0.83	0.17
O	52	25	51.1	22.0	0.70	0.30
Note:
Restricted to current ART users; n = 87

Region

adap.prop.region <- adapDF %>%
  group_by(ADAP, region) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(region) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

adap.prop.region %>% 
  tidyr::pivot_wider(names_from = ADAP,
              values_from = c(nobs,  n.wtd, prop.wtd)) %>%
  kable(col.names = c("Region", rep(c("no ADAP", "ADAP"), 3)),
        caption= paste("ADAP enrollment by region")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current ART users; n = ",
                 sum(adap.prop.region$nobs)))

ADAP enrollment by region

	Obs N		Wtd N		Wtd Prop
Region	no ADAP	ADAP	no ADAP	ADAP	no ADAP	ADAP
EasternWA	4	5	2.4	2.9	0.45	0.55
King	30	12	38.3	14.4	0.73	0.27
WesternWA	24	12	22.2	9.5	0.70	0.30
Note:
Restricted to current ART users; n = 87

Insurance

adap.prop.ins <- adapDF %>%
  group_by(ADAP, insurance) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(insurance) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

adap.prop.ins %>% 
  tidyr::pivot_wider(names_from = ADAP,
              values_from = c(nobs,  n.wtd, prop.wtd)) %>%
  kable(col.names = c("Insurance", rep(c("no ADAP", "ADAP"), 3)),
        caption= paste("ADAP enrollment by insurance")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current ART users; n = ",
                 sum(adap.prop.ins$nobs)))

ADAP enrollment by insurance

	Obs N		Wtd N		Wtd Prop
Insurance	no ADAP	ADAP	no ADAP	ADAP	no ADAP	ADAP
emplyr	45	6	48.7	5.8	0.89	0.11
indvdl_othr	1	4	1.2	3.0	0.29	0.71
medicare_caid	7	12	6.3	10.6	0.37	0.63
othr_gov	4	7	5.9	7.4	0.44	0.56
uninsrd_NA	1	NA	0.7	NA	1.00	NA
Note:
Restricted to current ART users; n = 87

Frequency of ART related check-ups

• Most of the respondents who are on ART visit clinics every 3-6 months.

kable(data.frame(table(cDF$ARTCHKFREQ_CURR)), 
      caption= "Frequency of ART related check-ups ") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste("Current ART users; n = ",
                 nrow(adapDF)))

Frequency of ART related check-ups

Var1	Freq
Every month	0
Every 3 months	18
Every 6 months	54
Every 9 months	2
Every 12 months	9
I returned only once	0
I return at some other interval of time	1
I did not return regularly	0
I dont know	0
I prefer not to answer	0
Note:
Current ART users; n = 87

---

PDAP

Considerations

There are two forms of PrEP drug assistance:

• the State program (PDAPA in the WHAMP survey)

• the Gilead programs (“Map” and “Cap” in the DOH PDAP dataset, “PDAPB” in the WHAMP survey)

• We track enrollment in the State program only.

The WHAMP survey shows that there is some overlap between these programs, but the enrollment is uncorrelated (the marginal enrollment probabilities are equal to the joint).

# with(cDF, table(PREP_CURRENT, statePDAP = PDAPA, useNA = "al"))
# with(cDF, table(PREP_CURRENT, gileadPDP = PDAPB, useNA = "al"))
# with(cDF, table(statePDAP = PDAPA, gileadPDP = PDAPB, useNA = "al"))

sjPlot::tab_xtab(var.row = cDF$PDAPA, 
                 var.col = cDF$PDAPB, 
                 title = "WHAMP survey PDAP enrollment",
                 var.labels = c("State PDAP", "Gilead PDAP"),
                 show.row.prc = T,
                 show.col.prc = T,
                 show.summary = F)

WHAMP survey PDAP enrollment

State PDAP	Gilead PDAP		Total
	Unchecked	Checked
Unchecked	48 36.1 % 82.8 %	85 63.9 % 83.3 %	133 100 % 83.1 %
Checked	10 37 % 17.2 %	17 63 % 16.7 %	27 100 % 16.9 %
Total	58 36.2 % 100 %	102 63.8 % 100 %	160 100 % 100 %

PrEP Use

• PrEP use is limited in the WHAMP survey, which limits the PDAP enrollment data we observe.

• DRW’s 2018 WHPP survey estimate of HIV- PrEP users is 18%.

# we exclude the dropouts is.na(prep.current)
hivnegDF <- cDF %>%
  filter(diag.status == 0 & !is.na(prep.current)) %>%
  mutate(prepCurrent = 
           factor(prep.current,
                  levels = c(0, 1), 
                  labels = c("not on PrEP", "on PrEP")))

tempDF <- hivnegDF %>%
  group_by(prepCurrent) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

tempDF %>% kable(caption= "Prevalence of PrEP use") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste0("Sample restricted to HIV-, n = ",
                 sum(tempDF$nobs),
                 "\nalso excludes ",
           sum(is.na(cDF$prep.current)), " dropouts"))

Prevalence of PrEP use

prepCurrent	nobs	n.wtd	prop.wtd
not on PrEP	594	581.3	0.78
on PrEP	161	167.5	0.22
Note:
Sample restricted to HIV-, n = 755 also excludes 83 dropouts

Enrolled in PDAP

We only use the overall estimate as a target. We show the breakdowns here, but the cell sizes are too small to be reliable.

Overall

Of HIV-

# We track enrollment in the State program only
pdapnegDF <- hivnegDF %>%
  mutate(PDAP = case_when(PDAPA == 1  ~ "PDAP",
                          TRUE ~ "no PDAP",),
         statePDAP = PDAPA,
         gileadPDAP = PDAPB)

pdap_prop_neg <- pdapnegDF %>%
  group_by(PDAP) %>%
  summarise(nobs = n(),
            n.wtd = sum(ego.wawt)) %>%
  mutate(prop.wtd = n.wtd / sum(n.wtd))

pdap_prop_neg %>% 
  kable(caption= "State PDAP enrollment among HIV-",
        digits = c(0,0,1,3)) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste("Sample restricted to HIV-, n =",
                 sum(pdap_prop_neg$nobs)))

State PDAP enrollment among HIV-

PDAP	nobs	n.wtd	prop.wtd
no PDAP	728	720.9	0.963
PDAP	27	27.9	0.037
Note:
Sample restricted to HIV-, n = 755

pdap.frac.neg <- pdap_prop_neg$prop.wtd[2]

Of current PrEP users

prepDF <- pdapnegDF %>%
  filter(prep.current == 1)

pdap_prop_prep <- prepDF %>%
  group_by(PDAP) %>%
  summarise(nobs = n(),
            n.wtd = sum(ego.wawt)) %>%
  mutate(prop.wtd = n.wtd / sum(n.wtd))

pdap_prop_prep %>% 
  kable(caption= "State PDAP enrollment among current PrEP users",
        digits = c(0,0,1,3)) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste("Sample restricted to current PrEP users, n =",
                 sum(pdap_prop_prep$nobs)))

State PDAP enrollment among current PrEP users

PDAP	nobs	n.wtd	prop.wtd
no PDAP	134	139.6	0.833
PDAP	27	27.9	0.167
Note:
Sample restricted to current PrEP users, n = 161

pdap.frac.prep <- pdap_prop_prep$prop.wtd[2]

Age

Table

pdap.prop.age <- prepDF %>%
  group_by(PDAP, age_group) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(age_group) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

pdap.prop.age %>% 
  tidyr::pivot_wider(names_from = PDAP,
              values_from = c(nobs, n.wtd, prop.wtd)) %>%
  kable(col.names = c("Age", rep(c("no PDAP", "PDAP"), 3)),
        caption= paste("PDAP enrollment by age")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current PrEP users; n = ",
                 sum(pdap_prop_prep$nobs)))

PDAP enrollment by age

	Obs N		Wtd N		Wtd Prop
Age	no PDAP	PDAP	no PDAP	PDAP	no PDAP	PDAP
15-24	6	1	8.1	1.8	0.82	0.18
25-34	45	10	42.2	10.1	0.81	0.19
35-44	30	8	36.5	8.1	0.82	0.18
45-54	39	2	39.1	1.9	0.95	0.05
55-65	14	6	13.6	6.0	0.69	0.31
Note:
Restricted to current PrEP users; n = 161

Plot

fig <- ggplot(pdap.prop.age, 
              aes(x = PDAP, y = prop.wtd, text = prop.wtd,
                  group = age_group, fill = age_group)) +
  geom_bar(stat="identity", position = "dodge", 
           color = "darkgrey", alpha = 0.7) +
  labs(x = "PDAP enrollment",
       y = "wtd. proportion") +
  theme(plot.margin = margin(t=20,r=5,b=5,l=0))

title.plotly <-
  list(text = paste0('PDAP enrollment by age group',
                     '<br>',
                     '<sup>',
                     'n = ',
                     sum(pdap_prop_prep$nobs),
                     '</sup>'))

ggplotly(fig , tooltip = "text") %>%
  layout(#hoverlabel = list(font=list(color="black")),
    title = title.plotly)

Race

Table

pdap.prop.race <- prepDF %>%
  group_by(PDAP, race) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(race) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

pdap.prop.race %>% 
  tidyr::pivot_wider(names_from = PDAP,
              values_from = c(nobs, n.wtd, prop.wtd)) %>%
  kable(col.names = c("Race", rep(c("no PDAP", "PDAP"), 3)),
        caption= paste("PDAP enrollment by race")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current PrEP users; n = ",
                 sum(pdap_prop_prep$nobs)))

PDAP enrollment by race

	Obs N		Wtd N		Wtd Prop
Race	no PDAP	PDAP	no PDAP	PDAP	no PDAP	PDAP
B	6	NA	17.1	NA	1.00	NA
H	11	2	9.5	1.9	0.83	0.17
O	117	25	113.0	26.0	0.81	0.19
Note:
Restricted to current PrEP users; n = 161

Plot

fig <- ggplot(pdap.prop.race, 
              aes(x = PDAP, y = prop.wtd, text = prop.wtd,
                  group = race, fill = race)) +
  geom_bar(stat="identity", position = "dodge", color = "darkgrey",
           alpha = 0.7, ) +
  labs(x = "PDAP enrollment",
       y = "wtd. proportion") +
  theme(plot.margin = margin(t=20,r=5,b=5,l=0)) 

title.plotly <-
  list(text = paste0('PDAP enrollment by race',
                     '<br>',
                     '<sup>',
                     'HIV-; n = ',
                     sum(pdap_prop_prep$nobs),
                     '</sup>'))

ggplotly(fig , tooltip = "text") %>%
  layout(#hoverlabel = list(font=list(color="black")),
    title = title.plotly)

Region

Table

pdap.prop.region <- prepDF %>%
  group_by(PDAP, region) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(region) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

pdap.prop.region %>% 
  tidyr::pivot_wider(names_from = PDAP,
              values_from = c(nobs, n.wtd, prop.wtd)) %>%
  kable(col.names = c("Region", rep(c("no PDAP", "PDAP"), 3)),
        caption= paste("PDAP enrollment by region")) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Restricted to current PrEP users; n = ",
                 sum(pdap_prop_prep$nobs)))

PDAP enrollment by region

	Obs N		Wtd N		Wtd Prop
Region	no PDAP	PDAP	no PDAP	PDAP	no PDAP	PDAP
EasternWA	11	NA	7.0	NA	1.00	NA
King	75	19	94.9	21.7	0.81	0.19
WesternWA	48	8	37.7	6.2	0.86	0.14
Note:
Restricted to current PrEP users; n = 161

Plot

fig <- ggplot(pdap.prop.region, 
              aes(x = PDAP, y = prop.wtd, text = prop.wtd,
                  group = region, fill = region)) +
  geom_bar(stat="identity", position = "dodge", color = "darkgrey",
           alpha = 0.7, ) +
  labs(x = "PDAP enrollment",
       y = "wtd. proportion") +
  theme(plot.margin = margin(t=20,r=5,b=5,l=0)) 

title.plotly <-
  list(text = paste0('PDAP enrollment by region',
                     '<br>',
                     '<sup>',
                     'HIV-; n = ',
                     sum(pdap_prop_prep$nobs),
                     '</sup>'))

ggplotly(fig , tooltip = "text") %>%
  layout(#hoverlabel = list(font=list(color="black")),
    title = title.plotly)

Insurance

Table

pdap.prop.ins <- prepDF %>%
  group_by(insurance, PDAP) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(insurance) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

pdap.prop.ins %>% 
  tidyr::pivot_wider(names_from = PDAP,
              values_from = c(nobs,  n.wtd, prop.wtd)) %>%
  kable(col.names = c("Region", rep(c("no PDAP", "PDAP"), 3)),
        caption= "PDAP enrollment by insurance") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Proportion of current PrEP users; n = ",
                 sum(pdap.prop.ins$nobs)))

PDAP enrollment by insurance

	Obs N		Wtd N		Wtd Prop
Region	no PDAP	PDAP	no PDAP	PDAP	no PDAP	PDAP
emplyr	104	16	111.3	16.7	0.87	0.13
indvdl_othr	6	7	5.8	7.7	0.43	0.57
medicare_caid	11	4	10.1	3.4	0.75	0.25
othr_gov	9	NA	8.9	NA	1.00	NA
uninsrd_NA	4	NA	3.4	NA	1.00	NA
Note:
Proportion of current PrEP users; n = 161

Plot

fig <- ggplot(pdap.prop.ins, 
              aes(x = PDAP, y = prop.wtd, text = prop.wtd,
                  group = insurance, fill = insurance)) +
  geom_bar(stat="identity", position = "dodge", color = "darkgrey",
           alpha = 0.7, ) +
  labs(x = "PDAP enrollment",
       y = "wtd. proportion") +
  theme(plot.margin = margin(t=20,r=5,b=5,l=0)) 

title.plotly <-
  list(text = paste0('PDAP enrollment by insurance',
                     '<br>',
                     '<sup>',
                     'HIV-; n = ',
                     sum(pdap.prop.ins$nobs),
                     '</sup>'))

ggplotly(fig , tooltip = "text") %>%
  layout(#hoverlabel = list(font=list(color="black")),
    title = title.plotly)

Income

Table

tempDF <- prepDF %>%
  group_by(income_cat, PDAP) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(income_cat) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

tempDF %>% 
  tidyr::pivot_wider(names_from = PDAP,
              values_from = c(nobs,  n.wtd, prop.wtd)) %>%
  kable(col.names = c("Income", rep(c("no PDAP", "PDAP"), 3)),
        caption= "PDAP enrollment by income") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Proportion of current PrEP users; n = ",
                 sum(tempDF$nobs)))

PDAP enrollment by income

	Obs N		Wtd N		Wtd Prop
Income	no PDAP	PDAP	no PDAP	PDAP	no PDAP	PDAP
fpl138	5	1	4.0	0.7	0.85	0.15
fpl138-300	19	4	18.4	3.9	0.83	0.17
fpl300-400	16	5	15.2	5.4	0.74	0.26
fpl400-500	10	NA	10.6	NA	1.00	NA
fpl500-600	9	2	7.5	2.1	0.78	0.22
fpl600+	72	15	80.6	15.7	0.84	0.16
NA	3	NA	3.3	NA	1.00	NA
Note:
Proportion of current PrEP users; n = 161

Plot

fig <- ggplot(tempDF, 
              aes(x = PDAP, y = prop.wtd, text = prop.wtd,
                  group = income_cat, fill = income_cat)) +
  geom_bar(stat="identity", position = "dodge", color = "darkgrey",
           alpha = 0.7, ) +
  labs(x = "PDAP enrollment",
       y = "wtd. proportion") +
  theme(plot.margin = margin(t=20,r=5,b=5,l=0)) 

title.plotly <-
  list(text = paste0('PDAP enrollment by FPL income category',
                     '<br>',
                     '<sup>',
                     'HIV-; n = ',
                     sum(tempDF$nobs),
                     '</sup>'))

ggplotly(fig , tooltip = "text") %>%
  layout(#hoverlabel = list(font=list(color="black")),
    title = title.plotly)

---

Interruptions

• For current PrEP users, ever interrupted, by PDAP participation.

Prevalence

• Interruption is more likely if not on PDAP

• But, we don’t know if the PDAP interruptions we observe happened when they were on PDAP

• We’ll also look at the length and reasons for the interruption.

Table

tempDF <- prepDF %>%
  group_by(PDAP, prep.int) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1)) %>%
  group_by(PDAP) %>%
  mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))

tempDF %>% 
  tidyr::pivot_wider(names_from = prep.int,
              values_from = c(nobs, n.wtd, prop.wtd)) %>%
  kable(col.names = c("PDAP", rep(c("no", "yes"), 3)),
        caption= "Ever interrupted PrEP by PDAP enrollment") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Current PrEP; n = ",
                 sum(tempDF$nobs)))

Ever interrupted PrEP by PDAP enrollment

	Obs N		Wtd N		Wtd Prop
PDAP	no	yes	no	yes	no	yes
no PDAP	103	31	112.3	27.2	0.81	0.19
PDAP	24	3	24.8	3.1	0.89	0.11
Note:
Current PrEP; n = 161

Plot

fig <- ggplot(tempDF %>% 
                mutate(prepint = ifelse(prep.int==0, "no", "yes")),
              aes(x = prepint, y = prop.wtd, text = n.wtd,
                  group = PDAP, fill = PDAP)) +
  geom_bar(stat="identity", position = "dodge", color = "darkgrey",
           alpha = 0.7, ) +
  labs(x = "ever interrupted PrEP",
       y = "wtd. proportion") +
  theme(plot.margin = margin(t=20,r=5,b=5,l=0)) 

title.plotly <-
  list(text = paste0('Ever interrupted by PDAP enrollment',
                     '<br>',
                     '<sup>',
                     'current PrEP users; n = ',
                     sum(tempDF$nobs),
                     '</sup>'))

ggplotly(fig , tooltip = "text") %>%
  layout(#hoverlabel = list(font=list(color="black")),
    title = title.plotly)

Interruption length

• Not much difference

tempDF <- prepDF %>%
  filter(prep.int == 1) %>%
  group_by(PDAP) %>%
  summarise(nobs = n(),
            n.wtd = round(sum(ego.wawt), 1),
            mean.wtd = round(Hmisc::wtd.mean(PREP_INTMOS, ego.wawt),1),
            median.wtd = Hmisc::wtd.quantile(PREP_INTMOS, ego.wawt,
                                             probs = 0.5)) 

t(tempDF)[-1,] %>% 
  kable(col.names = c("no", "yes"),
        caption= "PrEP interruption length by PDAP enrollment") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
   footnote(paste("Current PrEP users who ever interrupted; n = ",
                 sum(tempDF$nobs)))

PrEP interruption length by PDAP enrollment

	no	yes
nobs	31	3
n.wtd	27.2	3.1
mean.wtd	3.5	3.3
median.wtd	2.124291	2.144912
Note:
Current PrEP users who ever interrupted; n = 34

Reason

• REALLY small N’s here
  
• No risk is still either the largest or second largest reason
  
• Some indication that no HCP was important for PDAP participants
  
• No indication that cost was the issue, for either group.

tempDF <- prepDF %>%
  filter(prep.int == 1) %>%
  group_by(PDAP, prep.whyint) %>%
  dplyr::summarize(nobs = n(),
            n.wtd = sum(ego.wawt)) %>%
  group_by(PDAP) %>%
  mutate(prop.wtd = n.wtd/sum(n.wtd)) 

tempDF %>%
  tidyr::pivot_wider(names_from = PDAP,
              values_from = c(nobs, n.wtd, prop.wtd)) %>%
  kable(col.names = c("PDAP", rep(c("no PDAP", "PDAP"), 3)),
        caption= "Reason for interrupting by PDAP enrollment",
        digits = c(rep(0,3), rep(2,4))) %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  add_header_above(c("", "Obs N" = 2, "Wtd N" = 2, "Wtd Prop"= 2)) %>%
  footnote(paste("Current PrEP users who ever interrupted; n = ",
                 sum(tempDF$nobs)))

Reason for interrupting by PDAP enrollment

	Obs N		Wtd N		Wtd Prop
PDAP	no PDAP	PDAP	no PDAP	PDAP	no PDAP	PDAP
No HCP for Rx	6	1	6.78	1.20	0.25	0.39
No risk	14	2	11.74	1.86	0.43	0.61
Could not afford	4	NA	3.68	NA	0.14	NA
Side effects	1	NA	0.74	NA	0.03	NA
Other	6	NA	4.31	NA	0.16	NA
Note:
Current PrEP users who ever interrupted; n = 34

---

Frequency of PREP related check-ups

• Almost all respondents who are on PrEP visit clinics every 3 months.

kable(data.frame(table(cDF$PREPCHKFREQ_CURR)), 
      caption= "Frequency of PREP related check-ups ") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped")) %>%
  footnote(paste("Current PrEP users; n = ",
                 nrow(prepDF)))

Frequency of PREP related check-ups

Var1	Freq
Every month	1
Every 3 months	137
Every 6 months	11
Every 9 months	0
Every 12 months	5
I have returned only once	1
I return at some other interval of time	1
I do not return regularly	4
Im not sure	0
I prefer not to answer	0
Note:
Current PrEP users; n = 161

---

Construct and save output

The two overall enrollment fractions are saved as parameters.

# Parameters

adap.frac = list(art = adap.frac.art,
                 n.art = sum(adap_prop_art$nobs),
                 hivpos = adap.frac.pos,
                 n.pos = sum(adap_prop_pos$nobs))

pdap.frac = list(prep = pdap.frac.prep,
                 n.prep = sum(pdap_prop_prep$nobs),
                 hivneg = pdap.frac.neg,
                 n.neg = sum(pdap_prop_neg$nobs))

whampDAP <- list(adap.frac = adap.frac,
                 pdap.frac = pdap.frac,
                 makefile = "make_WhampDAP.Rmd")
descTable <- 
  tibble(Params = names(whampDAP), 
         Description = c("ADAP enrollment fractions & nobs", 
                         "PDAP enrollment fractions & nobs",
                         "source file"),
         Subset = c("ART current or HIV+",
                    "PrEP current or HIV-",
                    " "),
         Method = c(rep("wtd proportions", 2), " "),
         Levels = c(rep("overall", 2), " "),
         Notes = c("adap.frac art is starting value for calibration",
                   "pdap.frac prep is starting value for calibration",
                   "WHAMP survey estimates"))

whampDAP <- c(whampDAP, list(descTable = descTable))

saveRDS(whampDAP,
        file = here::here("Data", "Params", "WhampDAPparam.RDS"))

descTable %>%
  kable(caption= "DAP enrollment rate estimates from WHAMP survey") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped"))

DAP enrollment rate estimates from WHAMP survey

Params	Description	Subset	Method	Levels	Notes
adap.frac	ADAP enrollment fractions & nobs	ART current or HIV+	wtd proportions	overall	adap.frac art is starting value for calibration
pdap.frac	PDAP enrollment fractions & nobs	PrEP current or HIV-	wtd proportions	overall	pdap.frac prep is starting value for calibration
makefile	source file				WHAMP survey estimates

print("Structure of output object:")

## [1] "Structure of output object:"

str(whampDAP)

## List of 4
##  $ adap.frac:List of 4
##   ..$ art   : num 0.299
##   ..$ n.art : int 87
##   ..$ hivpos: num 0.282
##   ..$ n.pos : int 92
##  $ pdap.frac:List of 4
##   ..$ prep  : num 0.167
##   ..$ n.prep: int 161
##   ..$ hivneg: num 0.0372
##   ..$ n.neg : int 755
##  $ makefile : chr "make_WhampDAP.Rmd"
##  $ descTable: tibble[,6] [3 x 6] (S3: tbl_df/tbl/data.frame)
##   ..$ Params     : chr [1:3] "adap.frac" "pdap.frac" "makefile"
##   ..$ Description: chr [1:3] "ADAP enrollment fractions & nobs" "PDAP enrollment fractions & nobs" "source file"
##   ..$ Subset     : chr [1:3] "ART current or HIV+" "PrEP current or HIV-" " "
##   ..$ Method     : chr [1:3] "wtd proportions" "wtd proportions" " "
##   ..$ Levels     : chr [1:3] "overall" "overall" " "
##   ..$ Notes      : chr [1:3] "adap.frac art is starting value for calibration" "pdap.frac prep is starting value for calibration" "WHAMP survey estimates"

Stratified estimates (of those engaged in care) are saved as descriptives.

They should not be used for either targets or input parameters, due to small cell sizes.

The construction uses an external script – it simply calculates standard summary stats for the PDAP fraction of PrEP users and the ADAP fraction of ART users.

# Descriptives

source(here::here("MakeData", "MM", "Scripts",
                  "makeWhampDAPdesc.R"))

descOutput <- list(adap=adap, pdap=pdap)

descTable <- 
  tibble(Targets = names(descOutput), 
         Description = c("ADAP enrollment fraction", 
                         "PDAP enrollment fraction"),
         Subset = c("ART Users",
                    "PrEP Users"),
         Method = rep("wtd obs summaries", 2),
         Levels = rep("overall and by age, race, region & insurance", 2),
         Notes = rep("Sample size too small for target use", 2))

descOutput = c(descOutput, list(descTable = descTable))
saveRDS(descOutput,
        file = here::here("Data", "Descriptives", "WhampDAP.RDS"))

descTable %>%
  kable(caption= "DAP enrollment rate estimates from WHAMP survey") %>% 
  kable_styling(full_width=F, position="center", 
                bootstrap_options = c("striped"))

DAP enrollment rate estimates from WHAMP survey

Targets	Description	Subset	Method	Levels	Notes
adap	ADAP enrollment fraction	ART Users	wtd obs summaries	overall and by age, race, region & insurance	Sample size too small for target use
pdap	PDAP enrollment fraction	PrEP Users	wtd obs summaries	overall and by age, race, region & insurance	Sample size too small for target use

print("Structure of output object:")

## [1] "Structure of output object:"

str(descOutput)

## List of 3
##  $ adap     :List of 5
##   ..$ all   : tibble[,7] [1 x 7] (S3: tbl_df/tbl/data.frame)
##   .. ..$ nobs      : int 87
##   .. ..$ n.valid   : int 87
##   .. ..$ n.missing : int 0
##   .. ..$ wtd.n     : num 89.7
##   .. ..$ wtd.mean  : num 0.299
##   .. ..$ wtd.sd    : num 0.46
##   .. ..$ wtd.semean: num 0.0486
##   ..$ age   : tibble[,8] [4 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ age.grp   : num [1:4] 2 3 4 5
##   .. ..$ nobs      : int [1:4] 8 18 18 43
##   .. ..$ n.valid   : int [1:4] 8 18 18 43
##   .. ..$ n.missing : int [1:4] 0 0 0 0
##   .. ..$ wtd.n     : num [1:4] 7.23 20.47 19.61 42.37
##   .. ..$ wtd.mean  : num [1:4] 0.135 0.311 0.285 0.327
##   .. ..$ wtd.sd    : num [1:4] 0.368 0.474 0.463 0.475
##   .. ..$ wtd.semean: num [1:4] 0.137 0.105 0.105 0.073
##   ..$ race  : tibble[,8] [3 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ race      : chr [1:3] "B" "H" "O"
##   .. ..$ nobs      : int [1:3] 3 7 77
##   .. ..$ n.valid   : int [1:3] 3 7 77
##   .. ..$ n.missing : int [1:3] 0 0 0
##   .. ..$ wtd.n     : num [1:3] 5.34 11.29 73.06
##   .. ..$ wtd.mean  : num [1:3] 0.556 0.164 0.301
##   .. ..$ wtd.sd    : num [1:3] 0.551 0.388 0.462
##   .. ..$ wtd.semean: num [1:3] 0.239 0.116 0.054
##   ..$ region: tibble[,8] [3 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ region    : chr [1:3] "EasternWA" "King" "WesternWA"
##   .. ..$ nobs      : int [1:3] 9 42 36
##   .. ..$ n.valid   : int [1:3] 9 42 36
##   .. ..$ n.missing : int [1:3] 0 0 0
##   .. ..$ wtd.n     : num [1:3] 5.3 52.7 31.7
##   .. ..$ wtd.mean  : num [1:3] 0.542 0.273 0.3
##   .. ..$ wtd.sd    : num [1:3] 0.553 0.45 0.466
##   .. ..$ wtd.semean: num [1:3] 0.2402 0.062 0.0827
##   ..$ ins   : tibble[,8] [5 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ insurance : chr [1:5] "emplyr" "indvdl_othr" "medicare_caid" "othr_gov" ...
##   .. ..$ nobs      : int [1:5] 51 5 19 11 1
##   .. ..$ n.valid   : int [1:5] 51 5 19 11 1
##   .. ..$ n.missing : int [1:5] 0 0 0 0 0
##   .. ..$ wtd.n     : num [1:5] 54.493 4.267 16.881 13.302 0.741
##   .. ..$ wtd.mean  : num [1:5] 0.106 0.708 0.629 0.557 0
##   .. ..$ wtd.sd    : num [1:5] 0.31 0.52 0.498 0.517 0
##   .. ..$ wtd.semean: num [1:5] 0.042 0.252 0.121 0.142 0
##  $ pdap     :List of 5
##   ..$ all   : tibble[,7] [1 x 7] (S3: tbl_df/tbl/data.frame)
##   .. ..$ nobs      : int 161
##   .. ..$ n.valid   : int 161
##   .. ..$ n.missing : int 0
##   .. ..$ wtd.n     : num 167
##   .. ..$ wtd.mean  : num 0.167
##   .. ..$ wtd.sd    : num 0.374
##   .. ..$ wtd.semean: num 0.0289
##   ..$ age   : tibble[,8] [5 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ age.grp   : num [1:5] 1 2 3 4 5
##   .. ..$ nobs      : int [1:5] 7 55 38 41 20
##   .. ..$ n.valid   : int [1:5] 7 55 38 41 20
##   .. ..$ n.missing : int [1:5] 0 0 0 0 0
##   .. ..$ wtd.n     : num [1:5] 9.95 52.27 44.65 41 19.61
##   .. ..$ wtd.mean  : num [1:5] 0.183 0.193 0.182 0.046 0.304
##   .. ..$ wtd.sd    : num [1:5] 0.407 0.399 0.39 0.212 0.472
##   .. ..$ wtd.semean: num [1:5] 0.1292 0.0551 0.0584 0.0331 0.1067
##   ..$ race  : tibble[,8] [3 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ race      : chr [1:3] "B" "H" "O"
##   .. ..$ nobs      : int [1:3] 6 13 142
##   .. ..$ n.valid   : int [1:3] 6 13 142
##   .. ..$ n.missing : int [1:3] 0 0 0
##   .. ..$ wtd.n     : num [1:3] 17.1 11.4 139
##   .. ..$ wtd.mean  : num [1:3] 0 0.164 0.187
##   .. ..$ wtd.sd    : num [1:3] 0 0.388 0.392
##   .. ..$ wtd.semean: num [1:3] 0 0.115 0.0332
##   ..$ region: tibble[,8] [3 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ region    : chr [1:3] "EasternWA" "King" "WesternWA"
##   .. ..$ nobs      : int [1:3] 11 94 56
##   .. ..$ n.valid   : int [1:3] 11 94 56
##   .. ..$ n.missing : int [1:3] 0 0 0
##   .. ..$ wtd.n     : num [1:3] 6.99 116.63 43.86
##   .. ..$ wtd.mean  : num [1:3] 0 0.186 0.141
##   .. ..$ wtd.sd    : num [1:3] 0 0.391 0.352
##   .. ..$ wtd.semean: num [1:3] 0 0.0362 0.0531
##   ..$ ins   : tibble[,8] [5 x 8] (S3: tbl_df/tbl/data.frame)
##   .. ..$ insurance : chr [1:5] "emplyr" "indvdl_othr" "medicare_caid" "othr_gov" ...
##   .. ..$ nobs      : int [1:5] 120 13 15 9 4
##   .. ..$ n.valid   : int [1:5] 120 13 15 9 4
##   .. ..$ n.missing : int [1:5] 0 0 0 0 0
##   .. ..$ wtd.n     : num [1:5] 128.05 13.51 13.57 8.9 3.45
##   .. ..$ wtd.mean  : num [1:5] 0.131 0.571 0.253 0 0
##   .. ..$ wtd.sd    : num [1:5] 0.338 0.514 0.452 0 0
##   .. ..$ wtd.semean: num [1:5] 0.0299 0.1399 0.1226 0 0
##  $ descTable: tibble[,6] [2 x 6] (S3: tbl_df/tbl/data.frame)
##   ..$ Targets    : chr [1:2] "adap" "pdap"
##   ..$ Description: chr [1:2] "ADAP enrollment fraction" "PDAP enrollment fraction"
##   ..$ Subset     : chr [1:2] "ART Users" "PrEP Users"
##   ..$ Method     : chr [1:2] "wtd obs summaries" "wtd obs summaries"
##   ..$ Levels     : chr [1:2] "overall and by age, race, region & insurance" "overall and by age, race, region & insurance"
##   ..$ Notes      : chr [1:2] "Sample size too small for target use" "Sample size too small for target use"