Report_SARS-CoV-2_Target_Failures_and_Outliers.Rmd

---
#title: 'SARS-CoV-2 Report'
#author: "Brendan J. Kelly, MD, MS"
output: 
  html_document:
    self_contained: true
    #toc: true
    #toc_float: true
css: style_helvetica_double.css
---

```{r setup, include=FALSE, eval=TRUE, include=FALSE}

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, error = FALSE)
library(tidyverse)
library(readxl)
library(janitor)
library(stringr)
library(lubridate)
library(ggsci)
library(gt)
set.seed(16)



```



```{r read_process_data, echo=FALSE, eval=TRUE, include=FALSE, message=FALSE, error=FALSE}

#' possible instruments
instrument_search_string <- "8800|BDMax|Cepheid|Eplex|HUP54|LIAT|m2000|Roche|Simplexa|Thermo"

#' filepath to long-format Ct data
long_data_folder_path <- "../ct_data_reports_long_format/"

#' read and clean long-format Ct data
long_data <- list.files(path = long_data_folder_path, full.names = TRUE, pattern = "csv$")
dat <- map(.x = long_data, .f = ~ vroom::vroom(.x)) %>%
  bind_rows()


dat %>%
  # clean variable names and select variables
  janitor::clean_names() %>%
  select(accession_number, discrete_task_display, contains("date"), contains("order_"), contains("result")) %>%
  select(-birth_date) %>%
  select(-perform_result_ascii_text) %>%
  # recognize NA values
  mutate(numeric_result_value = replace(numeric_result_value, numeric_result_value > 100 | numeric_result_value == 88, NA)) %>%
  # format dates to permit temporal analysis
  mutate_at(.vars = vars(contains("date")), .funs = ~ as.Date(.x, format = "%Y/%m/%d")) %>%
  mutate(discrete_task_display = gsub("SARS-CoV-2 \\(|\\)","",discrete_task_display)) %>%
  # ensure platform and amplicon are clearly documented
  rename(amplicon = discrete_task_display,
         platform = text_result_value) %>%
  group_by(accession_number) %>%
  mutate(platform = stringr::str_extract(string = paste0(platform, collapse = " "), pattern = instrument_search_string),
         platform = replace(platform, is.na(platform), "Platform Not Recorded")) %>%
  ungroup() %>%
  # add year and week values
  mutate(year = lubridate::year(container_in_lab_date_time),
         week = lubridate::week(container_in_lab_date_time),
         year_week = glue::glue("{year}-{stringr::str_pad(week, width = 2, side = 'left', pad = '0')}")) %>%
  # filter to include only recorded Ct values
  filter(!is.na(numeric_result_value)) %>%
  distinct(accession_number, platform, amplicon, order_id, numeric_result_value, .keep_all = TRUE) %>%
  identity() -> ct

ct


#' patch Fluidigm platform label & filter to include only 8800, Fluidigm, Simplexa, Thermo
ct %>%
  mutate(platform = replace(platform, grepl("Ct N1|Ct N2",amplicon), "Fluidigm")) %>%
  # restrict platforms analyzed:
  filter(platform %in% c("8800", "Fluidigm", "Simplexa", "Thermo")) %>%
  filter(amplicon != "Ct Value") %>%
    distinct(accession_number, platform, amplicon, order_id, numeric_result_value, .keep_all = TRUE) %>%
  identity() -> ct



```
  

