explore_mix.Rmd

---
title: "Explore Mixtures Latest"
output: 
  bookdown::word_document2:
    fig_caption: yes
    toc: yes
editor_options: 
  chunk_output_type: console
---

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

library(tidyverse, quietly = TRUE)
library(rpart, quietly = TRUE)
library(party, quietly = TRUE)
library(partykit, quietly = TRUE)
library(rpartScore, quietly = TRUE)
library(readxl, quietly = TRUE)
library(flextable, quietly = TRUE)
# library(subselect) #sim annealing
library(leaps, quietly = TRUE) #regsubsets
library(glmnet, quietly = TRUE)
library(toxEval)
library(corrplot, quietly = TRUE)

library(here)
source(here("R/mixtures/mix_script.R"))
source(here("R/mixtures/prepare_mixture_data.R"))
source(file = "read_chemicalSummary.R")
source(file = "R/analyze/open_land_use.R")

n_sites <- 7
EAR_thresh <- 0.001

df <- get_final_mixtures(chemicalSummary,
                         EAR_thresh,
                         n_sites)

```

## Final mixture table

```{r joinedEVERYTHING}
df_flex <- flextable(df)
autofit(df_flex)

```


```{r openLand}
df_lu <- open_land_use()

```

```{r corLU, results="asis"}

big_enough <- function(x, thresh = 2){
  max(x, na.rm = TRUE) > thresh
}

df_lu_filtered <- df_lu %>% 
  select(-site) %>% 
  select_if(big_enough) %>% 
  bind_cols(select(df_lu, frac_2010, frac_2014))

M <- cor(df_lu_filtered)
corrplot(M, type="upper", tl.cex = 0.5)
```

Reduce the variables by those that have a corrolation coefficient > 0.9 or < -0.9.

```{r corLU_reduced, results="asis"}
lu_vars <- colnames(M)
lu_vars_rev <- rev(lu_vars)

exclude <- c()
exclude_rev <- c()

for(i in seq_along(lu_vars)){
  
  check_var <- lu_vars[i]
  
  if(check_var %in% exclude){
    next
  } else {
    check_cor <- names(which(M[,check_var] > 0.8 | 
                        M[,check_var] < -0.8))
    check_cor <- check_cor[check_cor != check_var]
    exclude <- c(exclude, check_cor)
  }

}

df_lu_filtered <- df_lu_filtered[,-which(names(df_lu_filtered) %in% exclude)]

M2 <- cor(df_lu_filtered)
corrplot(M2, type="upper", tl.cex = 0.5)
```

Now, we use the sites with detections and filter out variables that don't have a single site left that doesn't have 10% of that category.

```{r corLU_manual, results="asis"}

df_lu_filtered <- bind_cols(df_lu[,c("site","Urban","Crops")],
                            df_lu_filtered) %>% 
  select(-Basin_Area_mi2, -Population_Density)

all_chemicals_in_mixtures <- unique(unlist(df$Chemicals))

big_enough <- function(x, thresh = 10){
  max(x, na.rm = TRUE) > thresh
}

cs_mix <- chemicalSummary %>% 
  filter(chnm %in% all_chemicals_in_mixtures,
         EAR > 0) %>% 
  select(site) %>% 
  distinct() %>% 
  left_join(df_lu_filtered, by = "site") %>% 
  select(-site) %>% 
  select_if(big_enough)

auto_categories <- names(cs_mix)

```

If we were going to fine-tune the categories, we'd do it here!

```{r showCats}
auto_categories
```

```{r trees, results="asis"}

chemicalSummary <- chemicalSummary %>% 
  filter(site != "04010500") #half this watershed is in Canada, so the landuse stuff is/could be way off

for(i in seq_len(nrow(df))){
  
  chems <- unlist(df$Chemicals[i])
  mixture <- paste(chems, collapse = ",")
  endpoint <- df$Assay[i] 
  
  cat("\n")
  cat("\n##", mixture,"\n")
  
  mixture <- paste(chems, collapse = ",\n")
  
  sub_df <- chemicalSummary %>% 
    filter(endPoint == {{endpoint}},
           EAR > 0) %>% 
    group_by(site, shortName, date) %>% 
    summarize(sumEAR = sum(EAR)) %>% 
    group_by(site, shortName) %>% 
    summarize(maxEAR = max(sumEAR)) %>% 
    left_join(df_lu_filtered, by = "site") %>% 
    mutate(mix_st = mixture,
           lowEAR = maxEAR,
           highEAR = maxEAR) %>% 
    ungroup()
  
  cat("\n")
  
  cat("\n### Tree\n")
  form_bigger_lm <- formula(paste("maxEAR ~ ",
                               paste(auto_categories,
                                     collapse = " + ")))
  
  tree_return <- plot_trees(form_bigger_lm, sub_df, endpoint)
  cat("\n")
  
  cat("\n### Linear\n")
  
  new_form_lm <- get_formula(sub_df, auto_categories, 
                             sumEAR = "maxEAR",
                             lasso = FALSE, survival = FALSE)
  
  new_form_surv <- get_formula(sub_df, auto_categories, 
                               sumEAR = "maxEAR",
                               lasso = FALSE)
  
  form_surv <- reformulate(termlabels = attr(terms(new_form_surv),"term.labels"),
              response = 'survival::Surv(lowEAR,
                                      highEAR,
                                      type="interval2")')

 x_df <- plot_lm(new_form_lm, 
                form_surv, 
                sumEAR = "maxEAR",
                sub_df = sub_df)
  
  cat("\n")
  cat("\n#### Log\n")
  cat("\n")
  
  new_form_lm_log <- get_formula(sub_df, auto_categories, 
                             sumEAR = "maxEAR", log = TRUE,
                             lasso = FALSE, survival = FALSE)
  
  new_form_surv_log <- get_formula(sub_df, auto_categories, 
                               sumEAR = "maxEAR", log = TRUE,
                               lasso = FALSE)
  
  x_df_log <- plot_lm(new_form_lm_log, 
                      new_form_surv_log, 
                      sumEAR = "maxEAR",
                      sub_df = sub_df, log = TRUE)
  
  
}


```