---
title: "Van Wert et al 2023 Fig 1"
author: "Jacey C. Van Wert"
output: html_document
date: "2023-03-28"
---

##Set up and load libraries
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(here)
library(tidyverse)
library(fishualize)
library(GGally)
library(patchwork)
```


## Write SEM function
```{r}
se <- function(x, na.rm=FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x)/length(x))}
```


## Create plot styling
```{r}
plot_style_legend <- theme_classic() +
  theme(
    axis.text.y = element_text(size = 18, colour = "black"),
    axis.text.x = element_text(size = 18, colour = "black"),
    axis.line.y = element_line(colour = "black", size = 0.5),
    axis.line.x = element_line(colour = "black", size = 0.5),
    axis.ticks.y = element_line(size = 0.5),
    axis.ticks.x = element_line(size = 0),
    axis.line = element_line(colour = "black"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.title.y = element_text(size = 18),
    axis.title.x = element_blank(),
    strip.text = element_text(size = 18),
    plot.title = element_text(hjust = 0.5, size = 20),
    legend.title = element_text(size = 18),
    legend.text = element_text(size = 18)
  )

new_palette_5 <- c("#4ba1a4","#95A662","#DAAB18", "#EE7E09","#91d3c5")

```


##Load data
```{r}
nutrients<- read.csv(here::here("data","VanWert_etal_2023_fishfeces.csv"))
```


## Summarise data for stacked bar plot
```{r}
nutrients_sum_tro <- nutrients %>%
  filter(study == "ucsb_2018" | study == "cas_2019") %>%
  dplyr::select(fish_id, code, wt_mass_g, trophic_group, family, water_percent, protein_percent, carb_percent, lipid_percent, ash_percent) %>%
  dplyr::group_by(code, trophic_group) %>%
  dplyr::summarize(
    perc_water_avg = mean(water_percent, na.rm = TRUE),
    perc_protein_avg = mean(protein_percent, na.rm = TRUE),
    perc_carb_avg = mean(carb_percent, na.rm = TRUE),
    perc_lipid_avg = mean(lipid_percent, na.rm = TRUE),
    perc_ash_avg = mean(ash_percent, na.rm = TRUE),
    N_Value_water = length(water_percent[!is.na(water_percent)]),
    N_Value_protein = length(protein_percent[!is.na(protein_percent)]),
    N_Value_carb = length(carb_percent[!is.na(carb_percent)]),
    N_Value_lipid = length(lipid_percent[!is.na(lipid_percent)]),
    N_Value_ash = length(ash_percent[!is.na(ash_percent)])
  ) %>%
  filter(
    N_Value_water > 3,
    N_Value_carb > 3,
    N_Value_protein > 3,
    N_Value_lipid > 3,
    N_Value_ash > 3
  )

nutrients_sum_tro_norm <- nutrients_sum_tro %>%
  filter(
    N_Value_water > 3,
    N_Value_carb > 3,
    N_Value_protein > 3,
    N_Value_lipid > 3,
    N_Value_ash > 3
  ) %>%
  dplyr::group_by(trophic_group) %>%
  dplyr::summarise(
    sum = sum(perc_water_avg + perc_protein_avg + perc_carb_avg + perc_lipid_avg + perc_ash_avg),
    x100_water = perc_water_avg / sum * 100,
    x100_protein = perc_protein_avg / sum * 100,
    x100_carb = perc_carb_avg / sum * 100,
    x100_lipid = perc_lipid_avg / sum * 100,
    x100_ash = perc_ash_avg / sum * 100
  )

nutrients_sum_norm_long_tro <- nutrients_sum_tro_norm %>% 
  pivot_longer(cols = x100_water:x100_ash, names_to="nutrient", values_to = "perc") 
```

## Create stacked bar plot
```{r}
tro_dry <- nutrients_sum_norm_long_tro %>%
  mutate(percent_nutrient = case_when(
    nutrient %in% ("x100_ash") ~ "Ash",
    nutrient %in% ("x100_carb") ~ "Carbohydrates",
    nutrient %in% ("x100_lipid") ~ "Lipids",
    nutrient %in% ("x100_protein") ~ "Proteins",
    nutrient %in% ("x100_water") ~ "Water"
  )) %>%
  ggplot(aes(x = trophic_group, y = perc, fill = percent_nutrient)) +
  geom_bar(stat = "identity", width = .7) +
  labs(y = "Feces Weight (%)") +
  scale_fill_manual(values = new_palette_5, "Macronutrient") +
  plot_style_legend

```


## Summarise data for donut plots
```{r}
nutrients_sum_tro_norm_no_water <- nutrients_sum_tro %>%
  filter(
    N_Value_water > 3,
    N_Value_carb > 3,
    N_Value_protein > 3,
    N_Value_lipid > 3,
    N_Value_ash > 3
  ) %>%
  group_by(trophic_group) %>%
  summarise(
    sum = sum(perc_protein_avg + perc_carb_avg + perc_lipid_avg + perc_ash_avg),
    x100_protein = perc_protein_avg / sum * 100,
    x100_carb = perc_carb_avg / sum * 100,
    x100_lipid = perc_lipid_avg / sum * 100,
    x100_ash = perc_ash_avg / sum * 100
  )


nutrients_useful_long <- nutrients_sum_tro_norm_no_water %>% 
  pivot_longer(cols = x100_protein:x100_ash, names_to="nutrient", values_to = "perc") 

```

```{r}
# Compute the cumulative percentages (top of each rectangle)
nutrients_useful_long_coralli<- nutrients_useful_long %>%filter(trophic_group == "Corallivore")
nutrients_useful_long_herbi<- nutrients_useful_long %>%filter(trophic_group == "Herbivore")
nutrients_useful_long_detri<- nutrients_useful_long %>%filter(trophic_group == "Detritivore")
nutrients_useful_long_plankti<- nutrients_useful_long %>%filter(trophic_group == "Planktivore")
nutrients_useful_long_inverti<- nutrients_useful_long %>%filter(trophic_group == "Invertivore")


hsize <- 2 #size of donut 

coralli <- ggplot(nutrients_useful_long_coralli, aes(x = hsize, y = perc, fill = nutrient)) +
  geom_col() +
  coord_polar(theta = "y") +
  xlim(c(0.2, hsize + 0.5)) +
  scale_fill_manual(values = new_palette_5, "Macronutrient", guide = "none") +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.grid = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    text = element_text(size = 20)
  )


