#Rakowski and Leibold 2021
#Analysis of zooplankton data from pleid field mesocosm experiment,
 #beginning with conversion from length & count to biomass

#Load packages
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
library(lme4)
library(patchwork)

#Set working directory
setwd("~/Desktop")

#Import cleaned zoop data
pltaxamass0s <- read_csv("~/Downloads/processed zoop data pleid expt.csv")
View(pltaxamass0s)

#Ensure "tank" (mesocosm number) is a factor 
pltaxamass0s$tank <- as.factor(pltaxamass0s$tank)

#Reorder treatments for figures
pltaxamass0s$treatment <- factor(pltaxamass0s$treatment, 
                  levels=c("no zooplankton","zooplankton","40 pleids","80 pleids"))

#Find mean & median mass of each taxon overall & sort by dominance
sum.pl.zoop <- pltaxamass0s %>%
  dplyr::group_by(taxon) %>%
  dplyr::summarise(mn.mass=mean(mass, na.rm=TRUE), md.mass=median(mass, na.rm=TRUE)) %>%
  arrange(-md.mass, -mn.mass)
print(sum.pl.zoop, n=24)

#Separate copepods from non-copepod zoop
#Non-copepods
pl.NC0s <- subset(pltaxamass0s, ! taxon %in% c("Arctodiaptomus_dorsalis",
  "Microcyclops_varicans", "nauplius"))
unique(pl.NC0s$taxon)
#Copepods
pl.cop0s <- subset(pltaxamass0s, taxon %in% c("Arctodiaptomus_dorsalis",
  "Microcyclops_varicans", "nauplius"))
unique(pl.cop0s$taxon)

########### Analyze non-copepod zooplankton ############

#Sum non-copepod mass by sample
pl.sumNC0s <- pl.NC0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(mass=sum(mass, na.rm=TRUE))
#View(pl.sumNC0s)

#Calculate mean, SD, & SE over the two dates
sumNC <- pl.sumNC0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(NC.mn=mean(mass, na.rm=TRUE), NC.sd=sd(mass, na.rm=TRUE), NC.se=NC.sd/sqrt(2))
sumNC

#Prep data for zoop mass figure faceted by zoop group (Fig. 1)
sumNC$group <- "non-copepods"
sumNC2 <- sumNC %>%
  dplyr::rename(mean=NC.mn, SE=NC.se) %>%
  select(!NC.sd)
sumNC2

#Calculate mean, SD, & SE over all samples by treatment  ***NEW***
trt.sumNC <- pl.sumNC0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(NC.mn=mean(mass, na.rm=TRUE), NC.sd=sd(mass, na.rm=TRUE), NC.se=NC.sd/sqrt(10))
trt.sumNC

#Merge this treatment summary onto above 
trt.sumNC$group <- "non-copepods"
trt.sumNC$tank <- "all"
trt.sumNC2 <- trt.sumNC %>%
  dplyr::rename(mean=NC.mn, SE=NC.se) %>%
  select(!NC.sd)
trt.sumNC2
sumNC3 <- rbind(sumNC2, trt.sumNC2)
View(sumNC3)

#Add 0.01 to each sample to eliminate 0s
pl.sumNC0s$massplus.01 <- pl.sumNC0s$mass + .01

#Plot non-copepod zoop mass by treatment & date
ggplot(pl.sumNC0s, aes(x=treatment, y=massplus.01)) + geom_jitter(height=0, width=0.15) +
  scale_y_log10(breaks=c(1e-01,1e+00,1e+01,1e+02,1e+03), labels=c("0.1","1","10","100","1000")) +
  annotation_logticks(sides="l") + facet_wrap(~date)

#Pleids didn't appear to affect non-copepod zoop mass. Test this statistically:

#Model w/ Neoplea
NC.GLMM <- glmer(massplus.01 ~ zoop + Neo + (1|tank), family=Gamma(link="log"), 
                   data=pl.sumNC0s)
