run_model.Rmd

---
title: "Impact of antiviral treatment on influenza-associated hospitalizations"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = FALSE, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 10, fig.height = 6)
```

```{r}
pacman::p_load(here, tidyverse, scales, ggsci, data.table, patchwork, ggstats)

set.seed(123)
nsim <- 1000

cols <- c("#aa4499","#999933", "#88ccee")
cols_alt <- c("#ddcc77",  "#88ccee", "#44aa99","#117733", "#882255")
cols_alt2 <- c("#ddcc77",  "#88ccee", "#44aa99","#117733",  "#cc6677", "#882255", "#332288")
cols_risk <- c( "#cc6677", "#332288", "#44aa99")

# geom settings
boxpt   <- 21
boxsize <- 0.3

ptsize <- 2.5
dodgewdth <- 0.4
barsize <- 0.4
barwdth <- 0.7

## Load model functions
source("functions.R")
```


### Initial parameters

```{r}
# Age groups 
age.grps <- c("0-4", "5-17", "18-49", "50-64", "65+")
n.ages   <- length(age.grps)

# labels for age groups when plotting (\u2265 = greater than or equal to)
age_labels <- c(age.grps[1:4], '\u226565')

# Proportion high risk in each age group (Zimmerman 2010)
high.risk <- c(0.05, 0.10, 0.20, 0.35, 0.55)

# Calculate risk stratifications for parameters with population averages and relative ratios

# CHR
chr <- get_risk_stratifications(pop = 1/c(143.44, 364.71, 178.16, 94.30, 11.00),   
                                rr =  c(1, 1, 5.5, 8.9, 4.9),      # Matias 2017 
                                prop.high = high.risk)
 
# Care-seeking lower bound
seek.care.lower <- get_risk_stratifications(pop = c(0.25, 0.25, 0.35, 0.45, 0.45), # ONM
                                            rr =  rep(1.3, 5),    # Biggerstaff 2010       
                                            prop.high = high.risk) 

# Care-seeking upper bound
seek.care.upper <- get_risk_stratifications(pop = c(0.35, 0.35, 0.50, 0.55, 0.55), # ONM
                                            rr =  rep(1.3, 5),        # Biggerstaff 2010       
                                            prop.high = high.risk) 

# Table of baseline parameter ranges for each age group 
# (does not include parameters which are defined by low/int/high scenarios)
# lo/hi refer to lower / upper limits of assumed distribution
# different values in each vector c() are for the different age groups

# Low risk first
pars.table.lr <- 
    data.frame(
        risk = "low",
        age = c("0-4", "5-17", "18-49", "50-64", "65+"),  
        # population, 2022 cUS ensus
        population = c(18538353 + 3683113, 53912474, 140148894, 62892984, 57794852),
        # proportion low risk in each age group 
        proportion = 1 - high.risk,
        # clinical attack rate (CAR), 2022/23
        car.lo = c(0.11, 0.12, 0.06, 0.07, 0.02), 
        car.hi = c(0.34, 0.44, 0.14, 0.18, 0.09),
        # molecular/ rapid test sensitivity (for those tested)
        test.sens.lo = rep(0.8, n.ages),
        test.sens.hi = rep(1.0, n.ages),
        # care-seeking overall proportions 
        seek.care.lo = seek.care.lower$low, 
        seek.care.hi = seek.care.upper$low, 
        # early care-seeking proportion (within 2 days)
        seek.early.lo = c(0.45, 0.45, 0.45, 0.40, 0.35), 
        seek.early.hi = c(0.55, 0.55, 0.55, 0.50, 0.45), 
        # risk of hospitalization (assume constant for now, i.e. lo = hi)
        hosp.risk.lo = chr$low, 
        hosp.risk.hi = chr$low, 
        # % taking prescribed NAIs
        nai.compliance.lo = rep(0.4, n.ages), 
        nai.compliance.hi = rep(0.6, n.ages)
    )

# High risk next - most parameters the same so copy low-risk
pars.table.hr <- pars.table.lr %>%
                    mutate(
                       risk = "high",
                       proportion = high.risk,
                       # based on Biggerstaff et al 2010
                       seek.care.lo = seek.care.lower$high, 
                       seek.care.hi = seek.care.upper$high, 
                       # based on weighted average calcs above
                       hosp.risk.lo = chr$high,
                       hosp.risk.hi = chr$high
                    )

