general-statistics.Rmd

---
title: "Exposure and Statistics Overview"
author: "Bryan Mayer"
date: "`r Sys.Date()`"
output: workflowr::wflow_html
---

Here, we calculate some of the initial (pre-model) results from the infant cohort and exposure characteristcs. 
- Demographics
- Initial survial curves
- Exposure assessment

```{r,  warning = F, message = F, echo = F}
knitr::opts_chunk$set(
  comment = NA,
  fig.align = "center",
  tidy = FALSE
)
```

```{r, message = F, warning = F}
library(survival)
library(tidyverse)
library(conflicted)
library(kableExtra)
library(cowplot)
conflict_prefer("filter", "dplyr")
```


```{r setup}

theme_set(
  theme_bw() +
    theme(panel.grid.minor = element_blank(),
          legend.position = "top")
)

source("code/plot_labels.R")
source("code/processing_functions.R")

exposure_data = read_csv("data/exposure_data.csv")
exposure_data_long = read_csv("data/exposure_data_long.csv")

leg_plot = ggplot(data = exposure_data_long, aes(x = idpar, y = count, fill = factor(obs_infected))) +
  geom_tile() +
  scale_fill_manual("", values = infection_labels$colours, breaks = infection_labels$breaks,
                      labels = infection_labels$labels)
trans_legend = get_legend(leg_plot + theme(legend.position = "top"))

range_str = function(x, digits = 3) paste(round(range(x), digits), collapse = " - ")
IQR_range_str = function(x, digits = 3) paste(round(quantile(x, c(0.25, 0.75)), digits), collapse = " - ")
  
```

# Demographics

## Infant ages

```{r}
exposure_data %>%
  select(FamilyID, enrollment_age) %>%
  distinct() %>%
  summarize(
        N = n(),
        enroll_median_age_days = median(enrollment_age),
        IQR = paste(quantile(enrollment_age, c(0.25, 0.75)), collapse = ", "),
        range_days = paste(range(enrollment_age), collapse = ", ")
        ) %>%
  kable() %>% kable_styling(full_width = F)

```

## Mom HIV

```{r}
exposure_data %>% select(FamilyID, momhiv) %>% 
  distinct() %>%
  group_by(momhiv) %>%
  summarize(N = n()) %>% 
  kable() %>% kable_styling(full_width = F)
```


# Survival analysis

```{r}

exposure_data %>% 
  group_by(virus, FamilyID) %>% 
  summarize(obs_infected = max(infectious_1wk),
            is_infected = max(infected)) %>% group_by(virus) %>%
  summarize(
    total_infants = n_distinct(FamilyID),
    total_infected = sum(is_infected),
    total_outcome = sum(obs_infected)
    ) %>%
  kable() %>%
  kable_styling(full_width = F)

```


```{r}
surv_data = exposure_data %>% 
  group_by(FamilyID, virus, momhiv, final_infant_wk) %>%
  summarize(
    infected = max(infected)
  )

surv_fit = surv_data %>%
  group_by(virus) %>%
  nest() %>%
  mutate(
    surv_mod = map(data, ~survfit(Surv(final_infant_wk, infected) ~ 1, data = .)),
    surv_mod_hiv = map(data, ~survfit(Surv(final_infant_wk, infected) ~ momhiv, data = .)),
    logrank = map_dbl(data, ~coin::pvalue(coin::logrank_test(Surv(final_infant_wk, infected) ~ factor(momhiv),
                                                             data = ., distribution = "exact")))
    ) %>%
  select(-data)

surv_fit %>%
  select(virus, logrank) %>%
  rename(`Mother HIV Log-rank` = logrank) %>%
  kable() %>% kable_styling(full_width = F)


surv_res = pmap_df(surv_fit, function(virus, surv_mod, surv_mod_hiv, logrank){
   broom::tidy(surv_mod) %>%
    mutate(strata = "Pooled") %>%
    bind_rows(broom::tidy(surv_mod_hiv)) %>%
    mutate(
      virus = virus,
      momhiv = str_remove_all(strata, "momhiv=")
      ) %>%
    bind_rows(crossing(virus = virus, time = -1e-12, estimate = 1, momhiv = c("Pooled", "neg", "pos")))
})

```


