---
title: Tracer selection using design specific mixtures
author: Alex Koiter
---
  
## Load Libraries
```{r}  
library(tidyverse)
library(ggfortify)
library(patchwork)
```

## % classified function
```{r}
classify <- function(data, model, by = "site", type = "percent") {
  if(type == "count") x <- table(as.data.frame(data)[,"site"], model$class)
  if(type == "percent") x <- diag(prop.table(table(as.data.frame(data)[,"site"], model$class), 1)) * 100
  x
}
```

## Load in data
```{r}
geo_results <- read.csv(here::here("./notebooks/Geochemistry analysis - Copy 2.csv")) %>%
  pivot_longer(cols = Ag:Zr, names_to = "Fingerprint", values_to = "value") %>%
  filter(Fingerprint %in% c("Ag", "Al", "As","B","Ba","Be","Bi","Ca","Cd","Ce","Co", "Cr", "Cs", "Cu", "Fe", "Ga", "Hf", "Hg", "In", "K", "La", "Li", "Mg", "Mn", "Mo", "Nb", "Ni", "P", "Pb", "Rb", "S", "Sb", "Sc", "Se", "Sn", "Sr", "Te", "Th", "Tl", "U", "V", "Y", "Zn", "Zr")) %>% # excludes fingerprints that are below level of detection
  dplyr::select(-X) %>% # don't need this column
  filter(sample_design %in% c("Grid", "Transect", "Likely to erode")) 

col_results <- read.csv(here::here("./notebooks/final results revised.csv")) %>%
  pivot_longer(cols = X:B, names_to = "Fingerprint", values_to = "value") %>%
  dplyr::select(-X.1) %>% # don't need this column
  filter(sample_design %in% c("Grid", "Transect", "Likely to erode")) %>%
  mutate(Fingerprint = paste0(Fingerprint, "_col")) # appended _col as some of the colour coefficients have the same id eg Boron = B and Blue also = B

f_order <- c("Ag", "Al", "As","B","Ba","Be","Bi","Ca","Cd","Ce","Co", "Cr", "Cs", "Cu", "Fe", "Ga", "Hf", "Hg", "In", "K", "La", "Li", "Mg", "Mn", "Mo", "Nb", "Ni", "P", "Pb", "Rb", "S", "Sb", "Sc", "Se", "Sn", "Sr", "Te", "Th", "Tl", "U", "V", "Y", "Zn", "Zr", "R_col", "G_col", "B_col", "x_col", "y_col", "Y_col", "X_col", "Z_col","L_col", "a_col", "b_col", "u_col", "v_col", "c_col", "h_col")

# Bind data sets
results <- geo_results %>%
  bind_rows(col_results) 
```

## Virtual mixtures
```{r}
proportions <- seq(0, 1, 0.05)

mixtures <- results %>%
  group_by(Fingerprint, site, sample_design) %>%
  summarise(avg =  mean(value)) %>%
  pivot_wider(names_from = site, values_from = avg) %>%
  group_by(sample_design) %>%
  mutate(mix_ag = map(Agriculture, ~.x * (1 - proportions)),
         mix_forest = map(Forest, ~.x * proportions)) %>%
  group_by(Fingerprint, Agriculture, Forest, sample_design) %>%
  summarize(mix = map2(mix_ag, mix_forest, ~data.frame(mix = .x + .y, prop_forest = proportions))) %>%
  unnest(mix) %>%
  pivot_wider(id_cols = c("Fingerprint", "Agriculture", "Forest", "sample_design"), names_from = prop_forest, values_from = mix) 
```

## Range test
```{r}
range <- results %>%
  group_by(site, sample_design, Fingerprint) %>%
  summarise(upper_lim = quantile(value, 0.75),
            lower_lim = quantile(value, 0.25)) %>%
  group_by(Fingerprint, sample_design) %>%
  summarise(lower_lim = min(lower_lim),
            upper_lim = max(upper_lim)) %>%
  ungroup() %>%
  left_join(mixtures) %>%
  dplyr::select(-c(Agriculture, Forest)) %>%
  group_by(Fingerprint, sample_design) %>%
  mutate(pass = all(round(`0`:`1`, 5) >= round(lower_lim, 5)) & all(round(`0`:`1`, 5) <= round(upper_lim, 5)))  ## rounding as differences as small as 10^-18 were being removed otherwise

grid_range <- filter(range, pass == TRUE , sample_design == "Grid") %>%
  arrange(factor(Fingerprint, levels = f_order)) %>%
  pull(Fingerprint)

transect_range <- filter(range, pass == TRUE , sample_design == "Transect") %>%
  arrange(factor(Fingerprint, levels = f_order)) %>%
  pull(Fingerprint)

likely_range <- filter(range, pass == TRUE , sample_design == "Likely to erode") %>%
  arrange(factor(Fingerprint, levels = f_order)) %>%
  pull(Fingerprint)

range_results <- filter(range, pass == TRUE) %>%
  dplyr::select(Fingerprint, sample_design)

```