pars.table <- rbind(pars.table.lr, pars.table.hr)

# Get baseline parameters with ranges
pars.range <- get_parameters(nsamples = nsim, pars.table) %>% mutate(age = factor(age, levels = age.grps))
```

```{r, fig.height = 6, results='hide'}
# Plot select parameter distributions
pars_long <- pars.range %>% 
    gather(param, value, 
           clinical.attack.rate, hosp.risk, seek.care, seek.care.early, test.sens, nai.compliance) 

# add better parameter names for plotting
cairo_pdf(here("figs", "fig2.pdf"), width = 8, height = 5)
rename_pars(pars_long) %>%
    # hosp.risk doesn't vary --> don't need to plot distribution
    filter(param != "hosp.risk") %>%
    ggplot(aes(x = age, y = value, fill = risk)) + 
    geom_violin(alpha = 0.5, position = position_dodge(width = 0.75)) + 
    geom_boxplot(width = 0.1, position = position_dodge(width = 0.75)) +
    scale_fill_manual("Risk", values = cols_risk[2:3]) + ylim(c(0, 1)) +
    facet_wrap(~ parname, ncol = 3) + 
    get_theme(txt = 11, legend.position = c(0.85, 0.25)) +
    scale_x_discrete(labels = age_labels) +
    labs(x = "Age (years)", y = "Sampled value")
dev.off()
```

```{r}
# Define scenarios for remaining parameters

# Testing scenarios (% tested for early vs late care-seekers)
# Start with low risk
testing.pars.low <- data.frame(
                    test.level = paste0(c("low", "intermediate", "high"), " testing"),
                    perc.tested.early = c(0.05, 0.25, 0.5),
                    perc.tested.late  = c(0.05, 0.05, 0.1)
                )
 
testing.pars.high <- testing.pars.low
 
# % prescribed NAIs for early vs late care-seekers
prescribing.pars.low <- data.frame(
                        prescribe.level = paste0(c("low", "high"), " prescribing"),
                        # % prescribed NAIs, with +ve test
                        prescribed.nai.pos.early = c(0.25, 0.25), 
                        prescribed.nai.pos.late  = c(0.10, 0.10),
                        # % prescribed NAIs, with -ve test 
                        prescribed.nai.neg.early = c(0.05, 0.05), 
                        prescribed.nai.neg.late  = c(0.00, 0.00),
                        # % prescribed NAIs, without test
                        prescribed.nai.no.early = c(0.10, 0.10), 
                        prescribed.nai.no.late  = c(0.05, 0.05)
                )
 
prescribing.pars.high <- data.frame(
                        prescribe.level = paste0(c("low", "high"), " prescribing"),
                        # % prescribed NAIs, with +ve test
                        prescribed.nai.pos.early = c(0.25, 0.50), 
                        prescribed.nai.pos.late  = c(0.10, 0.20),
                        # % prescribed NAIs, with -ve test 
                        prescribed.nai.neg.early = c(0.05, 0.10), 
                        prescribed.nai.neg.late  = c(0.00, 0.05),
                        # % prescribed NAIs, without test
                        prescribed.nai.no.early = c(0.10, 0.20), 
                        prescribed.nai.no.late  = c(0.05, 0.10)
                )

# NAI effectiveness against hospitalization (early vs late care-seekers)
nai.pars.low <- data.frame(
                    effect.level = paste0(c("low", "intermediate", "high"), " NAI effect"),
                    nai.effect.early = c(0.2, 0.4, 0.70),
                    nai.effect.late  = c(0.0, 0.2, 0.35)
                )

nai.pars.high <- nai.pars.low

# join scenarios together
scenario.grid.low <- expand_grid(testing.pars.low, prescribing.pars.low, nai.pars.low,.name_repair = "unique") %>%
                        mutate(scenario = paste0(test.level,  ", ", prescribe.level,  ", ", effect.level)) %>% 
                        select(-contains("level")) 


scenario.grid.high <- expand_grid(testing.pars.high, prescribing.pars.high, nai.pars.high,.name_repair = "unique") %>%
                        mutate(scenario = paste0(test.level,  ", ", prescribe.level,  ", ", effect.level)) %>% 
                        select(-contains("level")) 

# join with other baseline parameters
pars.low  <- pars.range %>% filter(risk == "low")  %>% expand_grid(., scenario.grid.low) 
pars.high <- pars.range %>% filter(risk == "high") %>% expand_grid(., scenario.grid.high) 