summary(NC.GLMM) 
#              Estimate Std. Error t value Pr(>|z|)   
#(Intercept) -1.1836592  0.5586991  -2.119  0.03412 * 
#zoopy        2.4857501  0.7600263   3.271  0.00107 **
#Neo          0.0007758  0.0099557   0.078  0.93788   

#Model w/out Neoplea
NC.GLMM.noNeo <- glmer(massplus.01 ~ zoop + (1|tank), family=Gamma(link="log"),
                   data=pl.sumNC0s)
summary(NC.GLMM.noNeo)
#            Estimate Std. Error t value Pr(>|z|)    
#(Intercept)  -1.1836     0.5594  -2.116   0.0344 *  
#zoopy         2.5186     0.6457   3.900  9.6e-05 ***

#Null model
NC.GLMM.null <- glmer(massplus.01 ~ (1|tank), family=Gamma(link="log"), data=pl.sumNC0s)

#Compare models
anova(NC.GLMM.null, NC.GLMM.noNeo, NC.GLMM)
#                npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
#otherGLMM.null     3 191.62 196.69 -92.810   185.62                          
#otherGLMM.noNeo    4 182.48 189.23 -87.238   174.48 11.1423  1  0.0008438 ***
#otherGLMM          5 184.47 192.91 -87.235   174.47  0.0075  1  0.9310823 

 #Adding zooplankton significantly improves model, but not adding Neoplea


############### Copepods ################

#Sum copepod mass by sample
pl.sumcop0s <- pl.cop0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(mass=sum(mass, na.rm=TRUE))
#View(pl.sumcop0s)

#Calculate mean, SD, & SE over the two dates
sumcops <- pl.sumcop0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(copMN=mean(mass, na.rm=TRUE), copSD=sd(mass, na.rm=TRUE), copSE=copSD/sqrt(2))
sumcops

#Prep data for zoop mass figure faceted by zoop group (Fig. 1)
sumcops$group <- "copepods"
sumcops2 <- sumcops %>%
  dplyr::rename(mean=copMN, SE=copSE) %>%
  dplyr::select(!copSD)
sumcops2

#Calculate mean, SD, & SE over all samples by treatment  ***NEW*** 
trt.sumcops <- pl.sumcop0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(copMN=mean(mass, na.rm=TRUE), copSD=sd(mass, na.rm=TRUE), copSE=copSD/sqrt(10))
trt.sumcops

#Merge this treatment summary onto above 
trt.sumcops$group <- "copepods"
trt.sumcops$tank <- "all"
trt.sumcops2 <- trt.sumcops %>%
  dplyr::rename(mean=copMN, SE=copSE) %>%
  select(!copSD)
trt.sumcops2
sumcops3 <- rbind(sumcops2, trt.sumcops2)
View(sumcops3)

#Add 0.1 to each sample to eliminate 0s
pl.sumcop0s$massplus.1 <- pl.sumcop0s$mass + .1

#Plot copepod mass by treatment & date
ggplot(pl.sumcop0s, aes(x=treatment, y=massplus.1)) + geom_jitter(height=0, width=0.15) +
  scale_y_log10(breaks=c(1e-01,1e+00,1e+01,1e+02,1e+03,1e+04), labels=c("0.1","1","10","100","1000","10000")) +
  annotation_logticks(sides="l") + facet_wrap(~date)

#Pleids didn't appear to affect copepod mass. Test this statistically:

#Full model
copGLMM <- glmer(massplus.1 ~ zoop + Neo + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumcop0s)
summary(copGLMM) 
#            Estimate Std. Error t value Pr(>|z|)
#(Intercept)  0.38192    0.93010   0.411    0.681
#zoopy        1.84080    1.26247   1.458    0.145
#Neo          0.00488    0.01653   0.295    0.768 

#W/out Neoplea
copGLMM.noNeo <- glmer(massplus.1 ~ zoop + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumcop0s)
summary(copGLMM.noNeo)
#            Estimate Std. Error t value Pr(>|z|)  
#(Intercept)   0.3820     0.9319    0.41   0.6818  
#zoopy         2.0361     1.0772    1.89   0.0587 .

