---
title: "Hanlon_2023_Figure_3"
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(broom)
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)

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")
    )
}


x3 <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_3")
x3$Plot <- as.character(x3$Plot)
x3$Year <- as.character(x3$Year)
x3$Plant <- as.character(x3$Plant)

df <- x3 %>% 
  distinct() %>% 
  arrange(Element) %>% 
  pivot_wider(names_from = Element, values_from = Net) %>% 
  mutate(across(where(is.numeric), ~ as.numeric(scale(., scale = FALSE)))) %>% 
  drop_na
df2 <- df %>% dplyr::select(where(is.numeric))

mc <- Mclust(df2)  
summary(mc)
mc$modelName                # Optimal selected model ==> "VVE"
mc$G                        # Optimal number of cluster => 8

clust_plot <- fviz_mclust(mc, "classification", geom = "point", 
            pointsize = 2, 
            palette = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E68310", "#008695", "#CF1C90",  "#4B4B8F"),
            main = NULL,
            shape = 16,
            xlab = "Dimension 1 (28.5%)",
            ylab = "Dimension 2 (12.9%)",
            alpha = 0.5) +
  theme_mh() +
  theme(plot.title = element_blank(),
        plot.subtitle = element_blank(),
        legend.position = "top") +
    guides(color = guide_legend(ncol = 8, title.position = "left"))

clust_plot


# put cluster ID matched with p1_wide2
clusterID <- as_tibble(mc$classification) %>% 
  rename(Cluster = value)

p1_all <- bind_cols(clusterID, df) # see what variables influence cluster involvement

# look at cluster breakdown by treatment and genotype - maybe add in inbred/hybrid?
pc <- p1_all %>% 
  dplyr::select(Cluster, Plot, Genotype, Location, N_Trmt, Water_Trmt, Trmt, Plant, type)



clus_mod <- lm(Cluster ~  Location * Genotype * type* N_Trmt * Water_Trmt, data = pc)
summary(aov(clus_mod))
sa <- tidy(aov(clus_mod))
sa

pc$Cluster <- as.factor(pc$Cluster)

cp4 <- ggplot(pc)+
  geom_bar(aes(x = Genotype, fill = Cluster)) +
  scale_fill_manual(values= c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E68310", "#008695", "#CF1C90",  "#4B4B8F")) +
  theme_mh () +
  coord_flip() +
  theme(legend.position = "top",
        axis.title.x = element_blank()) +
  guides(fill = guide_legend(ncol = 4, title.position = "top"))
cp4

cp5 <- ggplot(pc)+
  geom_bar(aes(x = Cluster, fill = Trmt)) +
  scale_fill_manual(values= c("#1F78B4", "#A6CEE3", "#33A02C", "#B2DF8A"), name = "Treatment") +
  theme_mh () +
  guides(fill = guide_legend(ncol = 2, title.position = "top")) +
  coord_flip() +
  theme(legend.position = "top",
        axis.title.x = element_blank())
cp5

cp6 <- ggplot(pc)+
  geom_bar(aes(x = Cluster, fill = Location)) +
  scale_fill_manual(values= c("#E58606", "#5D69B1", "#52BCA3", "#A5AA99")) +
  guides(fill = guide_legend(ncol = 2, title.position = "top")) +
  theme_mh () +
  coord_flip() +
  theme(legend.position = "top",
        axis.title.x = element_blank(),
        axis.title.y = element_blank())
cp6



clust_plot + 
  theme(legend.position = "none") + 
  cp4 +
  theme(legend.position = "right",
        axis.title.y = element_blank(),
        legend.direction = "vertical",
        legend.box = "vertical") +
  guides(fill = guide_legend(direction='vertical')) +
  
  plot_layout(ncol = 2, width = c(4, 3.5))

ggsave("Figure3_AB.png", width = 8, height = 4, units = "in", dpi = 300)

cp5 + cp6 +
  plot_layout(ncol = 3, width = c(3, 3))


ggsave("Figure3_CD.png", width = 6, height = 4, units = "in", dpi = 300)
```
