5_Table_2_outcomes.Rmd

---
author: Leon Di Stefano
date: "`r Sys.Date()`"
output: 
    html_document:
      df_print: paged
params:
  fit_name: "main_fit"
  outcome_min: 28
  outcome_max: 35
title: "`r paste('Results tables, inc. secondary and safety outcomes', params$out_stub, params$fit_name, sep = '-')`"

---

```{r echo=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(here)
here::i_am(file.path("hcq_pooling_analysis", "5_Table_2_outcomes.Rmd"))
source(here("hcq_pooling_analysis", "common.R"))

require(gt)
require(MASS, exclude = "select") # clashes with dplyr
require(broom)
output_dir <- here("hcq_pooling_analysis", "output")

out_stub         <- paste(params$outcome_min, params$outcome_max, sep = '-')
output_dir       <- here("hcq_pooling_analysis", "output", out_stub)
output_model_dir <- file.path(output_dir, params$fit_name)
fit_file         <- file.path(output_model_dir, paste0(params$fit_name, ".rds"))

subgroup_dir <- file.path(output_model_dir, "subgroup_effects")
```

Load the main data table:

```{r}
patients <- 
  read_rds(file.path(output_dir, "patients.rds")) %>%
  filter(passed_screen == 1) %>%
  mutate(treat = factor(treat, levels = c("no_HCQ", "HCQ")))
```

## 6 effects

-   Primary outcome
    -   Plug-in
    -   Model-adjusted
    
(These have been computed in another script)
    
-   Secondary outcomes
    1.  Hospitalization length of stay
    2.  Duration of mechanical ventilation
    -   ~~(Worst NCOSS score before d28?) "To the extent possible, we will also examine the lowest NCOSS score achieved prior to day 28-30."~~
-   Safety outcomes
    3.  Number of AEs and SAEs overall and by category
    4.  All-cause mortality
    
### Ineligible sites

It doesn't make sense to compute these effects for WashU or Bassett (only one treatment value).

## Primary outcome 

### Plot (ORCHID-style)

```{r}
data_tbl <- read_rds(file.path(output_dir, "data_tbl.rds"))
```

```{r fig.height=2.5, fig.width=7}
orchid_plot_line_tbl <-
  data_tbl %>% 
      filter(!is.na(niaid_outcome)) %>% 
      count(niaid_outcome, treat) %>% 
      group_by(treat) %>% 
      arrange(rev(niaid_outcome)) %>% 
      mutate(
        cumulative_prop = cumsum(n)/sum(n), 
        treat = (as.numeric(treat) - 1.5)*.5 + 1.5) %>%
      filter(niaid_outcome != "1 - death")

orchid_style_barplot <-
  ggplot(
  data_tbl %>% filter(!is.na(niaid_outcome)), 
  aes(treat, fill = niaid_outcome)) + 
  geom_bar(width = .5, position = "fill",
           color = "white", size = .2) + 
  scale_fill_manual(values = niaid_colors) + 
  geom_line(
    data = orchid_plot_line_tbl, 
    aes(treat, cumulative_prop, group = niaid_outcome), 
    color = "grey50", 
    linetype = "dashed") + 
  coord_flip() + 
  scale_x_discrete(
    expand = expansion(.3),
    labels = c(
      HCQ = "HCQ or CQ",
      no_HCQ = "placebo or\nstandard care"
    )) + 
  scale_y_continuous(
    labels = scales::percent, 
    breaks = (0:5)/5, 
    expand = expansion(.01)) + 
  theme(axis.line.y = element_blank(), 
        axis.ticks.y = element_blank()) +
  xlab("") + ylab("")

orchid_style_barplot

ggsave(file.path(output_dir, "orchid_style_barplot.png"), orchid_style_barplot,
       height=2.5, width=7)
ggsave(file.path(output_dir, "orchid_style_barplot.svg"), orchid_style_barplot,
       height=2.5, width=7)
```

Could be good to add *missingness* to this.

