WHAMP Survey PrEP Analyses
Estimates input parameters and output targets for the PrEP Continuum
ZK, MM
2021-06-04
rm(list = ls())
library(tidyr)
library(stringr)
library(skimr)
library(dplyr)
library(broom)
library(Hmisc)
library(ggplot2)
library(ggpubr)
library(plotly)
library(gridExtra)
library(nnet)
library(kableExtra)
library(survey)
# Define a couple of functions for model comparison tables
# for formulas
get_formula <- function(prefix, n) {
forms <- eval(parse(text =
paste0("list(",
paste0(prefix,
1:n,
"$terms", collapse = ","),
")")))
out <- sapply(forms, FUN = format)
return(out)
}
# for aics, need to also specify model type because
# object structures returned from AIC() are different
get_aics <- function(prefix, n, type) {
mods <- eval(parse(text =
paste0("list(",
paste0(prefix,
1:n,
collapse = ","),
")")))
aics <- sapply(mods, FUN = AIC)
if(type == "mlm"){
out = aics
} else {
out = aics[2,]
}
return(out)
}
# Adjust annotations for facet_grid plots
# https://github.com/ropensci/plotly/issues/1224#issuecomment-727175492
# str(gp[['x']][['layout']][['annotations']])
layout_ggplotly <- function(
gg, x = -0.07, y = -0.065,
x_legend=1.05, y_legend=0.95,
mar=list(l=50, r=150)){
# The 1 and 2 contains the x and y axis labels, 9 is 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)
}
# Create the output object
save_out <- list()Back to Outline
Introduction
This report covers the data management and modeling associated with the PrEP continuum:
Input parameters for EpiModel – these will be read in via the
param_msmfunction inmod.initialize, and used in theprep.continuumfunction inmod.pdapOutput Dx targets – will be used in the PrEP output Dx report.
PrEP continuum
The PrEP continuum is represented as:
Aware -> Initiate -> Discontinue -> {Quit or Re-initiate}Parameters – Each flow is parameterized by a rate, and the rates are based on model predictions and back calculations. Model terms vary, and were selected on the basis of AICs (and the ability to use the terms in the simulation). The parameters from these are stored in data.tables. Backcalculations are often working with much smaller sample subsets, so breakdowns (if they exist) are only on one dimension. There is also one relative risk – for discontinuation based on SNAP. These are saved in an RDS file: WhampPrepStats.RDS. The full list of parameters is shown in Parameter description table at the end.
Data Validation Targets – The simulation output Dx data validation targets are a mix of current prevalence (e.g., of awareness, ever PrEP, current PrEP, etc) and durations (e.g., time since first started PrEP, interruption length). These are saved in an RDS file: prepTargets.RDS.
Adherence
WHAMP did not measure adherence, but Darcy’s WHPP study did. Roughly 95% of men were adherent (4+ pills/wk). Those data suggest lower adherence among younger MSM, Hispanic MSM, and MSM in Western WA, but there are not enough observations to draw firm conclusions.
Darcy’s WHPP study
From published paper in AIDS and Behavior
Background: Many state and local health departments now promote and support the use of HIV preexposure prophylaxis (PrEP), yet monitoring use of the intervention at the population level remains challenging.
Methods: We report the results of an online survey designed to measure PrEP use among men who have sex with men (MSM) in Washington State. Data on the proportion of men with indications for PrEP based on state guidelines and levels of awareness, interest, and use of PrEP are presented for 1080 cisgender male respondents who completed the survey between January 1 and February 28, 2017. We conducted bivariate and multivariable logistic regression to identify factors associated with current PrEP use. To examine patterns of discontinuation, we conducted Cox proportional hazards regression and fit a Kaplan-Meier curve to reported data on time on PrEP.
Results: Eighty percent of respondents had heard of PrEP, 19% reported current use, and 36% of men who had never used PrEP wanted to start taking it. Among MSM for whom state guidelines recommend PrEP, 31% were taking it. In multivariable analysis, current PrEP use was associated with older age, higher education, and meeting indications for PrEP use. Our data suggest that 20% of PrEP users discontinue within 12 months, and men with lower educational attainment were more likely to discontinue.
From her dissertation
Prevalence
The prevalence of PrEP use in 2017, based on data from the WHPP sample reweighted to the demographics of HIV-negative/undiagnosed men in Washington, is shown in the table below. For comparison and to inform plausible trends in uptake, we present the prevalence of PrEP in 2018/2019 from the second round of the interent survey, administered from November 2018 to January 2019.
df <- data.frame(Year = c("2017", "2018/19"),
Fraction = c(0.3331, 0.4555))
kable(df, caption = "Proportion of PrEP-eligible men using PrEP, by survey round") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote("Data reweighted to the demographics of HIV-negative/undiagnosed men in WA")| Year | Fraction |
|---|---|
| 2017 | 0.3331 |
| 2018/19 | 0.4555 |
| Note: | |
| Data reweighted to the demographics of HIV-negative/undiagnosed men in WA |
The predicted current usage by age and region is shown below.
Source: Darcy’s analysis of WHPP data
Adherence
The 2017 WHPP survey asked PrEP-users to self-report how many pills they took in the past 30 days.
df <- data.frame(Adherence = c("<2 pills", "2-3 pills", "4+ pills"),
Fraction = c(0.0235, 0.0192, 0.9573))
kable(df, caption = "Adherence over the past 30 days based on self-report") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote("2017 WHPP survey; pills per wk")| Adherence | Fraction |
|---|---|
| <2 pills | 0.0235 |
| 2-3 pills | 0.0192 |
| 4+ pills | 0.9573 |
| Note: | |
| 2017 WHPP survey; pills per wk |
Discontinuation
The plot below shows the Kaplan-Meier curve for all MSM in the original unweighted WHPP sample. The curve seems to be approaching an asymptote at around 70% still on PrEP after 2 yrs – but PrEP has only been widely used for about 2 years, so we don’t have any real info beyond that.
Source: Darcy’s analysis of WHPP data
Data
All of the data used in this report come from the WHAMP Survey. The 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. Final sample size was 927.
- Much of the PrEP analysis, which is restricted to HIV negatives, is based on about 800 cases. Filtering for specific conditional rate estimation, and including the SNAP variable in a model can reduce the sample to smaller sizes.
Note that while the STERGMs augmented the sample with the ARTnet WA cases, we do not do that here. PrEP dynamics are changing rapidly, so the 2017-18 field period of ARTnet was a barrier to inclusion.
# Either make or read in data
dataset <- "readin"
#dataset <- "makeit"
if(dataset == "makeit") {
source(here::here("MakeData", "MM", "Scripts",
"makeWideData.R")) # creates wDF
} else {
wDF <- readRDS(here::here("MakeData", "MM", "Data",
"wideDF.RDS"))
}
# Missing data ----
## assign missing insurance to the uninsured/NA category
## assign missing income to "Missing" category
wDF <- wDF %>%
mutate(insurance = ifelse(is.na(insurance) | insurance == "uninsrd",
"uninsrd_NA", insurance),
income_cat = forcats::fct_explicit_na(income_cat, "Missing"))
## complete case DF for attrs & diag.status -- only 4 cases removed
wDF$complete <- with(wDF, complete.cases(age_group, race, region, diag.status))
# table(wDF$complete)
cDF <- wDF %>% filter(complete == 1)
# Restrict to HIV- for all PrEP analyses
# create some char vars for tables
# and some other target vars
hivnegDF <- cDF %>%
filter(diag.status == 0) %>%
mutate(prepCurrent = ifelse(is.na(prep.current) | prep.current == 0,
"not on prep",
"on prep"),
prepAware = ifelse(prep.aware == 1, "aware", "unaware"),
prep.discont = if_else(prep.stop==1 | prep.int==1, 1, 0),
prep.reinit = case_when(prep.stop==1 ~ 0,
prep.int==1 ~ 1,
TRUE ~ NA_real_))
# Restrict to !is.na(prep.aware) for prep aware analyses
prepAwareDF <- hivnegDF %>% filter(!is.na(prep.aware))
# Restrict to current PrEP users for duration and interruption analyses
prepCurrentDF <- hivnegDF %>% filter(prepCurrent == "on prep")
# Restrict to ever PrEP users for discontinuation (= stop + interrupt)
prepDiscontDF <- hivnegDF %>% filter(prep.ever == 1)
# Restrict to discontinuers for re-initiation analysis
prepReinitDF <- prepDiscontDF %>% filter(prep.discont == 1)Missing data
For the main demographic attributes (including HIV Dx), the missing values are very few, so we use the complete case subset (excludes 4 cases).
Missing values for insurance and income_cat are more common, and can be found in this report. We group the missing insurance with the uninsured, based on exploration of their patterns. We use the missing income cases as a factor level in the analyses, so we don’t lose them. But income and insurance are not used in the models, so the impact is only on the descriptives we show in this report.
Missing values for sexual activity come mostly from attrtion. Attrition reduced the sample size to 771 by the SNAP variable, which was a key predictor for PrEP initiation and stopping rates.
PrEP aware HIV-
Descriptive statistics
Overall
Overall, very few HIV- are not aware of PrEP, but there is some important heterogeneity by race and region.
In the 2017 WHPP survey, 80% of HIV- MSM in WA had heard of PrEP.
# There are missing values in PREP_AWARE that have been removed
# here b/c the summary tables don't play nice with them
#with(hivnegDF, table(PREP_AWARE, useNA = "al")
prep.aware.all <- prepAwareDF %>%
group_by(prepAware) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.aware.all %>% kable(caption= "Distribution of PrEP awareness") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("HIV- and not missing PrEP aware; n = ",
nrow(prepAwareDF)))| prepAware | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| aware | 697 | 691.5 | 0.92 |
| unaware | 59 | 58.8 | 0.08 |
| Note: | |||
| HIV- and not missing PrEP aware; n = 756 |
Age
Table
prep.aware.age <- prepAwareDF %>%
group_by(age_group, prepAware) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(age_group) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.aware.age %>% select(-c(n.obs, n.wtd)) %>%
pivot_wider(names_from = age_group,
values_from = prop.wtd) %>%
kable(caption= "Distribution of PrEP awareness by age") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "Age group" = 5)) %>%
footnote(paste("HIV- not missing aware, wtd proportions; n = ",
nrow(prepAwareDF)))| prepAware | 15-24 | 25-34 | 35-44 | 45-54 | 55-65 |
|---|---|---|---|---|---|
| aware | 0.86 | 0.95 | 0.95 | 0.95 | 0.86 |
| unaware | 0.14 | 0.05 | 0.05 | 0.05 | 0.14 |
| Note: | |||||
| HIV- not missing aware, wtd proportions; n = 756 |
Plot
fig <- ggplot(prep.aware.age,
aes(x = prepAware, 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 = "PrEP awareness",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0)) +
scale_colour_brewer(palette = "YlOrRd")
title.plotly <-
list(text = paste0('PrEP awareness by age group',
'<br>',
'<sup>',
'n = ',
sum(prep.aware.age$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Race
Table
prep.aware.race <- prepAwareDF %>%
group_by(race, prepAware) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(race) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.aware.race %>% select(-c(n.obs, n.wtd)) %>%
pivot_wider(names_from = race,
values_from = prop.wtd) %>%
kable(caption= paste("Distribution of PrEP awareness by race")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "Race" = 3)) %>%
footnote(paste("HIV- not missing aware, wtd proportions; n = ",
sum(prep.aware.race$n.obs)))| prepAware | B | H | O |
|---|---|---|---|
| aware | 0.93 | 0.79 | 0.94 |
| unaware | 0.07 | 0.21 | 0.06 |
| Note: | |||
| HIV- not missing aware, wtd proportions; n = 756 |
Plot
fig <- ggplot(prep.aware.race,
aes(x = prepAware, y = prop.wtd, text = prop.wtd,
group = race, fill = race)) +
geom_bar(stat="identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
labs(x = "PrEP awareness",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('PrEP awareness by race',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.aware.race$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Region
Table
prep.aware.region <- prepAwareDF %>%
group_by(region, prepAware) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(region) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.aware.region %>% select(-c(n.obs, n.wtd)) %>%
pivot_wider(names_from = region,
values_from = prop.wtd) %>%
kable(caption= paste("Distribution of PrEP awareness by region")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "Region" = 3)) %>%
footnote(paste("HIV- not missing aware, wtd proportions; n = ",
sum(prep.aware.region$n.obs)))| prepAware | EasternWA | King | WesternWA |
|---|---|---|---|
| aware | 0.79 | 0.95 | 0.91 |
| unaware | 0.21 | 0.05 | 0.09 |
| Note: | |||
| HIV- not missing aware, wtd proportions; n = 756 |
Plot
fig <- ggplot(prep.aware.region,
aes(x = prepAware, y = prop.wtd, text = prop.wtd,
group = region, fill = region)) +
geom_bar(stat="identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
labs(x = "PrEP awareness",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('PrEP awareness by region',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.aware.region$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Insurance
Table
prep.aware.ins <- prepAwareDF %>%
group_by(insurance, prepAware) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(insurance) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.aware.ins %>% select(-c(n.obs, n.wtd)) %>%
pivot_wider(names_from = insurance,
values_from = prop.wtd) %>%
kable(caption= paste("Distribution of PrEP awareness by insurance")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "Insurance" = 5)) %>%
footnote(paste("HIV- not missing aware, wtd proportions; n = ",
sum(prep.aware.ins$n.obs)))| prepAware | emplyr | indvdl_othr | medicare_caid | othr_gov | uninsrd_NA |
|---|---|---|---|---|---|
| aware | 0.95 | 0.99 | 0.86 | 0.87 | 0.75 |
| unaware | 0.05 | 0.01 | 0.14 | 0.13 | 0.25 |
| Note: | |||||
| HIV- not missing aware, wtd proportions; n = 756 |
Plot
fig <- ggplot(prep.aware.ins,
aes(x = prepAware, y = prop.wtd, text = prop.wtd,
group = insurance, fill = insurance)) +
geom_bar(stat="identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
labs(x = "PrEP awareness",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('PrEP awareness by insurance',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.aware.ins$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Income
Table
prep.aware.income <- prepAwareDF %>%
group_by(income_cat, prepAware) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(income_cat) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.aware.income %>% select(-c(n.obs, n.wtd)) %>%
pivot_wider(names_from = income_cat,
values_from = prop.wtd) %>%
kable(caption= paste("Distribution of PrEP awareness by FPL Income Category")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "Income" = 7)) %>%
footnote(paste("HIV- not missing aware, wtd proportions; n = ",
sum(prep.aware.income$n.obs)))| prepAware | fpl138 | fpl138-300 | fpl300-400 | fpl400-500 | fpl500-600 | fpl600+ | Missing |
|---|---|---|---|---|---|---|---|
| aware | 0.76 | 0.92 | 0.95 | 0.94 | 0.96 | 0.96 | 0.86 |
| unaware | 0.24 | 0.08 | 0.05 | 0.06 | 0.04 | 0.04 | 0.14 |
| Note: | |||||||
| HIV- not missing aware, wtd proportions; n = 756 |
Plot
fig <- ggplot(prep.aware.income,
aes(x = prepAware, 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 = "PrEP awareness",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('PrEP awareness by FPL income category',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.aware.income$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Prediction
Logistic regression
Model 1
library(survey)
awarePrevDF <- prepAwareDF %>%
mutate(age.grp = factor(age.grp),
race = relevel(factor(race), ref = "O"),
region = relevel(factor(region), ref = "King"),
income_cat = relevel(factor(income_cat), ref = "fpl600+"),
insurance = relevel(factor(insurance), ref = "emplyr"))
design <- svydesign(id = ~ego.id,
data = awarePrevDF,
weights = ~ego.wawt)
logitm1 <- svyglm(prep.aware ~ age.grp,
data = awarePrevDF,
family = quasibinomial(),
design)
kable(tidy(logitm1, exponentiate = T),
caption= "Model 1", digits = 2) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 6.27 | 0.29 | 6.28 | 0.00 |
| age.grp2 | 3.08 | 0.43 | 2.63 | 0.01 |
| age.grp3 | 3.30 | 0.52 | 2.28 | 0.02 |
| age.grp4 | 2.87 | 0.46 | 2.31 | 0.02 |
| age.grp5 | 1.01 | 0.42 | 0.01 | 0.99 |
Model 2
logitm2 <- svyglm(prep.aware ~ race,
data = awarePrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm2, exponentiate = T),
caption= "Model 2", digits = 2) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 15.40 | 0.17 | 15.84 | 0.00 |
| raceB | 0.83 | 0.71 | -0.26 | 0.79 |
| raceH | 0.25 | 0.36 | -3.89 | 0.00 |
Model 3
logitm3 <- svyglm(prep.aware ~ region,
data = awarePrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm3, exponentiate = T),
caption= "Model 3", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 21.09 | 0.28 | 10.84 | 0.0000 |
| regionEasternWA | 0.18 | 0.41 | -4.20 | 0.0000 |
| regionWesternWA | 0.46 | 0.35 | -2.23 | 0.0264 |
Model 4
logitm4 <- svyglm(prep.aware ~ insurance,
data = awarePrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm4, exponentiate = T),
caption= "Model 4", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 19.52 | 0.23 | 13.01 | 0.0000 |
| insuranceindvdl_othr | 3.60 | 1.04 | 1.24 | 0.2171 |
| insurancemedicare_caid | 0.32 | 0.38 | -3.00 | 0.0028 |
| insuranceothr_gov | 0.36 | 0.57 | -1.82 | 0.0689 |
| insuranceuninsrd_NA | 0.15 | 0.41 | -4.61 | 0.0000 |
Model 5
logitm5 <- svyglm(prep.aware ~ race + region,
data = awarePrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm5, exponentiate = T),
caption= "Model 5", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 26.48 | 0.32 | 10.26 | 0.0000 |
| raceB | 0.75 | 0.71 | -0.41 | 0.6852 |
| raceH | 0.28 | 0.34 | -3.73 | 0.0002 |
| regionEasternWA | 0.20 | 0.39 | -4.06 | 0.0001 |
| regionWesternWA | 0.46 | 0.36 | -2.21 | 0.0277 |
Model 6
logitm6 <- svyglm(prep.aware ~ race + region + insurance,
data = awarePrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm6, exponentiate = T),
caption= "Model 6", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 33.04 | 0.34 | 10.41 | 0.0000 |
| raceB | 0.86 | 0.68 | -0.21 | 0.8304 |
| raceH | 0.37 | 0.40 | -2.47 | 0.0139 |
| regionEasternWA | 0.24 | 0.45 | -3.21 | 0.0014 |
| regionWesternWA | 0.54 | 0.38 | -1.64 | 0.1021 |
| insuranceindvdl_othr | 3.86 | 1.02 | 1.33 | 0.1855 |
| insurancemedicare_caid | 0.41 | 0.40 | -2.20 | 0.0284 |
| insuranceothr_gov | 0.54 | 0.71 | -0.87 | 0.3821 |
| insuranceuninsrd_NA | 0.23 | 0.42 | -3.53 | 0.0004 |
Comparison of AICs
models <- get_formula("logitm", 6)
aics <- get_aics("logitm", 6, "logit")
kable(data.frame(Model = models,
AIC = aics),
caption= "AIC comparison", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| Model | AIC |
|---|---|
| prep.aware ~ age.grp | 407.7 |
| prep.aware ~ race | 403.7 |
| prep.aware ~ region | 398.5 |
| prep.aware ~ insurance | 391.6 |
| prep.aware ~ race + region | 389.5 |
| prep.aware ~ race + region + insurance | 379.8 |
- Insurance is a strong predictor, but we can’t simulate with it, so we will use model 5 with race + region.
Predicted distribution
#### Predict
expand_df <- expand.grid(
race = levels(awarePrevDF$race),
region = levels(awarePrevDF$region))
pred_df <- cbind(expand_df,
prob = predict(logitm5,
newdata = expand_df,
type = "response",
npop = 1000))
pred.prep.aware <- pred_df %>%
select(-prob.SE) %>%
rename(prob = prob.response)%>%
data.table::setDT()
data.table::setkeyv(pred.prep.aware,
c("race", "region"))
## Save prediction
save_out$prep.aware.dt <- pred.prep.aware
## Plot predicted values
fig <- ggplot(data = pred.prep.aware,
aes(x = region, y = prob, text = round(prob,2),
group = race, fill = race)) +
geom_bar(stat = "identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
theme(plot.margin = margin(t=30,r=20,b=20,l=20),
legend.position = "bottom")
title.plotly <-
list(text = paste0('Predicted PrEP Awareness by race and region',
'<br>',
'<sup>',
'n = ',
dim(awarePrevDF)[[1]],
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(title = title.plotly)PrEP indicated
This is used as a validation target, and the overall level is used as an input parameter to restrict PrEP initiation to persons who meet this criterium.
The CDC and WADOH guidelines are a bit different, and the WADOH guidelines target more specifically:
We can not implement the full set of current WADOH guidelines because we are not simulating STI spread. So we use what we can, which corresponds more to the “Discussion” indication than the “Recommend”. See the GitHub discussion, and this comment in particular.
We implement 3 criteria:
prep.indic.hiv.part: HIV+ partner, regardless of VL statusprep.indic.cai: CAI outside of a mutually monogamous LT partnership- (ego snap > 1 or pconc >0) and at least one partner CAI OR
- 1 CAI partner but relationship < 12 mo (“not long term”)
prep.indic.snap2: SNAP > 1
And we combine these into a single indicator, prep.indic.any, which classifies someone as PrEP indicated if any of the 3 criteria above are met.
Indicator comparisons
HIV+ partner(s)
tDF <- hivnegDF %>%
group_by(prep.indic.hiv.part) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
tDF %>%
kable(caption= "PrEP indicated: Ongoing HIV+ partner(s)") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("HIV- : n =", nrow(hivnegDF)))| prep.indic.hiv.part | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 591 | 587.8 | 0.71 |
| 1 | 55 | 54.4 | 0.07 |
| NA | 185 | 186.6 | 0.23 |
| Note: | |||
| HIV- : n = 831 |
CAI
Outside of a mutually monogamous LT relationship
tDF <- hivnegDF %>%
group_by(prep.indic.cai) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
tDF %>%
kable(caption= "PrEP indicated: CAI") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("HIV- : n =", nrow(hivnegDF)))| prep.indic.cai | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 277 | 281.2 | 0.34 |
| 1 | 364 | 354.9 | 0.43 |
| NA | 190 | 192.6 | 0.23 |
| Note: | |||
| HIV- : n = 831 |
SNAP > 1
tDF <- hivnegDF %>%
group_by(prep.indic.snap2) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
tDF %>%
kable(caption= "PrEP indicated: SNAP > 1") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("HIV- : n =", nrow(hivnegDF)))| prep.indic.snap2 | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 314 | 308.0 | 0.37 |
| 1 | 332 | 334.2 | 0.40 |
| NA | 185 | 186.6 | 0.23 |
| Note: | |||
| HIV- : n = 831 |
Overlap
Here we filter out the cases missing SNAP (attrition + no anal sex since 1980).
tDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
group_by(prep.indic.snap2, prep.indic.hiv.part, prep.indic.cai) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
ungroup() %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
tDF %>% select(-n.wtd) %>%
kable(caption= "PrEP indicated",
digits = c(0,0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
pack_rows("Not indicated", 1, 1) %>%
pack_rows("SNAP < 2", 2, 5) %>%
pack_rows("SNAP = 2+", 6, 9) %>%
footnote(paste("HIV- n =", nrow(hivnegDF),
"; valid obs n =", sum(tDF$n.obs),
"; missing SNAP n =", table(is.na(hivnegDF$snap))[[2]]))| prep.indic.snap2 | prep.indic.hiv.part | prep.indic.cai | n.obs | prop.wtd |
|---|---|---|---|---|
| Not indicated | ||||
| 0 | 0 | 0 | 233 | 0.35 |
| SNAP < 2 | ||||
| 0 | 0 | 1 | 69 | 0.11 |
| 0 | 0 | NA | 5 | 0.01 |
| 0 | 1 | 0 | 5 | 0.01 |
| 0 | 1 | 1 | 2 | 0.00 |
| SNAP = 2+ | ||||
| 1 | 0 | 0 | 36 | 0.07 |
| 1 | 0 | 1 | 248 | 0.37 |
| 1 | 1 | 0 | 3 | 0.00 |
| 1 | 1 | 1 | 45 | 0.07 |
| Note: | ||||
| HIV- n = 831 ; valid obs n = 646 ; missing SNAP n = 185 | ||||
Indication by use
Use among non-indicated
If we want to use one or more of the PrEP indicators to restrict the population that can initiate PrEP, the non-indicated group must not be currently using PrEP.
As the table below shows, this is not the case in the WHAMP survey. Even if we take the broadest definition of indicated – indicated by at least one of our 3 criteria – a small fraction of the non-indicated are currently using PrEP.
However, the fraction of non-indicated on PrEP is quite low, so we can probably ignore it without much impact: those respondents are at very low risk, so they won’t contribute to transmission, and they will not be eligible for PDAP, so won’t contribute to costs.
sjPlot::tab_xtab(var.row = hivnegDF$prep.indic.any,
var.col = hivnegDF$prep.current,
weight.by = hivnegDF$ego.wawt,
title = "PrEP current use by indication",
show.row.prc = T,
show.col.prc = T,
show.summary = F)| prep.indic.any | prep.current | Total | |
|---|---|---|---|
| 0 | 1 | ||
| 0 |
219 96.9 % 46 % |
7 3.1 % 4.4 % |
226 100 % 35.6 % |
| 1 |
257 62.8 % 54 % |
152 37.2 % 95.6 % |
409 100 % 64.4 % |
| Total |
476 75 % 100 % |
159 25 % 100 % |
635 100 % 100 % |
Non use among indicated
The indication criteria are a mix of individual level behavior, like SNAP2, and network related risk, like the HIV status of a partner, and CAI outside of a mutually monogamous long-term partnership.
The table below shows that current PrEP use is much higher for those indicated by the individual level criteria (their own behavior) than for those indicated by network related risks. This strongly suggests that the decision to start PrEP is influenced more by a person’s assessment of their own risk behavior than their indirect risks.
aaa <- hivnegDF %>%
filter(!is.na(snap) & !is.na(prep.current)) %>%
mutate(snap2 = ifelse(snap>1, 2, snap)) %>%
filter(prep.indic.any==1)
sjPlot::tab_xtab(var.row = aaa$snap2,
var.col = aaa$prep.current,
weight.by = aaa$ego.wawt,
title = "For PrEP indicated: current use by SNAP2",
show.row.prc = T,
show.summary = F)| snap2 | prep.current | Total | |
|---|---|---|---|
| 0 | 1 | ||
| 1 |
68 88.3 % |
9 11.7 % |
77 100 % |
| 2 |
190 56.9 % |
144 43.1 % |
334 100 % |
| Total |
258 62.8 % |
153 37.2 % |
411 100 % |
Indication by attributes
These plots based on prep.indic.any.
Age
tDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
group_by(prep.indic.any, age_group) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(age_group) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
fig <- ggplot(tDF %>% filter(prep.indic.any==1),
aes(x = age_group, y = prop.wtd, text = prop.wtd)) +
geom_bar(stat="identity", position = "dodge",
fill = "blue", alpha = 0.5) +
labs(x = "PrEP indicated",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0)) +
scale_fill_brewer()
title.plotly <-
list(text = paste0('PrEP indicated by age: Any criteria',
'<br>',
'<sup>',
'n = ',
sum(tDF$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Race
tDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
group_by(prep.indic.any, race) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(race) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
fig <- ggplot(tDF %>% filter(prep.indic.any==1),
aes(x = race, y = prop.wtd, text = prop.wtd)) +
geom_bar(stat="identity", position = "dodge",
fill = "blue", alpha = 0.5) +
labs(x = "PrEP indicated",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0)) +
scale_fill_brewer()
title.plotly <-
list(text = paste0('PrEP indicated by race: Any criteria',
'<br>',
'<sup>',
'n = ',
sum(tDF$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Region
tDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
group_by(prep.indic.any, region) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(region) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
fig <- ggplot(tDF %>% filter(prep.indic.any==1),
aes(x = region, y = prop.wtd, text = prop.wtd)) +
geom_bar(stat="identity", position = "dodge",
fill = "blue", alpha = 0.5) +
labs(x = "PrEP indicated",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0)) +
scale_fill_brewer()
title.plotly <-
list(text = paste0('PrEP indicated by region: Any criteria',
'<br>',
'<sup>',
'n = ',
sum(tDF$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)SNAP2
tDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
mutate(snap2 = ifelse(snap >1, 2, snap)) %>%
group_by(prep.indic.any, snap2) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(snap2) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
fig <- ggplot(tDF %>% filter(prep.indic.any==1),
aes(x = snap2, y = prop.wtd, text = prop.wtd)) +
geom_bar(stat="identity", position = "dodge",
fill = "blue", alpha = 0.5) +
labs(x = "PrEP indicated",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0)) +
scale_fill_brewer()
title.plotly <-
list(text = paste0('PrEP indicated by SNAP2: Any criteria',
'<br>',
'<sup>',
'n = ',
sum(tDF$n.obs),
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)PrEP ever
This is used as a validation target.
PrEP ever - overall. Though the targets object will breakdown by attribute, the numbers are probably too small to be reliable.
PrEP ever x current. This is used to validate our intiation-discontinue-reinitiation dynamics.
NAs come from attrition
PrEP ever
hivnegDF %>%
group_by(prep.ever) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2)) %>%
kable(caption= "Distribution of ever PrEP use (HIV-)") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("HIV- : n =", nrow(hivnegDF)))| prep.ever | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 533 | 524.1 | 0.63 |
| 1 | 222 | 225.2 | 0.27 |
| NA | 76 | 79.5 | 0.10 |
| Note: | |||
| HIV- : n = 831 |
PrEP ever x current
sjPlot::tab_xtab(var.row = hivnegDF$prep.ever,
var.col = hivnegDF$prep.current,
weight.by = hivnegDF$ego.wawt,
title = "PrEP Ever by Current Use",
show.cell.prc = T,
show.summary = F)| prep.ever | prep.current | Total | |
|---|---|---|---|
| 0 | 1 | ||
| 0 |
524 70 % |
0 0 % |
524 70 % |
| 1 |
58 7.7 % |
167 22.3 % |
225 30 % |
| Total |
582 77.7 % |
167 22.3 % |
749 100 % |
Prevalence: Current PrEP use
This variable is used for two purposes:
Estimation of first time PrEP initiation rates (Incidence = Prevalence/Duration). Model based, restricted to current PrEP and never interrupted. Interruption identifies the sample of cyclers.
Validation target – Descriptives, for this purpose it only conditions on HIV- status. Comparison will also be to all HIV-
Variable:
prepCurrentfrom the originalprepCurrentIn the WHPP survey, the proportion of HIV- MSM who were current PrEP users was 19% in 2017, and among PrEP eligible HIV- was 33% in 2017 and 45% in 2018-19 (weighted to WA HIV- MSM).
Descriptive statistics
- Of the 227 respondents who report ever taking PrEP, r scales::percent(with(cDF, proportions(table(prep.ever, prepCurrent), 1))[[5]]) are currently on PrEP.
Overall
# all HIV-
# prepCurrent NA is a mix of not aware + never prep; we recode them to 0
#with(wDF, table(diag.status, prepCurrent, useNA ="al"))
prep.curr_all <- hivnegDF %>%
group_by(prepCurrent) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
mutate(prop.wtd = round(n.wtd / sum(n.wtd), 2))
prep.curr_all %>% kable(caption= "Distribution of current PrEP use (HIV-)") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("HIV- : n =", nrow(hivnegDF)))| prepCurrent | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| not on prep | 670 | 661.3 | 0.8 |
| on prep | 161 | 167.5 | 0.2 |
| Note: | |||
| HIV- : n = 831 |
Age
Table
prep.curr.age <- hivnegDF %>%
select(prepCurrent, age_group, ego.wawt) %>%
group_by(age_group, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(age_group) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, age_group, n.obs, n.wtd, prop.wtd)
prep.curr.age %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by Age",
col.names = c("Age", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing age_group; n = ",
table(is.na(hivnegDF$age_group))[[1]]))| Age | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| 15-24 | 124 | 10 | 0.06 |
| 25-34 | 251 | 52 | 0.23 |
| 35-44 | 166 | 45 | 0.26 |
| 45-54 | 155 | 41 | 0.28 |
| 55-65 | 135 | 19 | 0.15 |
| Note: | |||
| HIV- and not missing age_group; n = 831 |
Plot
fig <- ggplot(prep.curr.age,
aes(x = prepCurrent, 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 = "Current PrEP use",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0)) +
scale_fill_brewer()
title.plotly <-
list(text = paste0('Current PrEP use by age',
'<br>',
'<sup>',
'n = ',
sum(prep.curr.age$n.obs)/2,
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Race
Table
prep.curr.race <- hivnegDF %>%
select(prepCurrent, race, ego.wawt) %>%
group_by(race, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(race) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, race, n.obs, n.wtd, prop.wtd)
prep.curr.race %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by Race",
col.names = c("Race", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing race; n = ",
table(is.na(hivnegDF$snap))[[1]]))| Race | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| B | 27 | 17 | 0.31 |
| H | 89 | 11 | 0.12 |
| O | 715 | 139 | 0.20 |
| Note: | |||
| HIV- and not missing race; n = 646 |
Plot
fig <- ggplot(prep.curr.race,
aes(x = prepCurrent, y = prop.wtd, text = prop.wtd,
group = race, fill = race)) +
geom_bar(stat="identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
labs(x = "Current PrEP use",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('Current PrEP use by race',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.curr.race$n.obs)/2,
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Region
Table
prep.curr.region <- hivnegDF %>%
select(prepCurrent, region, ego.wawt) %>%
group_by(region, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(region) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, region, n.obs, n.wtd, prop.wtd)
prep.curr.region %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by Region",
col.names = c("Region", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing region; n = ",
table(is.na(hivnegDF$region))[[1]]))| Region | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| EasternWA | 120 | 7 | 0.08 |
| King | 376 | 117 | 0.25 |
| WesternWA | 335 | 44 | 0.16 |
| Note: | |||
| HIV- and not missing region; n = 831 |
Plot
fig <- ggplot(prep.curr.region,
aes(x = prepCurrent, y = prop.wtd, text = prop.wtd,
group = region, fill = region)) +
geom_bar(stat="identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
labs(x = "Current PrEP use",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('Current PrEP use by region',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.curr.region$n.obs)/2,
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Insurance
Table
prep.curr_ins <- hivnegDF %>%
select(prepCurrent, insurance, ego.wawt) %>%
group_by(insurance, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(insurance) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, insurance, n.obs, n.wtd, prop.wtd)
prep.curr_ins %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by Insurance",
col.names = c("Insurance", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing insurance; n = ",
table(is.na(hivnegDF$insurance))[[1]]))| Insurance | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| emplyr | 522 | 128 | 0.24 |
| indvdl_othr | 56 | 13 | 0.26 |
| medicare_caid | 113 | 14 | 0.14 |
| othr_gov | 43 | 9 | 0.21 |
| uninsrd_NA | 97 | 3 | 0.03 |
| Note: | |||
| HIV- and not missing insurance; n = 831 |
Plot
fig <- ggplot(prep.curr_ins,
aes(x = prepCurrent, y = prop.wtd, text = prop.wtd,
group = insurance, fill = insurance)) +
geom_bar(stat="identity", position = "dodge", color = "darkgrey",
alpha = 0.7, ) +
labs(x = "Current PrEP use",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('Current PrEP use by insurance',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.curr_ins$n.obs)/2,
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)Income
Table
prep.curr_income <- hivnegDF %>%
select(prepCurrent, income_cat, ego.wawt) %>%
group_by(income_cat, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(income_cat) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, income_cat, n.obs, n.wtd, prop.wtd)
prep.curr_income %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by FPL Income",
col.names = c("income_cat", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing income; n = ",
table(hivnegDF$income_cat=="Missing")[[1]]))| income_cat | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| fpl138 | 93 | 5 | 0.06 |
| fpl138-300 | 168 | 22 | 0.13 |
| fpl300-400 | 73 | 20 | 0.29 |
| fpl400-500 | 51 | 11 | 0.23 |
| fpl500-600 | 26 | 10 | 0.42 |
| fpl600+ | 341 | 97 | 0.27 |
| Missing | 79 | 3 | 0.04 |
| Note: | |||
| HIV- and not missing income; n = 752 |
Plot
fig <- ggplot(prep.curr_income,
aes(x = prepCurrent, 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 = "Current PrEP use",
y = "wtd. proportion") +
theme(plot.margin = margin(t=20,r=5,b=5,l=0))
title.plotly <-
list(text = paste0('Current PrEP use by FPL income category',
'<br>',
'<sup>',
'HIV-; n = ',
sum(prep.curr_income$n.obs)/2,
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(#hoverlabel = list(font=list(color="black")),
title = title.plotly)SNAP
Table (20)
tempDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
select(prepCurrent, snap, ego.wawt) %>%
mutate(snap20 = pmin(snap, 20)) %>%
group_by(snap20, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(snap20) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, snap20, n.obs, n.wtd, prop.wtd)
tempDF %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by SNAP",
col.names = c("SNAP", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing SNAP; n = ",
table(is.na(hivnegDF$snap))[[1]],
"; snap topcoded at 20+"))| SNAP | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| 0 | 104 | 3 | 0.03 |
| 1 | 210 | 13 | 0.06 |
| 2 | 62 | 9 | 0.14 |
| 3 | 46 | 22 | 0.45 |
| 4 | 37 | 14 | 0.41 |
| 5 | 28 | 8 | 0.26 |
| 6 | 20 | 11 | 0.55 |
| 7 | 6 | 1 | 0.20 |
| 8 | 16 | 6 | 0.38 |
| 9 | 5 | 4 | 0.80 |
| 10 | 23 | 8 | 0.38 |
| 12 | 12 | 4 | 0.36 |
| 14 | 1 | 1 | 1.00 |
| 15 | 16 | 15 | 0.94 |
| 16 | 1 | 1 | 1.00 |
| 18 | 2 | 1 | 0.33 |
| 20 | 54 | 39 | 0.68 |
| Note: | |||
| HIV- and not missing SNAP; n = 646 ; snap topcoded at 20+ |
Plot (20)
tempDF %>%
filter(prepCurrent == "on prep") %>%
ggplot(aes(x=snap20, y=prop.wtd, weight=n.obs)) +
ylim(c(0, 1)) +
labs(title = "Current PrEP use by SNAP",
x = "Number of anal sex partners last year",
y = "Proportion on PrEP") +
geom_line() +
geom_point(aes(size=n.obs), color = "blue", alpha = 0.5)+
geom_smooth(span = 1.2)
Table (5)
Our constraint here come from the fact that we will be using the prevalence as the numerator over duration to construct incidence. Duration is restaricted to on PrEP, only 161 cases, so we don’t have enough data to estimate the effect of the tail of SNAP.
prep.curr.snap5 <- hivnegDF %>%
filter(!is.na(snap)) %>%
select(prepCurrent, snap, ego.wawt) %>%
mutate(snap5 = pmin(snap, 5)) %>%
group_by(snap5, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(snap5) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, snap5, n.obs, n.wtd, prop.wtd)
prep.curr.snap5 %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by SNAP",
col.names = c("SNAP", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing SNAP; n = ",
table(is.na(hivnegDF$snap))[[1]],
"; snap topcoded at 5+"))| SNAP | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| 0 | 104 | 3 | 0.03 |
| 1 | 210 | 13 | 0.06 |
| 2 | 62 | 9 | 0.14 |
| 3 | 46 | 22 | 0.45 |
| 4 | 37 | 14 | 0.41 |
| 5 | 187 | 98 | 0.52 |
| Note: | |||
| HIV- and not missing SNAP; n = 646 ; snap topcoded at 5+ |
Plot (5)
prep.curr.snap5 %>%
filter(prepCurrent == "on prep") %>%
ggplot(aes(x=snap5, y=prop.wtd, weight=n.obs)) +
ylim(c(0, 1)) +
labs(title = "Current PrEP use by SNAP",
x = "Number of anal sex partners last year",
y = "Proportion on PrEP") +
geom_line() +
geom_point(aes(size=n.obs), color = "blue", alpha = 0.5)+
geom_smooth(span = 1.2)
Deg.tot
Table
tempDF <- hivnegDF %>%
filter(!is.na(deg.tot)) %>%
select(prepCurrent, deg.tot, ego.wawt) %>%
mutate(deg.tot3 = pmin(deg.tot, 3)) %>%
group_by(deg.tot3, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(deg.tot3) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, deg.tot3, n.obs, n.wtd, prop.wtd)
tempDF %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by Deg Tot (M+C)",
col.names = c("deg.tot", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "", "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing deg.tot; n = ",
table(is.na(hivnegDF$deg.tot))[[1]],
"; deg.tot topcoded at 3+"))| deg.tot | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| 0 | 266 | 24 | 0.10 |
| 1 | 272 | 46 | 0.17 |
| 2 | 74 | 37 | 0.51 |
| 3 | 82 | 52 | 0.57 |
| Note: | |||
| HIV- and not missing deg.tot; n = 694 ; deg.tot topcoded at 3+ |
Plot
tempDF %>%
filter(prepCurrent == "on prep") %>%
ggplot(aes(x=deg.tot3, y=prop.wtd, weight=n.obs)) +
ylim(c(0, 1)) +
labs(title = "Current PrEP use by total degree",
caption = paste("HIV- proportions; n = ",
sum(tempDF$n.wtd),
"; deg.tot topcoded at 3+"),
x = "Degree (Main + Casl)",
y = "Proportion on PrEP") +
geom_line() +
geom_point(aes(size=n.obs), color = "blue", alpha = 0.5)
Deg tot x SNAP interaction
Bottom line: SNAP is the driver, not active degree, so we need to use SNAP for the predictive model.
tempDF <- hivnegDF %>%
filter(!is.na(deg.tot)) %>%
select(prepCurrent, snap, deg.tot, ego.wawt) %>%
mutate(deg.tot3 = pmin(deg.tot, 3),
snap5 = pmin(snap, 5)) %>%
group_by(deg.tot3, snap5, prepCurrent) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(deg.tot3, snap5) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2)) %>%
select(prepCurrent, deg.tot3, snap5, n.obs, n.wtd, prop.wtd)Plot 1
tempDF %>%
filter(prepCurrent == "on prep") %>%
ggplot(aes(x=deg.tot3, y=prop.wtd, group=snap5, fill=snap5)) +
geom_bar(stat="identity", position = "dodge") +
labs(title = "PrEP use by SNAP, at each level of active degree",
x = "Currently active partnerships (active degree)",
y = "Wtd proportion on PrEP",
caption = paste0("HIV- and not missing deg.tot: n = ",
table(is.na(hivnegDF$deg.tot))[[1]],
"; Degree topcoded at 3+"))
Plot 2
tempDF <- hivnegDF %>%
select(snap, deg.tot, ego.wawt, prepCurrent) %>%
filter(!is.na(snap) & !is.na(deg.tot)) %>%
mutate(snap5 = pmin(snap, 5),
conc = ifelse(deg.tot < 2, 0, 1)) %>%
group_by(prepCurrent, snap5, conc) %>%
dplyr::summarize(n.obs = n(),
n.wtd = round(sum(ego.wawt))) %>%
group_by(snap5, conc) %>%
mutate(n.obs = sum(n.obs),
prop.wtd = round(n.wtd/sum(n.wtd), 2))
tempDF %>%
filter(prepCurrent == "on prep") %>%
ggplot(aes(x=snap5, y=prop.wtd,
group = factor(conc), color = factor(conc),
weight=n.obs)) +
ylim(c(0, 1)) +
labs(title = "Current PrEP use by SNAP & Concurrent Persistent P",
x = "Number of anal sex partners last year",
y = "Proportion on PrEP") +
geom_line() +
geom_point(aes(size=n.obs), alpha = 0.5)
Table
tempDF %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "Current PrEP use use by Deg Tot (M+C)",) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("Group" = 3, "On PrEP" = 2)) %>%
footnote(paste("HIV- and not missing deg.tot; n = ",
table(is.na(hivnegDF$deg.tot))[[1]],
"; deg.tot topcoded at 3+"))| snap5 | conc | n.obs | n.wtd | prop.wtd |
|---|---|---|---|---|
| 0 | 0 | 104 | 3 | 0.03 |
| 1 | 0 | 210 | 13 | 0.06 |
| 2 | 0 | 50 | 6 | 0.12 |
| 2 | 1 | 12 | 3 | 0.25 |
| 3 | 0 | 38 | 12 | 0.35 |
| 3 | 1 | 8 | 11 | 0.69 |
| 4 | 0 | 21 | 9 | 0.45 |
| 4 | 1 | 16 | 6 | 0.38 |
| 5 | 0 | 67 | 28 | 0.43 |
| 5 | 1 | 120 | 70 | 0.57 |
| Note: | ||||
| HIV- and not missing deg.tot; n = 694 ; deg.tot topcoded at 3+ |
PrEP indicated
Based on any indication.
Table
tempDF <- hivnegDF %>%
filter(!is.na(snap)) %>%
group_by(prep.indic.any, prepCurrent) %>%
summarise(n.obs = n(),
n.wtd = round(sum(ego.wawt), 1)) %>%
group_by(prep.indic.any) %>%
mutate(prop.wtd = round(n.wtd/sum(n.wtd), 2))
tempDF %>%
filter(prepCurrent == "on prep") %>%
select(-prepCurrent) %>%
kable(caption= "PrEP indication for Current Users",
col.names = c("Indicated", "n.obs", "Wtd N", "Wtd Prop")) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "On PrEP" = 3)) %>%
footnote(paste("HIV- and not missing SNAP; n = ",
sum(tempDF$n.obs)))| Indicated | n.obs | Wtd N | Wtd Prop |
|---|---|---|---|
| 0 | 7 | 6.7 | 0.03 |
| 1 | 145 | 152.3 | 0.37 |
| Note: | |||
| HIV- and not missing SNAP; n = 646 |
Plot
tempDF %>%
filter(prepCurrent == "on prep") %>%
ggplot(aes(x=prep.indic.any, y=prop.wtd, weight=n.obs)) +
ylim(c(0, 1)) +
labs(title = "Current PrEP use by PrEP Indication",
x = "Indicated",
y = "Proportion on PrEP") +
geom_bar(stat="identity", fill = "blue", alpha = 0.5)
Prediction
For this purpose we need to restrict to PrEP aware and never discontinued. Discontinuation is subject to different rates of re-initiation, and this variable will instead be used to estimate first initiation rates. We also need to filter out NA on SNAP.
Logistic regression
Model 1
# create filter var (tmp.interrupt) for prep interrupters
prepPrevDF <- hivnegDF %>%
mutate(tmp.interrupt = if_else(is.na(prep.int) | prep.int==0, 0, 1)) %>%
filter(prep.aware == 1 & tmp.interrupt == 0 & !is.na(snap)) %>%
mutate(age.grp = factor(age.grp),
race = relevel(factor(race), ref = "O"),
region = relevel(factor(region), ref = "King"),
income_cat = relevel(factor(income_cat), ref = "fpl600+"),
insurance = relevel(factor(insurance), ref = "emplyr"),
snap5 = pmin(snap, 5))
design <- svydesign(id = ~ego.id,
data = prepPrevDF,
weights = ~ego.wawt)
logitm1 <- svyglm(prep.current ~ age.grp,
data = prepPrevDF,
family = quasibinomial(),
design)
kable(tidy(logitm1, exponentiate = T),
caption= "Model 1", digits = 2) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.05 | 0.61 | -4.83 | 0.00 |
| age.grp2 | 5.69 | 0.63 | 2.74 | 0.01 |
| age.grp3 | 7.99 | 0.68 | 3.08 | 0.00 |
| age.grp4 | 8.50 | 0.65 | 3.29 | 0.00 |
| age.grp5 | 4.41 | 0.68 | 2.19 | 0.03 |
Model 2
logitm2 <- svyglm(prep.current ~ race,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm2, exponentiate = T),
caption= "Model 2", digits = 2) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.29 | 0.11 | -11.17 | 0.00 |
| raceB | 3.01 | 0.71 | 1.55 | 0.12 |
| raceH | 0.43 | 0.43 | -1.98 | 0.05 |
Model 3
logitm3 <- svyglm(prep.current ~ region,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm3, exponentiate = T),
caption= "Model 3", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.40 | 0.16 | -5.86 | 0.0000 |
| regionEasternWA | 0.20 | 0.47 | -3.46 | 0.0006 |
| regionWesternWA | 0.47 | 0.24 | -3.08 | 0.0022 |
Model 4
logitm4 <- svyglm(prep.current ~ insurance,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm4, exponentiate = T),
caption= "Model 4", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.35 | 0.14 | -7.49 | 0.0000 |
| insuranceindvdl_othr | 1.12 | 0.40 | 0.29 | 0.7755 |
| insurancemedicare_caid | 0.39 | 0.41 | -2.26 | 0.0243 |
| insuranceothr_gov | 0.41 | 0.64 | -1.37 | 0.1706 |
| insuranceuninsrd_NA | 0.00 | 0.27 | -57.69 | 0.0000 |
Model 5
logitm5 <- svyglm(prep.current ~ income_cat,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm5, exponentiate = T),
caption= "Model 5", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.43 | 0.17 | -4.98 | 0.0000 |
| income_catfpl138 | 0.07 | 0.75 | -3.56 | 0.0004 |
| income_catfpl138-300 | 0.33 | 0.34 | -3.30 | 0.0010 |
| income_catfpl300-400 | 0.94 | 0.38 | -0.17 | 0.8667 |
| income_catfpl400-500 | 0.58 | 0.47 | -1.17 | 0.2428 |
| income_catfpl500-600 | 1.68 | 0.49 | 1.05 | 0.2936 |
| income_catMissing | 0.00 | 0.27 | -58.02 | 0.0000 |
Model 6
logitm6 <- svyglm(prep.current ~ poly(snap5, 2),
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm6, exponentiate = T),
caption= "Model 6", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1.500000e-01 | 0.20 | -9.77 | 0.0000 |
| poly(snap5, 2)1 | 1.002303e+17 | 5.01 | 7.81 | 0.0000 |
| poly(snap5, 2)2 | 0.000000e+00 | 5.51 | -1.74 | 0.0817 |
Model 7
logitm7 <- svyglm(prep.current ~ prep.indic.any,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm7, exponentiate = T),
caption= "Model 7", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.02 | 0.46 | -8.11 | 0 |
| prep.indic.any | 21.75 | 0.48 | 6.47 | 0 |
Model 8
logitm8 <- svyglm(prep.current ~ age.grp + race + region,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm8, exponentiate = T),
caption= "Model 8", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.08 | 0.62 | -4.01 | 0.0001 |
| age.grp2 | 4.91 | 0.64 | 2.49 | 0.0132 |
| age.grp3 | 6.60 | 0.65 | 2.90 | 0.0039 |
| age.grp4 | 7.03 | 0.65 | 2.98 | 0.0030 |
| age.grp5 | 3.95 | 0.68 | 2.01 | 0.0450 |
| raceB | 2.19 | 0.66 | 1.18 | 0.2384 |
| raceH | 0.46 | 0.47 | -1.65 | 0.0998 |
| regionEasternWA | 0.24 | 0.46 | -3.12 | 0.0019 |
| regionWesternWA | 0.48 | 0.24 | -3.01 | 0.0027 |
Model 9
logitm9 <- svyglm(prep.current ~ age.grp + race + region +
poly(snap5, 2),
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm9, exponentiate = T),
caption= "Model 9", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 4.00000e-02 | 0.62 | -5.30 | 0.0000 |
| age.grp2 | 6.35000e+00 | 0.63 | 2.95 | 0.0033 |
| age.grp3 | 7.64000e+00 | 0.65 | 3.14 | 0.0018 |
| age.grp4 | 9.90000e+00 | 0.67 | 3.41 | 0.0007 |
| age.grp5 | 4.85000e+00 | 0.68 | 2.32 | 0.0208 |
| raceB | 1.31000e+00 | 0.65 | 0.42 | 0.6765 |
| raceH | 3.70000e-01 | 0.55 | -1.83 | 0.0680 |
| regionEasternWA | 1.70000e-01 | 0.49 | -3.69 | 0.0002 |
| regionWesternWA | 3.40000e-01 | 0.31 | -3.57 | 0.0004 |
| poly(snap5, 2)1 | 1.68751e+18 | 5.05 | 8.32 | 0.0000 |
| poly(snap5, 2)2 | 0.00000e+00 | 4.64 | -2.23 | 0.0264 |
Model 10
logitm10 <- svyglm(prep.current ~ age.grp + race + region +
prep.indic.any,
data = prepPrevDF,
family=quasibinomial(),
design)
kable(tidy(logitm10, exponentiate = T),
caption= "Model 10", digits = c(rep(2,4), 4)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.01 | 0.76 | -6.73 | 0.0000 |
| age.grp2 | 6.43 | 0.65 | 2.86 | 0.0044 |
| age.grp3 | 7.01 | 0.67 | 2.92 | 0.0036 |
| age.grp4 | 9.22 | 0.68 | 3.27 | 0.0012 |
| age.grp5 | 4.80 | 0.70 | 2.25 | 0.0249 |
| raceB | 1.26 | 0.70 | 0.33 | 0.7385 |
| raceH | 0.47 | 0.48 | -1.57 | 0.1176 |
| regionEasternWA | 0.21 | 0.47 | -3.30 | 0.0010 |
| regionWesternWA | 0.43 | 0.26 | -3.21 | 0.0014 |
| prep.indic.any | 22.85 | 0.47 | 6.67 | 0.0000 |
Comparison of AICs
It’s interesting to see that the SNAP variable performs better than the PrEP indication variable in predicting current PrEP use. One interpretation is that the SNAP specification – parametric up to 5+, with declining marginal impact – captures the tail impacts better than the binary PrEP indicated variable. Another is that the decision to initiate PrEP is influenced more by a person’s own behavior than by their partner-related risks, which is consistent with the pattern we identified above (see Non-use among indicated section).
models <- get_formula("logitm", 10)
aics <- get_aics("logitm", 10, "logit")
kable(data.frame(Model = models,
AIC = aics),
caption= "AIC comparison", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| Model | AIC |
|---|---|
| prep.current ~ age.grp | 596.1 |
| prep.current ~ race | 610.0 |
| prep.current ~ region | 598.3 |
| prep.current ~ insurance | 597.3 |
| prep.current ~ income_cat | 576.9 |
| prep.current ~ poly(snap5, 2) | 463.5 |
| prep.current ~ prep.indic.any | 513.7 |
| prep.current ~ age.grp + race + region | 581.7 |
| prep.current ~ age.grp + race + region + poly(snap5, 2) | 425.4 |
| prep.current ~ age.grp + race + region + prep.indic.any | 485.3 |
Predicted distribution
expand_df <- expand.grid(
age.grp = levels(prepPrevDF$age.grp),
race = levels(prepPrevDF$race),
region = levels(prepPrevDF$region),
snap5 = 0:5)
pred.prep.curr <- cbind(expand_df,
prob = predict(logitm9,
newdata = expand_df,
type = "response",
npop = 1000))
pred.prep.curr <- pred.prep.curr %>%
select(-prob.SE) %>%
rename(prob = prob.response)%>%
data.table::setDT()
data.table::setkeyv(pred.prep.curr,
c("age.grp", "race", "region", "snap5"))
## Save object
save_out$prep.current.dt <- pred.prep.curr
## Plot predicted values
fig <- ggplot(data = pred.prep.curr,
aes(x = snap5, y = prob, text = round(prob,2),
group = age.grp, color = age.grp)) +
geom_line(alpha = 0.7) +
geom_point() +
facet_grid(race ~ region) +
labs(x = "SNAP (topcoded at 5+)", y = "predicted probability") +
theme(plot.margin = margin(t=30,r=20,b=20,l=20),
legend.position = "bottom")
title.plotly <-
list(text = paste0('Predicted PrEP use by snap, age, race and region',
'<br>',
'<sup>',
'n = ',
dim(prepPrevDF)[[1]],
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(title = title.plotly)%>%
layout_ggplotly # fix marginsDuration of PrEP use
As with prevalence of current PrEP use, this variable will be used for two purposes:
Estimation of first-time PrEP initiation rates (Incidence = Prevalence/Duration). Model based, restricted to current PrEP and never interrupted. Interruption identifies the sample of cyclers.
Validation target – Descriptives, for this purpose it only conditions “current users” and measures the total interval from first initiation of PrEP to survey date.
Note: we can’t get an unbiased estimate of the duration of a current PrEP spell, because we don’t have the start/stop dates for interrupters.
Variable:
mos.frst.prep, filtering byprep.stopandprep.intWHPP estimates are reported as hazards with 30% discontinuing by about 2 yrs (see the Introduction).
Descriptive statistics
Overall
Table
# Uses the prepCurrentDF
prep.dur_all <- prepCurrentDF %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur_all %>%
kable(caption= "Months since first initiated PrEP", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current users, includes interrupters: n =",
nrow(prepCurrentDF)))| wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|
| 26.6 | 19.3 | 24.2 | 8 | 41 | 161 |
| Note: | |||||
| Current users, includes interrupters: n = 161 |
Plot
ggplot(prepCurrentDF, aes(x = mos.frst.prep)) +
geom_histogram(aes(y = ..density..),
color = "grey30", fill = "white") +
geom_density(alpha = .2, fill = "blue") +
geom_vline(xintercept = prep.dur_all["wtd.mean"][[1]],
color = "red") +
labs(title = "Months since first started PrEP",
x = "Months",
caption = paste("Current users, includes interrupters: n = ",
nrow(prepCurrentDF)))
Age
Too few observations in the breakdown to plot.
prep.dur.age <- prepCurrentDF %>%
group_by(age_group) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur.age %>%
kable(caption= "Months since first initiated PrEP by age",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current users, includes interrupters; n = ",
nrow(prepCurrentDF)))| age_group | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| 15-24 | 20.6 | 22.1 | 15.0 | 7.0 | 25.0 | 7 |
| 25-34 | 24.7 | 17.7 | 25.0 | 7.0 | 39.0 | 55 |
| 35-44 | 25.2 | 20.2 | 19.0 | 8.0 | 43.0 | 38 |
| 45-54 | 31.6 | 20.3 | 30.0 | 12.4 | 48.0 | 41 |
| 55-65 | 27.7 | 16.8 | 27.9 | 15.6 | 38.9 | 20 |
| Note: | ||||||
| Current users, includes interrupters; n = 161 |
Race
Table
prep.dur.race <- prepCurrentDF %>%
group_by(race) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur.race %>%
kable(caption= "Months since first initiated PrEP by race",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current users, includes interrupters; n = ",
nrow(prepCurrentDF)))| race | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| B | 18.2 | 16.1 | 8.0 | 8.0 | 36.6 | 6 |
| H | 36.6 | 17.3 | 41.3 | 28.1 | 52.2 | 13 |
| O | 26.8 | 19.4 | 23.1 | 9.5 | 41.0 | 142 |
| Note: | ||||||
| Current users, includes interrupters; n = 161 |
Plot
ggplot(prepCurrentDF,
aes(x = race, y = mos.frst.prep,
text = mos.frst.prep,
group = race, fill = race)) +
geom_boxplot(varwidth = T, notch = T) +
labs(title = "Months since first PrEP use",
caption = paste0("On PrEP: n = ", nrow(prepCurrentDF),
"; unweighted quantiles"),
y = "months")
Region
Table
prep.dur.region <- prepCurrentDF %>%
group_by(region) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur.region %>%
kable(caption= "Months since first initiated PrEP by region",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current users, includes interrupters; n = ",
nrow(prepCurrentDF)))| region | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| EasternWA | 26.0 | 19.4 | 27.8 | 19.9 | 41.0 | 11 |
| King | 28.7 | 20.3 | 27.0 | 8.0 | 43.0 | 94 |
| WesternWA | 21.0 | 15.3 | 19.3 | 9.7 | 35.9 | 56 |
| Note: | ||||||
| Current users, includes interrupters; n = 161 |
Plot
ggplot(prepCurrentDF,
aes(x = region, y = mos.frst.prep,
text = mos.frst.prep,
group = region, fill = region)) +
geom_boxplot(varwidth = T, notch = T) +
labs(title = "Months since first PrEP use",
caption = paste0("On PrEP: n = ", nrow(prepCurrentDF),
"; unweighted quantiles"),
y = "months")
Insurance
Table
Too few cases in some categories to plot.
prep.dur_ins <- prepCurrentDF %>%
group_by(insurance) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur_ins %>%
kable(caption= "Months since first initiated PrEP by ins",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current users, includes interrupters; n = ",
nrow(prepCurrentDF)))| insurance | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| emplyr | 28.0 | 19.8 | 27.0 | 8.0 | 43.0 | 120 |
| indvdl_othr | 26.2 | 15.7 | 31.3 | 18.0 | 36.3 | 13 |
| medicare_caid | 17.0 | 13.8 | 20.0 | 3.3 | 25.3 | 15 |
| othr_gov | 25.7 | 23.3 | 26.4 | 7.0 | 40.8 | 9 |
| uninsrd_NA | 16.4 | 15.5 | 20.4 | 12.5 | 36.6 | 4 |
| Note: | ||||||
| Current users, includes interrupters; n = 161 |
Income
Table
prep.dur_income <- prepCurrentDF %>%
group_by(income_cat) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur_income %>%
kable(caption= "Months since first initiated PrEP by FPL income",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current users, includes interrupters; n = ",
nrow(prepCurrentDF)))| income_cat | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| fpl138 | 22.2 | 19.0 | 25.8 | 18.5 | 43.4 | 6 |
| fpl138-300 | 25.7 | 15.2 | 25.2 | 14.3 | 38.9 | 23 |
| fpl300-400 | 25.9 | 20.7 | 24.7 | 6.9 | 45.8 | 21 |
| fpl400-500 | 25.9 | 22.5 | 26.2 | 7.0 | 46.2 | 10 |
| fpl500-600 | 26.8 | 14.2 | 25.6 | 21.0 | 34.0 | 11 |
| fpl600+ | 27.6 | 20.3 | 26.5 | 8.0 | 43.0 | 87 |
| Missing | 15.8 | 16.4 | 21.0 | 12.5 | 32.6 | 3 |
| Note: | ||||||
| Current users, includes interrupters; n = 161 |
Plot
ggplot(prepCurrentDF %>% filter(income_cat != "Missing"),
aes(x = income_cat, y = mos.frst.prep,
text = mos.frst.prep,
group = income_cat, fill = income_cat)) +
geom_boxplot(varwidth = T, notch = T) +
labs(title = "Months since first PrEP use",
caption = paste0("On PrEP and not missing income: n = ",
table(prepCurrentDF$income_cat=="Missing")[[1]],
"; unweighted quantiles"),
y = "months")
SNAP
Table (5)
prep.dur.snap5 <- prepCurrentDF %>%
mutate(snap5 = pmin(snap, 5)) %>%
group_by(snap5) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur.snap5 %>%
kable(caption= "Months since first initiated PrEP by SNAP",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste0("Current users, includes interrupters; n = ",
nrow(prepCurrentDF),
"; SNAP topcoded at 5+"))| snap5 | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| 0 | 8.1 | 2.6 | 10.3 | 8.1 | 11.0 | 3 |
| 1 | 19.6 | 14.7 | 19.5 | 9.6 | 35.8 | 13 |
| 2 | 26.1 | 20.8 | 21.9 | 12.0 | 48.8 | 10 |
| 3 | 15.9 | 14.8 | 8.0 | 6.7 | 29.2 | 15 |
| 4 | 31.6 | 22.9 | 38.4 | 17.4 | 54.0 | 15 |
| 5 | 29.7 | 19.3 | 27.0 | 15.0 | 43.0 | 96 |
| NA | 27.6 | 19.5 | 31.4 | 17.0 | 48.5 | 9 |
| Note: | ||||||
| Current users, includes interrupters; n = 161; SNAP topcoded at 5+ |
Plot
prepCurrentDF %>%
filter(snap != "Missing") %>%
mutate(snap5 = pmin(snap, 5)) %>%
ggplot(aes(x = snap5, y = mos.frst.prep,
text = mos.frst.prep,
group = snap5, fill = snap5)) +
geom_boxplot(varwidth = T, notch = T) +
labs(title = "Months since first PrEP use",
caption = paste0("On PrEP and not missing snap: n = ",
table(prepCurrentDF$snap=="Missing")[[1]],
"; unweighted quantiles"),
y = "months")
Deg.tot
Table
prep.dur_deg.tot3 <- prepCurrentDF %>%
mutate(deg.tot3 = pmin(deg.tot, 3)) %>%
group_by(deg.tot3) %>%
summarise(wtd.mean = wtd.mean(mos.frst.prep, weights=ego.wawt),
wtd.sd = sqrt(wtd.var(mos.frst.prep, weights=ego.wawt)),
wtd.med = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.5),
wtd.25 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.25),
wtd.75 = wtd.quantile(mos.frst.prep, weights=ego.wawt, probs=0.75),
n.obs = n())
prep.dur_deg.tot3 %>%
kable(caption= "Months since first initiated PrEP by active degree",
digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste0("Current users, includes interrupters; n = ",
nrow(prepCurrentDF),
"; deg.tot topcoded at 3+"))| deg.tot3 | wtd.mean | wtd.sd | wtd.med | wtd.25 | wtd.75 | n.obs |
|---|---|---|---|---|---|---|
| 0 | 19.7 | 19.4 | 10.6 | 6.0 | 31.9 | 27 |
| 1 | 28.0 | 20.3 | 21.5 | 11.1 | 45.8 | 45 |
| 2 | 28.6 | 16.9 | 31.0 | 20.0 | 41.0 | 38 |
| 3 | 27.0 | 19.8 | 24.6 | 8.0 | 41.7 | 42 |
| NA | 27.6 | 19.5 | 31.4 | 17.0 | 48.5 | 9 |
| Note: | ||||||
| Current users, includes interrupters; n = 161; deg.tot topcoded at 3+ |
Plot
prepCurrentDF %>%
filter(!is.na(deg.tot)) %>%
mutate(deg.tot3 = pmin(deg.tot, 3)) %>%
ggplot(aes(x = deg.tot3, y = mos.frst.prep,
text = mos.frst.prep,
group = deg.tot3, fill = deg.tot3)) +
geom_boxplot(varwidth = T, notch = T) +
labs(title = "Months since first PrEP use",
caption =
paste0("On PrEP and not missing deg.tot, deg.tot capped at 3+: n = ",
table(is.na(prepCurrentDF$deg.tot))[[1]],
"; unweighted quantiles"),
y = "months")
Prediction
Here we will just run the last model with all predictors. Sample excludes the interrupters, like the prevalence model, to get the group of first time PrEP users.
Full Model
# create filter var (tmp.interrupt) for prep interrupters
prepCurrentDF <- prepCurrentDF %>%
mutate(mop = if_else(mos.frst.prep==0, 0.5, mos.frst.prep),
tmp.interrupt = if_else(is.na(prep.int) | prep.int==0, 0, 1)) %>%
filter(tmp.interrupt == 0 & !is.na(snap)) %>%
mutate(age.grp = factor(age.grp),
race = relevel(factor(race), ref = "O"),
region = relevel(factor(region), ref = "King"),
snap5 = pmin(snap, 5))
glmdur <- glm(mop ~ age.grp + race + region + snap5,
family = poisson(link = "log"),
data = prepCurrentDF,
weights = ego.wawt)
kable(tidy(glmdur),
caption= "Full Model", digits = c(rep(2,4), 3)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 2.43 | 0.12 | 20.79 | 0.000 |
| age.grp2 | -0.08 | 0.09 | -0.83 | 0.409 |
| age.grp3 | 0.00 | 0.10 | 0.05 | 0.963 |
| age.grp4 | 0.33 | 0.09 | 3.54 | 0.000 |
| age.grp5 | 0.18 | 0.10 | 1.75 | 0.080 |
| raceB | -0.43 | 0.06 | -6.79 | 0.000 |
| raceH | 0.65 | 0.07 | 8.85 | 0.000 |
| regionEasternWA | -0.13 | 0.10 | -1.31 | 0.189 |
| regionWesternWA | -0.73 | 0.05 | -14.12 | 0.000 |
| snap5 | 0.21 | 0.02 | 13.14 | 0.000 |
Predicted distribution
expand_df <- expand.grid(
age.grp = levels(prepPrevDF$age.grp),
race = levels(prepPrevDF$race),
region = levels(prepPrevDF$region),
snap5 = 0:5)
pred.prep.dur <- cbind(expand_df,
dur = predict(glmdur,
newdata = expand_df,
type = "response",
npop = 1000)) %>%
data.table::setDT()
data.table::setkeyv(pred.prep.dur,
c("age.grp", "race", "region", "snap5"))
## Save objects
save_out$prep.dur.dt <- pred.prep.dur
## Plot predicted values
fig <- ggplot(data = pred.prep.dur,
aes(x = snap5, y = dur, text = round(dur, 1),
group = age.grp, color = age.grp)) +
geom_line(alpha = 0.7) +
geom_point() +
facet_grid(race ~ region) +
labs(x = "SNAP (topcoded at 5+)", y = "months on PrEP") +
theme(plot.margin = margin(t=30,r=20,b=20,l=20),
legend.position = "bottom")
title.plotly <-
list(text = paste0('Months on PrEP by snap, age, race and region',
'<br>',
'<sup>',
'Non-interrupters only, n = ',
dim(prepCurrentDF)[[1]],
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(title = title.plotly) %>%
layout_ggplotly # fix marginsInitiation rates (first time users only)
We estimate the rate of PrEP initiation from the ratio of the predicted values for prevalence and duration.
Te resulting data table will replace the current prep.rx 3-vector.
Estimated PrEP Incidence (prep.rx)
# this calculates a *weekly* incidence rate
pred.prep.newrx.wk <- pred.prep.dur %>%
left_join(pred.prep.curr) %>%
rename(prev = prob) %>%
mutate(inci.unadj = prev/(dur*4),
inci.adj = (prev/(1-prev))/(dur*4),
age.grp = as.numeric(age.grp)) %>%
data.table::setDT()
data.table::setkeyv(pred.prep.newrx.wk,
c("age.grp", "race", "region", "snap5"))
# save out predicted values
save_out$prep.newrx.wk.dt <- pred.prep.newrx.wkUnadjusted
IR = Prev / Dur
fig <- ggplot(data = pred.prep.newrx.wk,
aes(x = snap5, y = inci.unadj,
text = round(inci.unadj,2),
group = factor(age.grp), color = factor(age.grp))) +
geom_line(alpha = 0.7) +
geom_point() +
#ylim(c(0, 0.05)) +
facet_grid(race ~ region) +
labs(x = "SNAP (topcoded at 5+)", y = "rate per week") +
theme(plot.margin = margin(t=30,r=20,b=20,l=20),
legend.position = "bottom")
title.plotly <-
list(text = paste0('Unadjusted PrEP initiation rates/wk by snap, age, race and region',
'<br>',
'<sup>',
'numerator n = ',
dim(prepPrevDF)[[1]],
'; denominator n =',
dim(prepCurrentDF)[[1]],
'; first time PrEP users only',
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(title = title.plotly) %>%
layout_ggplotly # fix marginsAdjusted
IR = (Prev/(1-Prev)) / Dur
fig <- ggplot(data = pred.prep.newrx.wk,
aes(x = snap5, y = inci.adj,
text = round(inci.adj,2),
group = factor(age.grp), color = factor(age.grp))) +
geom_line(alpha = 0.7) +
geom_point() +
#ylim(c(0, 0.05)) +
facet_grid(race ~ region) +
labs(x = "SNAP (topcoded at 5+)", y = "rate per week") +
theme(plot.margin = margin(t=30,r=20,b=20,l=20),
legend.position = "bottom")
title.plotly <-
list(text = paste0('Adjusted PrEP initiation rates/wk by snap, age, race and region',
'<br>',
'<sup>',
'numerator n = ',
dim(prepPrevDF)[[1]],
'; denominator n =',
dim(prepCurrentDF)[[1]],
'; first time PrEP users only',
'</sup>'))
ggplotly(fig , tooltip = "text") %>%
layout(title = title.plotly) %>%
layout_ggplotly # fix marginsComparison
#fig <- #
ggplot(data = pred.prep.newrx.wk,
aes(x = inci.unadj, y = inci.adj,
text = round(inci.adj,2))) +
geom_point(aes(color = factor(snap5))) +
geom_abline(color = "grey") +
labs(title = "Incidence rate estimate comparison",
x = "Unadjusted IR = P/D", y = "Adjusted IR = P(1-P)/D") +
scale_color_brewer(palette = "Blues") +
theme_bw()
Discontinuation
The remaining dynamics of PrEP use are driven by quitting and cycling.
Components
In the WHAMP survey, we do not directly observe quitters, and we only observe a subsample of the cyclers. What we observe in the WHAMP survey is:
- Stopped – Respondents who have used PrEP (prep.ever) but are not using now (prep.stop). Note that this is a mix of those who quit, and those who cycle (will reinitiate at some time in the future).
- Interrupted – Respondents currently on PrEP (prepCurrent) who have ever interrupted (prep.int). This is the subsample of cyclers.
- The length of the interruption.
- For both stopping and interruption we ask the reason why. Comparison of the reason provides some insight into the quitters and cyclers.
Descriptives for each follow.
Stopped
For ever users: probability and reasons for stopping.
There are not enough observations to break down by race (e.g., no cases of stopping PrEP among the 6 black MSM reporting ever use), but the probability for Hispanics is the same as Others.
prepStopDF <- hivnegDF %>%
filter(prep.ever == 1) %>%
group_by(prep.stop) %>%
dplyr::summarize(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prepStopDF %>%
kable(caption= "Stopped PrEP", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever PrEP use: n = ", sum(prepStopDF$n.obs)))| prep.stop | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 161 | 167.5 | 0.7 |
| 1 | 61 | 57.8 | 0.3 |
| Note: | |||
| Ever PrEP use: n = 222 |
prep.stop_why <- hivnegDF %>%
filter(prep.stop == 1) %>%
group_by(prep.whystop) %>%
dplyr::summarize(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.stop_why %>%
mutate(prop.wtd = n.wtd/sum(n.wtd)) %>%
kable(caption= "Reason for stopping PrEP", digits = c(0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Stopped PrEP: n = ", sum(prep.stop_why$n.obs)))| prep.whystop | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| No HCP for Rx | 2 | 1.8 | 0.03 |
| No risk | 13 | 11.2 | 0.19 |
| Could not afford | 10 | 8.8 | 0.15 |
| Side effects | 10 | 10.2 | 0.18 |
| Take a break | 7 | 6.7 | 0.12 |
| Other | 11 | 10.9 | 0.19 |
| NA | 8 | 8.1 | 0.14 |
| Note: | |||
| Stopped PrEP: n = 61 |
Interrupted
Among current users, probability and reasons for interrupting.
prepIntDF <- hivnegDF %>%
filter(prepCurrent == "on prep") %>%
group_by(prep.int) %>%
dplyr::summarize(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prepIntDF %>%
kable(caption= "Interrupted PrEP", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Current PrEP users: n = ",
sum(prepIntDF$n.obs)))| prep.int | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 127 | 137.1 | 0.8 |
| 1 | 34 | 30.3 | 0.2 |
| Note: | |||
| Current PrEP users: n = 161 |
prep.int_why <- hivnegDF %>%
filter(prep.int == 1) %>%
group_by(prep.whyint) %>%
dplyr::summarize(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.int_why %>%
kable(caption= "Reason for interrupting PrEP",
digits = c(0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Interrupted PrEP: n = ", sum(prep.int_why$n.obs)))| prep.whyint | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| No HCP for Rx | 7 | 7.9 | 0.26 |
| No risk | 16 | 13.7 | 0.45 |
| Could not afford | 4 | 3.7 | 0.12 |
| Side effects | 1 | 0.8 | 0.03 |
| Other | 6 | 4.3 | 0.14 |
| Note: | |||
| Interrupted PrEP: n = 34 |
Plot comparisons
Frequency
tempDF <- hivnegDF %>%
filter(prep.ever == 1) %>%
mutate(type = case_when(prep.int==1 ~ "interrupted",
prep.stop==1 ~ "currently stopped",
TRUE ~ "other")) %>%
group_by(type) %>%
dplyr::summarize(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd)) %>%
select(type, prop.wtd) %>%
filter(type != "other")
p <- ggplot(tempDF,
aes(x=type, y=prop.wtd, text = round(prop.wtd,2))) +
geom_bar(stat="identity", fill="blue", alpha=.7) +
labs(title = "Prevalence of PrEP interrupting and stopping",
x = "Of those who have ever used PrEP")
ggplotly(p, tooltip = "text")Reason
compDF <- prep.stop_why %>%
left_join(prep.int_why, by = c("prep.whystop" = "prep.whyint")) %>%
select(reason = prep.whystop,
stop = prop.wtd.x,
interrupt = prop.wtd.y) %>%
pivot_longer(-reason,
names_to = "var",
values_to = "prop") %>%
mutate(reason = if_else(is.na(reason), "NA", as.character(reason)),
reason = factor(reason),
reason = forcats::fct_relevel(reason, "NA", after = Inf))
p <- ggplot(compDF,
aes(x = reason, y = prop, text = round(prop,2),
group = var, fill = var)) +
geom_bar(stat="identity", position = "dodge2", alpha = 0.7) +
labs(title = "Reasons for stopping or interrupting PrEP",
y = "weighted prop.") +
theme(axis.text.x = element_text(size = 8, angle = 50,
vjust = 0.6))
ggplotly(p, tooltip = "text") Interruption length
hivnegDF %>% filter(prep.int == 1) %>%
mutate(PREP_INTMOS = ifelse(PREP_INTMOS == 1.5, 2, PREP_INTMOS)) %>%
group_by(PREP_INTMOS) %>%
dplyr::summarize(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd)) %>%
kable(caption= "Months of PrEP interruption",
digits = c(0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Interrupted PrEP: n = ",
table(hivnegDF$prep.int)[[2]]))| PREP_INTMOS | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 1 | 11 | 9.7 | 0.32 |
| 2 | 7 | 7.2 | 0.24 |
| 3 | 2 | 1.6 | 0.05 |
| 5 | 3 | 2.0 | 0.06 |
| 6 | 5 | 3.7 | 0.12 |
| 7 | 3 | 3.6 | 0.12 |
| 8 | 3 | 2.6 | 0.08 |
| Note: | |||
| Interrupted PrEP: n = 34 |
hivnegDF %>% filter(prep.int==1) %>%
select(PREP_INTMOS) %>%
summarytools::descr(stats = "common", style = "rmarkdown") %>%
kable(caption= "Months of PrEP interruption", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| PREP_INTMOS | |
|---|---|
| Mean | 3.5 |
| Std.Dev | 2.6 |
| Min | 1.0 |
| Median | 2.0 |
| Max | 8.0 |
| N.Valid | 34.0 |
| Pct.Valid | 100.0 |
mos.prep.int.summary <- hivnegDF %>%
filter(prep.int==1) %>%
with(., summary(PREP_INTMOS))
mos.prep.int.summary.snap2 <- hivnegDF %>%
filter(prep.int==1 & !is.na(snap)) %>%
mutate(snap2 = if_else(snap < 2, 0, 1)) %>%
group_by(snap2) %>%
dplyr::summarize(mean.mos = mean(PREP_INTMOS, na.rm = T))Discontinuation rate
This combines the (currently) stopped with the ever interrupted to estimate the overall discontinuation rate among ever users. There are too few cases to do a full analysis by demographic attributes and SNAP.
Note this discontinuation rate will be used for both first time users and those who are cycling on and off (interrupters).
The currently stopped are a mix of quitters and cyclers.
Marginally, the SNAP effects are very large, and we should consider including these as a simplified relative risk (see below). The region effects are the only other moderately large effect, so we show those as a 3-element vector of rates, but will not use them for simulation.
We observe the prevalence of ever discontinued at the time of the survey. To translate this into a weekly rate of discontinuation we use the expression:
\[ p(\textsf{discontinue by t}) = 1 - (1-p(\textsf{discontinue per timestep}))^t \]
where:
\(p(\textsf{discontinue by t})\) is the observed prevalence of stopping or having interrupted PrEP, and
\(t\) is the observed mean number of months since first started PrEP.
We can then solve for the per timestep probability.
Overall
Prevalence
prep.discont_all <- prepDiscontDF %>%
group_by(prep.discont) %>%
summarise(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.discont_all %>%
kable(caption= "PrEP discontinuation rate",
digits = c(0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever PrEP: n = ", sum(prep.discont_all$n.obs)))| prep.discont | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 127 | 137.1 | 0.61 |
| 1 | 95 | 88.1 | 0.39 |
| Note: | |||
| Ever PrEP: n = 222 |
Rate (per week)
# Summary stats for months ago first started prep
hivnegDF %>% filter(prep.ever==1) %>%
select(mos.frst.prep) %>%
summarytools::descr(stats = "common", style = "rmarkdown") %>%
kable(caption = "Months since started first PrEP", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| mos.frst.prep | |
|---|---|
| Mean | 27.9 |
| Std.Dev | 18.9 |
| Min | 0.0 |
| Median | 27.0 |
| Max | 70.0 |
| N.Valid | 220.0 |
| Pct.Valid | 99.1 |
mos.frst.prep.summary <- hivnegDF %>%
filter(prep.ever==1) %>%
with(., summary(mos.frst.prep))
cDp <- prep.discont_all$prop.wtd[prep.discont_all$prep.discont==1]
mt <- mos.frst.prep.summary["Mean"]
# weekly rate
pred.prep.discont.wk.all <- 1 - (1-cDp)^(1/(4*mt))
print(paste("Weekly discontinuation rate of: ",
round(pred.prep.discont.wk.all,4)))## [1] "Weekly discontinuation rate of: 0.0044"# save out parameter
save_out$prep.discont.wk.all <- pred.prep.discont.wk.allRace
Too few cases to support a breakdown by race.
prep.discont.race <- prepDiscontDF %>%
group_by(race, prep.discont) %>%
summarise(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
group_by(race) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.discont.race %>%
kable(caption= "PrEP discontinuation rate by race",
digits = c(0,0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever PrEP & not NA for discontinuation: n = ",
sum(prep.discont.race$n.obs)))| race | prep.discont | n.obs | n.wtd | prop.wtd |
|---|---|---|---|---|
| B | 0 | 4 | 15.5 | 0.91 |
| B | 1 | 2 | 1.6 | 0.09 |
| H | 0 | 8 | 6.9 | 0.43 |
| H | 1 | 10 | 9.1 | 0.57 |
| O | 0 | 115 | 114.7 | 0.60 |
| O | 1 | 83 | 77.4 | 0.40 |
| Note: | ||||
| Ever PrEP & not NA for discontinuation: n = 222 |
Region
Moderately large effects here.
Prevalence
prep.discont.region <- prepDiscontDF %>%
group_by(region, prep.discont) %>%
summarise(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
group_by(region) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.discont.region %>%
kable(caption= "PrEP discontinuation rate by region",
digits = c(0,0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever PrEP & not NA for discontinuation: n = ",
sum(prep.discont.region$n.obs)))| region | prep.discont | n.obs | n.wtd | prop.wtd |
|---|---|---|---|---|
| EasternWA | 0 | 6 | 3.7 | 0.39 |
| EasternWA | 1 | 9 | 5.8 | 0.61 |
| King | 0 | 83 | 103.8 | 0.68 |
| King | 1 | 43 | 48.4 | 0.32 |
| WesternWA | 0 | 38 | 29.6 | 0.47 |
| WesternWA | 1 | 43 | 33.9 | 0.53 |
| Note: | ||||
| Ever PrEP & not NA for discontinuation: n = 222 |
Rate (weekly)
regiontab <- hivnegDF %>% filter(prep.ever==1) %>%
select(mos.frst.prep, region) %>%
with(., summarytools::stby(data=mos.frst.prep, INDICES = region,
FUN = summarytools::descr,
stats = c("common"),
style = "rmarkdown")) %>%
summarytools::tb(drop.var.col=TRUE)
regiontab %>%
kable(caption = "Months since started first PrEP", digits = 1) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| region | variable | mean | sd | min | med | max | n.valid | pct.valid |
|---|---|---|---|---|---|---|---|---|
| EasternWA | value | 29.7 | 21.0 | 0 | 28 | 61 | 15 | 100.0 |
| King | value | 31.3 | 19.8 | 0 | 30 | 70 | 125 | 99.2 |
| WesternWA | value | 22.2 | 15.6 | 1 | 20 | 68 | 80 | 98.8 |
cDp <- prep.discont.region$prop.wtd[prep.discont_all$prep.discont==1]
mt <- regiontab$mean
# weekly rate calculated in case we decide to use it
pred.prep.discont.wk.region <- 1 - (1-cDp)^(1/(4*mt))
print(paste("Weekly discontinuation rate of: ",
round(pred.prep.discont.wk.region,4)))## [1] "Weekly discontinuation rate of: 0.008"
## [2] "Weekly discontinuation rate of: 0.0031"
## [3] "Weekly discontinuation rate of: 0.0085"# save out parameter
save_out$prep.discont.wk.region <- pred.prep.discont.wk.regionSNAP
“No risk” is the single most important reason for interrupting PrEP, and also important among those who have stopped taking PrEP. Unlike race and region, SNAP is dynamic, so we can only see its impact from the currently stopped population, not from those who are on PrEP now, but interrupted at some prior period. So:
- Prevalence is estimated here from the stop variable,
- SNAP is binarized to {0,1} vs 2+
- How to estimate the weekly rate is unclear
- So we may want to parameterize this as a relative risk
This is exploratory, so should not be used unless the overall weekly probability does not produce a satisfactory model. In that event, we will need to verify the use of the relative risk approach.
Descriptive Table
prepDiscontDF %>%
mutate(snap3 = if_else(snap >2, 3, snap)) %>%
group_by(snap3, prep.stop) %>%
dplyr::summarize(n=n()) %>%
mutate(prop = n/sum(n)) %>%
pivot_wider(names_from = prep.stop,
values_from = c(n, prop)) %>%
kable(caption = "Stopped PrEP by SNAP",
col.names = c("SNAP3", rep(c("no", "yes"), 2)),
digits = c(0,0,0,2,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
add_header_above(c("", "N"=2, "Prop"=2)) %>%
add_header_above(c("", "Stopped PrEP"=4)) %>%
footnote(paste("Ever PrEP and not missing snap: n = ",
table(is.na(prepDiscontDF$snap))[[1]]))| SNAP3 | no | yes | no | yes |
|---|---|---|---|---|
| 0 | 3 | 2 | 0.60 | 0.40 |
| 1 | 13 | 17 | 0.43 | 0.57 |
| 2 | 10 | 4 | 0.71 | 0.29 |
| 3 | 126 | 31 | 0.80 | 0.20 |
| NA | 9 | 7 | 0.56 | 0.44 |
| Note: | ||||
| Ever PrEP and not missing snap: n = 206 |
Relative risk
prep.discont.snap2 <- prepDiscontDF %>%
filter(!is.na(snap)) %>%
mutate(snap2 = if_else(snap < 2, 0, 1)) %>%
group_by(snap2, prep.stop) %>%
summarise(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
group_by(snap2) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.discont.snap2 %>%
kable(caption= "PrEP stopped rate by snap2",
digits = c(0,0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever PrEP and not missing SNAP: n = ",
sum(prep.discont.snap2$n.obs)))| snap2 | prep.stop | n.obs | n.wtd | prop.wtd |
|---|---|---|---|---|
| 0 | 0 | 16 | 15.3 | 0.45 |
| 0 | 1 | 19 | 18.5 | 0.55 |
| 1 | 0 | 136 | 143.6 | 0.82 |
| 1 | 1 | 35 | 32.0 | 0.18 |
| Note: | ||||
| Ever PrEP and not missing SNAP: n = 206 |
## This is the relative risk if discontinuing PrEP if SNAP < 2
pred.prep.discont.snap2.rr <-
prep.discont.snap2$prop.wtd[prep.discont.snap2$prep.stop==1 &
prep.discont.snap2$snap2==1] /
prep.discont.snap2$prop.wtd[prep.discont.snap2$prep.stop==1 &
prep.discont.snap2$snap2==0]
print(paste("Relative risk of stopping PrEP for SNAP > 2 = ",
round(pred.prep.discont.snap2.rr, 2)))## [1] "Relative risk of stopping PrEP for SNAP > 2 = 0.33"# save out parameter
save_out$prep.discont.snap2.rr <- pred.prep.discont.snap2.rrQuitting
The answers given to the “why did you stop taking PrEP” question suggest some people are unlikely to start again, while others may cycle on and off as their behavioral risk changes. Some of the reasons are more common for stoppers, others for interrupters.
Side effects are the only reason reported almost exclusively by stoppers: 18% of stoppers vs. only 2% of interrupters (see the comparison plot in the Discontinuation section above).
So, if we need an estimate of the quitter fraction, 16% of those who discontinue the first time might be a good starting place.
prep.quit.frac <- 0.16
save_out$prep.quit.frac <- prep.quit.fracRe-initiation rate
This will be based on an even smaller sample of data: the “interrupter” fraction of those who ever discontinued PrEP. The rate of re-initiation will be based on the average length of interruption, using the same approach as for the discontinuation rate (see above).
It is very likely that re-initiation also reflects risk behavior, but we have a very sample size (34 re-initiators) and the quex didn’t ask about activity at the time of re-initiation. Still, the results are unambiguous.
Prevalence
Overall
prep.reinit_all <- prepReinitDF %>%
group_by(prep.reinit) %>%
summarise(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.reinit_all %>%
kable(caption= "PrEP re-initiation rate",
digits = c(0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever discontinued PrEP: n = ",
sum(prep.reinit_all$n.obs)))| prep.reinit | n.obs | n.wtd | prop.wtd |
|---|---|---|---|
| 0 | 61 | 57.8 | 0.66 |
| 1 | 34 | 30.3 | 0.34 |
| Note: | |||
| Ever discontinued PrEP: n = 95 |
SNAP
prep.reinit_snap2 <- prepReinitDF %>%
filter(!is.na(snap)) %>%
mutate(snap2 = if_else(snap < 2, 0, 1)) %>%
group_by(snap2, prep.reinit) %>%
summarise(n.obs = n(),
n.wtd = sum(ego.wawt)) %>%
mutate(prop.wtd = n.wtd/sum(n.wtd))
prep.reinit_snap2 %>%
kable(caption= "PrEP re-initiation rate by SNAP",
digits = c(0,0,0,1,2)) %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped")) %>%
footnote(paste("Ever discontinued PrEP: n = ",
sum(prep.reinit_all$n.obs)))| snap2 | prep.reinit | n.obs | n.wtd | prop.wtd |
|---|---|---|---|---|
| 0 | 0 | 19 | 18.5 | 0.84 |
| 0 | 1 | 4 | 3.5 | 0.16 |
| 1 | 0 | 35 | 32.0 | 0.56 |
| 1 | 1 | 28 | 25.1 | 0.44 |
| Note: | ||||
| Ever discontinued PrEP: n = 95 |
Rate (per week)
# Summary stats for calculating weekly rate
## HACK: we calculate the rate for all folks, but apply it only to SNAP > 1
cDp <- prep.reinit_all$prop.wtd[prep.reinit_all$prep.reinit==1]
cDp.snap2 <- prep.reinit_snap2$prop.wtd[prep.reinit_snap2$snap2==1 &
prep.reinit_snap2$prep.reinit==1]
mt <- as.numeric(mos.prep.int.summary["Mean"])
mt.snap2 <- mos.prep.int.summary.snap2$mean.mos[mos.prep.int.summary.snap2$snap2==1]
# weekly rate
pred.prep.reinit.rate.wk.all <- 1 - (1-cDp)^(1/(4*mt))
pred.prep.reinit.rate.wk.snap2 <- 1 - (1-cDp.snap2)^(1/(4*mt.snap2))
print(paste("All: Weekly PrEP re-initiation rate of: ",
round(pred.prep.reinit.rate.wk.all,4)))## [1] "All: Weekly PrEP re-initiation rate of: 0.0293"print(paste("SNAP2: Weekly PrEP re-initiation rate of: ",
round(pred.prep.reinit.rate.wk.snap2,4)))## [1] "SNAP2: Weekly PrEP re-initiation rate of: 0.037"# save out parameter
save_out$prep.reinit.wk <- list(all = pred.prep.reinit.rate.wk.all,
snap2 = pred.prep.reinit.rate.wk.snap2)Save
Parameters
descTable <-
tibble(Params = names(save_out),
Description = c("Awareness",
"Prevalence",
"Duration (mo)",
"Initiation rate (newrx, wk)",
"Discontinuation rate (wk)",
"Discontinuation rate (wk)",
"RR of discontinuation",
"Quitting fraction",
"Re-initiation rate (wk)"),
Subset = c("HIV-",
"1st time users",
"1st time users",
"1st time users",
"1st time users",
"1st time users",
"Stopped",
"1st time users",
"Interrupted"),
Method = c("logistic regr",
"logistic regr",
"Poisson glm",
"Ratio of Prevalence:Duration",
"Backcalculation",
"Backcalculation",
"Ratio of SNAP2+:SNAP{0,1}",
"Side effects comparison",
"Backcalculation"),
Levels = c("race x region",
rep("age x race x region x SNAP5",3),
"overall",
"region",
"SNAP2+",
"overall",
"overall and SNAP2+"))
descTable %>%
kable(caption = "WHAMP PrEP Parameters") %>%
kable_styling(full_width=F, position="center",
bootstrap_options = c("striped"))| Params | Description | Subset | Method | Levels |
|---|---|---|---|---|
| prep.aware.dt | Awareness | HIV- | logistic regr | race x region |
| prep.current.dt | Prevalence | 1st time users | logistic regr | age x race x region x SNAP5 |
| prep.dur.dt | Duration (mo) | 1st time users | Poisson glm | age x race x region x SNAP5 |
| prep.newrx.wk.dt | Initiation rate (newrx, wk) | 1st time users | Ratio of Prevalence:Duration | age x race x region x SNAP5 |
| prep.discont.wk.all | Discontinuation rate (wk) | 1st time users | Backcalculation | overall |
| prep.discont.wk.region | Discontinuation rate (wk) | 1st time users | Backcalculation | region |
| prep.discont.snap2.rr | RR of discontinuation | Stopped | Ratio of SNAP2+:SNAP{0,1} | SNAP2+ |
| prep.quit.frac | Quitting fraction | 1st time users | Side effects comparison | overall |
| prep.reinit.wk | Re-initiation rate (wk) | Interrupted | Backcalculation | overall and SNAP2+ |
save_out$descTable <- descTable # description table
print("Structure of parameter list object:")## [1] "Structure of parameter list object:"for(i in 1:length(names(save_out))){
str(save_out[i], vec.len=0, give.attr=F)}## List of 1
## $ prep.aware.dt:Classes 'data.table' and 'data.frame': 9 obs. of 3 variables:
## ..$ race : Factor w/ 3 levels "O","B","H": NULL ...
## ..$ region: Factor w/ 3 levels "King","EasternWA",..: NULL ...
## ..$ prob : num [1:9] NULL ...
## List of 1
## $ prep.current.dt:Classes 'data.table' and 'data.frame': 270 obs. of 5 variables:
## ..$ age.grp: Factor w/ 5 levels "1","2","3","4",..: NULL ...
## ..$ race : Factor w/ 3 levels "O","B","H": NULL ...
## ..$ region : Factor w/ 3 levels "King","EasternWA",..: NULL ...
## ..$ snap5 : int [1:270] NULL ...
## ..$ prob : num [1:270] NULL ...
## List of 1
## $ prep.dur.dt:Classes 'data.table' and 'data.frame': 270 obs. of 5 variables:
## ..$ age.grp: Factor w/ 5 levels "1","2","3","4",..: NULL ...
## ..$ race : Factor w/ 3 levels "O","B","H": NULL ...
## ..$ region : Factor w/ 3 levels "King","EasternWA",..: NULL ...
## ..$ snap5 : int [1:270] NULL ...
## ..$ dur : num [1:270] NULL ...
## List of 1
## $ prep.newrx.wk.dt:Classes 'data.table' and 'data.frame': 270 obs. of 8 variables:
## ..$ age.grp : num [1:270] NULL ...
## ..$ race : Factor w/ 3 levels "O","B","H": NULL ...
## ..$ region : Factor w/ 3 levels "King","EasternWA",..: NULL ...
## ..$ snap5 : int [1:270] NULL ...
## ..$ dur : num [1:270] NULL ...
## ..$ prev : num [1:270] NULL ...
## ..$ inci.unadj: num [1:270] NULL ...
## ..$ inci.adj : num [1:270] NULL ...
## List of 1
## $ prep.discont.wk.all: 'table' Named num NULL ...
## List of 1
## $ prep.discont.wk.region: num [1:3] NULL ...
## List of 1
## $ prep.discont.snap2.rr: num NULL ...
## List of 1
## $ prep.quit.frac: num NULL ...
## List of 1
## $ prep.reinit.wk:List of 2
## ..$ all : num NULL ...
## ..$ snap2: num NULL ...
## List of 1
## $ descTable: tibble [9 x 5] (S3: tbl_df/tbl/data.frame)
## ..$ Params : chr [1:9] ...
## ..$ Description: chr [1:9] ...
## ..$ Subset : chr [1:9] ...
## ..$ Method : chr [1:9] ...
## ..$ Levels : chr [1:9] ...saveRDS(save_out, here::here("Data", "Params", "WhampPrepParam.RDS"))Targets
Calls the external script makePrepTargets.R. Targets are the prevalence of PrEP status variables.
source(here::here("MakeData", "MM", "Scripts", "makePrepTargets.R"),
echo = T)##
## > library(dplyr)
##
## > library(ggpubr)
##
## > dataset <- "makeit"
##
## > if (dataset == "makeit") {
## + wideDF <- readRDS(here::here("MakeData", "MM", "Data", "wideDF.RDS"))
## + hivnegDF <- wideDF %>% filter(diag.stat .... [TRUNCATED]
##
## > summaries <- list(nobs = ~n(), n.valid = ~sum(!is.na(.x)),
## + n.missing = ~sum(is.na(.x)), wtd.n = ~sum(ego.wawt[!is.na(.x)]),
## + wtd.mean = .... [TRUNCATED]
##
## > all <- hivnegDF %>% summarise(across(prep.aware:prep.int.mos,
## + summaries)) %>% data.frame() %>% tidyr::pivot_longer(everything(),
## + names .... [TRUNCATED]
##
## > age <- hivnegDF %>% group_by(age.grp) %>% summarise(across(prep.aware:prep.int.mos,
## + summaries)) %>% data.frame() %>% tidyr::pivot_longer(!age .... [TRUNCATED]
##
## > race <- hivnegDF %>% group_by(race) %>% summarise(across(prep.aware:prep.int.mos,
## + summaries)) %>% data.frame() %>% tidyr::pivot_longer(!race, .... [TRUNCATED]
##
## > region <- hivnegDF %>% group_by(region) %>% summarise(across(prep.aware:prep.int.mos,
## + summaries)) %>% data.frame() %>% tidyr::pivot_longer(!r .... [TRUNCATED]
##
## > snap5 <- hivnegDF %>% filter(!is.na(snap5)) %>% group_by(snap5) %>%
## + summarise(across(prep.aware:prep.int.mos, summaries)) %>%
## + data.fra .... [TRUNCATED]
##
## > xtab <- hivnegDF %>% filter(!is.na(prep.ever)) %>%
## + group_by(prep.ever, prep.current) %>% summarise(nobs = n(),
## + n.wtd = round(sum(ego.w .... [TRUNCATED]
##
## > targets_out <- list(all = all, age = age, race = race,
## + region = region, snap5 = snap5, everxcurrent = xtab, makefile = "makePrepTargets, make ..." ... [TRUNCATED]print("Structure of target list object:")## [1] "Structure of target list object:"str(targets_out)## List of 7
## $ all : tibble [8 x 9] (S3: tbl_df/tbl/data.frame)
## ..$ var : chr [1:8] "prep.aware" "prep.ever" "prep.indic.any" "prep.current" ...
## ..$ nobs : int [1:8] 831 831 831 831 831 831 831 831
## ..$ n.valid : int [1:8] 756 755 641 755 222 222 161 34
## ..$ n.missing : int [1:8] 75 76 190 76 609 609 670 797
## ..$ wtd.n : num [1:8] 750 749 636 749 225 ...
## ..$ wtd.mean : num [1:8] 0.922 0.301 0.645 0.223 0.391 ...
## ..$ wtd.sd : num [1:8] 0.269 0.459 0.479 0.417 0.489 ...
## ..$ wtd.semean: num [1:8] 0.00934 0.01594 0.01663 0.01448 0.01699 ...
## ..$ wtd.median: Named num [1:8] 1 0 1 0 0 0 0 2
## .. ..- attr(*, "names")= chr [1:8] "50%" "50%" "50%" "50%" ...
## $ age : tibble [40 x 10] (S3: tbl_df/tbl/data.frame)
## ..$ age.grp : num [1:40] 1 1 1 1 1 1 1 1 2 2 ...
## ..$ var : chr [1:40] "prep.aware" "prep.ever" "prep.indic.any" "prep.current" ...
## ..$ nobs : int [1:40] 124 124 124 124 124 124 124 124 251 251 ...
## ..$ n.valid : int [1:40] 112 112 86 112 14 14 7 4 232 231 ...
## ..$ n.missing : int [1:40] 12 12 38 12 110 110 117 120 19 20 ...
## ..$ wtd.n : num [1:40] 135 135 104.9 135 17.1 ...
## ..$ wtd.mean : num [1:40] 0.863 0.126 0.623 0.073 0.728 ...
## ..$ wtd.sd : num [1:40] 0.346 0.334 0.487 0.261 0.459 ...
## ..$ wtd.semean: num [1:40] 0.0275 0.0266 0.0388 0.0208 0.0365 ...
## ..$ wtd.median: Named num [1:40] 1 0 1 0 1 ...
## .. ..- attr(*, "names")= chr [1:40] "50%" "50%" "50%" "50%" ...
## $ race : tibble [24 x 10] (S3: tbl_df/tbl/data.frame)
## ..$ race : chr [1:24] "B" "B" "B" "B" ...
## ..$ var : chr [1:24] "prep.aware" "prep.ever" "prep.indic.any" "prep.current" ...
## ..$ nobs : int [1:24] 27 27 27 27 27 27 27 27 89 89 ...
## ..$ n.valid : int [1:24] 24 24 21 24 6 6 6 2 82 82 ...
## ..$ n.missing : int [1:24] 3 3 6 3 21 21 21 25 7 7 ...
## ..$ wtd.n : num [1:24] 43.9 43.9 37.6 43.9 17.1 ...
## ..$ wtd.mean : num [1:24] 0.9273 0.3901 0.9312 0.3901 0.0922 ...
## ..$ wtd.sd : num [1:24] 0.263 0.493 0.256 0.493 0.298 ...
## ..$ wtd.semean: num [1:24] 0.0354 0.0665 0.0346 0.0665 0.0402 ...
## ..$ wtd.median: Named num [1:24] 1 0 1 0 0 0 0 2 1 0 ...
## .. ..- attr(*, "names")= chr [1:24] "50%" "50%" "50%" "50%" ...
## $ region : tibble [24 x 10] (S3: tbl_df/tbl/data.frame)
## ..$ region : chr [1:24] "EasternWA" "EasternWA" "EasternWA" "EasternWA" ...
## ..$ var : chr [1:24] "prep.aware" "prep.ever" "prep.indic.any" "prep.current" ...
## ..$ nobs : int [1:24] 120 120 120 120 120 120 120 120 376 376 ...
## ..$ n.valid : int [1:24] 105 105 88 105 15 15 11 5 344 343 ...
## ..$ n.missing : int [1:24] 15 15 32 15 105 105 109 115 32 33 ...
## ..$ wtd.n : num [1:24] 75.4 75.4 62.6 75.4 9.5 ...
## ..$ wtd.mean : num [1:24] 0.7871 0.126 0.578 0.0919 0.6131 ...
## ..$ wtd.sd : num [1:24] 0.412 0.334 0.498 0.291 0.515 ...
## ..$ wtd.semean: num [1:24] 0.0446 0.0361 0.0539 0.0314 0.0557 ...
## ..$ wtd.median: Named num [1:24] 1 0 1 0 1 ...
## .. ..- attr(*, "names")= chr [1:24] "50%" "50%" "50%" "50%" ...
## $ snap5 : tibble [48 x 10] (S3: tbl_df/tbl/data.frame)
## ..$ snap5 : num [1:48] 0 0 0 0 0 0 0 0 1 1 ...
## ..$ var : chr [1:48] "prep.aware" "prep.ever" "prep.indic.any" "prep.current" ...
## ..$ nobs : int [1:48] 104 104 104 104 104 104 104 104 210 210 ...
## ..$ n.valid : int [1:48] 104 104 104 104 5 5 3 2 210 210 ...
## ..$ n.missing : int [1:48] 0 0 0 0 99 99 101 102 0 0 ...
## ..$ wtd.n : num [1:48] 93.98 93.98 93.98 93.98 4.97 ...
## ..$ wtd.mean : num [1:48] 0.9056 0.0529 0 0.0288 0.784 ...
## ..$ wtd.sd : num [1:48] 0.294 0.225 0 0.168 0.46 ...
## ..$ wtd.semean: num [1:48] 0.0303 0.0232 0 0.0173 0.0475 ...
## ..$ wtd.median: Named num [1:48] 1 0 0 0 1 ...
## .. ..- attr(*, "names")= chr [1:48] "50%" "50%" "50%" "50%" ...
## $ everxcurrent: tibble [3 x 5] (S3: tbl_df/tbl/data.frame)
## ..$ prep.ever : num [1:3] 0 1 1
## ..$ prep.current: num [1:3] 0 0 1
## ..$ nobs : int [1:3] 533 61 161
## ..$ n.wtd : num [1:3] 524.1 57.8 167.5
## ..$ prop.wtd : num [1:3] 0.7 0.08 0.22
## $ makefile : chr "makePrepTargets, makePrepDynamics"saveRDS(targets_out,
file = here::here("Data", "Targets", "prepTargets.RDS"))