pars <- rbind(pars.low, pars.high) 
```

```{r}
out.initial <- get_outcomes(pars) %>% 
                    mutate(scenario.all = scenario) %>%
                    separate(scenario.all, c("testing", "prescribing", "effectiveness"), ", ") 
```


### Initial scenarios vs baseline

Simulate the initial scenarios above relative to a theoretical baseline scenario in which antivirals ('NAIs') are *not* used (can be modeled by setting NAI effectiveness to 0)

```{r}
pars.worst <- pars %>% filter(scenario == "high testing, high prescribing, high NAI effect") %>%
                # setting nai effectiveness to 0 means nothing else matters -- all boxes have same hospitalization risk
                mutate(nai.effect.early = 0, 
                       nai.effect.late  = 0,
                       scenario = "worst")  
```

```{r}
# simulate baseline scenario
out.baseline <- get_outcomes(pars.worst) 

# join with output from baseline scenarios
tmp.baseline <- bind_rows(out.initial, out.baseline) %>% 
                    select(sim, risk, age, scenario, hosp.total) %>%
                    spread(risk, hosp.total) %>% mutate(all = high + low) %>%
                    gather(risk, hosp.total, high, low, all) %>% spread(age, hosp.total) %>%
                    mutate(all = `0-4`+ `5-17`+ `18-49`+ `50-64`+ `65+`) %>%
                    gather(age, hosp.total, `0-4`:all) %>%
                    spread(scenario, hosp.total) 

# calculate changes relative to the baseline scenario (absolute and percent)
change.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x))) 
change.perc.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x) / worst * 100)) 
```

Total hospitalizations in the baseline scenario without antivirals

```{r}
change.worst %>% filter(risk == "all", age == "all") %>% do(., get_summary(., var = "worst", groups = NULL))
```


```{r}
# Total hospitalizations averted and number needed to treat --------
total_averted <- change.worst %>% filter(risk == "all", age == "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ") %>%
    do(., get_summary(., var = "change", groups = c("testing", "prescribing", "effectiveness")))

total_perc_averted <- change.perc.worst %>% filter(risk == "all", age == "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ") %>%
    do(., get_summary(., var = "change", groups = c("testing", "prescribing", "effectiveness")))

nai.used <- out.initial %>% group_by(sim, testing, prescribing, effectiveness) %>%
    summarize(nai.total = sum(nai.total), .groups = "drop") 
```

Total averted hospitalizations in each initial scenario

```{r}
total_averted %>% 
    mutate(effectiveness = gsub(effectiveness, pattern = " NAI effect", replacement = ""))
```

Total percent averted hospitalizations in each initial scenario

```{r}
total_perc_averted %>% 
    mutate(effectiveness = gsub(effectiveness, pattern = " NAI effect", replacement = ""))
```

Total antivirals used 

```{r}
nai.used %>% filter(effectiveness == "high NAI effect") %>%
    do(., get_summary(., var = "nai.total", groups = c("testing", "prescribing"))) 
```

NNT

```{r}
change.worst %>% filter(risk == "all", age == "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ") %>%
    left_join(nai.used) %>%
    mutate(nnt = nai.total/change) %>%
    do(., get_summary(., var = "nnt", groups = c("effectiveness", "testing", "prescribing"))) %>%
    mutate(effectiveness = gsub(effectiveness, pattern = " NAI effect", replacement = ""))
```


```{r, fig.height = 10, fig.width = 10}
# plot total averted
p_total <- total_averted %>% 
    mutate(effectiveness = gsub(effectiveness, pattern = " NAI effect", replacement = "")) %>%
    ggplot(aes(x = effectiveness, y = mean, color = effectiveness)) + 
    geom_point(size = ptsize) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), width = barsize - 0.2, linewidth = barwdth) +
    facet_grid(testing ~ prescribing) +
    get_theme(txt =11) + 
    scale_y_continuous(labels = comma) +
    scale_color_manual(guide = "none", values = cols) +
    labs(x = "Antiviral effectiveness", y = "Number of hospitalizations averted relative to baseline")