Adding counts to the bars:

```{r fig.height=2.5, fig.width=7}
bar_label_tbl <-
  data_tbl %>% 
  filter(!is.na(niaid_outcome)) %>%
  count(treat, niaid_outcome) %>%
  add_count(treat, wt = n, name = "n_treated") %>%
  mutate(prop = n/n_treated,
         pct = paste0(round(100 * prop), "%"),
         label = paste0(pct, "\n(n=", n, ")"))

bar_treat_label_tbl <-
  data_tbl %>%
  mutate(missing_outcome = is.na(niaid_outcome)) %>%
  count(treat, missing_outcome, name = "n_missing") %>%
  add_count(treat, wt = n_missing, name = "n_treated") %>%
  filter(missing_outcome) %>%
  mutate(treat_label = paste0(
    ifelse(treat == "HCQ", "HCQ/CQ", "placebo or\nstandard care"),
    "\n",
    "n=", n_treated,
    "\n(", n_missing," missing)"
  ))

treat_labels <- bar_treat_label_tbl$treat_label
names(treat_labels) <- bar_treat_label_tbl$treat

bar_line_tbl <-
  data_tbl %>% 
      filter(!is.na(niaid_outcome)) %>% 
      count(niaid_outcome, treat) %>% 
      group_by(treat) %>% 
      arrange(rev(niaid_outcome)) %>% 
      mutate(
        cumulative_prop = cumsum(n)/sum(n), 
        treat = (as.numeric(treat) - 1.5)*.5 + 1.5) %>%
      filter(niaid_outcome != "1 - death")

orchid_style_barplot_annotated <-
  ggplot(
  bar_label_tbl, 
  aes(treat, prop, fill = niaid_outcome)) + 
  geom_col(width = .5, position = "fill",
           color = "white", size = .2) + 
  scale_fill_manual(values = niaid_colors) + 
  geom_line(
    data = orchid_plot_line_tbl,
    aes(treat, cumulative_prop, group = niaid_outcome),
    color = "grey50",
    linetype = "dashed") +
  coord_flip() +
  scale_x_discrete(
    expand = expansion(.3),
    labels = treat_labels) + 
  scale_y_continuous(
    labels = scales::percent,
    breaks = (0:5)/5,
    expand = expansion(.01)) +
  theme(axis.line.y = element_blank(), 
        axis.ticks.y = element_blank()) +
  xlab("") + ylab("") + 
  geom_text(aes(label = label), 
            data = bar_label_tbl, 
            position = position_stack(vjust = .5),
            size = 4)

orchid_style_barplot_annotated

ggsave(file.path(output_dir, "orchid_style_barplot_annotated.png"), orchid_style_barplot_annotated,
       height=2.5, width=7)
ggsave(file.path(output_dir, "orchid_style_barplot_annotated.svg"), orchid_style_barplot_annotated,
       height=2.5, width=7)
```

### Table

```{r}
outcome_table <- 
  #read_rds(file.path(output_dir, "standardized_effects/overall_analysis.rds")$summary_table
  read_rds(file.path(subgroup_dir, paste0(params$fit_name, "overall_analysis.rds")))$summary_table
  
```

### 1) Hospitalization length of stay

Missing values:

```{r}
sum(is.na(patients$hosp_days_calc))
```

Exploratory plot:

```{r fig.height=9, fig.width=6}
patients %>% 
  ggplot(aes(hosp_days_calc)) +
  facet_rep_grid(siteid~treat, scales = "free_y") +
  geom_histogram(binwidth = 1) +
  scale_y_continuous(expand = expansion(0))
```

Table:

```{r}
table(patients$hosp_days_calc, patients$siteid)
```

This suggests

-   Using median/nonparametric test
-   Capping at 28d for uniformity/compliance:

```{r}
patients$hosp_days_calc_cap <-
  pmin(patients$hosp_days_calc, 28)
```

