09_exp_electivity
2023-08-07
This code calculates electivity indices (Chesson’s alpha) for the experiments. Since not all taxa were present in all experiments, the top nine taxa (in abudnance) from each experiment were chosen in order for alpha values to be comparable. For calculations, the zero-jelly tanks were treated as starting densities, and the highest jelly tanks were treated as ending densities.
Report is here: https://rpubs.com/HailaSchultz/exp_electivity
load packages
library(electivity)
library(reshape2)
library(dplyr)
library(ggplot2)
library(forcats)Set up data
Read database into R. If current database changes, just change the path
Database <- read.csv("/Users/hailaschultz/Dropbox/Other studies/Aurelia project/Data Analysis/data/current_data/Final_Aurelia_Database_Jan11_2023.csv")Add control tanks
First, read in the file that designates different samples to tanks
with jellies other, tanks with zero jellies sampled at the
beginning of the experiment remove, and tanks with zero
jellies sampled at the end of the experiment zero
Control_tanks <- read.csv("/Users/hailaschultz/Dropbox/Other studies/Aurelia project/Data Analysis/data/current_data/Control_Tanks.csv")Add control tank data to the database file
Database$Control_info <- Control_tanks$c1[match(Database$Sample.Code, Control_tanks$Sample.Code)]Subset data
Subset to trials with large jellies where the number of jellies in the tanks was manipulated
Exp_numb_trials<- subset(Database, Trial.Type=='Number')
Exp_numb_large<-subset(Exp_numb_trials,Jelly.Size=="Large")remove aug 22 2019 - the jellies we used for this experiment were younger/smaller than the other experiments, so I decided to omit it from analysis
Exp_sub<-subset(Exp_numb_large,Sample.Date!="08/22/2019")
Exp_sub<-subset(Exp_sub,Sample.Date!="08/13/2019")remove tanks sampled at the beginning of experiment - we decided not to use these tanks because they were different than the tanks with zero jellies sampled at the end of the experiment and didn’t add anything beneficial
Exp_sub<-subset(Exp_sub,Control_info!="remove")remove miscellaneous taxa
Exp_sub<-subset(Exp_sub,Genus.species!="Diatom-Centric")recode acartia as calanoida
Exp_sub <- Exp_sub %>%
mutate(Genus.species = fct_recode(Genus.species,
"CALANOIDA" = "ACARTIA"))Calculate indices for each experiment date
9/23/19
subset date
Sep23_19<-subset(Exp_sub,Sample.Date=="09/23/2019")group by tank
Sep23_19<-Sep23_19 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Sep23_19_ave<-Sep23_19 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Sep23_19_ave$Number.of.Jellies = as.factor(Sep23_19_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Sep23_19_ave<-subset(Sep23_19_ave,Number.of.Jellies!=c("3"))reshape data
Sep23_19_ave<-dcast(Sep23_19_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Sep23_19_ave<-Sep23_19_ave[complete.cases(Sep23_19_ave), ]rename columns
Sep23_19_ave<-Sep23_19_ave %>%
rename(
"p" = "0",
"remaining"="6")calculate amount consumed
Sep23_19_ave$r<-Sep23_19_ave$p-Sep23_19_ave$remainingkeep top nine taxa
Sep23_19_top <- top_n(Sep23_19_ave, 9, p)get selectivity indices
chesson<-chesson_alpha(Sep23_19_top$r,Sep23_19_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Sep23_19_alpha <- cbind(Sep23_19_top, chesson)08/29/19
subset date
Aug29_19<-subset(Exp_sub,Sample.Date=="08/29/2019")group by tank
Aug29_19<-Aug29_19 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Aug29_19_ave<-Aug29_19 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Aug29_19_ave$Number.of.Jellies = as.factor(Aug29_19_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Aug29_19_ave<-subset(Aug29_19_ave,Number.of.Jellies!=c("3"))reshape data
Aug29_19_ave<-dcast(Aug29_19_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Aug29_19_ave<-Aug29_19_ave[complete.cases(Aug29_19_ave), ]rename columns
Aug29_19_ave<-Aug29_19_ave %>%
rename(
"p" = "0",
"remaining"="6")calculate amount consumed
Aug29_19_ave$r<-Aug29_19_ave$p-Aug29_19_ave$remainingkeep top nine taxa
Aug29_19_top <- top_n(Aug29_19_ave, 9, p)get selectivity indices
chesson<-chesson_alpha(Aug29_19_top$r,Aug29_19_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Aug29_19_alpha <- cbind(Aug29_19_top, chesson)09/05/2019
subset date
Sep5_19<-subset(Exp_sub,Sample.Date=="09/05/2019")group by tank
Sep5_19<-Sep5_19 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Sep5_19_ave<-Sep5_19 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Sep5_19_ave$Number.of.Jellies = as.factor(Sep5_19_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Sep5_19_ave<-subset(Sep5_19_ave,Number.of.Jellies!=c("3"))reshape data
Sep5_19_ave<-dcast(Sep5_19_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Sep5_19_ave<-Sep5_19_ave[complete.cases(Sep5_19_ave), ]rename columns
Sep5_19_ave<-Sep5_19_ave %>%
rename(
"p" = "0",
"remaining"="6")calculate amount consumed
Sep5_19_ave$r<-Sep5_19_ave$p-Sep5_19_ave$remainingkeep top nine taxa
Sep5_19_top <- top_n(Sep5_19_ave, 9, p)replace negative number with 1
Sep5_19_top[Sep5_19_top < 0] <- 1 get selectivity indices
chesson<-chesson_alpha(Sep5_19_top$r,Sep5_19_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Sep5_19_alpha <- cbind(Sep5_19_top, chesson)09/12/2019
subset date
Sep12_19<-subset(Exp_sub,Sample.Date=="09/12/2019")group by tank
Sep12_19<-Sep12_19 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Sep12_19_ave<-Sep12_19 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Sep12_19_ave$Number.of.Jellies = as.factor(Sep12_19_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Sep12_19_ave<-subset(Sep12_19_ave,Number.of.Jellies!=c("3"))reshape data
Sep12_19_ave<-dcast(Sep12_19_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Sep12_19_ave<-Sep12_19_ave[complete.cases(Sep12_19_ave), ]rename columns
Sep12_19_ave<-Sep12_19_ave %>%
rename(
"p" = "0",
"remaining"="6")calculate amount consumed
Sep12_19_ave$r<-Sep12_19_ave$p-Sep12_19_ave$remainingkeep top nine taxa
Sep12_19_top <- top_n(Sep12_19_ave, 9, p)replace negative number with 1
Sep12_19_top[Sep12_19_top < 0] <- 1 get selectivity indices
chesson<-chesson_alpha(Sep12_19_top$r,Sep12_19_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Sep12_19_alpha <- cbind(Sep12_19_top, chesson)08/24/2020
subset date
Aug24_20<-subset(Exp_sub,Sample.Date=="08/24/2020")group by tank
Aug24_20<-Aug24_20 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Aug24_20_ave<-Aug24_20 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Aug24_20_ave$Number.of.Jellies = as.factor(Aug24_20_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Aug24_20_ave<-subset(Aug24_20_ave,Number.of.Jellies!=c("2"))reshape data
Aug24_20_ave<-dcast(Aug24_20_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Aug24_20_ave<-Aug24_20_ave[complete.cases(Aug24_20_ave), ]rename columns
Aug24_20_ave<-Aug24_20_ave %>%
rename(
"p" = "0",
"remaining"="4")calculate amount consumed
Aug24_20_ave$r<-Aug24_20_ave$p-Aug24_20_ave$remainingkeep top nine taxa
Aug24_20_top <- top_n(Aug24_20_ave, 9, p)replace negative number with 1
Aug24_20_top[Aug24_20_top < 0] <- 1 get selectivity indices
chesson<-chesson_alpha(Aug24_20_top$r,Aug24_20_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Aug24_20_alpha <- cbind(Aug24_20_top, chesson)08/17/2020
subset date
Aug17_20<-subset(Exp_sub,Sample.Date=="08/17/2020")group by tank
Aug17_20<-Aug17_20 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Aug17_20_ave<-Aug17_20 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Aug17_20_ave$Number.of.Jellies = as.factor(Aug17_20_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Aug17_20_ave<-subset(Aug17_20_ave,Number.of.Jellies!=c("2"))reshape data
Aug17_20_ave<-dcast(Aug17_20_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Aug17_20_ave<-Aug17_20_ave[complete.cases(Aug17_20_ave), ]rename columns
Aug17_20_ave<-Aug17_20_ave %>%
rename(
"p" = "0",
"remaining"="4")calculate amount consumed
Aug17_20_ave$r<-Aug17_20_ave$p-Aug17_20_ave$remainingkeep top nine taxa
Aug17_20_top <- top_n(Aug17_20_ave, 9, p)get selectivity indices
chesson<-chesson_alpha(Aug17_20_top$r,Aug17_20_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Aug17_20_alpha <- cbind(Aug17_20_top, chesson)08/14/2020
subset date
Aug14_20<-subset(Exp_sub,Sample.Date=="08/14/2020")group by tank
Aug14_20<-Aug14_20 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Aug14_20_ave<-Aug14_20 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Aug14_20_ave$Number.of.Jellies = as.factor(Aug14_20_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Aug14_20_ave<-subset(Aug14_20_ave,Number.of.Jellies!=c("2"))reshape data
Aug14_20_ave<-dcast(Aug14_20_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Aug14_20_ave<-Aug14_20_ave[complete.cases(Aug14_20_ave), ]rename columns
Aug14_20_ave<-Aug14_20_ave %>%
rename(
"p" = "0",
"remaining"="4")calculate amount consumed
Aug14_20_ave$r<-Aug14_20_ave$p-Aug14_20_ave$remainingkeep top nine taxa
Aug14_20_top <- top_n(Aug14_20_ave, 9, p)get selectivity indices
chesson<-chesson_alpha(Aug14_20_top$r,Aug14_20_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Aug14_20_alpha <- cbind(Aug14_20_top, chesson)08/13/2020
subset date
Aug13_20<-subset(Exp_sub,Sample.Date=="08/13/2020")group by tank
Aug13_20<-Aug13_20 %>%
group_by(Genus.species,Number.of.Jellies,Station,Sample.Date,Sample.Year) %>%
summarise(
sum = sum(Density....m3.))average by number of jellies
Aug13_20_ave<-Aug13_20 %>%
group_by(Genus.species,Number.of.Jellies,Sample.Date,Sample.Year) %>%
summarise(
mean = mean(sum))
Aug13_20_ave$Number.of.Jellies = as.factor(Aug13_20_ave$Number.of.Jellies)subset to zero jelly and highest jelly tanks
Aug13_20_ave<-subset(Aug13_20_ave,Number.of.Jellies!=c("1"))
Aug13_20_ave<-subset(Aug13_20_ave,Number.of.Jellies!=c("2"))reshape data
Aug13_20_ave<-dcast(Aug13_20_ave,Genus.species+Sample.Date+Sample.Year ~ Number.of.Jellies, value.var = "mean")remove incomplete cases
Aug13_20_ave<-Aug13_20_ave[complete.cases(Aug13_20_ave), ]rename columns
Aug13_20_ave<-Aug13_20_ave %>%
rename(
"p" = "0",
"remaining"="4")calculate amount consumed
Aug13_20_ave$r<-Aug13_20_ave$p-Aug13_20_ave$remainingkeep top nine taxa
Aug13_20_top <- top_n(Aug13_20_ave, 9, p)get selectivity indices
chesson<-chesson_alpha(Aug13_20_top$r,Aug13_20_top$p, na.rm = TRUE)
chesson<-as.data.frame(chesson)add dataframes together
Aug13_20_alpha <- cbind(Aug13_20_top, chesson)Add all dates together
AllDates <- rbind(Sep23_19_alpha, Aug29_19_alpha, Sep5_19_alpha, Sep12_19_alpha, Aug24_20_alpha, Aug17_20_alpha, Aug14_20_alpha, Aug13_20_alpha)Make Figures
remove taxa with only one experiment
Electivity<-subset(AllDates,Genus.species!="HYPERIIDEA")
Electivity<-subset(Electivity,Genus.species!="LITTORINA")Recode taxa names
unique(Electivity$Genus.species)## [1] CALANOIDA BIVALVIA
## [3] BRYOZOA DITRICHOCORYCAEUS ANGLICUS
## [5] GASTROPODA HARPACTICOIDA
## [7] OITHONA SHRIMP
## [9] BARNACLES COPEPODA
## [11] CRABS Chaet/Euphaus Egg
## [13] POLYCHAETA
## 31 Levels: CALANOIDA AETIDEUS BARNACLES BIVALVIA BRYOZOA ... SHRIMPlevels(Electivity$Genus.species) <- list("Medium Calanoids" = "CALANOIDA",
"Larval Bivalves" = "BIVALVIA",
"Bryozoans"="BRYOZOA",
"Ditrichocorycaeus anglicus"="DITRICHOCORYCAEUS ANGLICUS",
"Larval Gastropods"="GASTROPODA",
"Barnacle Nauplii"="BARNACLES",
"Copepod Nauplii"="COPEPODA",
"Larval Shrimp"="SHRIMP",
"Harpacticoids/Cyclopoids"="HARPACTICOIDA",
"Oithonids"="OITHONA",
"Larval Crabs"="CRABS",
"Chaetognath or Euphausiid Eggs"="Chaet/Euphaus Egg",
"Larval Polychaetes"="POLYCHAETA")reorder factors
Electivity$Genus.species <- factor(Electivity$Genus.species,
levels = c('Medium Calanoids','Ditrichocorycaeus anglicus',
'Oithonids',
'Harpacticoids/Cyclopoids','Copepod Nauplii',
'Larval Shrimp','Larval Crabs',
'Barnacle Nauplii',
'Chaetognath or Euphausiid Eggs',
'Larval Bivalves','Larval Gastropods','Bryozoans',
'Larval Polychaetes'))electivity_plot<-ggplot(Electivity, aes(x = Genus.species, y = chesson)) + xlab("Taxa")+ylab(expression(alpha))+
geom_boxplot()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=7,colour="black"),
axis.text.y = element_text(size=7,colour="black"),
axis.title.x = element_text(size=7,colour="black"),
axis.title.y=element_text(size=7,colour="black"))+
geom_hline(yintercept=0.11,color="red",lty="dashed")
electivity_plot
save figure
setwd("/Users/hailaschultz/Dropbox/Other studies/Aurelia project/Data Analysis/figures/Final_Figures")
ggsave(filename = "Fig3.png", plot = electivity_plot, width = 80, height = 80, units="mm", device='png', dpi=600)
ggsave(filename = "Fig3.tif", plot = electivity_plot, width = 80, height = 80, units="mm", device='tiff', dpi=600)