```{r drop_data, echo=FALSE, eval=TRUE, include=FALSE, message=FALSE, error=FALSE}

# data on dropped values
ct %>%
  mutate(amplicon_code = case_when(amplicon == "Ct ORF1ab" ~ "target_1",
                                   amplicon == "Ct S gene" ~ "target_2",
                                   amplicon == "Ct N gene" ~ "target_3",
                                   amplicon == "Ct N1" ~ "target_4",
                                   amplicon == "Ct N2" ~ "target_5",
                                   amplicon == "Ct E gene" ~ "target_6"
                                   )) %>%
  select(accession_number, container_in_lab_date_time, year, week, year_week, platform, amplicon_code, numeric_result_value) %>%
  distinct() %>%
  filter(complete.cases(.)) %>%
  pivot_wider(id_cols = c(accession_number, container_in_lab_date_time, year, week, year_week, platform),
              names_from = amplicon_code,
              values_from = numeric_result_value) %>%
  group_by(accession_number) %>%
  mutate(dropout = (platform == "8800" & ((!is.na(target_1) & is.na(target_6)) | (!is.na(target_6) & is.na(target_1)))) | (platform == "Fluidigm" & ((!is.na(target_4) & is.na(target_5)) | (!is.na(target_5) & is.na(target_4)))) | (platform == "LIAT" & ((!is.na(target_1) & is.na(target_2)) | (!is.na(target_2) & is.na(target_1)))) | (platform == "m2000" & ((!is.na(target_1) & is.na(target_2)) | (!is.na(target_2) & is.na(target_1)))) | (platform == "Roche" & ((!is.na(target_1) & is.na(target_2)) | (!is.na(target_2) & is.na(target_1)))) | (platform == "Simplexa" & ((!is.na(target_1) & is.na(target_2)) | (!is.na(target_2) & is.na(target_1)))) | (platform == "Thermo" & ((!is.na(target_1) & !is.na(target_2) & is.na(target_3)) | (!is.na(target_2) & !is.na(target_3) & is.na(target_1)) | (!is.na(target_3) & !is.na(target_1) & is.na(target_2))))) %>%
  mutate(s_dropout = (platform == "LIAT" & ((!is.na(target_1) & is.na(target_2)))) | (platform == "m2000" & ((!is.na(target_1) & is.na(target_2)))) | (platform == "Roche" & ((!is.na(target_1) & is.na(target_2)))) | (platform == "Simplexa" & ((!is.na(target_1) & is.na(target_2)))) | (platform == "Thermo" & (!is.na(target_3) & !is.na(target_1) & is.na(target_2)))) %>%
  mutate(n1_dropout = (platform == "Fluidigm" & (!is.na(target_5) & is.na(target_4)))) %>%
  mutate(n2_dropout = (platform == "Fluidigm" & (!is.na(target_4) & is.na(target_5)))) %>%
  ungroup() %>%
  identity() -> ct_drop




```


  
## Clinical Pathology SARS-CoV-2 Test Report `r format(Sys.Date(), '%d %B %Y')`  {.tabset .tabset-fade .tabset-pills}  
  
### Specimens    
  
This report summarizes `r length(unique(ct$accession_number))` specimens submitted between <u>`r min(ct$container_in_lab_date_time, na.rm = TRUE)`</u> and <u>`r max(ct$container_in_lab_date_time, na.rm = TRUE)`</u>:    
  
```{r platform_specimens, echo=FALSE, eval=FALSE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

library(downloadthis)

ct %>%
  select(accession_number, platform, amplicon, container_in_lab_date_time, year_week, numeric_result_value) %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_all_specimens_",Sys.Date()),
    output_extension = ".xlsx",
    button_label = "Download Data as XLSX",
    button_type = "success",
    has_icon = TRUE,
    icon = "fa fa-save"
  )

ct %>%
  select(accession_number, platform, amplicon, container_in_lab_date_time, year_week, numeric_result_value) %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_all_specimens_",Sys.Date()),
    output_extension = ".csv",
    button_label = "Download Data as CSV",
    button_type = "info",
    has_icon = TRUE,
    icon = "fa fa-save",
    csv2 = FALSE
  )


```
  
<br>  
  
```{r platform_summary, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

ct %>%
  group_by(platform, amplicon) %>%
  summarise(n = n_distinct(accession_number)) %>%
  ungroup() %>%
  filter(grepl("Ct",amplicon)) %>%
  gt() %>%
  gt::cols_label(platform = "Platform", amplicon = "Target Description", n = "# Specimens with Measured Ct Values") %>%
  gt::opt_table_lines(extent = "all")


```
  
-----  
\newpage
<P style = "page-break-before: always">
  
### Ct Values    

