This script is to analyze the energy density of the categories of zooplankton used in out SEM food web model.
library(tidyverse)
library(ggplot2)
source("./breakout_zoop_energy_groups.R")
All the files are generated with a excel macro so their file names all start the same way. I use map to read them all into a list
file_name_start = "./zoop_energy_data/zoo_energy_"
catagories = na.omit(read.csv(paste0(file_name_start, "catagories.csv"))) %>%
mutate_if(is.character, str_trim)
file_names = as.list(c( "dry", "wet", "digestibility", "ratio"))
data = file_names %>%
map(~na.omit(read.csv(paste0(file_name_start, .x, ".csv"))) %>%
mutate_if(is.character, str_trim)) %>%
set_names(file_names)
If there are groups that are larger then our categories, make entries for each sub category containing that data the rows for each one
higher_groups = catagories %>%
na.omit()
# break out the catagories
data_brokenout = data %>%
map(~break_out_groups(.x, higher_groups, "group"))
data_ratio = data_brokenout$ratio %>%
group_by(group) %>%
summarize(dry_to_wet_ratio = mean(dry_to_wet_ratio, na.rm = TRUE))
# convert wet to dry if data available
converted_dry = data_ratio %>%
na.omit() %>%
right_join(data_brokenout$wet, by = "group") %>%
mutate(energy_density_j_per_g_dry_mass =
energy_density_j_per_g_wet_mass/dry_to_wet_ratio) %>%
select(-energy_density_j_per_g_wet_mass, -dry_to_wet_ratio) %>%
bind_rows(data_brokenout$dry)
converted_wet = data_ratio %>%
na.omit() %>%
right_join(data_brokenout$dry, by = "group") %>%
mutate(energy_density_j_per_g_wet_mass =
energy_density_j_per_g_dry_mass*dry_to_wet_ratio) %>%
select(-energy_density_j_per_g_dry_mass, -dry_to_wet_ratio) %>%
bind_rows(data_brokenout$wet)
Convert all to wet weight and print a table
dry_summary = converted_wet %>%
group_by(group) %>%
summarise(mean_energy_density = mean(energy_density_j_per_g_wet_mass),
energy_density_sd = sd(energy_density_j_per_g_wet_mass))
print(dry_summary)
## # A tibble: 7 x 3
## group mean_energy_density energy_density_sd
## <chr> <dbl> <dbl>
## 1 Amphipods 4096. 616.
## 2 Cladoceran 1975. 1087.
## 3 Cyclopoid 2885. 146.
## 4 Herbivorous Calanoid 3126. 375.
## 5 Mysids 3301. 482.
## 6 Predatory Calanoid 2947. 211.
## 7 Rotifers 1588. 76.4
Print a box plot of the data
density_plot = ggplot(data = converted_wet,
aes(x = fct_reorder(group, energy_density_j_per_g_wet_mass),
y = energy_density_j_per_g_wet_mass)) +
theme_classic(base_size = 20) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(x = "Group",
y = "Energy Density (j/g wet mass)" ) +
geom_boxplot() +
geom_dotplot(binaxis='y', stackdir='center', dotsize=2, binwidth = 30, alpha = 0.2)+
stat_summary(fun=mean, geom="point", color = "red", shape=23, size=4, stroke = 2)
print(density_plot)