```{r survival-plot}

# NA censors recoded to 0 (occurs at time ~ 0)

surv_res %>%
  arrange(virus, momhiv, time) %>%
  mutate(n.censor = if_else(is.na(n.censor), 0, n.censor)) %>%
  ungroup() %>%
  mutate(virus = factor(virus, levels = c("HHV-6", "CMV"))) %>%
  ggplot(aes(x = time, y = 1 - estimate, colour = momhiv)) + 
  geom_step() +
  geom_point(aes(shape = n.censor > 0)) +
  scale_shape_manual(guide = F, values = c(-1, 3)) +
  scale_x_continuous("Weeks after infant birth", breaks = 0:10 * 10) +
  scale_y_continuous("Proportion infected", expand = c(0.01, 0)) +
  geom_vline(xintercept = 52, colour = "black", linetype = "dashed") +
  scale_color_discrete("", breaks = c("neg", "pos", "Pooled"), 
                     labels = c("Mother HIV-", "Mother HIV+", "Pooled")) +
  facet_wrap(~virus) +
  geom_text(data = mutate(ungroup(surv_fit), virus = factor(virus, levels = c("HHV-6", "CMV"))), 
            aes(label = str_c("p = ", round(logrank, 2))), 
            x = Inf, y = 0, colour = "black", vjust = -.2, hjust = 1.2) +
  theme(legend.position = "top")

```

# Exposure Analysis

## Overall

```{r}

exposure_data_long %>% 
  subset(idpar != "HH") %>%
  group_by(virus, idpar) %>% 
  summarize(
    total = n(),
    total_observed = total - sum(interpolated),
    total_interpolate = stat_paste(sum(interpolated),  100*mean(interpolated), digits = 1)
  ) %>%
  mutate(`Exposure source` = factor(idpar, levels = c("S", "M"),
      labels = c("Secondary Children", "Mother"))) %>%
  select(virus, `Exposure source`, everything(), -idpar) %>%
  arrange(desc(virus), `Exposure source`) %>%
  write_csv("output/results-tables/supp_table1.csv") %>%
  kable() %>%
  kable_styling(full_width = F)


```


```{r exposure-profile, fig.height=8}

exposure_data_summary = exposure_data_long %>%
  mutate(
    pos_count = count > 0
    ) %>%
  subset(!interpolated) %>%
  group_by(virus, FamilyID, obs_infected, idpar) %>%
  mutate(
    total_pos = sum(pos_count),
    pct_pos = 100 * mean(pos_count)
    ) %>%
  group_by(pct_pos, total_pos, add = T) %>%
  summarise_at(vars(count), list(~n(), mean = mean, median = median, maximum = max)) %>%
  rename(N = n)

plot_labels = exposure_data_summary %>%
  gather(stat, estimate, mean, maximum, pct_pos) %>%
  group_by(stat) %>%
  summarize(min_lim = min(estimate), max_lim = ceiling(max(estimate))) %>%
  left_join(tibble(stat = c("pct_pos", "mean", "maximum"), 
                   out_lab = c("Percent~Positive", "Mean~Log[10]~VL", "Max.~Log[10]~VL"))) %>%
  mutate(stat = factor(stat, levels = c("pct_pos", "mean", "maximum"))) %>%
  arrange(stat) %>%
  ungroup() %>%
  mutate(letter_code = 1:3)

pl_exposure = map(plot_labels$stat %>% levels(), function(s){
  exposure_data_summary = exposure_data_summary %>%
    ungroup() %>%
    mutate(virus = factor(virus, levels = c("HHV-6", "CMV")))
  
  
  tmp_theme = theme(
    legend.position = "none",
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 9)
  )
    
  pl_lab = subset(plot_labels, stat == s)
  out_label = pl_lab$out_lab[1]
  lower_limit = pl_lab$min_lim
  upper_limit = pl_lab$max_lim
    
  pl1 = exposure_data_summary %>%
    gather(stat, estimate, mean, maximum, pct_pos) %>%
    select(-median,-total_pos,-N) %>%
    spread(idpar, estimate) %>%
    filter(stat == s) %>%
    ggplot(aes(
      x = S,
      y = M,
      colour = factor(obs_infected)
    )) +
    geom_point() +
    geom_abline() +
    geom_point() +
    scale_colour_manual(values = c("#4393C3", "#D6604D")) +
    scale_y_continuous(parse(text = paste0("Mother~", out_label)),
                       limits = c(lower_limit, upper_limit)) +
    scale_x_continuous(parse(text = paste0("Secondary~Children~", out_label)),
                       limits = c(lower_limit, upper_limit)) +
    facet_wrap( ~ virus) +
    tmp_theme
  
  pl2 = exposure_data_summary %>%
    ungroup() %>%
    mutate(
      obs_infected = factor(obs_infected),
      idpar = fct_recode(
        fct_rev(idpar),
        "Secondary\nChildren" = "S",
        "Mother" = "M",
        "Household\nSum" = "HH"
      )
    ) %>%
    ggplot(aes_string(x = "idpar", y = s, colour = "obs_infected")) +
    geom_boxplot() +
    scale_colour_manual(values = c("#4393C3", "#D6604D")) +
    geom_point(position = position_dodge(width = 0.75)) +
    scale_y_continuous(parse(text = out_label), limits = c(lower_limit, upper_limit)) +
    xlab(parse(text = c(""))) +
    facet_wrap( ~ virus) +
    tmp_theme + theme(axis.text.x = element_text(size = 7))
  
    plot_grid(pl1, pl2, nrow = 1, labels = paste0("(",letters[c(pl_lab$letter_code, pl_lab$letter_code + 3)], ")"))
    
})

plot_grid(plot_grid(plotlist = pl_exposure, nrow = 3), trans_legend, nrow = 2, rel_heights = c(11, 1))

```

