#Rakowski and Leibold 2021 - pleid experiment: analysis of phytoplankton data

#Load packages
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
library(nlme)

#Load & view pleid experiment phyto counts
phyto <- read_csv("~/Downloads/processed phyto data pleid expt.csv")
View(phyto)

#Make sure tank is a factor (not continuous/numeric)
phyto$tank <- as.factor(phyto$tank)

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

#Reorder phyto taxon levels for figure 3
sum.Augphyto$taxon <- factor(sum.Augphyto$taxon, levels=c("big_ovoids", "ovoids", "Oocystis",
                        "pennate diatom", "picoplankton", "Chlorella", "Selenastrum"))

#Bar chart showing phyto composition across treatments (Fig. S2)

#no legend
figS2 <- ggplot(sum.Augphyto, aes(fill=taxon, y=mn.biovol, x=treatment)) +
  geom_bar(position="stack", stat="identity") + theme_bw(10) +
  labs(y=expression(paste("Mean biovolume (",mu,m^3,"/nL)")), x="Treatment") +
  scale_y_continuous(breaks=c(0, 1000, 2000, 3000), labels=c("0","1","2","3")) +
  scale_x_discrete(breaks=c("A", "DZ", "40P", "80P"),
    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", "#0072B2", "#E69F00", "#56B4E9", "#CC79A7", "#009E73", "#F0E442"),
    breaks=c("big_ovoids", "ovoids", "Oocystis", "pennate diatom", "picoplankton", "Chlorella",  
             "Selenastrum"), name="Phytoplankton\nmorphospecies",
    labels=c("large ovoid chlorophyte", "small ovoid chlorophyte", expression(italic("Oocystis")), "pennate diatom",
             "green picoplankton", expression(italic("Chlorella")), expression(italic("Selenastrum")))) +
  theme(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"))

#w/ legend
pdf("figure S2 phyto bar chart legend.pdf", width=5.25, height=3.75)
ggplot(sum.Augphyto, aes(fill=taxon, y=mn.biovol, x=treatment)) +
  geom_bar(position="stack", stat="identity") + theme_bw(10) +
  labs(y=expression(paste("Mean biovolume (",mu,m^3,"/nL)")), x="Treatment") +
  scale_y_continuous(breaks=c(0, 1000, 2000, 3000), labels=c("0","1","2","3")) +
  scale_x_discrete(breaks=c("A", "DZ", "40P", "80P"),
    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", "#0072B2", "#E69F00", "#56B4E9", "#CC79A7", "#009E73", "#F0E442"),
    breaks=c("big_ovoids", "ovoids", "Oocystis", "pennate diatom", "picoplankton", "Chlorella",  
             "Selenastrum"), name="Phytoplankton\nmorphospecies",
    labels=c("large ovoid chlorophyte", "small ovoid chlorophyte", expression(italic("Oocystis")), "pennate diatom",
             "green picoplankton", expression(italic("Chlorella")), expression(italic("Selenastrum")))) +
  theme(legend.title=element_text(size=11), 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()


#Sum total phyto mass by sample
sumphyto0s <- phyto %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(bv=sum(bv, na.rm=TRUE))
sumphyto0s

#Calculate mean, SD, & SE over the two dates
sumphyt <- sumphyto0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(mn.bv=mean(bv, na.rm=TRUE), sd.bv=sd(bv, na.rm=TRUE), se.bv=sd.bv/sqrt(2),
            hi=mn.bv+se.bv, lo=mn.bv-se.bv)
sumphyt$group <- "total"
sumphyt

#Calculate mean, SD, & SE over all samples by treatment  ***NEW***
trt.sumphyt <- sumphyto0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(mn.bv=mean(bv, na.rm=TRUE), sd.bv=sd(bv, na.rm=TRUE), se.bv=sd.bv/sqrt(10),
            hi=mn.bv+se.bv, lo=mn.bv-se.bv)
#Prep this treatment summary for the figure (new Fig. 2)
trt.sumphyt$group <- "total"
trt.sumphyt$tank <- "all"
trt.sumphyt

#dot plot 
ggplot(sumphyt, aes(x=treatment, y=mn.bv, group=tank)) + 
  geom_errorbar(aes(ymin=mn.bv-se.bv, ymax=mn.bv+se.bv), width=0.4,
                 position=position_dodge(width=0.4)) +
  geom_point(position=position_dodge(width=0.4)) + theme_bw(10) +
  scale_y_log10(breaks=c(1e+02,1e+03,1e+04,1e+05),
                labels=c("0.1","1","10","100")) + annotation_logticks(sides="l") +
  labs(y=expression(paste("Mean phytoplankton biovolume (",mu,m^3,"/nL)")), 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"))


########### Compare treatments, looking at all subsets #############

###### Total phyto #######

#Full model
all.gam <- glmer(bv ~ zoop + Neo + (1|tank), family=Gamma(link="log"), 
                 data=sumphyto0s)
summary(all.gam)

#W/out Neoplea
all.gam.noNeo <- glmer(bv ~ zoop + (1|tank), family=Gamma(link="log"), data=sumphyto0s)
summary(all.gam.noNeo)
#            Estimate Std. Error t value Pr(>|z|)    
#(Intercept)   9.0003     0.5939  15.155   <2e-16 ***
#zoopy        -0.7078     0.6855  -1.033    0.302 

#Null model
all.gam.null <- glmer(bv ~ (1|tank), family=Gamma(link="log"), data=sumphyto0s)

#Compare models
anova(all.gam.null, all.gam.noNeo, all.gam)
#              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
#all.gam.null     3 792.53 797.60 -393.27   786.53                     
#all.gam.noNeo    4 793.50 800.26 -392.75   785.50 1.0317  1     0.3098
#all.gam          5 793.90 802.34 -391.95   783.90 1.6016  1     0.2057

#Neither zoop nor Neoplea significantly affected total phyto biovolume


####### Larger phyto (large ovoids, Oocystis, & pennate diatoms) #######

#Subset the large phyto
bigphyto0s <- subset(phyto, taxon %in% c("big_ovoids", "Oocystis", "pennate diatom"))
bigphyto0s

#Sum large phyto mass by sample
sumbigphyto0s <- bigphyto0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(bv=sum(bv, na.rm=TRUE))
sumbigphyto0s$logplus10.bv <- log(sumbigphyto0s$bv + 10)
sumbigphyto0s$bvplus10 <- sumbigphyto0s$bv + 10
sumbigphyto0s

#Calculate mean, SD, & SE over the two dates
sumbigphyt <- sumbigphyto0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(mn.bv=mean(bv, na.rm=TRUE), sd.bv=sd(bv, na.rm=TRUE), se.bv=sd.bv/sqrt(2),
         mn.bvplus10=mn.bv+10, hi=mn.bv+se.bv, lo=mn.bv-se.bv, hi.plus10=hi+10, lo.plus10=lo+10)
sumbigphyt

#Prep data for phyto mass figure faceted by size grouping (Fig. 2)
sumbigphyt$group <- "large"
sumbig2 <- sumbigphyt %>%
  select(!c(mn.bv, hi, lo)) %>%
  dplyr::rename(mn.bv=mn.bvplus10, hi=hi.plus10, lo=lo.plus10) 
sumbig2

#Calculate mean, SD, & SE over all samples by treatment  ***NEW***
trt.sumbigphyt <- sumbigphyto0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(mn.bv=mean(bv, na.rm=TRUE), sd.bv=sd(bv, na.rm=TRUE), se.bv=sd.bv/sqrt(10),
            hi=mn.bv+se.bv, lo=mn.bv-se.bv)
#Prep this treatment summary for the figure (new Fig. 2)
trt.sumbigphyt$group <- "large"
trt.sumbigphyt$tank <- "all"
trt.sumbigphyt

#dot plot
ggplot(sumbigphyt, aes(x=treatment, y=mn.bvplus10, group=tank)) + 
  geom_errorbar(aes(ymin=lo.plus10, ymax=hi.plus10), width=0.4,
                 position=position_dodge(width=0.4)) +
  geom_point(position=position_dodge(width=0.4)) + theme_bw(10) +
  scale_y_log10(breaks=c(1e+01,1e+02,1e+03,1e+04,1e+05),
                labels=c("0.01","0.1","1","10","100")) + annotation_logticks(sides="l") +
  labs(y=expression(paste("Mean larger phytoplankton biovolume (",mu,m^3,"/nL)")), 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"))

#Analyze big phyto biovolume

#Full model
big.gam <- glmer(bvplus10 ~ zoop + Neo + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=sumbigphyto0s)
summary(big.gam)
#             Estimate Std. Error t value Pr(>|z|)    
#(Intercept)  8.153787   0.637273  12.795   <2e-16 ***
#zoopy       -1.901844   0.853543  -2.228   0.0259 *  
#Neo          0.009605   0.011003   0.873   0.3827  

#W/out Neoplea
big.gam.noNeo <- glmer(bvplus10 ~ zoop + (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=sumbigphyto0s)
summary(big.gam.noNeo)
#            Estimate Std. Error t value Pr(>|z|)    
#(Intercept)   8.1465     0.6502   12.53   <2e-16 ***
#zoopy        -1.5107     0.7478   -2.02   0.0434 *

#Null model
big.gam.null <- glmer(bvplus10 ~ (1|tank), family=Gamma(link="log"), nAGQ=10,
                 data=sumbigphyto0s)

#Compare models
anova(big.gam.null, big.gam.noNeo, big.gam)
#              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#big.gam.null     3 107.79 112.86 -50.896  101.793                       
#big.gam.noNeo    4 106.06 112.81 -49.029   98.059 3.7341  1    0.05331 .
#big.gam          5 107.32 115.76 -48.659   97.318 0.7409  1    0.38937  

#zoop marginally significantly reduced large phyto (p = 0.0533)


####### Small phyto #######

#Subset small phyto
smallphyto0s <- subset(phyto, taxon %in% c("ovoids", "Chlorella", "picoplankton",
                                                "Selenastrum"))

#Sum small phyto mass by sample
sumsmallphyto0s <- smallphyto0s %>%
  dplyr::group_by(tank, treatment, zoop, Neo, date) %>%
  dplyr::summarise(bv=sum(bv, na.rm=TRUE))
sumsmallphyto0s$logplus10.bv <- log(sumsmallphyto0s$bv)
sumsmallphyto0s

#Calculate mean, SD, & SE over the two dates
sumsmallphyt <- sumsmallphyto0s %>%
  dplyr::group_by(tank, treatment) %>%
  dplyr::summarise(mn.bv=mean(bv, na.rm=TRUE), sd.bv=sd(bv, na.rm=TRUE), se.bv=sd.bv/sqrt(2),
            hi=mn.bv+se.bv, lo=mn.bv-se.bv)
sumsmallphyt$group <- "small"
sumsmallphyt

#Calculate mean, SD, & SE over all samples by treatment  ***NEW***
trt.sumsmallphyt <- sumsmallphyto0s %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(mn.bv=mean(bv, na.rm=TRUE), sd.bv=sd(bv, na.rm=TRUE), se.bv=sd.bv/sqrt(10),
            hi=mn.bv+se.bv, lo=mn.bv-se.bv)
#Prep this treatment summary for the figure (new Fig. 2)
trt.sumsmallphyt$group <- "small"
trt.sumsmallphyt$tank <- "all"
trt.sumsmallphyt

#dot plot
ggplot(sumsmallphyt, aes(x=treatment, y=mn.bv, group=tank)) + 
  geom_errorbar(aes(ymin=lo, ymax=hi), width=0.4,
                 position=position_dodge(width=0.4)) +
  geom_point(position=position_dodge(width=0.4)) + theme_bw(10) +
  scale_y_log10(breaks=c(1e+01,1e+02,1e+03,1e+04,1e+05),
                labels=c("0.01","0.1","1","10","100")) + annotation_logticks(sides="l") +
  labs(y=expression(paste("Mean smaller phytoplankton biovolume (",mu,m^3,"/nL)")), 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"))

#Analyze small phyto biovolume

#Full model
small.gam <- glmer(bv ~ zoop + Neo + (1|tank), family=Gamma(link="log"), data=sumsmallphyto0s)
summary(small.gam)
#             Estimate Std. Error t value Pr(>|z|)    
#(Intercept)  8.269154   0.646907  12.783   <2e-16 ***
#zoopy       -0.792739   0.872469  -0.909    0.364    
#Neo          0.007051   0.011324   0.623    0.534  

#W/out Neoplea
small.noNeo <- glmer(bv ~ zoop + (1|tank), family=Gamma(link="log"), data=sumsmallphyto0s)
summary(small.noNeo)
#            Estimate Std. Error t value Pr(>|z|)    
#(Intercept)   8.2692     0.6539  12.645   <2e-16 ***
#zoopy        -0.5117     0.7543  -0.678    0.498

#Null model
small.null <- glmer(bv ~ (1|tank), family=Gamma(link="log"), data=sumsmallphyto0s)

#Compare models
anova(small.null, small.noNeo, small.gam)
#            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
#small.null     3 755.71 760.78 -374.85   749.71                     
#small.noNeo    4 757.25 764.01 -374.63   749.25 0.4544  1     0.5003
#small.gam      5 758.87 767.31 -374.43   748.87 0.3850  1     0.5349

#no sig. effects of zoop or Neoplea on small phyto biovolume


#Combine total, large, & small phyto summary subsets
phytsum <- rbind(sumphyt, sumbig2, sumsmallphyt); phytsum
trt.phytsum <- rbind(trt.sumphyt, trt.sumbigphyt, trt.sumsmallphyt); trt.phytsum


#Reorder phyto groups for figure
phytsum$group <- factor(phytsum$group, levels=c("total", "large", "small"))
trt.phytsum$group <- factor(trt.phytsum$group, levels=c("total", "large", "small"))

#Figure 2: faceted dot plot of total, large, & small phyto biovolume

pdf("figure 2 phyto faceted dot plot.pdf", width=7.5, height=3.5)
ggplot(phytsum, aes(x=treatment, y=mn.bv, group=tank)) + 
  geom_jitter(width=.2, color="grey") + 
  geom_point(data=trt.phytsum, aes(x=treatment, y=mn.bv)) +
  geom_errorbar(data=trt.phytsum, aes(x=treatment, ymin=lo, ymax=hi), width=0.4) +
  facet_wrap(~group) + theme_bw(10) +
  scale_y_log10(breaks=c(1e+01,1e+02,1e+03,1e+04,1e+05),
                labels=c("0.01","0.1","1","10","100")) + annotation_logticks(sides="l") +
  labs(y=expression(paste("Mean phytoplankton biovolume (",mu,m^3,"/nL)")), 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()