---
title: "Hanlon_2023_Supplemental_Figure_8"
author: "Meredith Hanlon"
date: "2023-04-14"
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")
    )
}
```



```{r}

pall2 <- read_excel("Hanlon_2023_data.xlsx", sheet = "SFigure_6-8_1") %>% 
  filter(Location == "ARL")

cores <-  read_excel("Hanlon_2023_data.xlsx", sheet = "SFigure_6-8_2") %>% 
  filter(Site == "ARL")

```


```{r}
cores$Plot <- as.character(cores$Plot)

all <- cores %>% 
  dplyr::select(Plot, Core, Trmt, Genotype, Segment, Length, SA, Weight) 

all$Length <- as.numeric(all$Length)
all$SA<- as.numeric(all$SA)

tots <- all %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  summarise(total_length = sum(Length, na.rm = TRUE))

deep40 <- all %>% group_by( Genotype, Plot, Trmt, Core) %>% 
  filter(Segment %in% c(5,6)) %>% 
  summarise(deep40 = sum(Length, na.rm = TRUE))
deep30 <- all %>% group_by( Genotype, Plot, Trmt, Core) %>% 
  filter(Segment %in% c(4,5,6)) %>% 
  summarise(deep30 = sum(Length, na.rm = TRUE))
shallow30 <- all %>% group_by( Genotype, Plot, Trmt, Core) %>% 
  filter(Segment %in% c(1,2,3)) %>% 
  summarise(shallow30 = sum(Length, na.rm = TRUE))
shallow20 <- all %>% group_by( Genotype, Plot, Trmt, Core) %>% 
  filter(Segment %in% c(1,2)) %>% 
  summarise(shallow20 = sum(Length, na.rm = TRUE))

ls <- left_join(tots, deep40) %>% 
  left_join(., deep30) %>% 
  left_join(., shallow30) %>% 
  left_join(., shallow20)

percs <- ls %>% 
  group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(perc_deep_30 = 100*(deep30/total_length),
         perc_deep_40 = 100*(deep40/total_length))

# change the deep lengths from NA to 0 - as they should be
percs[is.na(percs)] <- 0

dvals <- all %>% group_by(Genotype, Plot, Trmt, Core) %>%
  summarise(d95_l = 0.95*sum(Length, na.rm = TRUE),
            d90_l = 0.90*sum(Length, na.rm = TRUE),
            d80_l = 0.80*sum(Length, na.rm = TRUE),
            d50_l = 0.50*sum(Length, na.rm = TRUE),
            d20_l = 0.20*sum(Length, na.rm = TRUE))

dcalcs <- left_join(all, dvals)
dcalcs <- dcalcs %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  replace_na(list(Length = 0)) %>% 
  mutate(cumul_length = cumsum(Length))

