---
title: "Hanlon_2023_Figure_4"
author: "Meredith Hanlon"
date: "2023-04-11"
output: html_document
---


```{r global_options, echo = FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE, echo=FALSE, error = TRUE)
```

```{r}
library(tidyverse)
library(rstatix)
library(ggpubr)
library(readxl)
library(stats)
library(viridis)
library(patchwork)
library(ggpmisc)
library(ggthemes)
library(ggrepel)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(mclust)
library(broom)
library(RColorBrewer)
library(ggh4x)

stderror <- function(x) sd(x, na.rm = TRUE) / sqrt(length(na.omit(x)))

theme_mh <- function(){ 
    font <- "sans"   #assign font family up front
    theme(
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      panel.border = element_rect(color = "black", fill = NA),
      axis.ticks = element_line(color = "black"),   
      axis.line = element_line(color = "black"),
      strip.background = element_blank(),
      legend.key = element_rect(fill = NA),
      text = element_text(family = font, size = 10, color = "black"),
      axis.text = element_text(family = font, size = 8, color = "black")
    )
}

pall <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_1")
pa_cores <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_2")

pp <- pall %>% 
  group_by(Plot) %>% 
  summarise(across(Cl:Ba, ~mean(., na.rm = TRUE)))

cx <- left_join(pa_cores, pp) %>% 
  mutate(Core = as.numeric(Core))

cxs <- cx %>% 
  pivot_longer(Cl:Ba, names_to = "Element", values_to = "Count") %>% 
  dplyr::select(Genotype, Trmt, Plot, total_length, shallow20, deep30, deep40, perc_deep_30, perc_deep_40, d50, d80, d95, Element, Count) %>% 
  group_by(Genotype, Trmt, Element) %>% 
  distinct()

pa_cores$Core <- as.numeric(pa_cores$Core)
po <- pa_cores %>% 
  filter(!is.na(Core)) %>% 
  unite(cp, Plot, Core) %>% 
  separate(cp, into = c("Plot", "Core")) 

x2 <- cxs %>% 
  distinct() %>% 
  pivot_wider(names_from = Element, values_from = Count) %>% 
  dplyr::select(Genotype, Trmt, Plot, Cl:Ba) %>%
  distinct() %>% 
  dplyr::select(-Co) %>% 
  drop_na()

x2 %>% 
  pivot_longer(Cl:Ba, names_to = "Element", values_to = "count") %>% 
  mutate(count2 = count/10000) %>% 
  filter(Element %in% c("Ca", "Br", "K")) %>% 
  ggplot(aes(x = Trmt, y = count2)) +
  geom_boxplot(outlier.shape = NA) +
  geom_quasirandom(aes(color = Genotype, alpha = 0.7, shape = Trmt)) +
  theme_mh() +
  stat_compare_means(method = "t.test", aes(label = paste0("p = ", after_stat(p.format)))) +
  facet_wrap(.~Element) +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99"))  +
  theme(legend.position = "none",
        strip.text = element_text(size = 9)) +
  xlab("Nitrogen Treatment") +
  ylab("XRF Count x 10,000") +
  scale_y_continuous(expand = expansion(mult = c(0.00, 0.1)))

ggsave("Figure4A.png", dpi = 300, height = 2, width = 3.5)

# show deep 40 split with middle excluded
trl <- po %>% 
  dplyr::select(Genotype, Trmt, Plot, Core, deep40) %>% 
  group_by(Genotype, Trmt, Plot) %>% 
  summarise(deep40 = mean(deep40, na.rm = TRUE)) %>% 
  ungroup()
trl_q <- trl %>% 
  ungroup() %>% 
  summarise(deep40 = quantile(deep40, probs = seq(0, 1, 1/5)), q = c("k", "k", "d", "d","k", "k")) %>% 
  filter(q == "d")
trl2 <- trl %>% 
  filter(deep40 <= min(trl_q$deep40) | deep40 >= max(trl_q$deep40)) 
mv <- mean(trl2$deep40)
trl2 <- trl2 %>% 
  mutate(group = case_when(deep40 >= mv ~ "High",
                           TRUE ~ "Low"))

trl2 %>% 
  ggplot(aes(x = group, y = deep40)) +
  geom_boxplot(outlier.shape = NA) +
  geom_quasirandom(aes(color = Genotype, alpha = 0.7, shape = Trmt)) +
  theme_mh() +
  stat_compare_means(method = "t.test", aes(label = paste0("p = ", after_stat(p.format)))) +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99"))  +
  scale_shape_manual(values = c(16, 17), name = "Nitrogen Treatment") + 
  theme(legend.position = "none",
        strip.text = element_text(size = 9)) +
  xlab("Root length >40cm group") +
  ylab("Root length >40 cm (cm)") +
  scale_y_continuous(expand = expansion(mult = c(0.00, 0.1)))

ggsave("Figure4B.png", dpi = 300, height = 2, width = 3.5)


# Four group
trl_HN <- po %>% 
  dplyr::select(Genotype, Trmt, Plot, Core, deep40) %>% 
  filter(Trmt == "HN") %>% 
  group_by(Genotype, Trmt, Plot) %>% 
  summarise(deep40 = mean(deep40, na.rm = TRUE)) 

trl_q_HN <- trl_HN %>% 
  ungroup() %>% 
  summarise(deep40 = quantile(deep40, probs = seq(0, 1, 1/5)), q = c("k", "k", "d", "d","k", "k")) %>% 
  filter(q == "d")

trl2_HN <- trl_HN %>% 
  filter(deep40 <= min(trl_q_HN$deep40) | deep40 >= max(trl_q_HN$deep40)) 

mv <- mean(trl2_HN$deep40)

trl2_HN2 <- trl2_HN %>% 
  mutate(group = case_when(deep40 >= mv ~ "High_HN",
                           TRUE ~ "Low_HN"))

trl_LN <- po %>% 
  dplyr::select(Genotype, Trmt, Plot, Core, deep40) %>% 
  filter(Trmt == "LN") %>% 
  group_by(Genotype, Trmt, Plot) %>% 
  summarise(deep40 = mean(deep40, na.rm = TRUE)) 

trl_q_LN <- trl_LN %>% 
  ungroup() %>% 
  summarise(deep40 = quantile(deep40, probs = seq(0, 1, 1/5)), q = c("k", "k", "d", "d","k", "k")) %>% 
  filter(q == "d")

trl2_LN <- trl_LN %>% 
  filter(deep40 <= min(trl_q_LN$deep40) | deep40 >= max(trl_q_LN$deep40)) 

mv <- mean(trl2_LN$deep40)

trl2_LN2 <- trl2_LN %>% 
  mutate(group = case_when(deep40 >= mv ~ "High_LN",
                           TRUE ~ "Low_LN"))

grp <- bind_rows(trl2_HN2, trl2_LN2)
grp2 <- grp
grp2$group <- factor(grp$group,
               levels = c("High_HN", "High_LN", "Low_HN", "Low_LN"),
               labels = c("HL, HN", "HL, LN", "LL, HN", "LL, LN"))

stat.test <- grp2 %>% 
  ungroup() %>% 
  dunn_test(deep40 ~ group) 

st <- stat.test %>% add_y_position() %>% 
  filter(p.adj.signif != "ns")


grp2 %>% 
  ggplot(aes(x = group, y = deep40)) +
  geom_boxplot(outlier.shape = NA) +
  geom_quasirandom(aes(color = Genotype, shape = Trmt), alpha = 0.7) +
  theme_mh() +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99"), guide = guide_legend(ncol = 3, nrow = 4, title = NULL))  +
  scale_shape_manual(values = c(16, 17), guide = "none") + 
  theme(legend.position = "bottom",
        strip.text = element_text(size = 9),
        legend.justification = c(1,0),
        legend.text = element_text(size = 7)) +
  xlab("Root length >40cm group, Nitrogen treatment") +
  ylab("Root length >40 cm (cm)") +
  stat_pvalue_manual(st)

ggsave("Figure4C.png", dpi = 300, height = 4, width = 3.5)
```