```{r ct_values_aggregate_fig, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

# figures
ct %>%
  mutate(platform = replace(platform, platform == "Platform Not Recorded", "Platform<br>Not Recorded")) %>%
  ggplot(data = .) +
  geom_boxplot(aes(x = amplicon, y = numeric_result_value, fill = amplicon), outlier.size = 0.5, outlier.shape = 21) +
  facet_wrap(facets = ~ platform, scales = "free") +
  theme_bw() +
  theme(axis.text.y = element_text(color = "black", size = 6),
        axis.text.x = element_text(color = "black", size = 6, angle = 325, hjust = 0),
        strip.text = ggtext::element_markdown(color = "black", size = 6),
        strip.background = element_blank()) +
  ggsci::scale_fill_nejm() +
  labs(x = "", y = "Ct Value", fill = "RT-PCR\nTarget")


```

<br>

```{r ct_values_aggregate_tab, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

# table: aggregate central tendency/spread of Ct values
ct %>%
  group_by(platform, amplicon) %>%
  summarise(n = n_distinct(accession_number),
            median_ct = median(numeric_result_value, na.rm = TRUE),
            mean_ct = mean(numeric_result_value, na.rm = TRUE),
            sd_ct = sd(numeric_result_value, na.rm = TRUE),
            first_quartile_ct = quantile(numeric_result_value, probs = 0.25, na.rm = TRUE),
            third_quartile_ct = quantile(numeric_result_value, probs = 0.75, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(mean_sd = glue::glue("{round(mean_ct,1)} ({round(sd_ct,1)})"),
         median_iqr = glue::glue("{round(median_ct,1)} ({round(first_quartile_ct,1)} - {round(third_quartile_ct,1)})")) %>%
  select(platform, amplicon, n, median_iqr, mean_sd) %>%
  gt() %>%
  gt::cols_label(platform = "Platform",
                 amplicon = "Amplicon",
                 n = "Positive Specimens",
                 median_iqr = "Median (IQR) Ct",
                 mean_sd = "Mean (SD) Ct") %>%
  gt::opt_table_lines(extent = "all")



```


  
  
-----  
\newpage
<P style = "page-break-before: always">
  
### Temporal Trends      

Restricted to platforms with >= 2 weeks of recorded data.  


```{r ct_over_time, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

# weekly changes in central tendency/spread of Ct values
ct %>%
    mutate(platform = replace(platform, platform == "Platform Not Recorded", "Platform<br>Not Recorded")) %>%
  # only include platforms with at least 2 weeks of data
  group_by(platform) %>%
  filter(n_distinct(year_week) >= 2) %>%
  ungroup() %>%
  ggplot(data = .) +
  geom_boxplot(aes(x = year_week, y = numeric_result_value, fill = amplicon), outlier.size = 0.5, outlier.shape = 21) +
  facet_grid(rows = vars(amplicon), cols = vars(platform), scales = "free_y") +
  theme_bw() +
  theme(axis.text.y = element_text(color = "black", size = 6),
        axis.text.x = element_text(color = "black", size = 6, angle = 325, hjust = 0),
        strip.text = ggtext::element_markdown(color = "black", size = 6),
        strip.background = element_blank()) +
  ggsci::scale_fill_nejm() +
  labs(x = "Year-Week", y = "Ct Value", fill = "RT-PCR\nTarget")


```
    
<br>  
  
```{r ct_over_time_table, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

# weekly positive specimen totals
ct %>%
  mutate(platform = replace(platform, platform == "Platform Not Recorded", "Platform<br>Not Recorded")) %>%
  # only include platforms with at least 2 weeks of data
  group_by(platform) %>%
  filter(n_distinct(year_week) >= 2) %>%
  ungroup() %>%
  group_by(platform, year_week) %>%
  summarise(number_pos_specimens = n_distinct(accession_number)) %>%
  ungroup() %>%
  pivot_wider(id_cols = year_week, names_from = platform, values_from = number_pos_specimens, values_fill = 0) %>%
  arrange(year_week) %>%
  gt() %>%
  gt::cols_label(year_week = "Year - Week",
                 ) %>%
  gt::opt_table_lines(extent = "all") %>%
  gt::tab_spanner(label = "Positive Specimen Count", columns = 2:tidyselect::last_col())


```

  
-----  
\newpage
<P style = "page-break-before: always">
  
### Target Failures      
  
Target failures are defined by failure to detect a SARS-CoV-2 target where other another target is detected from the same specimen.  
  