herbi <- ggplot(nutrients_useful_long_herbi, aes(x = hsize, y = perc, fill = nutrient)) +
  geom_col() +
  coord_polar(theta = "y") +
  xlim(c(0.2, hsize + 0.5)) +
  scale_fill_manual(values = new_palette_5, "Macronutrient", guide = "none") +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.grid = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    text = element_text(size = 20)
  )


inverti <- ggplot(nutrients_useful_long_inverti, aes(x = hsize, y = perc, fill = nutrient)) +
  geom_col() +
  coord_polar(theta = "y") +
  xlim(c(0.2, hsize + 0.5)) +
  scale_fill_manual(values = new_palette_5, "Macronutrient", guide = "none") +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.grid = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    text = element_text(size = 20)
  )


detri <- ggplot(nutrients_useful_long_detri, aes(x = hsize, y = perc, fill = nutrient)) +
  geom_col() +
  coord_polar(theta = "y") +
  xlim(c(0.2, hsize + 0.5)) +
  scale_fill_manual(values = new_palette_5, "Macronutrient", guide = "none") +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.grid = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    text = element_text(size = 20)
  )


plankti <- ggplot(nutrients_useful_long_plankti, aes(x = hsize, y = perc, fill = nutrient)) +
  geom_col() +
  coord_polar(theta = "y") +
  xlim(c(0.2, hsize + 0.5)) +
  scale_fill_manual(values = new_palette_5, "Macronutrient", guide = "none") +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.grid = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    text = element_text(size = 20)
  )

```


## Compute percent in donut for each trophic group
```{r}
nutrients_useful_long_coral <- nutrients_useful_long %>%
  group_by(nutrient, trophic_group)%>%
  summarise(total = round(sum(perc),1))

```


## Combine and save figure
```{r}
layout <-
  "
FFFFF#
ABCDE#
"
wrap_plots(A = coralli, B = detri, C = herbi, D = inverti, E = plankti, F = tro_dry, design = layout)
ggsave(("fig1.png"), width = 12, height = 8)

```