# Figure 4D

LDAs run separately using this code, organized in a separate spreadsheet

This was done for each trait, each location, and each treatment combination individually

LDAs were run using all data or just the top2 vs bottom2 quantiles. The accuracy and pvalue were extracted and recorded.

This is an example with PA19 deep40 data.

```{r}
x4 <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_4")

# Training set LDA (70/30 split)
index = createDataPartition(y=x4$group, p=0.7, list=FALSE)

xtrain2 = x4[index,]
xtest2 = x4[-index,]

lda.fit <- train(group ~ ., data = xtrain2, method="lda")

pred.group <- predict(lda.fit, xtrain2)
tab <- table(pred.group, xtrain2$group)
confusionMatrix(tab)

# Test set LDA
pred.group <- predict(lda.fit, xtest2)
tab <- table(pred.group, xtest2$group)
confusionMatrix(tab)

# LOOCV on whole set
x5 <- x4 

lda.fit2 <- lda(group ~ ., data = x5, CV = T)
conf2 <- table(list(predicted=lda.fit2$class, observed=x5$group))
caret::confusionMatrix(conf2)
```


```{r}
set.seed(123)

lds <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_3")

lds2 <- lds %>% 
  pivot_longer(cols = c(Training, Training_p, Testing, Testing_p, LOOCV, LOOCV_p), names_to = c("Set"), values_to = c("value")) %>% 
  separate(Set, into = c("Test", "p")) %>% 
  pivot_wider(names_from = p, values_from = value) %>% 
  rename(accuracy = `NA`)

lds2$Data <- factor(lds2$Data,
                    levels = c("All", "top_bottom"),
                    labels = c("All Data", "Top 40%/Bottom 40%"))
lds2$Experiment <- factor(lds2$Experiment,
                           levels = c("ARL19", "CO19", "HAN19", "PA19", "PA20"),
                           labels = c("ARL WI", "CO", "HAN WI", "PA 2019", "PA 2020"))

lds2$Trait <- factor(lds2$Trait, 
         levels = c("total", "shallow20", "deep30", "deep40", "d50", "d80", "d95"),
         labels = c("Total", "<20cm",">30cm", ">40cm", "d50", "d80", "d95"))
lds2$Split <- factor(lds2$Split,
                    levels = c("none", "HN_LN", "WW_WS"),
                    labels = c("All", "Nitrogen", "Water"))

lds2 %>% 
  filter(Data != "All Data") %>% 
  mutate(sig = case_when(p < 0.05 ~ "sig1",
                         p >= 0.05 & p < 0.1 ~ "sig2")) %>% 
  mutate(acc2 = case_when(!is.na(sig) ~ accuracy)) %>% 
  ggplot(aes(x = Trait, y = Test, fill = accuracy)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(limits = c(0,1), expand = c(0,0), name = expression(Accuracy), na.value = "white") +
  geom_text(aes(label = round(acc2, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  theme_mh() +
  scale_x_discrete(position = "top") +
  facet_nested(Experiment+Split ~ ., nest_line = element_line(linetype = 1)) +
  theme(legend.position = "top",
        panel.spacing.y = unit(0.2, "lines"),
        legend.box.margin = margin(b = -10),
        legend.margin = margin(0, 0, 0, 0),
        axis.ticks.x = element_blank(),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = margin(l = 10, r = 10),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = 2)),
        strip.placement = "outside",
        strip.text = element_text(size = 9)) 

ggsave("Figure4D.png", width = 4, height = 8, units = "in")
```