```{r ct_dropout_download, echo=FALSE, eval=FALSE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

library(downloadthis)

ct_drop %>%
  filter(dropout == TRUE) %>%
  select(accession_number, year_week, platform, dropout, s_dropout) %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_dropouts_",Sys.Date()),
    output_extension = ".xlsx",
    button_label = "Download Data as XLSX",
    button_type = "success",
    has_icon = TRUE,
    icon = "fa fa-save"
  )

ct_drop %>%
  filter(dropout == TRUE) %>%
  select(accession_number, year_week, platform, dropout, s_dropout) %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_dropouts_",Sys.Date()),
    output_extension = ".csv",
    button_label = "Download Data as CSV",
    button_type = "info",
    has_icon = TRUE,
    icon = "fa fa-save",
    csv2 = FALSE
  )



```
  
<br>  
  
```{r ct_dropout_fig, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

# figure
ct_drop %>%
  filter(platform %in% c("8800","Fluidigm","LIAT","m2000","Roche","Simplexa","Thermo")) %>%
  group_by(year_week) %>%
  summarise(n_two_target = n_distinct(accession_number),
            total_dropout = sum(dropout, na.rm = TRUE),
            s_gene_dropout = sum(s_dropout, na.rm = TRUE),
            n1_gene_dropout = sum(n1_dropout, na.rm = TRUE),
            n2_gene_dropout = sum(n2_dropout, na.rm = TRUE),
            tdo_prop = sum(dropout, na.rm = TRUE) / n_distinct(accession_number),
            sgd_prop = sum(s_dropout, na.rm = TRUE) / n_distinct(accession_number),
            n1gd_prop = sum(n1_dropout, na.rm = TRUE) / n_distinct(accession_number),
            n2gd_prop = sum(n2_dropout, na.rm = TRUE) / n_distinct(accession_number),
            total_report = glue::glue("{total_dropout} ({round(tdo_prop*100,1)}%)"),
            s_gene_report = glue::glue("{s_gene_dropout} ({round(sgd_prop*100,1)}%)"),
            n1_gene_report = glue::glue("{n1_gene_dropout} ({round(n1gd_prop*100,1)}%)"),
            n2_gene_report = glue::glue("{n2_gene_dropout} ({round(n2gd_prop*100,1)}%)")) %>%
  ungroup() %>%
  select(year_week, contains("prop")) %>%
  pivot_longer(cols = contains("prop"), names_to = "target_type", values_to = "proportion_dropout") %>%
  mutate(target_type = case_when(grepl("tdo",target_type) ~ "Any Dropout",
                                 grepl("sgd",target_type) ~ "S Gene Dropout",
                                 grepl("n1gd",target_type) ~ "N1 Gene Dropout",
                                 grepl("n2gd",target_type) ~ "N2 Gene Dropout")) %>%
  identity() %>%
  mutate(target_type = factor(target_type, levels = c("Any Dropout", "S Gene Dropout", "N1 Gene Dropout", "N2 Gene Dropout"))) %>%
  ggplot(data = .) +
  geom_line(aes(x = year_week, y = proportion_dropout, group = target_type, color = target_type), alpha = 0.8) +
  geom_point(aes(x = year_week, y = proportion_dropout, color = target_type), alpha = 0.8) +
  theme_bw() +
  theme(axis.text = element_text(color = "black", size = 5),
        strip.background = element_blank(),
        legend.position = "top") +
  ggsci::scale_color_nejm() +
  labs(x = "Year-Week", y = "Proportion Dropout", color = "")


```

<br>  
  