# plot NNT: need to combine # hosps averted with # NAIs used
p_nnt <- change.worst %>% filter(risk == "all", age == "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ") %>%
    left_join(nai.used) %>%
    mutate(nnt = nai.total/change) %>%
    do(., get_summary(., var = "nnt", groups = c("effectiveness", "testing", "prescribing"))) %>%
    mutate(effectiveness = gsub(effectiveness, pattern = " NAI effect", replacement = "")) %>%
    ggplot(aes(x = effectiveness, y = mean, color = effectiveness)) + 
    geom_point(size = ptsize) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), width = barsize - 0.2, linewidth = barwdth) +
    facet_grid(testing ~ prescribing) +
    get_theme(txt = 11) + 
    scale_y_continuous(labels = comma) +
    scale_color_manual(guide = "none", values = cols) +
    labs(x = "Antiviral effectiveness", y = "Number of prescriptions per averted hospitalization")

p_total / p_nnt + plot_annotation(tag_levels = 'A')
ggsave(here("figs", "figS2.pdf"), width = 8, height = 10)
```

Partitioned by age and risk group:

Total hospitalizations averted 

```{r}
tmp <- change.worst %>% filter(age != "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ") %>%
    filter(effectiveness == "high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "testing", "prescribing", "effectiveness"))) %>%
    mutate(age = factor(age, levels = age.grps)) %>% arrange(age)

tmp %>% filter(risk == "high", testing == "high testing", prescribing == "high prescribing") 
tmp %>% filter(risk == "high", testing == "low testing", prescribing == "low prescribing") 
```

```{r, results='hide'}
cairo_pdf(here("figs", "fig3.pdf"), width = 8, height = 7)
tmp %>%
    ggplot(aes(x = age, y = mean, color = risk)) + 
    geom_point(position = position_dodge(width = dodgewdth), size = ptsize) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth), width = barsize, linewidth = barwdth) +
    facet_grid(testing ~ prescribing) +
    get_theme(txt = 11, legend.position = "top") + 
    scale_y_log10(labels = comma) + scale_x_discrete(labels = age_labels) +
    scale_color_manual("Risk", values = cols_risk) +
    labs(x = "Age (years)", y = "Number of hospitalizations averted relative to baseline")
dev.off()
```

Total percent of hospitalizations averted and NNT

```{r}
tmp2 <- change.perc.worst %>% filter(age != "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ")  %>%
    filter(effectiveness == "high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "testing", "prescribing", "effectiveness"))) %>%
    mutate(age = factor(age, levels = age.grps))

tmp2 %>% arrange(age) %>% filter(risk == "high", testing == "high testing", prescribing == "high prescribing") 
```

```{r, results='hide'}
cairo_pdf(here("figs", "figS3.pdf"), width = 8, height = 7)
tmp2 %>%
    ggplot(aes(x = age, y = mean, color = risk)) + 
    geom_point(position = position_dodge(width = dodgewdth), size = ptsize) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth), width = barsize, linewidth = barwdth) +
    facet_grid(testing ~ prescribing) +
    get_theme(txt = 11, legend.position = "top") + 
    scale_y_continuous(labels = comma) + scale_x_discrete(labels = age_labels) +
    scale_color_manual("Risk", values = cols_risk) +
    labs(x = "Age (years)", y = "Percent of hospitalizations averted relative to baseline")
dev.off()

# NNT
nai.used.age.risk0 <- out.initial %>% 
    select(sim, risk, age, scenario, nai.total) %>% 
    spread(risk, nai.total) %>% mutate(all = high + low) %>%
    gather(risk, nai.total, high, low, all) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ")  %>%
    group_by(sim, age, risk, testing, prescribing, effectiveness) %>%
    summarize(nai.total = sum(nai.total), .groups = "drop") 

# quoted values
tmp <- change.worst %>% filter(age != "all") %>%
    gather(scenario, change, contains("testing")) %>%
    separate(scenario, c("testing", "prescribing", "effectiveness"), ", ") %>%
    left_join(nai.used.age.risk0) %>%
    filter(effectiveness == "high NAI effect") %>%
    mutate(nnt = nai.total/change,
           age = factor(age, levels = age.grps)) %>%
    do(., get_summary(., var = "nnt", groups = c("risk", "age", "testing", "prescribing", "effectiveness"))) 

tmp %>% filter(risk == "high", age %in% c("0-4", "65+") )