```{r}
wilcox.test(patients$hosp_days_calc_cap ~ patients$treat)
```

```{r}
hosp_days_tests <-
  bind_rows(
    patients,
    patients %>% mutate(siteid = "OVERALL")
  ) %>%
  group_by(siteid) %>%
  summarise(
    test = list(possibly(wilcox.test, NA)(hosp_days_calc_cap ~ treat))) %>% 
  mutate(
    p_val = map_dbl(
      test,
      (function(x) {
        if (class(x) == "htest")  { x$p.value } else { NA }})
      ))

hosp_days_medians <-
  bind_rows(
    patients,
    patients %>% mutate(siteid = "OVERALL")
  ) %>%
  group_by(siteid, treat) %>%
  summarise(
    median_hosp_days = median(hosp_days_calc_cap, na.rm = TRUE),
    n = n(),
    # n_missing = sum(is.na(hosp_days_calc_cap))
    )

hosp_days_missing <-
  bind_rows(
    patients,
    patients %>% mutate(siteid = "OVERALL")
  ) %>%
  group_by(siteid) %>%
  summarise(n_missing = sum(is.na(hosp_days_calc_cap)))

(hosp_days_summary_table <-
  hosp_days_medians %>% 
  pivot_wider(
    names_from = treat, 
    values_from = c(median_hosp_days, n)) %>%
  left_join(hosp_days_tests %>% select(-test)) %>%
  left_join(hosp_days_missing))
```

### For table 2

```{r}
t2_hosp_days <-
  hosp_days_medians %>%
  ungroup() %>%
  mutate(col = paste(siteid, treat),
         measure = "Days of hospitalization between enrollment and d28 (median)",
         value = median_hosp_days
         ) %>%
  select(col, measure, value) %>%
  pivot_wider(names_from = col, values_from = value) %>%
  mutate(n_missing = sum(is.na(patients$hosp_days_calc_cap)))

t2_hosp_days
```

### 2) Duration of mechanical ventilation

```{r}
sum(is.na(patients$mechvent_days_calc))
```

```{r fig.height=9, fig.width=6}
patients %>% 
  ggplot(aes(mechvent_days_calc)) +
  facet_rep_grid(siteid~treat, scales = "free_y") +
  geom_histogram(binwidth = 1) +
  scale_y_continuous(expand = expansion(0))
```

Looks like it's really a question of zero vs. nonzero:

```{r}
table(patients$mechvent_days_calc, patients$siteid)
```

```{r fig.height=9, fig.width=6}
mech_vent_summary_long <- patients %>% 
  bind_rows(patients %>% mutate(siteid = "OVERALL")) %>%
  mutate(mech_vent = (mechvent_days_calc > 0)) %>%
  add_count(siteid, treat, name = "n_treated") %>%
  add_count(siteid, treat, mech_vent) %>%
  group_by(siteid, treat) %>%
  summarise(
    n_mech_vent = sum(mech_vent, na.rm = TRUE),
    n_treated      = sum(!is.na(mech_vent)),
    p_mech_vent = n_mech_vent/n_treated,
    n_missing      = sum(is.na(mech_vent))
  )

mech_vent_summary_long
```

```{r}
(mech_vent_ests <-
  bind_rows(
    patients,
    patients %>% mutate(siteid = "OVERALL")
  ) %>%
  group_by(siteid) %>%
  summarise(
    estimate = list(possibly(tidy, NA)(glm(I(mechvent_days_calc %in% c(0,1))~treat, 
                                       family = binomial()),
                    conf.int   = TRUE,
                    conf.level = .95))) %>%
  unnest(estimate) %>%
  select(-estimate) %>%
  filter(term %in% c("treatHCQ", NA)) %>%
  select(-term, -std.error) %>%
  mutate(or_hat = exp(statistic),
         or_hat_025 = exp(conf.low),
         or_hat_975 = exp(conf.high)) %>%
  select(-statistic, -conf.low, -conf.high))
```

### For Table 2