```{r ct_dropout_tab, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

# table  
ct_drop %>%
  filter(platform %in% c("8800","Fluidigm","LIAT","m2000","Roche","Simplexa","Thermo")) %>%
  group_by(year_week) %>%
  summarise(n_two_target = n_distinct(accession_number),
            total_dropout = sum(dropout, na.rm = TRUE),
            s_gene_dropout = sum(s_dropout, na.rm = TRUE),
            n1_gene_dropout = sum(n1_dropout, na.rm = TRUE),
            n2_gene_dropout = sum(n2_dropout, na.rm = TRUE),
            tdo_prop = sum(dropout, na.rm = TRUE) / n_distinct(accession_number),
            sgd_prop = sum(s_dropout, na.rm = TRUE) / n_distinct(accession_number),
            n1gd_prop = sum(n1_dropout, na.rm = TRUE) / n_distinct(accession_number),
            n2gd_prop = sum(n2_dropout, na.rm = TRUE) / n_distinct(accession_number),
            total_report = glue::glue("{total_dropout} ({round(tdo_prop*100,1)}%)"),
            s_gene_report = glue::glue("{s_gene_dropout} ({round(sgd_prop*100,1)}%)"),
            n1_gene_report = glue::glue("{n1_gene_dropout} ({round(n1gd_prop*100,1)}%)"),
            n2_gene_report = glue::glue("{n2_gene_dropout} ({round(n2gd_prop*100,1)}%)")) %>%
  ungroup() %>%
  select(year_week, n_two_target, contains("report")) %>%
  gt() %>%
  gt::cols_label(total_report = "Any Dropout (%)",
                 s_gene_report = "S Gene Dropout (%)",
                 n1_gene_report = "N1 Gene Dropout (%)",
                 n2_gene_report = "N2 Gene Dropout (%)",
                 year_week = "Year - Week",
                 n_two_target = "Number of Two-Target Specimens") %>%
  gt::opt_table_lines(extent = "all")



```
  
<br>
  
```{r ct_dropout_specimens, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

ct_drop %>%
  arrange(container_in_lab_date_time) %>%
  filter(dropout == TRUE) %>%
  select(accession_number, year_week, platform, contains("dropout")) %>%
  mutate(which_dropout = case_when(s_dropout == TRUE ~ "S Gene",
                               n1_dropout == TRUE ~ "N1 Gene",
                               n2_dropout == TRUE ~ "N2 Gene",
                               s_dropout == FALSE & n1_dropout == FALSE & n2_dropout == FALSE ~ "Other")) %>%
  gt() %>%
  tab_style(
    style = list(
      cell_fill(color = ggsci::pal_nejm(alpha = 0.5)(4)[1])#,
        #cell_text(weight = "bold")
        ),
    locations = cells_body(
        columns = "which_dropout",
        rows = which_dropout == "Other")
    ) %>%
  tab_style(
    style = list(
      cell_fill(color = ggsci::pal_nejm(alpha = 0.5)(4)[2])#,
        #cell_text(weight = "bold")
        ),
    locations = cells_body(
        columns = "which_dropout",
        rows = which_dropout == "S Gene")
    ) %>%
  tab_style(
    style = list(
      cell_fill(color = ggsci::pal_nejm(alpha = 0.5)(4)[3])#,
        #cell_text(weight = "bold")
        ),
    locations = cells_body(
        columns = "which_dropout",
        rows = which_dropout == "N1 Gene")
    ) %>%
    tab_style(
    style = list(
      cell_fill(color = ggsci::pal_nejm(alpha = 0.5)(4)[4])#,
        #cell_text(weight = "bold")
        ),
    locations = cells_body(
        columns = "which_dropout",
        rows = which_dropout == "N2 Gene")
    ) %>%
    gt::cols_label(accession_number = "Accession Number",
                 year_week = "Year - Week",
                 platform = "Platform",
                 which_dropout = "Failed Target") %>%
  gt::cols_hide(columns = c("s_dropout","n1_dropout","n2_dropout","dropout")) %>%
  gt::opt_table_lines(extent = "all")



```

 
  
-----  
\newpage
<P style = "page-break-before: always">
  
### Outlier Rules      
  
Download rules for calling outlier specimens based on modeling performed between February and May 2021:    
  