```{r exp-corr}

exposure_data_summary %>%
  subset(idpar != "HH") %>%
  group_by(virus, FamilyID, idpar) %>%
  gather(outcome, value, pct_pos, mean, maximum) %>%
  select(virus, FamilyID, idpar, outcome, value) %>%
  spread(idpar, value) %>%
  ungroup() %>%
  dplyr::mutate(
    fid = factor(FamilyID)
  ) %>%
  group_by(virus, outcome) %>%
  nest() %>%
  mutate(
    cor_test = map(data, ~cor.test(.x$M, .x$S), method = "Spearmen"),
    cor_res = map(cor_test, broom::tidy)
  ) %>%
  unnest(cor_res) %>%
  mutate(
    Endpoint = factor(outcome, levels = c("pct_pos", "mean", "maximum"),
                      labels = c("Pct. Positive", "Mean VL", "Maximum VL")),
    Estimate = stat_paste(estimate, conf.low, conf.high, digits = 2),
    p.value = round(p.value, 3)
  ) %>%
  ungroup() %>%
  select(virus, Endpoint, Estimate, p.value) %>%
  arrange(desc(virus), Endpoint) %>%
  write_csv("output/results-tables/supp_table2.csv") %>%
  kable(digits = 3) %>%
  kable_styling(full_width = F) %>%
  collapse_rows(1)

```

```{r}

exposure_data_tests = 
  exposure_data_summary %>% 
  gather(stat, estimate, mean, maximum, pct_pos) %>%
  group_by(virus, stat, idpar) %>%
  mutate(obs_infected = factor(obs_infected)) %>%
  summarize(
    wilcox_pvalue =  coin::pvalue(coin::wilcox_test(estimate ~ obs_infected, distribution = "exact"))
  )

overall_summary = exposure_data_summary %>% 
  gather(stat, estimate, mean, maximum, pct_pos) %>%
  group_by(virus, stat, idpar, obs_infected) %>%
  summarize(
    median = median(estimate)
  ) %>%
  mutate(obs_infected = recode_factor(obs_infected, `1` = "Transmission", `0` = "No transmission", .ordered = T)) %>%
  spread(obs_infected, median) %>%
  mutate(Difference = `Transmission`-`No transmission`) %>%
  left_join(exposure_data_tests, by = c("virus", "stat", "idpar"))

overall_summary %>%
  ungroup() %>%
  mutate(
    wilcox_pvalue = as.character(round_away_0(wilcox_pvalue, 3, T)),
    Endpoint = factor(
      stat,
      levels = c("pct_pos", "mean", "maximum"),
      labels = c("Pct. Positive", "Mean VL", "Maximum VL")
    ),
    `Exposure source` = factor(idpar, levels = c("S", "M", "HH"),
      labels = c("Secondary child", "Mother", "Household Sum"))
    ) %>%
  select(virus, Endpoint, `Exposure source`, everything(), -idpar, -stat) %>%
  mutate_if(is.numeric, round, 2) %>% 
  arrange(desc(virus), Endpoint, `Exposure source`) %>%
  write_csv("output/results-tables/supp_table3.csv") %>%
  kable() %>%
  kable_styling(full_width = F) %>%
  collapse_rows(columns = 1:2)

```