## Mann-W test
```{r}
wilcox_test <- range_results %>%
  inner_join(results) %>%
  nest(data = c(-sample_design, -Fingerprint)) %>%
  mutate(test = map(data, ~wilcox.test(value ~ site, data = ., exact = FALSE))) %>%
  mutate(p_value = map_dbl(test, ~ .$p.value),
         statistic =  map_dbl(test, ~ .$statistic))

grid_wilcox <- filter(wilcox_test, p_value < 0.01 , sample_design == "Grid") %>%
  arrange(factor(Fingerprint, levels = f_order))
grid_wilcox
pull(grid_wilcox, Fingerprint)

transect_wilcox <- filter(wilcox_test, p_value < 0.01 , sample_design == "Transect") %>%
  arrange(factor(Fingerprint, levels = f_order))
transect_wilcox
pull(transect_wilcox, Fingerprint)
  
likely_wilcox <- filter(wilcox_test, p_value < 0.01 , sample_design == "Likely to erode") %>%
  arrange(factor(Fingerprint, levels = f_order))
likely_wilcox
pull(likely_wilcox, Fingerprint)

wilcox_results <- filter(wilcox_test, p_value < 0.01) %>%
  dplyr::select(Fingerprint, sample_design)
```

## DFA
```{r} 
dfa_test <- wilcox_results %>%
  inner_join(results) %>%
  ungroup() %>%
  nest(data = c(-sample_design)) %>%
  mutate(newdata = map(data, ~pivot_wider(., names_from = Fingerprint, values_from = value))) %>%
  mutate(newdata = map(newdata, ~dplyr::select(., -sample_number))) %>%
  mutate(test = map(newdata, ~ klaR::greedy.wilks(site ~ ., data =. , niveau = 0.1, na.action = "na.omit"))) %>%
  mutate(Fingerprint = map(test, ~ .$results$vars)) %>%
  mutate(wilks = map(test, ~.x$results$Wilks.lambda))

grid_dfa <- dfa_test %>%
  unnest(Fingerprint) %>%
  filter(sample_design == "Grid") %>%
  pull(Fingerprint)
grid_dfa

transect_dfa <- 
  dfa_test %>%
  unnest(Fingerprint) %>%
  filter(sample_design == "Transect") %>%
  pull(Fingerprint)
transect_dfa

likely_dfa <- 
  dfa_test %>%
  unnest(Fingerprint) %>%
  filter(sample_design == "Likely to erode") %>%
  pull(Fingerprint)
likely_dfa
```


## PCA
```{r}
df_grid <- results %>% 
  filter(sample_design == "Grid") %>%
  filter(Fingerprint %in% grid_dfa) %>%
  pivot_wider(names_from = Fingerprint, values_from = value) %>%
  dplyr::select(-sample_design, -sample_number)

df_transect<- results %>%
  filter(sample_design == "Transect") %>%
  filter(Fingerprint %in% transect_dfa) %>%
  pivot_wider(names_from = Fingerprint, values_from = value) %>%
  dplyr::select(-sample_design, -sample_number)

df_likely <- results %>%
  filter(sample_design == "Likely to erode") %>%
  filter(Fingerprint %in% likely_dfa) %>%
  pivot_wider(names_from = Fingerprint, values_from = value) %>%
  dplyr::select(-sample_design, -sample_number)

pca_grid <- prcomp(select(df_grid, 2:ncol(df_grid)), scale. = TRUE)

autoplot(pca_grid)
p1 <- autoplot(pca_grid, 
               data = df_grid, 
               colour = 'site',
               frame = TRUE, 
               frame.type = 'norm') +
  theme_bw() +
  scale_colour_viridis_d(begin = 0, end = 0.6) +
  scale_fill_viridis_d(begin = 0, end = 0.6) +
  theme(legend.title = element_blank()) +
  labs(title = "Grid")
p1

pca_transect <- prcomp(select(df_transect, 2:ncol(df_transect)) , scale. = TRUE)
autoplot(pca_transect)
p2 <- autoplot(pca_transect, 
               data = df_transect, 
               colour = 'site', 
               frame = TRUE, 
               frame.type = 'norm') +
  theme_bw() +
  scale_colour_viridis_d(begin = 0, end = 0.6) +
  scale_fill_viridis_d(begin = 0, end = 0.6) +
  theme(legend.title = element_blank()) +
  labs(title = "Transect")
p2

pca_likely <- prcomp(select(df_likely, 2:ncol(df_likely)) , scale. = TRUE)
autoplot(pca_likely)
p3 <- autoplot(pca_likely, 
               data = df_likely, 
               colour = 'site', 
               frame = TRUE, 
               frame.type = 'norm') +
  theme_bw() +
  scale_colour_viridis_d(begin = 0, end = 0.6) +
  scale_fill_viridis_d(begin = 0, end = 0.6) +
  theme(legend.title = element_blank())+
  labs(title = "Likely to erode")
p3


p4 <- p1 + p2 + p3+
  plot_layout(guides = 'collect') & theme(legend.position = "bottom")
p4

#ggsave(filename = "Figures R1/design_pca.png", plot = p4, width = 174, height = 100, units = "mm", dpi = 600)
```