```{r ct_diff_model_rules, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}


#' ################################
#' add posterior quantiles to individual data points
#' - from Matthew Kay: https://mjskay.github.io/tidybayes/articles/tidybayes-residuals.html
#' ################################

library(tidybayes)
library(downloadthis)

m_dif <- read_rds(file = "./models/ct_dif_by_platform_brms_model.rds.bz2")

ct_drop %>%
  arrange(container_in_lab_date_time) %>%
  # # impute missing (target failure / dropout) values at 40 cycles
  # replace_na(list("target_1" = 40, "target_2" = 40, "target_3" = 40, "target_4" = 40, "target_5" = 40, "target_6" = 40)) %>%
  group_by(accession_number) %>%
  mutate(target_dif = case_when(platform == "8800" ~ target_6 - target_1,
                                platform == "Fluidigm" ~ target_5 - target_4,
                                platform == "LIAT" ~ target_2 - target_1,
                                platform == "m2000" ~ target_2 - target_1,
                                platform == "Roche" ~ target_2 - target_1,
                                platform == "Simplexa" ~ target_2 - target_1,
                                platform == "Thermo" ~ target_2 - target_1)) %>%
  nest() %>%
  mutate(min_ct = map(.x = data, .f = ~ min(c(.x[["target_1"]], .x[["target_2"]], .x[["target_3"]], .x[["target_4"]], .x[["target_5"]], .x[["target_6"]]), na.rm = TRUE))) %>%
  unnest(cols = c("data","min_ct")) %>%
  ungroup() %>%
  filter(!is.na(target_dif) & !is.na(min_ct)) %>%
  group_by(platform) %>%
  filter(n_distinct(accession_number) >= 10) %>%
  ungroup() %>%
  identity() -> m_dat


m_dat %>%
  add_predicted_draws(m_dif) %>%
  select(accession_number, platform, target_dif, min_ct, .row, .prediction) %>%
  summarise(p_residual = mean(.prediction < target_dif),
            p_residual_corrected = ifelse(p_residual < 0.5, 1 - p_residual, p_residual),
            z_residual = qnorm(p_residual),
            z_residual_corrected = qnorm(p_residual_corrected)) %>%
  mutate(p_residual_corrected_999 = p_residual_corrected > 0.999,
         p_residual_corrected_99 = p_residual_corrected > 0.99,
         p_residual_corrected_95 = p_residual_corrected > 0.95,
         ) %>%
  ungroup() %>%
    select(accession_number, year_week, platform, target_dif, min_ct, p_residual_corrected, p_residual_corrected_95, p_residual_corrected_99) %>%
  identity() -> m_predictive_residuals


m_dat %>%
  select(platform, min_ct) %>%
  expand(platform = unique(platform),
         #min_ct = modelr::seq_range(min_ct, n = 100),
         min_ct = 5:40) %>%
  add_predicted_draws(m_dif) %>%
  mutate(min_ct = round(min_ct)) %>%
  group_by(platform, min_ct) %>%
  # choose threshold for calling outliers (e.g., 99% or 95% or 90% or 80%)
  tidybayes::median_hdi(.prediction, .width = 0.95) %>%
  #
  ungroup() %>%
  select(platform, min_ct, .prediction, .upper, .lower, .width) %>%
  rename_at(.vars = vars(.prediction, .upper, .lower), .funs = ~ paste0("target_ct_dif",.x)) %>%
  # annotate platforms with sign of Ct value comparison
  mutate(platform = case_when(platform == "8800" ~ "8800 (E Gene - ORF1ab)",
                              platform == "Fluidigm" ~ "Fluidigm (N2 Gene - N1 Gene)",
                              platform == "Simplexa" ~ "Simplexa (S Gene - ORF1ab)",
                              platform == "Thermo" ~ "Thermo (S Gene - ORF1ab)")) %>%
  arrange(platform, min_ct) %>%
  identity() -> m_rules


m_rules %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_outlier_rules_",Sys.Date()),
    output_extension = ".xlsx",
    button_label = "Download Outlier Rules as XLSX",
    button_type = "success",
    has_icon = TRUE,
    icon = "fa fa-save"
  )


m_rules %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_outlier_rules_",Sys.Date()),
    output_extension = ".csv",
    button_label = "Download Outlier Rules as CSV",
    button_type = "info",
    has_icon = TRUE,
    icon = "fa fa-save",
    csv2 = FALSE
  )






```
  
  
-----  
\newpage
<P style = "page-break-before: always">
  
### Export Target Failures & Outliers      
  