# plot
cairo_pdf(here("figs", "figS4.pdf"), width = 8, height = 7)
tmp %>%
    ggplot(aes(x = age, y = mean, color = risk)) + 
    geom_point(position = position_dodge(width = dodgewdth), size = ptsize) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth), width = barsize, linewidth = barwdth) +
    facet_grid(testing ~ prescribing) +
    get_theme(txt = 11, legend.position = "top") + 
    scale_y_log10(labels = comma) + scale_x_discrete(labels = age_labels) +
    scale_color_manual("Risk", values = cols_risk) +
    labs(x = "Age (years)", y = "Number of prescriptions per averted hospitalization")
dev.off()
```


### Alternative scenarios: increased individual-level parameters

Increase parameters that represent 

* Increased % seeking care
* Increased % of care-seekers doing so within 2d
* Increased antiviral compliance 

```{r}
tmp.pars <- pars %>% filter(scenario == "high testing, high prescribing, high NAI effect") 

# alter input parameters: assume increase of 50% for each
pars.alt.cs1 <- tmp.pars %>% mutate(seek.care = seek.care * 1.5)
pars.alt.cs2 <- tmp.pars %>% mutate(seek.care.early = seek.care.early * 1.5)
pars.alt.cs3 <- tmp.pars %>% mutate(nai.compliance = nai.compliance * 1.5)
pars.alt.cs4 <- tmp.pars %>% mutate(seek.care = seek.care * 1.5, 
                                seek.care.early = seek.care.early * 1.5, 
                                nai.compliance = nai.compliance * 1.5)

ind_levels <- c("initial", "seek care", "seek care early", "comply", "seek care, seek care early, and comply")

# join all new scenarios together
pars.alt.cs <- bind_rows(tmp.pars, pars.alt.cs1, pars.alt.cs2, pars.alt.cs3, pars.alt.cs4) %>% 
                mutate(alternative = rep(ind_levels, each = nrow(tmp.pars))) %>%
    # make sure new inputs don't exceed 1
    mutate(across(c(seek.care, seek.care.early, nai.compliance), function(x) pmin(x, 1)))

# simulate new outputs
out.alt.cs <- get_outcomes(pars.alt.cs) 
```

Plot new parameter distributions

```{r, fig.height = 6, results='hide'}
pars_long <- 
    pars.alt.cs %>% distinct(sim, risk, age, seek.care, seek.care.early, nai.compliance, alternative) %>%
    gather(param, value, seek.care, seek.care.early, nai.compliance) %>%
    filter( alternative == "initial" | (param == "seek.care" & alternative ==  "seek care") | 
            (param == "seek.care.early" & alternative ==  "seek care early") | 
                (param == "nai.compliance" & alternative ==  "comply")) %>%
    mutate(scenario = ifelse(alternative == "initial", alternative, "50% increase"))

# add better parameter names for plotting
cairo_pdf(here("figs", "figS1.pdf"), width = 8, height = 6)
rename_pars(pars_long) %>%
    mutate(risk = ifelse(risk == "high", "High risk", "Low risk")) %>%
    ggplot(aes(x = age, y = value, fill = scenario)) + 
    geom_violin(alpha = 0.5, position = position_dodge(width = 0.75)) + 
    geom_boxplot(alpha = 0.75, width = 0.1, position = position_dodge(width = 0.75)) +
    scale_fill_manual("Value", values = cols_alt) +
    ylim(c(0, 1)) +
    facet_grid(parname ~ risk) + get_theme(txt = 11) +
    scale_x_discrete(labels = age_labels) +
    labs(x = "Age (years)", y = "Sampled value")
dev.off()
```

Compare to baseline scenario, assuming NAI effectiveness is at its highest value.

```{r}
n <- length(unique(pars.alt.cs$alternative))

# join output with output for baseline scenarios
tmp.baseline <- out.baseline %>% uncount(n) %>% 
    mutate(alternative = rep(unique(pars.alt.cs$alternative), times = nrow(out.baseline))) %>%
    bind_rows(., out.alt.cs) %>% select(sim, risk, age, scenario, alternative, hosp.total) %>%
    spread(risk, hosp.total) %>% mutate(all = high + low) %>%
    gather(risk, hosp.total, high, low, all) %>% spread(scenario, hosp.total) 