```{r}
t2_mech_vent <-
  mech_vent_summary_long %>%
  ungroup() %>%
  mutate(col = paste(siteid, treat),
         measure = "Patients on mechanical ventilation between enrollment and day 28 (count (percent))",
         value = str_c(n_mech_vent, " (", scales::percent(p_mech_vent, accuracy = 1), ")"),
         n_missing = sum(is.na(patients$mechvent_days_calc))) %>%
  select(col, measure, value, n_missing) %>%
  pivot_wider(names_from = col, values_from = value)

t2_mech_vent
```

### 3) Number of #AEs, #SAEs, and by category, per person

```{r}
patients %>%
  group_by(siteid, treat) %>%
  summarise(tot_sae = sum(num_sae))
```

```{r}
patients %>%
  group_by(siteid, treat) %>%
  summarise(tot_sae = sum(num_sae)) %>%
  ggplot(aes(siteid, tot_sae, fill = treat)) +
  geom_col(position = "dodge") +
  coord_flip()
```

Because of the T/C imbalance across groups, need to make this per person, i.e. include an offset in the Poisson regression:

```{r}
(tot_saes_fits <- 
  patients %>%
  bind_rows(patients %>% mutate(siteid = "OVERALL")) %>%
  group_by(siteid, treat) %>%
  summarise(tot_sae = sum(num_sae),
            n_patients = n()) %>%
  group_by(siteid) %>%
  summarise(
    fit = list(
      possibly(tidy, NA)(glm(tot_sae ~ offset(log(n_patients)) + treat, family = poisson()),
                         conf.int = TRUE,
                         conf.level = .95))
  ) %>%
  unnest(fit) %>%
  filter(term %in% c("treatHCQ", NA)) %>%
  mutate(rr_hat = exp(estimate),
         rr_hat_025 = exp(conf.low),
         rr_hat_975 = exp(conf.high)) %>%
   select(-estimate, -conf.low, -conf.high, -statistic, -std.error, -fit))
```

### For Table 2, all in one

```{r}
ae_summary_function <- function(x) {
  num_per_patient <- mean(x, na.rm = TRUE)
  count           <- sum(x, na.rm = TRUE)
  
  if(all(is.na(x))) {return(NA)}
  
  if(count == 0) {return("0")}
  
  return(paste0(count, " (", scales::number(num_per_patient, accuracy = .01), ")"))
}

mean_aes <- 
  patients %>%
  bind_rows(
    patients %>% mutate(siteid = "OVERALL")
  ) %>% 
  select(siteid, treat, matches("num_s?ae.*")) %>%
  group_by(siteid, treat) %>%
  summarise(across(matches("num_s?ae.*"), ae_summary_function)) %>%
  ungroup() %>%
  mutate(col = paste(siteid, treat)) %>%
  select(-c(siteid, treat)) %>%
  pivot_longer(-col, names_to = "measure") %>%
  pivot_wider(names_from = col, values_from = value)

mean_aes
```

```{r}
missing_aes <- 
  patients %>%
  summarise_at(vars(matches("num_s?ae.*")), 
               (function(x) sum(is.na(x)))) %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("measure") %>%
  rename(n_missing = V1)

missing_aes
```

```{r}
t2_aes <- inner_join(mean_aes, missing_aes) %>%
  mutate(
    measure = case_when(
      measure == "num_ae"              ~ "Adverse events (AEs) (count (per patient))",
      measure == "num_sae"             ~ "Serious adverse events (SAEs) (count (per patient))",
      measure == "num_ae_qtc_prlngtn"  ~ "QTc prolongation AEs (count (per patient))",
      measure == "num_sae_qtc_prlngtn" ~ "QTc prolongation SAEs (count (per patient))",
      measure == "num_ae_elvtd_lfts"   ~ "Elevated LFT AEs (count (per patient))",
      measure == "num_sae_elvtd_lfts"  ~ "Elevated LFT SAEs (count (per patient))",
      measure == "num_ae_arrhythmia"   ~ "Arrhythmia AEs (count (per patient))",
      measure == "num_sae_arrhythmia"  ~ "Arrhythmia SAEs (count (per patient))"
    )
  )

t2_aes
```