## Household composition

Household sum composition was determined and reported by taking the mean proportion over measurements within a household. The summary across the households uses median and IQR to match the box plot statistics.

```{r hh-composition, fig.width = 4}

hh_summary = exposure_data %>%
  filter(HH > 0) %>%
  mutate(   
    S_pctHH = if_else(HH == 0, 0, 100 * (10^S/10^HH)),
    M_pctHH = if_else(HH == 0, 0, 100 * (10^M/10^HH))
    ) %>%
  group_by(virus, FamilyID, obs_infected) %>%
  summarise_at(vars(S_pctHH, M_pctHH), list(mean = mean, median = median))

hh_summary %>%
  ungroup() %>%
  select(-obs_infected) %>%
  gather(stat, est, -virus, -FamilyID) %>%
  group_by(virus, stat) %>%
  summarize_if(is.double, list(mean = mean, median = median, IQR = IQR_range_str, range = range_str)) %>%
  ungroup() %>%
  mutate_if(is.double, round, digits = 2) %>%
  filter(str_detect(stat, "mean")) %>%
  mutate(
    stat = substr(stat, 1, 1)
  ) %>%
  rename(`Household member` = stat) %>%
  rename_at(vars(-IQR), list(str_to_title)) %>%
  kable(caption = "Percent household composition") %>% kable_styling(full_width = F)
  
hh_summary %>%
  ungroup() %>%
  mutate(virus = factor(virus, levels = c("HHV-6", "CMV"))) %>%
  ggplot(aes(x = virus, y = S_pctHH_mean)) +
  geom_boxplot() +
  geom_point(aes(colour = factor(obs_infected))) +
  xlab("") +
  scale_colour_manual("", values = infection_labels$colours, breaks = infection_labels$breaks,
                      labels = infection_labels$labels) +
  ylab("Percent of household shedding attributable to secondary children") 

```


## Interpolated exposure

In the sensitivity analysis, we'd like to assess the interpolation. All of the summary statistical analysis is limited to the observed exposures. Assuming the interpolation was unbiased, the precision in tests could still be artificially inflated without some correction. Because linear interpolation is used, the mean and maximum exposure estimates should be largely unaffected. However, the percent positive estimate may not be precise without a larger sample size.

```{r interpolate-plot}

exposure_interpolated_summary = exposure_data_long %>%
  mutate(
    pos_count = count > 0
    ) %>%
  group_by(virus, FamilyID, obs_infected, idpar) %>%
  mutate(
    interpolated_pct = 100*mean(interpolated),
    total_pos = sum(pos_count),
    pct_pos = 100 * mean(pos_count)
    ) %>%
  group_by(pct_pos, total_pos, interpolated_pct, add = T) %>%
  dplyr::summarise_at(vars(count), list(~n(), mean = mean, median = median, maximum = max))

exposure_interpolated_summary %>%
  ungroup() %>%
  mutate(
    idpar = fct_recode(fct_rev(idpar), 
                       "Secondary\nChildren" = "S", "Mother" = "M", "Household\nSum" = "HH")
    ) %>%
  ggplot(aes(x = idpar, y = interpolated_pct)) +
  geom_boxplot() +
  geom_point() +
  xlab("") +
  ylab("Percent of weekly exposures interpolated") +
  facet_wrap(~virus)

```

```{r interpolate-exposure}

exposure_interpolated_summary %>% 
  ungroup() %>%
  subset(idpar != "HH") %>%
  gather(stat, estimate, mean, maximum, pct_pos) %>%
  mutate(
    idpar = fct_recode(fct_rev(idpar), 
                       "Secondary\nChildren" = "S", "Mother" = "M"),
    stat = fct_recode(fct_rev(stat), 
                       "% positive (weekly)" = "pct_pos", "Mean VL" = "mean", "Maximum VL" = "maximum"),
    ) %>%
  ggplot(aes(x = interpolated_pct, y = estimate, colour = factor(obs_infected))) +
  geom_point() +
  ylab("") +
  xlab("Percent of weekly exposures interpolated") +
  facet_grid(stat ~ virus+idpar, scales = "free_y", switch = "y") +
  scale_colour_manual("", values = infection_labels$colours, breaks = infection_labels$breaks,
                      labels = infection_labels$labels) +
  theme(strip.placement = "outside")

```