# calculate change relative to baseline scenarios
change.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x))) 
change.perc.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x) / worst * 100)) 
```

```{r, fig.width = 10, fig.height = 8, results='hide'}
# Plot compared to worst case (absolute differences)
p_alt_num <- 
    change.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    # add formal age group labels
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = ind_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.35), size = ptsize - 0.3) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.35), width = barsize + 0.3, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11, legend.position = "top") + 
    scale_y_log10(labels = comma) +
    scale_color_manual(NULL, values = cols_alt) +
    labs(x = "Risk stratification", y = "Number of hospitalizations\naverted relative to baseline")

# Plot compared to worst case (% differences)
p_alt_perc <- 
    change.perc.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = ind_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.35), size = ptsize - 0.3) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.35), width = barsize + 0.3, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11) + 
    scale_y_continuous(labels = comma) + 
    scale_color_manual(guide = "none", values = cols_alt) +
    labs(x = "Risk stratification", y = "Percent of hospitalizations\naverted relative to baseline")

# Plot NNT
nai.used.age.risk <- out.alt.cs %>% 
    select(sim, alternative, risk, age, nai.total) %>% 
    spread(risk, nai.total) %>% mutate(all = high + low) %>%
    gather(risk, nai.total, high, low, all) %>%
    group_by(sim, age, risk, alternative) %>%
    summarize(nai.total = sum(nai.total), .groups = "drop") 

p_alt_nnt <- 
    change.worst %>% filter(age != "all") %>%
    rename(change = 'high testing, high prescribing, high NAI effect') %>%
    left_join(nai.used.age.risk) %>%
    mutate(nnt = nai.total/change) %>%
    do(., get_summary(., var = "nnt", groups = c("risk", "age", "alternative")) ) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = ind_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.35), size = ptsize - 0.3) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.35), width = barsize + 0.3, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11) + 
    scale_y_log10(labels = comma) + 
    scale_color_manual(guide = "none", values = cols_alt) +
    labs(x = "Risk stratification", y = "Number of prescriptions\nper averted hospitalization")

cairo_pdf(file = here("figs", "fig4.pdf"), width = 8.25, height = 9)
p_alt_num / p_alt_perc / p_alt_nnt + plot_annotation(tag_levels = 'A') + plot_layout(heights = c(1, 0.9, 0.9))
dev.off()
```

Total averted 

```{r}
change.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(alternative == "seek care, seek care early, and comply", risk == "high")
```

Percent averted 

```{r}
change.perc.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(alternative == "seek care, seek care early, and comply", risk == "high")
```


### Alternative scenarios: increased clinical-level parameters

* Increased testing (regardless of timing)
* Increased prescribing among care-seekers with a positive test (regardless of timing)

```{r}
tmp.pars <- pars %>% filter(scenario == "high testing, high prescribing, high NAI effect") 

pars.alt1 <- tmp.pars %>% mutate(prescribed.nai.pos.early = ifelse(risk == "high", 0.75, 0.37), 
                                 prescribed.nai.pos.late  = ifelse(risk == "high", 0.30, 0.15))
pars.alt2 <- tmp.pars %>% mutate(perc.tested.early = 0.75, perc.tested.late = 0.15)
pars.alt3 <- tmp.pars %>% mutate(prescribed.nai.pos.early = ifelse(risk == "high", 0.75, 0.37), 
                                 prescribed.nai.pos.late  = ifelse(risk == "high", 0.30, 0.15),
                                 perc.tested.early = 0.75, perc.tested.late = 0.15)

clinical_levels <- c("initial", "prescribing", "testing", "testing and prescribing")

pars.alt <- bind_rows(tmp.pars, pars.alt1, pars.alt2, pars.alt3) %>% 
                mutate(alternative = rep(clinical_levels, each = nrow(tmp.pars)))

out.alt <- get_outcomes(pars.alt) 
```

Compare to best and worst-case scenarios, assuming NAI effectiveness is at its highest value.

```{r}
n <- length(unique(pars.alt$alternative))

tmp.baseline <- out.baseline %>% uncount(n) %>% 
    mutate(alternative = rep(unique(pars.alt$alternative), times = nrow(out.baseline))) %>%
    bind_rows(., out.alt) %>% select(sim, risk, age, scenario, alternative, hosp.total) %>%
    spread(risk, hosp.total) %>% mutate(all = high + low) %>%
    gather(risk, hosp.total, high, low, all) %>% spread(scenario, hosp.total) 