### 4) All-cause mortality

```{r}
patients %>% count(died, died_days_after_enrdt)
```

### Checking consistency with the outcome data

Need to 

1. **exclude those who died after day 28**
2. **check this against the outcome data**

*Computing from the outcome data*:

```{r}
outcomes <- read_rds(file.path(output_dir, "outcomes_all.rds"))

outcomes_died_before_d28 <-
  outcomes %>% 
  group_by(siteid, patient_id) %>%
  arrange(days_post_enrollment) %>%
  summarise(
    dead_before_d28_outcomes = any((covid_scale == "1 - death") & (days_post_enrollment <= 28)),
    dead_after_d28_outcomes  = any((covid_scale == "1 - death") & (days_post_enrollment >  28))
  ) %>%
  ungroup()
```

Compare the three variables:

```{r}
# patients %>%
#   mutate(died_before_d28_patient_tbl = died & (died_days_after_enrdt <= 28),
#          died_after_d28_patient_tbl  = died & (died_days_after_enrdt >  28)) %>%
#   select(patient_id, died, died_before_d28_patient_tbl, died_after_d28_patient_tbl) %>%
#   left_join(
#     outcomes_died_before_d28 %>% select(patient_id, dead_before_d28_outcomes, dead_after_d28_outcomes),
#     by = "patient_id"
#   ) %>%
#   count(died_before_d28_patient_tbl, dead_before_d28_outcomes, died_after_d28_patient_tbl, dead_after_d28_outcomes)
```

Consistency checks:

-   If the patient table says the patient died on day x, and the outcome table has measures on days >=x, then they're all "1 - died".
-   If the outcome table says the patient had died before day y, then the patient table records them dying someday <= y.

```{r}
patients %>% 
  select(patient_id, died, died_days_after_enrdt) %>% 
  left_join(outcomes) %>%
  group_by(patient_id) %>%
  summarise(
    inconsistency_1 = any(died & (days_post_enrollment > died_days_after_enrdt) & (covid_scale != "1 - death")),
    inconsistency_2 = any((days_post_enrollment < died_days_after_enrdt) & (covid_scale == "1 - death"))
  ) %>%
  ungroup() %>%
  select(-patient_id) %>%
  count(inconsistency_1, inconsistency_2)
```

```{r}
patients %>% 
  select(patient_id, died, died_days_after_enrdt) %>% 
  left_join(outcomes) %>%
  group_by(patient_id) %>%
  mutate(
    inconsistency_1 = died & (days_post_enrollment > died_days_after_enrdt) & (covid_scale != "1 - death"),
    inconsistency_2 = (days_post_enrollment < died_days_after_enrdt) & (covid_scale == "1 - death"),
    inconsistency_3 = (days_post_enrollment <= 28) & (covid_scale == "1 - death") & (!died)
  ) %>%
  filter(inconsistency_1 | inconsistency_2 | inconsistency_3)
```

### Compiling results for Table 2

```{r}
patients <-
  patients %>%
  mutate(
    died_before_d28 = 
      ifelse(died, died_days_after_enrdt <= 28, FALSE))

mortality_summary_long <-
  patients %>%
  bind_rows(
    patients %>% mutate(siteid = "OVERALL")
  ) %>%
  add_count(siteid, treat, name = "n_treated") %>%
  add_count(siteid, treat, died_before_d28) %>%
  group_by(siteid, treat) %>%
  summarise(
    n_died_before_d28 = sum(died_before_d28, na.rm = TRUE),
    n_treated      = sum(!is.na(died_before_d28)),
    p_died_before_d28 = n_died_before_d28/n_treated,
    n_missing      = sum(is.na(died_before_d28))
  )

mortality_summary_long
```