```{r ct_diff_model_results, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

ct_drop %>%
  filter(dropout == TRUE) %>%
  select(accession_number, year_week, platform, contains("target")) %>%
  pivot_longer(cols = contains("target")) %>%
  group_by(accession_number, year_week, platform) %>%
  summarise(min_ct = min(value, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(specimen_class = "Target Failure") %>%
  identity() -> target_failures



m_predictive_residuals %>% 
  # choose threshold for calling outliers (e.g., 99% or 95% or 90% or 80%)
  filter(p_residual_corrected >= 0.95) %>%
  #
  select(accession_number, year_week, platform, target_dif, min_ct, p_residual_corrected) %>%
  mutate(specimen_class = "Outlier") %>%
  bind_rows(target_failures) %>%
  # annotate platforms with sign of Ct value comparison
  mutate(platform = case_when(platform == "8800" ~ "8800 (E Gene - ORF1ab)",
                              platform == "Fluidigm" ~ "Fluidigm (N2 Gene - N1 Gene)",
                              platform == "Simplexa" ~ "Simplexa (S Gene - ORF1ab)",
                              platform == "Thermo" ~ "Thermo (S Gene - ORF1ab)")) %>%
  arrange(year_week, accession_number) %>%
  identity() -> outliers_and_target_failures


outliers_and_target_failures %>%
  filter((min_ct <= 28 & (grepl("8800",platform) | grepl("Simplexa",platform) | grepl("Thermo",platform))) | (min_ct <= 20 & grepl("Fluidigm",platform))) %>%
  identity() -> outliers_and_target_failures_low_Ct



```
  
  
Export all (`r outliers_and_target_failures %>% filter(specimen_class == "Target Failure") %>% nrow(.)`) target failures and (`r outliers_and_target_failures %>% filter(specimen_class == "Outlier") %>% nrow(.)`) outliers:  
  
```{r ct_diff_model_export_all, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}


outliers_and_target_failures %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_target_failures_and_outliers_all_",Sys.Date()),
    output_extension = ".xlsx",
    button_label = "Download Data as XLSX",
    button_type = "success",
    has_icon = TRUE,
    icon = "fa fa-save"
  )


outliers_and_target_failures %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_target_failures_and_outliers_all_",Sys.Date()),
    output_extension = ".csv",
    button_label = "Download Data as CSV",
    button_type = "info",
    has_icon = TRUE,
    icon = "fa fa-save",
    csv2 = FALSE
  )



```
  
<br>    
  
Export (`r outliers_and_target_failures_low_Ct %>% filter(specimen_class == "Target Failure") %>% nrow(.)`) target failures and (`r outliers_and_target_failures_low_Ct %>% filter(specimen_class == "Outlier") %>% nrow(.)`) outliers with low Ct values (< 28 for 8800, Simplexa, and Thermo; < 20 for Fluidigm):  
  
```{r ct_diff_model_export_low_Ct, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

outliers_and_target_failures_low_Ct %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_target_failures_and_outliers_low_Ct_",Sys.Date()),
    output_extension = ".xlsx",
    button_label = "Download Data as XLSX",
    button_type = "success",
    has_icon = TRUE,
    icon = "fa fa-save"
  )


outliers_and_target_failures_low_Ct %>%
  download_this(
    output_name = paste0("SARS-CoV-2_Ct_target_failures_and_outliers_low_Ct_",Sys.Date()),
    output_extension = ".csv",
    button_label = "Download Data as CSV",
    button_type = "info",
    has_icon = TRUE,
    icon = "fa fa-save",
    csv2 = FALSE
  )



```

  
<br>    
  
  
```{r ct_diff_model_export_table, echo=FALSE, eval=TRUE, include=TRUE, message=FALSE, error=FALSE, fig.align='center', out.width="120%"}

outliers_and_target_failures_low_Ct %>%
  gt() %>%
  gt::cols_label(accession_number = "Accession Number",
                 year_week = "Year - Week",
                 platform = "Platform",
                 min_ct = "Minimum Ct",
                 target_dif = "Target Difference (Ct)",
                 specimen_class = "Class") %>%
  gt::opt_table_lines(extent = "all") %>%
  gt::cols_hide("p_residual_corrected") %>%
  gt::fmt_missing(columns = 4:tidyselect::last_col())




```
  
  
-----  
\newpage
<P style = "page-break-before: always">