## Assess accuracy of discrimination 
```{r}
grid.lda <- MASS::lda(site ~ ., data = df_grid, na.action = "na.omit", CV = T)
transect.lda <- MASS::lda(site ~ ., data = df_transect, na.action = "na.omit", CV = T)
likely.lda <- MASS::lda(site ~ ., data = df_likely, na.action = "na.omit", CV = T)

#' ## percent correct for each category 
#' ### Grid 
table(as.data.frame(df_grid)[,"site"], grid.lda$class)
diag(prop.table(table(as.data.frame(df_grid)[,"site"], grid.lda$class), 1))

#' ### Transect
table(as.data.frame(df_transect)[,"site"], transect.lda$class)
diag(prop.table(table(as.data.frame(df_transect)[,"site"], transect.lda$class), 1))

#' ### Likely to erode
table(as.data.frame(df_likely)[,"site"], likely.lda$class)
diag(prop.table(table(as.data.frame(df_likely)[,"site"], likely.lda$class), 1))


#' ##  total percent correct  
#' ### Grid 
sum(diag(prop.table(table(as.data.frame(df_grid)[,"site"], grid.lda$class))))

#' ### Transect
sum(diag(prop.table(table(as.data.frame(df_transect)[,"site"], transect.lda$class))))

#' ### Likely to erode
sum(diag(prop.table(table(as.data.frame(df_likely)[,"site"], likely.lda$class))))
```

## Percent Classified correctly and Wilks Lamda
### Grid design
```{r}
grid_classify <- tibble(n = 1:length(grid_dfa),
            params = map(n, ~grid_dfa[1:.x])) %>%
  mutate(params = map_chr(params, ~paste0(.x, collapse = " + ")),
         f = map(params, ~as.formula(paste0("site ~ ", .x))),
         m = map(f, ~MASS::lda(.x, data = df_grid, na.action = "na.omit", CV = TRUE)),
         c = map_df(m, ~classify(data = df_grid, model = .x))) %>%
  select(-n, -f, -m) %>%
  mutate(wilks = dfa_test$wilks[dfa_test$sample_design == "Grid"][[1]])

grid_classify
```
### Transect design
```{r}
transect_classify <- tibble(n = 1:length(transect_dfa),
                        params = map(n, ~transect_dfa[1:.x])) %>%
  mutate(params = map_chr(params, ~paste0(.x, collapse = " + ")),
         f = map(params, ~as.formula(paste0("site ~ ", .x))),
         m = map(f, ~MASS::lda(.x, data = df_transect, na.action = "na.omit", CV = TRUE)),
         c = map_df(m, ~classify(data = df_transect, model = .x))) %>%
  select(-n, -f, -m) %>%
  mutate(wilks = dfa_test$wilks[dfa_test$sample_design == "Transect"][[1]])

transect_classify
```

### Likely to erode design
```{r}
likely_classify <- tibble(n = 1:length(likely_dfa),
                            params = map(n, ~likely_dfa[1:.x])) %>%
  mutate(params = map_chr(params, ~paste0(.x, collapse = " + ")),
         f = map(params, ~as.formula(paste0("site ~ ", .x))),
         m = map(f, ~MASS::lda(.x, data = df_likely, na.action = "na.omit", CV = TRUE)),
         c = map_df(m, ~classify(data = df_likely, model = .x))) %>%
  select(-n, -f, -m) %>%
  mutate(wilks = dfa_test$wilks[dfa_test$sample_design == "Likely to erode"][[1]])

likely_classify
```

## Plotting
```{r}
plotting3 <- results %>%
  filter(Fingerprint %in% unique(c(grid_dfa, likely_dfa, transect_dfa))) %>%
  group_by(Fingerprint, sample_design, site) %>%
  summarise(avg = mean(value),
            sd = sd(value)) %>%
  rename("Source" = site) %>%
  mutate(Fingerprint = as.factor(Fingerprint)) %>%
  mutate(Fingerprint = fct_recode(Fingerprint, "italic(`a*`)" = "a_col", "italic(`b*`)" = "b_col", "italic(`c*`)" = "c_col", "italic(`h*`)" = "h_col", "italic(`x*`)" = "x_col", "`Ca (%)`" = "Ca", "`Fe (%)`" = "Fe"))



p1 <- ggplot(data = plotting3, aes(x = sample_design, y = avg, colour = Source)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(data = plotting3, aes(ymax = avg + sd, ymin = avg -sd, x = sample_design), position = position_dodge(width = 0.5), width = 0.5) +
  scale_colour_viridis_d(begin = 0, end = 0.6) +
  theme_bw() +
  labs(y = "Concentation (ppm) or value") +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom",
        legend.title = element_blank()) +
  facet_wrap(~Fingerprint, ncol = 4, scales = "free_y", labeller = "label_parsed")

p1

#ggsave(plot = p1, filename = "Figures R1/design_means_sds.png", height = 150, width = 175, units = "mm", dpi = 600)
```