```{r}
t2_mortality <-
  mortality_summary_long %>%
  ungroup() %>%
  mutate(col = paste(siteid, treat),
         measure = "All cause mortality at or before day 28 (count (percent))",
         value = str_c(scales::number(n_died_before_d28, accuracy = 1), 
                       " (", 
                       scales::percent(p_died_before_d28, accuracy = 1), 
                       ")"),
         n_missing = sum(is.na(patients$died_before_d28))) %>%
  select(col, measure, value, n_missing) %>%
  pivot_wider(names_from = col, values_from = value)

t2_mortality
```

# Compiling Table 2

Counts to order and format the table:


```{r}
site_arm_counts <-
  patients %>% 
  bind_rows(patients %>% mutate(siteid = "OVERALL")) %>%
  count(siteid, treat, name = "n_group") %>%
  mutate(group = paste0(siteid, " ", treat, " (n = ", n_group, ")")) %>% 
  group_by(siteid) %>%
  mutate(n_site = sum(n_group)) %>%
  mutate(siteid = paste0(siteid, " (n = ", n_site, ")")) %>%
  arrange(desc(n_site), treat) %>%
  mutate(n_group_string = paste0(ifelse(treat == "HCQ", "HCQ/CQ", "control"), " (n = ", n_group, ")"))

site_arm_counts

write_tsv(site_arm_counts, file.path(output_dir, "Table_2_site_arm_counts.tsv"))
```

Pull in the primary outcome data:

```{r}
primary_outcome_overall <- 
  #read_rds(file.path(output_dir, "standardized_effects/overall_analysis.rds"))$summary_table %>%
  read_rds(file.path(subgroup_dir, paste0(params$fit_name, "overall_analysis.rds")))$summary_table %>%
  mutate(value = ci_fmt(exp(logor_hat), exp(logor_95_lower), exp(logor_95_upper))) %>%
  select(overall_fct, analysis, value) %>%
  filter(analysis %in% c("superpop_posterior", "maximum likelihood")) %>%
  rename(siteid = overall_fct)
           

primary_outcome_sites <- 
  #read_rds(file.path(output_dir, "standardized_effects/site_analysis.rds"))$summary_table %>%
  read_rds(file.path(subgroup_dir, paste0(params$fit_name, "_site_analysis.rds")))$summary_table %>%
  mutate(value = ci_fmt(exp(logor_hat), exp(logor_95_lower), exp(logor_95_upper))) %>%
  select(siteid, analysis, value) %>%
  filter(analysis %in% c("superpop_posterior", "maximum likelihood"))

t2_primary_outcome <-
  bind_rows(primary_outcome_overall, primary_outcome_sites) %>%
  pivot_wider(names_from = siteid, values_from = value)

write_tsv(t2_primary_outcome, file.path(output_dir, "Table_2_primary.tsv"))

t2_primary_outcome
```

Compile the other outcomes, adding "categories" and formatting the numbers nicely:

```{r}
rate_formatter <- function(x) as.character(x)
  # Can't get this to work smoothly (acknowledging "true zeroes")
  # function(x) {ifelse(x == 0, " 0 ", paste0(" ", scales::number(x, accuracy = .01), " "))}

t2_nonprimary <-
  bind_rows(
    t2_mortality %>% mutate(category = "secondary outcome"),
    t2_hosp_days %>% mutate_at(vars(-measure, -n_missing), rate_formatter) %>% mutate(category = "secondary outcome"),
    t2_mech_vent %>% mutate(category = "secondary outcome"),
    t2_aes %>% mutate_at(vars(-measure, -n_missing), rate_formatter) %>% mutate(category = "safety outcome")
  )

t2_nonprimary
```

```{r}
knitr::kable(t2_nonprimary)
```

```{r}
write_tsv(t2_nonprimary, file.path(output_dir, "Table_2_nonprimary.tsv"))
```



```{r}
sessionInfo()
```



```{r}
Sys.time()
```