d95 <- dcalcs %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(segment_length = Segment *10) %>% #only use if you have segment numbers versus depths
  mutate(d95_cal = (10*(d95_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d95_cal2 = replace_na(d95_cal, 0)) %>% 
  filter(d95_cal2>=0 & d95_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d95 = (ifelse(d95_cal2 == 0, (d95_l/(cumul_length/10)), d95_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d95)

d90 <- dcalcs %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d90_cal = (10*(d90_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d90_cal2 = replace_na(d90_cal, 0)) %>% 
  filter(d90_cal2>=0 & d90_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d90 = (ifelse(d90_cal2 == 0, (d90_l/(cumul_length/10)), d90_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d90)

d80 <- dcalcs %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d80_cal = (10*(d80_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d80_cal2 = replace_na(d80_cal, 0)) %>% 
  filter(d80_cal2>=0 & d80_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d80 = (ifelse(d80_cal2 == 0, (d80_l/(cumul_length/10)), d80_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d80)

d50 <- dcalcs %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d50_cal = (10*(d50_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d50_cal2 = replace_na(d50_cal, 0)) %>% 
  filter(d50_cal2>=0 & d50_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d50 = (ifelse(d50_cal2 == 0, (d50_l/(cumul_length/10)), d50_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d50)

d20 <- dcalcs %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d20_cal = (10*(d20_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d20_cal2 = replace_na(d20_cal, 0)) %>% 
  filter(d20_cal2>=0 & d20_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d20 = (ifelse(d20_cal2 == 0, (d20_l/(cumul_length/10)), d20_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d20)

dmes1 <- left_join(d95, d90)
dmes2 <- left_join(dmes1, d80)
dmes3 <- left_join(dmes2, d50)
dmes4 <- left_join(dmes3, d20)

d90d20 <-  dmes4 %>% group_by(Genotype, Plot, Trmt, Core) %>% 
  mutate(d90_to_d20 = d90/d20)

depths <- left_join(dmes4, d90d20)

depths$Core[depths$Core == 0] <- NA
percs$Core[percs$Core == 0] <- NA
all_coresumm <- left_join(depths, percs, by = c("Plot", "Genotype", "Trmt", "Core"))
all_coresumm$Plot <- as.factor(all_coresumm$Plot)

#remove cores with no length
pa_cores <- all_coresumm %>% 
  filter(total_length !=0)

# add inbred/hybrid designation
pa_cores$type <- ifelse(grepl("X", pa_cores$Genotype, ignore.case = T), "hybrid", "inbred")

# add a column indicating class from 2018
cores_plant <- pa_cores %>% 
  mutate(Pre_class = ifelse(Genotype %in% c("B101", "PHG50", "PHZ51", "DKPB80 X 3IIH6","PHG50 X PHG72", "PHZ51 X LH82"), "deep18", "shallow18"))
```


```{r}

all <- cores %>% 
  dplyr::select(Plot, Core, Trmt, Genotype, Segment, Length, SA, Weight)

all$Length <- as.numeric(all$Length)
all$SA<- as.numeric(all$SA)

tots <- all %>% group_by(Genotype, Plot, Trmt) %>% 
  summarise(total_length = sum(Length, na.rm = TRUE))

deep40 <- all %>% group_by( Genotype, Plot, Trmt) %>% 
  filter(Segment %in% c(5,6)) %>% 
  summarise(deep40 = sum(Length, na.rm = TRUE))
deep30 <- all %>% group_by( Genotype, Plot, Trmt) %>% 
  filter(Segment %in% c(4,5,6)) %>% 
  summarise(deep30 = sum(Length, na.rm = TRUE))
shallow30 <- all %>% group_by( Genotype, Plot, Trmt) %>% 
  filter(Segment %in% c(1,2,3)) %>% 
  summarise(shallow30 = sum(Length, na.rm = TRUE))
shallow20 <- all %>% group_by( Genotype, Plot, Trmt) %>% 
  filter(Segment %in% c(1,2)) %>% 
  summarise(shallow20 = sum(Length, na.rm = TRUE))

ls <- left_join(tots, deep40) %>% 
  left_join(., deep30) %>% 
  left_join(., shallow30) %>% 
  left_join(., shallow20)

percs <- ls %>% 
  group_by(Genotype, Plot, Trmt) %>% 
  mutate(perc_deep_30 = 100*(deep30/total_length),
         perc_deep_40 = 100*(deep40/total_length))

# change the deep lengths from NA to 0 - as they should be
percs[is.na(percs)] <- 0

dvals <- all %>% group_by(Genotype, Plot, Trmt) %>%
  summarise(d95_l = 0.95*sum(Length, na.rm = TRUE),
            d90_l = 0.90*sum(Length, na.rm = TRUE),
            d80_l = 0.80*sum(Length, na.rm = TRUE),
            d50_l = 0.50*sum(Length, na.rm = TRUE),
            d20_l = 0.20*sum(Length, na.rm = TRUE))

dcalcs <- left_join(all, dvals)
dcalcs <- dcalcs %>% 
  dplyr::select(-Core, -SA, -Weight) %>% 
  group_by(Genotype, Plot, Trmt, Segment) %>% 
  replace_na(list(Length = 0)) %>% 
  summarise(Length = sum(Length)) %>% 
  mutate(cumul_length = cumsum(Length)) %>% 
  left_join(dvals)

d95 <- dcalcs %>% group_by(Genotype, Plot, Trmt) %>% 
  mutate(segment_length = Segment *10) %>% #only use if you have segment numbers versus depths
  mutate(d95_cal = (10*(d95_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d95_cal2 = replace_na(d95_cal, 0)) %>% 
  filter(d95_cal2>=0 & d95_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d95 = (ifelse(d95_cal2 == 0, (d95_l/(cumul_length/10)), d95_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d95)

d90 <- dcalcs %>% group_by(Genotype, Plot, Trmt) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d90_cal = (10*(d90_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d90_cal2 = replace_na(d90_cal, 0)) %>% 
  filter(d90_cal2>=0 & d90_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d90 = (ifelse(d90_cal2 == 0, (d90_l/(cumul_length/10)), d90_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d90)

d80 <- dcalcs %>% group_by(Genotype, Plot, Trmt) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d80_cal = (10*(d80_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d80_cal2 = replace_na(d80_cal, 0)) %>% 
  filter(d80_cal2>=0 & d80_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d80 = (ifelse(d80_cal2 == 0, (d80_l/(cumul_length/10)), d80_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d80)

d50 <- dcalcs %>% group_by(Genotype, Plot, Trmt) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d50_cal = (10*(d50_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d50_cal2 = replace_na(d50_cal, 0)) %>% 
  filter(d50_cal2>=0 & d50_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d50 = (ifelse(d50_cal2 == 0, (d50_l/(cumul_length/10)), d50_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d50)

d20 <- dcalcs %>% group_by(Genotype, Plot, Trmt) %>% 
  mutate(segment_length = Segment *10) %>% 
  mutate(d20_cal = (10*(d20_l-lag(cumul_length))/(cumul_length-lag(cumul_length)))) %>% 
  mutate(d20_cal2 = replace_na(d20_cal, 0)) %>% 
  filter(d20_cal2>=0 & d20_cal2 <10) %>%
  slice(which.max(segment_length)) %>% 
  mutate(d20 = (ifelse(d20_cal2 == 0, (d20_l/(cumul_length/10)), d20_cal2 + segment_length - 10))) %>% 
  dplyr::select(Plot, d20)

dmes1 <- left_join(d95, d90)
dmes2 <- left_join(dmes1, d80)
dmes3 <- left_join(dmes2, d50)
dmes4 <- left_join(dmes3, d20)

d90d20 <-  dmes4 %>% group_by(Genotype, Plot, Trmt) %>% 
  mutate(d90_to_d20 = d90/d20)

depths <- left_join(dmes4, d90d20)

all_coresumm <- left_join(depths, percs, by = c("Plot", "Genotype", "Trmt"))
all_coresumm$Plot <- as.factor(all_coresumm$Plot)

#remove cores with no length
pa_cores <- all_coresumm %>% 
  filter(total_length !=0)

# add inbred/hybrid designation
pa_cores$type <- ifelse(grepl("X", pa_cores$Genotype, ignore.case = T), "hybrid", "inbred")

# add a column indicating class from 2018
cores_plot <- pa_cores %>% 
  mutate(Pre_class = ifelse(Genotype %in% c("B101", "PHG50", "PHZ51", "DKPB80 X 3IIH6","PHG50 X PHG72", "PHZ51 X LH82"), "deep18", "shallow18"))

```


# Process

```{r}
cs <- all %>% 
  dplyr::select(Trmt, Genotype, Segment, Length) %>% 
  group_by(Genotype, Trmt, Segment) %>% 
  summarise(SE = stderror(Length),
            root_length = mean(Length, na.rm = TRUE))

cs2 <- all_coresumm %>%   
  dplyr::select(Genotype, Trmt, d95) %>% 
  pivot_longer(d95, names_to = "trait", values_to = "dtr") %>%
  group_by(Genotype, Trmt, trait) %>% 
  summarise(SE = stderror(dtr),
            dtrs = mean(dtr, na.rm = TRUE))

cs$Segment <- as.numeric(cs$Segment) * 10
```


```{r, fig.width = 8, fig.height = 4}
ggplot() +
  geom_vline(data = cs2, mapping = aes(xintercept = dtrs, color = Trmt), linetype = 2, show.legend = F) +
  geom_rect(data = cs2, aes(xmin = dtrs-SE, xmax = dtrs+SE, ymin = 0, ymax = 200, fill = Trmt), alpha = 0.07) +
  geom_point(data = cs, mapping = aes(x = Segment, y = root_length, color = Trmt), alpha = 0.5, size = 2) +
  geom_line(data = cs, mapping = aes(x = Segment, y = root_length, color = Trmt), show.legend = F) + 
  geom_errorbar(data = cs, mapping = aes(x = Segment, y = root_length, ymin = root_length - SE, 
                    ymax = root_length + SE, color = Trmt), 
                alpha = 0.5, width = 0, show.legend = F) + 
  scale_x_reverse() +
  facet_wrap(.~Genotype, ncol = 6) +
  coord_flip() +
  scale_color_manual(values = c("#1F78B4", "#FF7F00"), name = "Treatment") +
  scale_fill_manual(values = c("#1F78B4", "#FF7F00"), guide = "none") +
  xlab("Core Segment Depth (cm)") +
  ylab("Root Length (cm)") +
  theme_mh() +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -10),
        legend.key.size = unit(0, 'lines'))
```

# XRF correlations

```{r}
#get plot level xrf data
pp <- pall2 %>% 
  group_by(Plot) %>% 
  summarise(across(Cl:S, ~mean(., na.rm = TRUE))) %>% 
  mutate(Plot = as.factor(Plot))

cx <- left_join(pa_cores, pp)

cxs <- cx %>% 
  pivot_longer(Cl:S, 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()

cor_all <- cxs %>%
  pivot_longer(cols = total_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element) %>% 
  dplyr::summarize(cos = cor(Count, root_value))%>% 
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) 
```

```{r, fig.width = 7, fig.height = 5}
cor_all2 <- cor_all %>% 
  filter(root_trait %in% c("total_length", "shallow20", "deep30", "deep40", "d50", "d80", "d95"))

cor_all2$root_trait <- factor(cor_all2$root_trait, 
                          levels = c("total_length", "shallow20", "deep30", "deep40", "d50", "d80", "d95"),
                          labels = c("Total", "<20cm",">30cm", ">40cm", "d50", "d80", "d95"))
```

```{r}
cp <- cor_all2 %>% 
  ggplot(aes(x = root_trait, y = Element, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(limits = c(0,1), expand = c(0,0), 
                     name = expression(R^2)) +
  geom_text(aes(label = r2), color = "black", size = 3) +
  theme_mh() +
  xlab(NULL) +
  theme(legend.position = "top",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -5),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2))) 
cp
```

# genotype by treatment means

```{r}
#get plot level xrf data
pp <- pall2 %>% 
  group_by(Plot) %>% 
  summarise(across(Cl:S, ~mean(., na.rm = TRUE))) %>% 
  mutate(Plot = as.factor(Plot))

cx <- left_join(pa_cores, pp)

cxs <- cx %>% 
  ungroup() %>% 
  dplyr::select(-Plot) %>% 
  group_by(Genotype, Trmt) %>% 
  summarise(across(everything(), ~mean(., na.rm = TRUE))) %>% 
  pivot_longer(Cl:S, names_to = "Element", values_to = "Count") %>% 
  dplyr::select(Genotype, Trmt, total_length, shallow20, deep30, deep40, perc_deep_30, perc_deep_40, d50, d80, d95, Element, Count) %>%
  pivot_longer(cols = total_length:d95, names_to = 'root_trait', values_to = "root_value") 

cor_all <- cxs  %>% 
  drop_na() %>% 
  group_by(root_trait, Element) %>% 
  dplyr::summarize(cos = cor(Count, root_value))%>% 
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) 
```

```{r, fig.width = 7, fig.height = 5}
cor_all2 <- cor_all %>% 
  filter(root_trait %in% c("total_length", "shallow20", "deep30", "deep40", "d50", "d80", "d95"))

cor_all2$root_trait <- factor(cor_all2$root_trait, 
                          levels = c("total_length", "shallow20", "deep30", "deep40", "d50", "d80", "d95"),
                          labels = c("Total", "<20cm",">30cm", ">40cm", "d50", "d80", "d95"))
```

```{r}
cp <- cor_all2 %>% 
  ggplot(aes(x = root_trait, y = Element, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(limits = c(0,1), expand = c(0,0), 
                     name = expression(R^2)) +
  geom_text(aes(label = r2), color = "black", size = 3) +
  theme_mh() +
  xlab(NULL) +
  theme(legend.position = "top",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -5),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2))) 
cp
```

# PCA mean

```{r}
x2 <- cxs %>% 
  distinct() %>% 
  pivot_wider(names_from = Element, values_from = Count) %>% 
  dplyr::select(Genotype:S) %>% 
  drop_na()

pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

# get variations
summary(pca_fit)

```

```{r}

pca_fit %>%
  augment(x2) %>% 
  ggplot(aes(.fittedPC1, .fittedPC2, color = Genotype, shape = Trmt)) + 
  geom_point(size = 3) +
  # scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#A5AA99")) +
  xlab("PC 1 (31.7%)") +
  ylab("PC 2 (29.1%)") +
  theme_mh()
```



```{r}
arrow_style <- arrow(angle = 20, ends = "first", type = "closed", length = grid::unit(8, "pt"))

rots <- pca_fit %>%
  tidy(matrix = "rotation") %>% 
  filter(PC <= 2) %>% 
  rename(Element = column)

rots %>% 
    pivot_wider(names_from = "PC", names_prefix = "PC", values_from = "value") %>%
    ggplot(aes(PC1, PC2)) +
    geom_segment(xend = 0, yend = 0, arrow = arrow_style, color = "grey60") +
    geom_text(aes(label = Element), hjust = 1, color = "#904C2F") +
  theme_mh()
  
```

```{r}

pcd <- pca_fit %>% augment(x2)
rd <- rots %>% 
    # need to scale rots to match
  mutate(value = value * 15) %>% 
  pivot_wider(names_from = "PC", names_prefix = "PC", values_from = "value") 
  

bips1 <- ggplot() +
  geom_segment(data = rd, aes(PC1, PC2), xend = 0, yend = 0, arrow = arrow_style, color = "grey60") +
  geom_point(data = pcd, aes(.fittedPC1, .fittedPC2, color = Genotype, shape = Trmt), size = 3, alpha = 0.1) +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99")) +
  xlab("PC 1 (31.7%)") +
  ylab("PC 2 (29.1%)") +
  geom_text_repel(data = rd, aes(PC1, PC2, label = Element), color = "black", box.padding = 0.1) +
  theme_mh() +
  theme(legend.position = "top",
        legend.box = "vertical",
        legend.margin = margin(t = -5, b = -5))

bips1
```

Check correlation with root traits

 PC1
 

```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC1, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)

ge_pc1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```

PC2

```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC2, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m2 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```

PC3

```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC3, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```

PC1 + PC2


```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC1, pm_pc2 = .fittedPC2, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```



PC1 + PC3


```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC1, pm_pc2 = .fittedPC3, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```


PC2 + PC3


```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC2, pm_pc2 = .fittedPC3, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```


# PCA raw


```{r}

cx <- left_join(pa_cores, pp)

cx$Plot <- as.character(cx$Plot)

cxa <- cx %>% 
  ungroup() %>% 
  dplyr::select(-type) %>% 
  group_by(Genotype, Trmt, Plot) %>% 
  summarise(across(everything(), ~mean(., na.rm = TRUE))) %>% 
  pivot_longer(Cl:S, names_to = "Element", values_to = "Count") %>% 
  dplyr::select(Genotype, Trmt, total_length, shallow20, deep30, deep40, perc_deep_30, perc_deep_40, d50, d80, d95, Element, Count) %>%
  pivot_longer(cols = total_length:d95, names_to = 'root_trait', values_to = "root_value") 

x2 <- cxa %>% 
  pivot_wider(names_from = Element, values_from = Count) %>% 
  dplyr::select(Genotype:S) %>% 
  drop_na()


pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

# get variations
summary(pca_fit)

```

```{r}

pca_fit %>%
  augment(x2) %>% 
  ggplot(aes(.fittedPC1, .fittedPC2, color = Genotype, shape = Trmt)) + 
  geom_point(size = 3) +
  # scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#A5AA99")) +
  xlab("PC 1 (27.2%)") +
  ylab("PC 2 (20.2%)") +
  theme_mh()
```



```{r}
arrow_style <- arrow(angle = 20, ends = "first", type = "closed", length = grid::unit(8, "pt"))

rots <- pca_fit %>%
  tidy(matrix = "rotation") %>% 
  filter(PC <= 2) %>% 
  rename(Element = column)

rots %>% 
    pivot_wider(names_from = "PC", names_prefix = "PC", values_from = "value") %>%
    ggplot(aes(PC1, PC2)) +
    geom_segment(xend = 0, yend = 0, arrow = arrow_style, color = "grey60") +
    geom_text(aes(label = Element), hjust = 1, color = "#904C2F") +
  theme_mh()
  
```

```{r}

pcd <- pca_fit %>% augment(x2)
rd <- rots %>% 
    # need to scale rots to match
  mutate(value = value * 15) %>% 
  pivot_wider(names_from = "PC", names_prefix = "PC", values_from = "value") 
  

bips2 <- ggplot() +
  geom_segment(data = rd, aes(PC1, PC2), xend = 0, yend = 0, arrow = arrow_style, color = "grey60") +
  geom_point(data = pcd, aes(.fittedPC1, .fittedPC2, color = Genotype), size = 3, alpha = 0.1) +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99")) +
  xlab("PC 1 (27.2%)") +
  ylab("PC 2 (20.2%)") +
  geom_text_repel(data = rd, aes(PC1, PC2, label = Element), color = "black", box.padding = 0.1) +
  theme_mh() +
  theme(legend.position = "top",
        legend.box = "vertical",
        legend.margin = margin(t = -5, b = -5))
bips2
```

Check correlation with root traits

PC1

```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC1, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)

all_pc1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```

PC2

```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC2, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m2 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```

PC3

```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC3, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```

PC1 + PC2


```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC1, pm_pc2 = .fittedPC2, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, pm_pc2, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```



PC1 + PC3


```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC1, pm_pc2 = .fittedPC3, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, pm_pc2, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```


PC2 + PC3


```{r}
pca_fit <- x2 %>%
  ungroup() %>% 
  select(Cl:S) %>%
  drop_na() %>% 
  prcomp(center = TRUE, scale = TRUE)

pmpc <- pca_fit %>% 
  augment(x2) %>% 
  pivot_wider(names_from = root_trait, values_from = root_value) %>% 
  dplyr::select(Trmt, Genotype, pm_pc1 = .fittedPC2, pm_pc2 = .fittedPC3, total_length:d95) 

models <- pmpc %>% 
  dplyr::select(Trmt, Genotype, pm_pc1, pm_pc2, total_length:d95) %>% 
  pivot_longer(total_length:d95, names_to = "trait", values_to = "trait_value") %>% 
  dplyr::select(-Trmt, -Genotype) %>% 
  ungroup() %>% 
  nest(data = -trait) %>% 
  mutate(fit = purrr::map(data, ~lm(trait_value ~ ., data = .x)),
         tidied = purrr::map(fit, tidy),
    glanced = purrr::map(fit, glance),
    augmented = purrr::map(fit, augment)) 

models %>% 
  unnest(glanced) %>% 
  paged_table()

m1 <- models %>% unnest(glanced) %>% 
  dplyr::select(trait, r.squared:nobs)
```



# Make the figure

layout will include plot data, two pc biplots, and correlations

panel A is root traits - color by genotype the same as bips and then that will serve as the legend


```{r, fig.width = 8, fig.height = 2.5}
ggplot() +
  geom_vline(data = cs2, mapping = aes(xintercept = dtrs, color = Genotype, linetype = Trmt)) +
  geom_rect(data = cs2, aes(xmin = dtrs-SE, xmax = dtrs+SE, ymin = 0, ymax = 500, fill = Genotype), alpha = 0.07) +
  geom_point(data = cs, mapping = aes(x = Segment, y = root_length, color = Genotype, shape = Trmt), alpha = 0.5, size = 2) +
  geom_line(data = cs, mapping = aes(x = Segment, y = root_length, color = Genotype, linetype = Trmt), show.legend = F) + 
  geom_errorbar(data = cs, mapping = aes(x = Segment, y = root_length, ymin = root_length - SE, 
                    ymax = root_length + SE, color = Genotype), 
                alpha = 0.5, width = 0, show.legend = F) + 
  scale_x_reverse() +
  facet_wrap(.~Genotype, ncol = 6) +
  coord_flip() +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99"), name = "Treatment", guide = "none") +
  scale_fill_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#F97B72", "#4B4B8F", "#A5AA99"), guide = "none") +
  xlab("Core Segment Depth (cm)") +
  ylab("Root Length (cm)") +
  theme_mh() +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -10),
        legend.key.size = unit(0, 'lines'))

ggsave("Supplemental_Figure_8A.png", width = 8, height = 2.5, dpi = 300, units = "in")
```


panel b and c are biplots - raw and plot means - pc1 vs pc2

```{r}
bips2 + bips1 +
plot_layout(width = c(3,3), nrow = 1, ncol = 2, guides = "collect") & theme(legend.position = "none")

ggsave("Supplemental_Figure_8BC.png", width = 5.6, height = 3, units = "in")
```


panel d is correlations

```{r}
pd <- all_pc1 %>% 
  dplyr::select(trait, adj.r.squared, p.value) %>% 
  mutate(type = "all") %>% 
  bind_rows(., ge_pc1 %>% dplyr::select(trait, adj.r.squared, p.value) %>% mutate(type = "avg")) %>% 
  mutate(adj.r.squared = case_when(adj.r.squared <0 ~ 0,
                                   TRUE ~ adj.r.squared))

pd2 <- pd %>% 
  dplyr::select(trait, adjr2 = adj.r.squared, type, p = p.value) %>% 
  filter(trait %in% c("total_length", "shallow20", "deep30", "deep40", "d50", "d80", "d95"))

pd2$trait <- factor(pd2$trait, 
         levels = c("total_length", "shallow20", "deep30", "deep40", "d50", "d80", "d95"),
         labels = c("Total", "<20cm",">30cm", ">40cm", "d50", "d80", "d95"))

pd2$type <- factor(pd2$type, levels = c("all", "avg"), labels = c("All", "Mean"))

pd2 %>% 
  mutate(sig = case_when(p < 0.05 ~ "sig1",
                         p >= 0.05 & p < 0.1 ~ "sig2")) %>% 
  mutate(r3 = case_when(!is.na(sig) ~ adjr2)) %>% 
  ggplot(aes(x = type, y = trait, fill = adjr2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(limits = c(0,1), expand = c(0,0), name = expression(R^2), na.value = "white") +
  geom_text(aes(label = round(r3, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  theme_mh() +
  scale_y_discrete(limits = rev) +
  xlab(NULL) +
  theme(legend.position = "top",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -5),
        axis.ticks.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2))) 


ggsave("Supplemental_Figure_8D.png", width = 2, height = 3, units = "in")

```