change.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x))) 
change.perc.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x) / worst * 100)) 
```

```{r, fig.height = 8, fig.width = 10, results='hide'}
# Compare worst-case to baseline scenarios (absolute difference)
p_alt2_num <- 
    change.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = clinical_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.35), size = ptsize - 0.25) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.35), width = barsize + 0.3, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11, legend.position = "top") + 
    scale_y_log10(labels = comma) +
    scale_color_manual("Scenario", values = cols_alt[c(1, 2, 4,5)]) +
    labs(x = "Risk stratification", y = "Number of hospitalizations\naverted relative to baseline")

# Compare worst-case to baseline scenarios (% difference)
p_alt2_perc <- 
    change.perc.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = clinical_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.35), size = ptsize - 0.25) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.35), width = barsize + 0.3, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11) + 
    scale_y_continuous(labels = comma) + 
    scale_color_manual(guide = "none", values = cols_alt[c(1, 2, 4,5)]) +
    labs(x = "Risk stratification", y = "Percent of hospitalizations\naverted relative to baseline")

# Plot NNT
nai.used.age.risk2 <- out.alt %>% 
    select(sim, alternative, risk, age, nai.total) %>% 
    spread(risk, nai.total) %>% mutate(all = high + low) %>%
    gather(risk, nai.total, high, low, all) %>%
    group_by(sim, age, risk, alternative) %>%
    summarize(nai.total = sum(nai.total), .groups = "drop") 

p_alt2_nnt <- 
    change.worst %>% filter(age != "all") %>%
    rename(change = 'high testing, high prescribing, high NAI effect') %>%
    left_join(nai.used.age.risk2) %>%
    mutate(nnt = nai.total/change) %>%
    do(., get_summary(., var = "nnt", groups = c("risk", "age", "alternative")) ) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = clinical_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.35), size = ptsize - 0.25) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.35), width = barsize + 0.3, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11) + 
    scale_y_log10(labels = comma) + 
    scale_color_manual(guide = "none", values = cols_alt[c(1, 2, 4,5)]) +
    labs(x = "Risk stratification", y = "Number of prescriptions\nper averted hospitalization")

cairo_pdf(here("figs", "figS5.pdf"), width = 8.25, height = 9)
p_alt2_num / p_alt2_perc / p_alt2_nnt + plot_annotation(tag_levels = 'A') + plot_layout(heights = c(1, 0.9, 0.9))
dev.off()
```

Total averted

```{r}
change.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(risk == "high", alternative == "testing and prescribing")
``` 
    
Percent averted

```{r}
change.perc.worst %>% 
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(risk == "high", alternative == "testing and prescribing")
``` 
 
### Alternative scenarios: everything

```{r}
pars0 <- pars %>% filter(scenario == "high testing, high prescribing, high NAI effect") 

pars.all <- pars0 %>% mutate(prescribed.nai.pos.early = ifelse(risk == "high", 0.75, 0.37), 
                             prescribed.nai.pos.late  = ifelse(risk == "high", 0.30, 0.15),
                             perc.tested.early = 0.75, 
                             perc.tested.late = 0.15,
                             seek.care = seek.care * 1.5, 
                             seek.care.early = seek.care.early * 1.5, 
                             nai.compliance = nai.compliance * 1.5)

all_levels <- c("initial", "increase in all 5 parameters")

pars.all <- bind_rows(pars0, pars.all) %>% 
            mutate(alternative = rep(all_levels, each = nrow(pars0))) %>%
            # make sure new inputs don't exceed 1
            mutate(across(c(seek.care, seek.care.early, nai.compliance), function(x) pmin(x, 1)))

out.all <- get_outcomes(pars.all) 

n <- length(unique(pars.all$alternative))

tmp.baseline <- out.baseline %>% uncount(n) %>% 
    mutate(alternative = rep(unique(pars.all$alternative), times = nrow(out.baseline))) %>%
    bind_rows(., out.all) %>% select(sim, risk, age, scenario, alternative, hosp.total) %>%
    spread(risk, hosp.total) %>% mutate(all = high + low) %>%
    gather(risk, hosp.total, high, low, all) %>% 
    spread(age, hosp.total) %>%
    mutate(all = `0-4`+ `5-17`+ `18-49`+ `50-64`+ `65+`) %>%
    gather(age, hosp.total, `0-4`:all) %>%
    spread(scenario, hosp.total) 