#Null model
copGLMM.null <- glmer(massplus.1 ~ (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumcop0s)
#Compare models
anova(copGLMM.null, copGLMM.noNeo, copGLMM)
#              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#copGLMM.null     3 86.450 91.516 -40.225   80.450                       
#copGLMM.noNeo    4 85.195 91.951 -38.598   77.195 3.2546  1    0.07122 .
#copGLMM          5 87.108 95.553 -38.554   77.108 0.0870  1    0.76802  


######### Cladocerans & ostracods (mostly Scaph) ###########

#Subset these taxa
pl.CladOst0s <- subset(pltaxamass0s, taxon %in% c("Scapholeberis_embryo", "cladoceran_embryo",
                                               "Scapholeberis_kingi", "ostracod"))
unique(pl.CladOst0s$taxon)

#Sum mass by sample
pl.sumCladOst0s <- pl.CladOst0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(mass=sum(mass, na.rm=TRUE))
View(pl.sumCladOst0s)

#Calculate mean, SD, & SE over the two dates
sumCladOst <- pl.sumCladOst0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(CladOstMN=mean(mass, na.rm=TRUE), CladOstSD=sd(mass, na.rm=TRUE),
            CladOstSE=CladOstSD/sqrt(2))
sumCladOst

#Prep data for zoop mass figure faceted by group (Fig. 1)
sumCladOst$group <- "cladocerans + ostracods"
sumCladOst2 <- sumCladOst %>%
  dplyr::rename(mean=CladOstMN, SE=CladOstSE) %>%
  dplyr::select(!CladOstSD)
sumCladOst2

#Calculate mean, SD, & SE over all samples by treatment  ***NEW*** 
trt.sumCladOst <- pl.sumCladOst0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(CladOstMN=mean(mass, na.rm=TRUE), CladOstSD=sd(mass, na.rm=TRUE), CladOstSE=CladOstSD/sqrt(10))
trt.sumCladOst

#Merge this treatment summary onto above 
trt.sumCladOst$group <- "cladocerans + ostracods"
trt.sumCladOst$tank <- "all"
trt.sumCladOst2 <- trt.sumCladOst %>%
  dplyr::rename(mean=CladOstMN, SE=CladOstSE) %>%
  select(!CladOstSD)
trt.sumCladOst2
sumCladOst3 <- rbind(sumCladOst2, trt.sumCladOst2)
View(sumCladOst3)

#Add 0.1 to eliminate 0s
pl.sumCladOst0s$massplus.1 <- pl.sumCladOst0s$mass + .1

#Plot mass by treatment & date
ggplot(pl.sumCladOst0s, aes(x=treatment, y=massplus.1)) + geom_jitter(height=0, width=0.15) +
  scale_y_log10(breaks=c(1e-01,1e+00,1e+01,1e+02,1e+03,1e+04), labels=c("0.1","1","10","100","1000","10000")) +
  annotation_logticks(sides="l") + facet_wrap(~date)

#No apparent effect of pleids on Scaph. Test statistically:

#Full model
CladOstGLMM <- glmer(massplus.1 ~ zoop + Neo + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumCladOst0s)
summary(CladOstGLMM) 
#            Estimate Std. Error t value Pr(>|z|)  
#(Intercept) -1.34675    0.70153  -1.920   0.0549 .
#zoopy        1.66659    0.96035   1.735   0.0827 .
#Neo          0.00830    0.01263   0.657   0.5109  

#W/out Neoplea
CladOstGLMM.noNeo <- glmer(massplus.1 ~ zoop + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumCladOst0s)
summary(CladOstGLMM.noNeo)
#            Estimate Std. Error t value Pr(>|z|)  
#(Intercept)  -1.3464     0.7044  -1.912   0.0559 .
#zoopy         2.0032     0.8162   2.454   0.0141 *

#Null model
CladOstGLMM.null <- glmer(massplus.1 ~ (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumCladOst0s)
#Compare models
anova(CladOstGLMM.null, CladOstGLMM.noNeo, CladOstGLMM)
#                  npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#CladOstGLMM.null     3 87.668 92.734 -40.834   81.668                       
#CladOstGLMM.noNeo    4 84.539 91.294 -38.269   76.539 5.1293  1    0.02353 *
#CladOstGLMM          5 86.106 94.550 -38.053   76.106 0.4325  1    0.51077  


########## Spirostomum ###########

#Subset
pl.Spiro0s <- subset(pltaxamass0s, taxon %in% c("Spirostomum_sp_1", "Spirostomum_sp_2"))
unique(pl.Spiro0s$taxon)

#Sum Spiro mass by sample
pl.sumSpiro0s <- pl.Spiro0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(mass=sum(mass, na.rm=TRUE))
View(pl.sumSpiro0s)

#Calculate mean, SD, & SE over the two dates
sumSpiro <- pl.sumSpiro0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(SpiroMN=mean(mass, na.rm=TRUE), SpiroSD=sd(mass, na.rm=TRUE),
            SpiroSE=SpiroSD/sqrt(2))
sumSpiro

#Prep data for zoop mass figure faceted by group
sumSpiro$group <- "Spirostomum"
sumSpiro2 <- sumSpiro %>%
  dplyr::rename(mean=SpiroMN, SE=SpiroSE) %>%
  dplyr::select(!SpiroSD)
sumSpiro2

#Calculate mean, SD, & SE over all samples by treatment  ***NEW*** 
trt.sumSpiro <- pl.sumSpiro0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(SpiroMN=mean(mass, na.rm=TRUE), SpiroSD=sd(mass, na.rm=TRUE), SpiroSE=SpiroSD/sqrt(10))
trt.sumSpiro

#Merge this treatment summary onto above 
trt.sumSpiro$group <- "Spirostomum"
trt.sumSpiro$tank <- "all"
trt.sumSpiro2 <- trt.sumSpiro %>%
  dplyr::rename(mean=SpiroMN, SE=SpiroSE) %>%
  select(!SpiroSD)
trt.sumSpiro2
sumSpiro3 <- rbind(sumSpiro2, trt.sumSpiro2)
View(sumSpiro3)

#Add 0.1 to eliminate 0s
pl.sumSpiro0s$massplus.1 <- pl.sumSpiro0s$mass + .1

#Analyze: full model
SpiroGLMM <- glmer(massplus.1 ~ zoop + Neo + (1|tank), family=Gamma(link="log"), nAGQ=10,
                  data=pl.sumSpiro0s)
summary(SpiroGLMM) 
#             Estimate Std. Error t value Pr(>|z|)   
#(Intercept) -1.257345   0.457544  -2.748    0.006 **
#zoopy        0.948931   0.632187   1.501    0.133   
#Neo          0.001766   0.008502   0.208    0.835 

#W/out Neoplea
SpiroGLMM.noNeo <- glmer(massplus.1 ~ zoop + (1|tank), family=Gamma(link="log"), nAGQ=10,
                  data=pl.sumSpiro0s)
summary(SpiroGLMM.noNeo)
#            Estimate Std. Error t value Pr(>|z|)   
#(Intercept)  -1.2569     0.4587  -2.740  0.00613 **
#zoopy         1.0182     0.5383   1.892  0.05854 .

#Null model
SpiroGLMM.null <- glmer(massplus.1 ~ (1|tank), family=Gamma(link="log"), nAGQ=10,
                  data=pl.sumSpiro0s)

#Compare models
anova(SpiroGLMM.null, SpiroGLMM.noNeo, SpiroGLMM)
#                npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#SpiroGLMM.null     3 58.082 63.148 -26.041   52.082                       
#SpiroGLMM.noNeo    4 57.017 63.773 -24.509   49.017 3.0644  1    0.08003 .
#SpiroGLMM          5 58.974 67.419 -24.487   48.974 0.0430  1    0.83567  


######## Rotifers ########

#Subset rotifers
pl.rot0s <- subset(pltaxamass0s, taxon %in% c("small_Euchlanis", "Monostyla_quadridentata",
             "Lecane_luna", "Monostyla_lunaris", "bdelloid", "Monostyla_copeis",
             "Monostyla_closterocerca", "ploimid", "Monostyla_bulla",
             "smaller_Monostyla_bulla", "small_ploimid", "Monostyla_cornuta",
             "Euchlanis", "Anuraeopsis", "Lecane_sp"))
unique(pl.rot0s$taxon)

#Sum mass by sample
pl.sumrot0s <- pl.rot0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(mass=sum(mass, na.rm=TRUE))
View(pl.sumrot0s)

#Calculate mean, SD, & SE over the two dates
sumrot <- pl.sumrot0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(rotMN=mean(mass, na.rm=TRUE), rotSD=sd(mass, na.rm=TRUE), rotSE=rotSD/sqrt(2))
sumrot

#Prep data for zoop mass figure faceted by group
sumrot$group <- "rotifers"
sumrot2 <- sumrot %>%
  dplyr::rename(mean=rotMN, SE=rotSE) %>%
  dplyr::select(!rotSD)
sumrot2

#Calculate mean, SD, & SE over all samples by treatment  ***NEW*** 
trt.sumrot <- pl.sumrot0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(rotMN=mean(mass, na.rm=TRUE), rotSD=sd(mass, na.rm=TRUE), rotSE=rotSD/sqrt(10))
trt.sumrot

#Merge this treatment summary onto above 
trt.sumrot$group <- "rotifers"
trt.sumrot$tank <- "all"
trt.sumrot2 <- trt.sumrot %>%
  dplyr::rename(mean=rotMN, SE=rotSE) %>%
  select(!rotSD)
trt.sumrot2
sumrot3 <- rbind(sumrot2, trt.sumrot2)
View(sumrot3)

#Calculate mean rotifer mass by treatment
sumrot.trt <- sumrot2 %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(mn.mass=mean(mean), md.mass=median(mean))
sumrot.trt
# treatment         mn.mass md.mass
# no zooplankton    0.0713  0.0259
# zooplankton       0.936   0.533 
# zoop + 40 pleids  0.814   0.135 
# zoop + 80 pleids  0.494   0.527
100*(0.814-0.936)/0.936 #40 Neoplea reduced mean rotifer mass by 13%...
100*(0.135-0.533)/0.533 #...and median rotifer mass by 75%.
100*(0.494-0.936)/0.936 #80 Neoplea reduced mean rotifer mass by 47%...
100*(0.527-0.533)/0.533 #...but median rotifer mass only by 1%.

#Add 0.001 to eliminate 0s
pl.sumrot0s$massplus.001 <- pl.sumrot0s$mass + .001

#Fit & compare models

#Full model
rotGLMM <- glmer(massplus.001 ~ zoop + Neo + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumrot0s)
summary(rotGLMM)
#             Estimate Std. Error t value Pr(>|z|)    
#(Intercept) -3.155534   0.501809  -6.288 3.21e-10 ***
#zoopy        2.840588   0.671734   4.229 2.35e-05 ***
#Neo         -0.006806   0.008541  -0.797    0.425

#W/out Neoplea
rotGLMM.noNeo <- glmer(massplus.001 ~ zoop + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumrot0s)
summary(rotGLMM.noNeo)
#            Estimate Std. Error t value Pr(>|z|)    
#(Intercept)  -3.1591     0.5098  -6.197 5.77e-10 ***
#zoopy         2.5722     0.5879   4.376 1.21e-05 ***

#Null model
rotGLMM.null <- glmer(massplus.001 ~ (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=pl.sumrot0s)

#Compare models
anova(rotGLMM.null, rotGLMM.noNeo, rotGLMM)
#              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
#rotGLMM.null     3 91.087 96.154 -42.544   85.087                         
#rotGLMM.noNeo    4 79.911 86.666 -35.955   71.911 13.177  1  0.0002835 ***
#rotGLMM          5 81.287 89.731 -35.643   71.287  0.624  1  0.4295669   


#Merge zoop group summaries

require(purrr)
pl.sumZ <- list(sumNC, sumcops, sumCladOst, sumSpiro, sumrot) %>% 
 reduce(full_join, by=c("tank", "treatment"))
View(pl.sumZ)
#add total zoop mass variable
pl.sumZ$totZ <- pl.sumZ$NC.mn + pl.sumZ$copMN

#Stacked bar plot of zoop composition by treatment (Fig. S1)

#Calculate proportions each zoop group comprises of total zoop, for each tank (averaged over days)
pl.sumZ$copprop <- pl.sumZ$copMN / pl.sumZ$totZ
pl.sumZ$CladOstprop <- pl.sumZ$CladOstMN / pl.sumZ$totZ
pl.sumZ$Spiroprop <- pl.sumZ$SpiroMN / pl.sumZ$totZ
pl.sumZ$rotprop <- pl.sumZ$rotMN / pl.sumZ$totZ

#Average these props by treatment
props <- pl.sumZ %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(mn.copprop=mean(copprop), mn.CladOstprop=mean(CladOstprop),
            mn.Spiroprop=mean(Spiroprop), mn.rotprop=mean(rotprop))
#props$total <- props$mn.copprop + props$mn.CladOstprop + props$mn.Spiroprop + props$mn.rotprop
props

#Put into long form
require(reshape2)
propslong <- melt(props, id.vars="treatment") %>%
  dplyr::rename(group=variable, prop=value)
propslong

#Plot these mean props as stacked percentage bar chart
figS1B <- ggplot(propslong, aes(fill=group, y=prop, x=treatment)) +
  geom_bar(position="fill", stat="identity") + theme_bw(10) +
  labs(y="Proportion of zooplankton mass", x="Treatment") +
  scale_x_discrete(breaks=c("no zooplankton", "zooplankton", "40 pleids", "80 pleids"),
     labels=c("no zoop", expression(paste("no ", italic(Neoplea))),
      expression(paste("40 ", italic(Neoplea))), expression(paste("80 ", italic(Neoplea))))) +
  scale_fill_manual(name="Zooplankton group", values=c("#D55E00", "#56B4E9", "#CC79A7", "#009E73"),
    breaks=c("mn.copprop", "mn.CladOstprop", "mn.Spiroprop", "mn.rotprop"),
    labels=c("copepods", "other Crustacea", expression(italic("Spirostomum")), "rotifers")) +
  theme(legend.position="none", #remove legend
    panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
    axis.title.x=element_text(vjust=.2),
    axis.text.x=element_text(color="black"), axis.text.y=element_text(color="black"))

#Stacked bar chart showing relative average masses rather than proportions

#Average zoop group masses by treatment
aves <- pl.sumZ %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(mn.cop=mean(copMN), mn.CladOst=mean(CladOstMN),
            mn.Spiro=mean(SpiroMN), mn.rot=mean(rotMN))
aves

#Long form
aveslong <- melt(aves, id.vars="treatment") %>%
  dplyr::rename(group=variable, mean.mass=value)
aveslong

#Bar chart
figS1A <- ggplot(aveslong, aes(fill=group, y=mean.mass, x=treatment)) +
  geom_bar(position="stack", stat="identity") + theme_bw(10) +
  labs(y=expression(paste("Mean mass (", mu, "g/L)")), x="Treatment") +
  scale_x_discrete(breaks=c("no zooplankton", "zooplankton", "40 pleids", "80 pleids"),
    labels=c("no zoop", expression(paste("no ", italic(Neoplea))),
      expression(paste("40 ", italic(Neoplea))), expression(paste("80 ", italic(Neoplea))))) +
  scale_fill_manual(values=c("#D55E00", "#56B4E9", "#CC79A7", "#009E73"),
    breaks=c("mn.cop", "mn.CladOst", "mn.Spiro", "mn.rot"), name="Zooplankton group",
    labels=c("copepods", "other Crustacea", expression(italic("Spirostomum")), "rotifers")) +
  theme( #legend.title=element_text(size=11), legend.text.align=0, 
        legend.position="none", #remove legend
        panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text.x=element_text(color="black"), axis.text.y=element_text(color="black"))

#With legend
pdf("figure S1 zoop composition legend.pdf", width=5, height=3.5)
ggplot(aveslong, aes(fill=group, y=mean.mass, x=treatment)) +
  geom_bar(position="stack", stat="identity") + theme_bw(10) +
  labs(y=expression(paste("Mean mass (", mu, "g/L)")), x="Treatment") +
  scale_x_discrete(breaks=c("no zooplankton", "zooplankton", "40 pleids", "80 pleids"),
    labels=c("no zoop", expression(paste("no ", italic(Neoplea))),
      expression(paste("40 ", italic(Neoplea))), expression(paste("80 ", italic(Neoplea))))) +
  scale_fill_manual(values=c("#D55E00", "#56B4E9", "#CC79A7", "#009E73"),
    breaks=c("mn.cop", "mn.CladOst", "mn.Spiro", "mn.rot"), name="Zooplankton group",
    labels=c("copepods", "cladocerans + ostracods", expression(italic("Spirostomum")), "rotifers")) +
  theme(legend.text.align=0, 
        panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text.x=element_text(color="black"), axis.text.y=element_text(color="black"))
dev.off()

#Put fig panels together
pdf("figure S1 zoop composition bar charts.pdf", width=7.5, height=3.5)
figS1A + figS1B + plot_layout(ncol=2)
dev.off()


#Faceted figure showing biomass of different zoop groups by treatment (Fig. 1)

#First, merge summarized zoop group data
sumZlong <- rbind(sumcops3, sumNC3, sumCladOst3, sumSpiro3, sumrot3)
#View(sumZlong)
#Calculate high & low values of error bars & add a little bit to all values to eliminate 0s
sumZlong <- mutate(sumZlong, hi=mean+SE+0.001, lo=mean-SE+0.001, meanplus.001=mean+0.001)
View(sumZlong)
#Reorder zoop groups for figure
sumZlong$group <- factor(sumZlong$group, levels=c("copepods", "non-copepods",
                                    "cladocerans + ostracods", "Spirostomum", "rotifers"))
#Now, the plot
pdf("figure 1 faceted zoop group mass v2.pdf", width=7.5, height=5)
ggplot(subset(sumZlong, tank!="all"), aes(x=treatment, y=meanplus.001)) + theme_bw(10) +
  geom_jitter(width=.2, color="grey") +
  geom_point(data=subset(sumZlong, tank=="all"), aes(x=treatment, y=meanplus.001)) +
  geom_errorbar(data=subset(sumZlong, tank=="all"), aes(x=treatment, ymin=lo, ymax=hi), width=0.4) +
  facet_wrap(~group) + 
  scale_y_log10(breaks=c(1e-03,1e-02,1e-01,1e+00,1e+01,1e+02,1e+03),
        labels=c("0.001","0.01","0.1","1","10","100","1000")) + annotation_logticks(sides="l") +
  labs(y=expression(paste("Mean mass (", mu, "g/L)")), x="Treatment") +
  scale_x_discrete(breaks=c("no zooplankton", "zooplankton", "40 pleids", "80 pleids"),
    labels=c("no zoop", expression(paste("no ", italic(Neo.))),
      expression(paste("40 ", italic(Neo.))), expression(paste("80 ", italic(Neo.))))) +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text.x=element_text(color="black"), axis.text.y=element_text(color="black"))
dev.off()