change.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x))) 
change.perc.worst <- tmp.baseline %>% mutate(across(contains("testing"), function(x) (worst - x) / worst * 100)) 
```

```{r, fig.height = 8, fig.width = 10, results='hide'}
# Compare absolute difference
p_alt3_num <- 
    change.worst %>% filter(age != "all") %>%
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = all_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.1), size = ptsize - 0.25) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.1), width = barsize + 0.1, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11, legend.position = "top") + 
    scale_y_log10(labels = comma) + 
    scale_color_manual("Scenario", values = cols_alt[c(1, 5)]) + 
    labs(x = "Risk stratification", y = "Number of hospitalizations\naverted relative to baseline")

# Compare % difference
p_alt3_perc <- change.perc.worst %>% filter(age != "all") %>%
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = all_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.1), size = ptsize - 0.25) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.1), width = barsize + 0.1, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11) + 
    scale_y_continuous(labels = comma) + 
    scale_color_manual(guide = "none", values = cols_alt[c(1, 5)]) + 
    labs(x = "Risk stratification", y = "Percent of hospitalizations\naverted relative to baseline")

# Plot NNT
nai.used.age.risk3 <- out.all %>% 
    select(sim, alternative, risk, age, nai.total) %>% 
    spread(risk, nai.total) %>% mutate(all = high + low) %>%
    gather(risk, nai.total, high, low, all) %>%
    group_by(sim, age, risk, alternative) %>%
    summarize(nai.total = sum(nai.total), .groups = "drop") 

p_alt3_nnt <- 
    change.worst %>% filter(age != "all") %>%
    rename(change = 'high testing, high prescribing, high NAI effect') %>%
    left_join(nai.used.age.risk3) %>%
    mutate(nnt = nai.total/change) %>%
    do(., get_summary(., var = "nnt", groups = c("risk", "age", "alternative")) ) %>%
    mutate(age_names = factor(paste(age_labels[match(age, age.grps)], "years"), levels = paste(age_labels, "years")),
           alternative = factor(alternative, levels = all_levels)) %>%
    ggplot(aes(x = risk, y = mean, color = alternative)) + 
    geom_point(position = position_dodge(width = dodgewdth + 0.1), size = ptsize - 0.25) + 
    geom_errorbar(aes(ymin = lower, ymax = upper), 
                  position = position_dodge(width = dodgewdth + 0.1), width = barsize + 0.1, linewidth = barwdth - 0.1) +
    facet_grid(~ age_names) + get_theme(txt = 11) + 
    scale_y_log10(labels = comma) + 
    scale_color_manual(guide = "none", values = cols_alt[c(1, 5)]) +
    labs(x = "Risk stratification", y = "Number of prescriptions\nper averted hospitalization")

cairo_pdf(here("figs", "figS6.pdf"), width = 8.25, height = 9)
p_alt3_num / p_alt3_perc / p_alt3_nnt + plot_annotation(tag_levels = 'A') + plot_layout(heights = c(1, 0.9, 0.9))
dev.off()
```

Total averted 
```{r}
change.worst %>% filter(age != "all") %>%
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(risk == "high", alternative == "increase in all 5 parameters")
```

Percent averted 
```{r}
change.perc.worst %>% filter(age != "all") %>%
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(risk == "high", alternative == "increase in all 5 parameters")
```

Across all age groups and risk stratifications:

Total and percent averted

```{r}
change.worst %>% filter(age == "all") %>%
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(risk == "all", alternative == "increase in all 5 parameters")

change.perc.worst %>% filter(age == "all") %>%
    rename(change = "high testing, high prescribing, high NAI effect") %>%
    do(., get_summary(., var = "change", groups = c("risk", "age", "alternative"))) %>%
    filter(risk == "all", alternative == "increase in all 5 parameters")
```

NNT
 
```{r}
 nai.used.all <- out.all %>% 
    select(sim, alternative, risk, age, nai.total) %>% 
    spread(risk, nai.total) %>% mutate(all = high + low) %>%
    gather(risk, nai.total, high, low, all) %>%
    group_by(sim, risk, alternative) %>%
    summarize(nai.total = sum(nai.total), .groups = "drop") %>%
    filter(risk == "all")

change.worst %>% filter(age == "all", risk == "all") %>%
    rename(change = 'high testing, high prescribing, high NAI effect') %>%
    left_join(nai.used.all) %>%
     mutate(nnt = nai.total/change) %>%
    do(., get_summary(., var = "nnt", groups = c("risk", "age", "alternative")) )
```