library(piecewiseSEM)
library(robustbase)
library(ciTools)
library(doBy)
library(ggplot2)
library(gridExtra)


####Zooplankton Abundance Data Files
n<-read.csv("mesocosmcounts_peerj.csv") #abundance data
#calculating total densities per mesocosm
colnames(n)[8:39]<-c("D.pl", "D.am", "D.pr", "D.mn", "D.ct", "D.db", "E.tb", "E.ln", "B.fr", "E.cr", "C.sc", "Cp", "Np", "C.lc", "C.sp", "D.th", "L.mn", "L.as", "M.rb", "S.vt", "A.vr","S.mc","S.or","M.ed", "D.br","S.Cr","E.Lc", "Rt", "P.pd", "Ch","Hl","By")
#summing up all of the individuals per species per mesocosm and then dividing by the total volume samples per mesocosm to get the density
x<-summaryBy(D.pl + D.am + D.pr + D.mn + D.ct + D.db + E.tb + E.ln + B.fr + E.cr + C.sc + Cp + Np + C.lc + C.sp + D.th + L.mn + L.as + M.rb + S.vt + A.vr + S.mc + S.or + M.ed + D.br + S.Cr + E.Lc + Rt + P.pd + Ch + Hl + By ~ Mesocosm + Origin + Dispersal + Bythotrephes + Week, data = n , FUN = (function(x) sum(x)/229.6125)) #whole mesocosm density
x1<-summaryBy(D.pl + D.am + D.pr + D.mn + D.ct + D.db + E.tb + E.ln + B.fr + E.cr + C.sc + Cp + Np + C.lc + C.sp + D.th + L.mn + L.as + M.rb + S.vt + A.vr + S.mc + S.or + M.ed + D.br + S.Cr + E.Lc + Rt + P.pd + Ch + Hl + By ~ Mesocosm + Origin + Dispersal + Bythotrephes + Week, data = n , FUN = (function(x)sum(x))) #full sample counts
colnames(x)<-c("Mesocosm", "Origin", "Dispersal", "Bythotrephes", "Week", "D.pl", "D.am", "D.pr", "D.mn", "D.ct", "D.db", "E.tb", "E.ln", "B.fr", "E.cr", "C.sc", "Cp", "Np", "C.lc", "C.sp", "D.th", "L.mn", "L.as", "M.rb", "S.vt", "A.vr","S.mc","S.or","M.ed", "D.br","S.Cr","E.Lc", "Rt", "P.pd", "Ch","Hl","By")
colnames(x1)<-c("Mesocosm", "Origin", "Dispersal", "Bythotrephes", "Week", "D.pl", "D.am", "D.pr", "D.mn", "D.ct", "D.db", "E.tb", "E.ln", "B.fr", "E.cr", "C.sc", "Cp", "Np", "C.lc", "C.sp", "D.th", "L.mn", "L.as", "M.rb", "S.vt", "A.vr","S.mc","S.or","M.ed", "D.br","S.Cr","E.Lc", "Rt", "P.pd", "Ch","Hl","By")
m1<-x

#####Minimum Detection Density#####
#determining the minimum detection density across the whole mesocosm
#minimum detection count would be 1 individuals per mesocosm
#total volume of water sampled in a mesocosm is 229.6125L
#Minimum number of trays per sample 7
#three layers per sample 
Tray<-3*7 #Minimum number of trays counted per sample
voltray<-6.657 #ml
volmes<-100*3 #each layer sample was  diluted to 100ml, 3 layers per mesocosm
totalvol<-229.6125 #Total volume sampled from each mesocosm
tottray<-Tray*voltray #total volume processed ml per mesocosm
propml<-tottray/volmes #proportion of volume processed 
volproc<-propml*totalvol #volume processed
Min<-1/volproc #0.0093
MinL<-Min/3 #the minimum detection density divided for the three layers so that the probability of a single individual in the water column of occuring in this layer is equal. 
#Want to add a number lower than the minimum detection limit
#0.009 will be added to all of the zeroes
#minimum detection count
propproc<-propml/totalvol 
MinC<-propproc*1 #Minimum detection count


Min2<-round(Min, 3) #rounding minimum detection density per mesocosm to three digits
MinL2<-round(MinL, 3) #rounding minimum detection density per layer to three digits
m2<-cbind(m1[,c(1:5)],(m1[,-c(1:5)] + Min2)) #dataset with final densities
str(m2)


####Data for Vertical position#####
layer<-read.csv("layerdensity_peerj.csv") #zooplankton density calculated by volume of layer
str(layer)
#abbreviating column names 
colnames(layer)[8:39]<-c("D.pl", "D.am", "D.pr", "D.mn", "D.ct", "D.db", 
	"E.tb", "E.ln", "B.fr", "E.cr", "C.sc", "Cp", "Np", "C.lc", "C.sp", 
	"D.th", "L.mn", "L.as", "M.rb", "S.vt", "A.vr","S.mc","S.or","M.ed", 
	"D.br","S.Cr","E.Lc", "Rt", "P.pd", "Ch","Hl","By")

#creating separate dataset for hypolimnetic layer
Epi0<-subset(layer, Layer == "Epi" & Week == 0)[,-c(1:7)] #Epi week 0
Meta0<-subset(layer, Layer == "Meta" & Week == 0)[,-c(1:7)] #Meta week 0
Hypo0<-subset(layer, Layer == "Hypo" & Week == 0)[,-c(1:7)] #Hypo week 0
Epi3<-subset(layer, Layer == "Epi" & Week == 3)[,-c(1:7)] #Epi week 3
Meta3<-subset(layer, Layer == "Meta" & Week == 3)[,-c(1:7)] #Meta week 3
Hypo3<-subset(layer, Layer == "Hypo" & Week == 3)[,-c(1:7)] #Hypo week 3

#creating a new dataset with minimum detectable densities added to 0s
Epi0.1<-Epi0+MinL2
Meta0.1<-Meta0+MinL2
Hypo0.1<-Hypo0+MinL2
Epi3.1<-Epi3+MinL2
Meta3.1<-Meta3+MinL2
Hypo3.1<-Hypo3+MinL2
#remove detection density for mean depth 
treatments<-rbind(subset(layer, Layer == "Epi" & Week == 0)[,c(1:7)], subset(layer, Layer == "Meta" & Week == 0)[,c(1:7)], subset(layer, Layer == "Hypo" & Week == 0)[,c(1:7)], subset(layer, Layer == "Epi" & Week == 3)[,c(1:7)], subset(layer, Layer == "Meta" & Week == 3)[,c(1:7)],subset(layer, Layer == "Hypo" & Week == 3)[,c(1:7)])
Dens<-rbind(Epi0.1, Meta0.1, Hypo0.1, Epi3.1, Meta3.1, Hypo3.1)
layer2<-cbind(treatments, Dens)

#Proportion Hypolimnetic individuals
Phypo0<-cbind(subset(layer, Layer == "Hypo" & Week == 3)[,c(1:6)],(Hypo0.1/(Epi0.1+Meta0.1+Hypo0.1))) 
#proportion of hypolimnetic individuals that are hypolimnetic in Week 0
Phypo3<-cbind(subset(layer, Layer == "Hypo" & Week == 3)[,c(1:6)],(Hypo3.1/(Epi3.1+Meta3.1+Hypo3.1))) 
#proportion of hypolimnetic individuals that are hypolimnetic in Week 3


#Creating dataset with proportion of individuals that are hypolimnetic
Vert1<-cbind(rbind(subset(layer, Layer == "Hypo" & Week == 0)[,c(1:7)],subset(layer, Layer == "Hypo" & Week == 3)[,c(1:7)]), rbind(Phypo0, Phypo3))

#Separating count data by week 0 and week 3
m0<-subset(x1, Week == 0) #count data
m3<-subset(x1, Week == 3) #count data

#Aggregating count data by groups
#Determining density for each group
daphinitial<-((m0[,6])+(m0[,7])+(m0[,8])+ (m0[,9])+(m0[,10]) +(m0[,11]))/229.6125 + Min2
daphfinal<-((m3[,6])+(m3[,9])+(m3[,10]) + (m3[,11])+(m3[,12]) +(m3[,13]))/229.6125 + Min2
smallcinitial<-((m0[,12]) + (m0[,13]) + (m0[,14]) + (m0[,19]) + (m0[,20]) + (m0[,30]))/229.6125 + Min2
smallcfinal<-((m3[,12] ) + (m3[,13] ) + (m3[,14]) + (m3[,19])+ (m3[,20]) + (m3[,30]))/229.6125 + Min2
largecinitial<-((m0[,31]) + (m0[,36]))/229.6125 + Min2 
largecfinal<-((m3[,31]) + (m3[,36]))/229.6125 + Min2
copepodidinitial<-(m0[,17])/229.6125 + Min2
copepodidfinal<-(m3[,17])/229.6125 + Min2
naupliiinitial<-(m0[,18])/229.6125 + Min2
naupliifinal<-(m3[,18])/229.6125 + Min2
calinitial<-((m0[,22]) + (m0[,30]))/229.6125 + Min2
calfinal<-((m3[,22]) + (m3[,30]))/229.6125 + Min2
cycinitial<-((m0[,18]) + (m0[,23]) + (m0[,31]))/229.6125 + Min2
cycfinal<-((m3[,18]) + (m3[,23]) + (m3[,31]))/229.6125 + Min2
juvenileinitial <-(m0[,17] + m0[,18])/229.6125 + Min2
juvenilefinal <- (m3[,17] + m3[,18])/229.6125 + Min2

#plot of starting densities
pdf("InitialDensities.pdf")
par(mfrow = c(3,2))
boxplot(daphinitial ~ Bythotrephes, data = m0, main = "Daphnia",  xlab = "Bythotrephes Treatment", ylab = "Density (Individuals/L)")
boxplot(smallcinitial ~ Bythotrephes, data = m0, main = "Small Cladocerans", xlab = "Bythotrephes Treatment", ylab = "Density (Individuals/L)")
boxplot(largecinitial ~ Bythotrephes, data = m0, main = "Large Cladocerans", xlab = "Bythotrephes Treatment", ylab = "Density (Individuals/L)" )
boxplot(calinitial ~ Bythotrephes, data = m0, main = "Calanoids", xlab = "Bythotrephes Treatment", ylab = "Density (Individuals/L)" ) 
boxplot(cycinitial ~ Bythotrephes, data = m0, main = "Cyclopoids", xlab = "Bythotrephes Treatment", ylab = "Density (Individuals/L)" ) 
boxplot(juvenileinitial ~ Bythotrephes, data = m0, main = "Juveniles", xlab = "Bythotrephes Treatment", ylab = "Density (Individuals/L)" ) 
dev.off()


#Change in density for groups
daphnia<-((daphfinal)/(daphinitial))
smallclad<-(smallcfinal/smallcinitial)
largeclad<-(largecfinal/largecinitial)
copepodid<-(copepodidfinal/copepodidinitial)
nauplii<-(naupliifinal/naupliiinitial)
calanoids<-(calfinal/calinitial)
cyclopoids<-(cycfinal/cycinitial)
juvenile <-(juvenilefinal/juvenileinitial)
#Dataset with per capita change in density for major taxonomic groups
mg<-cbind(m0[,c(1:5)],daphnia, smallclad, largeclad, calanoids, cyclopoids, copepodid, nauplii, juvenile)
#Creating datasets with per capita change in density at the species level between week 0 and week 3
Diff<-((subset(m1, Week == 3)[,-c(1:5)] + Min2)/(subset(m1, Week == 0)[,-c(1:5)] + Min2))
Diff2<-cbind(subset(m1, Week == 0)[,c(1:5)], Diff)

#####Calculating proportion of hypolimnetic Daphnia species in each mesocosm#####
#Calculating proportion hypolimnetic Daphnia
hypo1<-Phypo0[order(Phypo0$Mesocosm),] #ordering data by mesocosm 
Daphnia0<-((rowSums(Hypo0.1[,c(1:6)])/(rowSums(Epi0.1[,c(1:6)]+Meta0.1[,c(1:6)]+Hypo0.1[,c(1:6)]))))
Daphnia0.1<-data.frame(as.numeric(Daphnia0), Phypo0$Bythotrephes, Phypo0$Mesocosm)
colnames(Daphnia0.1)<-c("PropDaphHypo", "Byth", "Mesocosm")
Daphnia0.2<-Daphnia0.1[order(Daphnia0.1$PropDaphHypo),]

#putting vertical position data in same order as per capita growth data
#Proportion Hypolimnetic Daphnia ordered
Daphnia0.3<-Daphnia0.1[order(Daphnia0.1$Mesocosm),]

#####Merging Daphnia vertical position and per capita growth data#####
#combining difference in per capita growth data for major zooplankton taxonomic groups with the proportion of hypolimnetic Daphnia
mg1<-cbind(mg, Daphnia0.3)
#combining difference in per capita growth for each zooplankton species with the proportion of hypolimnetic Daphnia
Diff3<-cbind(Diff2, Daphnia0.3)

##### Loading and Preparing Algae dataset#####
algae<-read.csv("Algae.csv")
colnames(algae)<-c("Lake", "Mesocosm", "Origin", "Bythotrephes", "Week", "SampleDate", "Layer", "Date/Time", "Comment", "Total", "Green", "BlueGreen", "Diatoms", "Cryptophyta", "Yellowsubstances")
al3<-summaryBy(Total + Green + BlueGreen + Diatoms + Cryptophyta ~ Mesocosm + Origin + Bythotrephes + Week, data = subset(algae, Week == 3) , FUN = (function(x) sum(x))) #full water column algae
colnames(al3)<-c("Mesocosm", "Origin", "Bythotrephes", "Week", "Total", "Green", "BlueGreen", "Diatoms", "Cryptophyta")
al3.1<-al3[,-c(1:4)] + 0.01 #density in week 3


####Merging Algae and Zooplankton datasets
#merging algae and zooplankton groups datasets
mg1.1 <-mg1[,c(-15,-16)] #removing duplicate columns for Bythotrephes presence and absence and mesocosm number
zoop_algae <- merge(mg1.1, al3,  by = "Mesocosm")
#merging algae and zooplankton species datasets
mgdiff<-Diff3[,-c(39:40)] #removing duplicate columns for Bythotrephes presence and absence and mesocosm number
zoop_species_algae <- merge(mgdiff, al3, by = "Mesocosm")
#after checking the merge, removing double columns for Origin, Week and Bythotrephes
#also adding final densities for all daphnia, small cladocerans, large cladocerans, calanoids, and juveniles 
#I will recreate SEMs with total algae and all algal groups fit with final densities
zoop_algae_1 <-cbind(zoop_algae[,-c(15:22)], al3.1$Total, al3.1$Green, al3.1$BlueGreen, al3.1$Diatoms, al3.1$Cryptophyta)
zoop_species_algae1<-cbind(zoop_species_algae[,-c(39:46)],al3.1$Total, al3.1$Green, al3.1$BlueGreen, al3.1$Diatoms, al3.1$Cryptophyta)
str(zoop_algae_1)
#correcting names for Origin, Bythotrephes, Week, TotalAlgae
colnames(zoop_algae_1)[c(2,4:5,15:19)] <- c("Origin", "Bythotrephes", "Week", "FinalTotalAlgae","FinalGreen", "FinalBlueGreen", "FinalDiatoms", "FinalCryptophyta")
colnames(zoop_species_algae1)[c(2,4:5,39:43)] <- c("Origin", "Bythotrephes", "Week", "FinalTotalAlgae","FinalGreen", "FinalBlueGreen", "FinalDiatoms", "FinalCryptophyta")

#####Building Structural Equation Models#####
#creating a new variable where Bythotrephes is converted into an ordinal variable
zoop_algae_1$Bythotrephes1 <- as.numeric(zoop_algae_1$Bythotrephes =="Y")
zoop_species_algae1$Bythotrephes1 <- as.numeric(zoop_species_algae1$Bythotrephes == "Y")

#using piecewise SEMs because of low sample size
all_daphnia_sem<- psem(
 glm(daphnia ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(smallclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_algae_1),
 glm(juvenile ~ Bythotrephes1 + PropDaphHypo,  family = gaussian(link = log), data = zoop_algae_1),
 glm(calanoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(largeclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(cyclopoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(FinalTotalAlgae ~ Bythotrephes1 + PropDaphHypo + daphnia + smallclad + largeclad + calanoids + juvenile, family = gaussian(link = log), data = zoop_algae_1),
data = zoop_algae_1
)

#examining output and specifiying correlated errors
summary(update(all_daphnia_sem, smallclad %~~% daphnia, juvenile %~~% daphnia, largeclad %~~% daphnia, calanoids %~~% daphnia, 
	cyclopoids %~~% daphnia, juvenile %~~% smallclad, largeclad %~~% smallclad, calanoids %~~% smallclad,  cyclopoids %~~% smallclad, 
	largeclad %~~% juvenile, calanoids %~~% largeclad, cyclopoids %~~% largeclad, cyclopoids %~~% calanoids, FinalTotalAlgae %~~% cyclopoids))

#####Fitting separate SEM models for different algal groups
#Green Algae
all_daphnia_greenalgae_sem<- psem(
 glm(daphnia ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(smallclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_algae_1),
 glm(juvenile ~ Bythotrephes1 + PropDaphHypo,  family = gaussian(link = log), data = zoop_algae_1),
 glm(calanoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(largeclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(cyclopoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(FinalGreen ~ Bythotrephes1 + PropDaphHypo + daphnia + smallclad + largeclad + calanoids + juvenile, family = gaussian(link = log), data = zoop_algae_1),
data = zoop_algae_1
)
summary(update(all_daphnia_greenalgae_sem, smallclad %~~% daphnia, juvenile %~~% daphnia, largeclad %~~% daphnia, juvenile %~~% daphnia, 
	cyclopoids %~~% daphnia, juvenile %~~% smallclad, largeclad %~~% smallclad, calanoids %~~% smallclad, cyclopoids %~~% smallclad, 
	largeclad %~~% juvenile, calanoids %~~% largeclad, cyclopoids %~~% largeclad, cyclopoids %~~% calanoids, FinalGreen %~~% cyclopoids))


#Blue-Green Algae - Cyanobacteria
all_daphnia_bluegreenalgae_sem<- psem(
 glm(daphnia ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(smallclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_algae_1),
 glm(juvenile ~ Bythotrephes1 + PropDaphHypo,  family = gaussian(link = log), data = zoop_algae_1),
 glm(calanoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(largeclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(cyclopoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(FinalBlueGreen ~ Bythotrephes1 + PropDaphHypo + daphnia + smallclad + largeclad + calanoids + juvenile, family = gaussian(link = log), data = zoop_algae_1),
data = zoop_algae_1
)
summary(update(all_daphnia_bluegreenalgae_sem, smallclad %~~% daphnia, juvenile %~~% daphnia, largeclad %~~% daphnia, juvenile %~~% daphnia, 
	cyclopoids %~~% daphnia, juvenile %~~% smallclad, largeclad %~~% smallclad, calanoids %~~% smallclad, cyclopoids %~~% smallclad, 
	largeclad %~~% juvenile, calanoids %~~% largeclad, cyclopoids %~~% largeclad, cyclopoids %~~% calanoids, FinalBlueGreen %~~% cyclopoids))

#Diatoms
all_daphnia_diatoms_sem<- psem(
 glm(daphnia ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(smallclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_algae_1),
 glm(juvenile ~ Bythotrephes1 + PropDaphHypo,  family = gaussian(link = log), data = zoop_algae_1),
 glm(calanoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(largeclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(cyclopoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(FinalDiatoms ~ Bythotrephes1 + PropDaphHypo + daphnia + smallclad + largeclad + calanoids + juvenile, family = gaussian(link = log), data = zoop_algae_1),
data = zoop_algae_1
)
summary(update(all_daphnia_diatoms_sem, smallclad %~~% daphnia, juvenile %~~% daphnia, largeclad %~~% daphnia, calanoids %~~% daphnia, 
	cyclopoids %~~% daphnia, juvenile %~~% smallclad, largeclad %~~% smallclad, calanoids %~~% smallclad,  cyclopoids %~~% smallclad, 
	largeclad %~~% juvenile, calanoids %~~% largeclad, cyclopoids %~~% largeclad, cyclopoids %~~% calanoids, FinalDiatoms %~~% cyclopoids))

#Cryptophyta
all_daphnia_crytophyta_sem2<- psem(
 glm(daphnia ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(smallclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_algae_1),
 glm(juvenile ~ Bythotrephes1 + PropDaphHypo,  family = gaussian(link = log), data = zoop_algae_1),
 glm(calanoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(largeclad ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_algae_1),
 glm(cyclopoids ~ Bythotrephes1 + PropDaphHypo + juvenile, family = gaussian(link = log), data = zoop_algae_1),
 glm(FinalCryptophyta ~ Bythotrephes1 + PropDaphHypo + daphnia + smallclad + largeclad + calanoids + juvenile, family = gaussian(link = log), data = zoop_algae_1),
data = zoop_algae_1
)
x<-summary(update(all_daphnia_crytophyta_sem2, smallclad %~~% daphnia, juvenile %~~% daphnia, largeclad %~~% daphnia, calanoids %~~% daphnia, 
	cyclopoids %~~% daphnia, juvenile %~~% smallclad, largeclad %~~% smallclad, calanoids %~~% smallclad,  cyclopoids %~~% smallclad, 
	largeclad %~~% juvenile, calanoids %~~% largeclad, cyclopoids %~~% largeclad, cyclopoids %~~% calanoids, FinalCryptophyta %~~% cyclopoids))

####SEMS for the most common species within each zooplankton group
common_species_sem<- psem(
 glm(D.mn ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(D.ct ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(B.fr ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_species_algae1),
 glm(E.tb ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(E.ln ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(S.or ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(C.sc ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
data = zoop_species_algae1
) 
z<-summary(update(common_species_sem,
	B.fr %~~% D.mn, D.ct %~~% D.mn, B.fr %~~% D.ct,
	E.tb %~~% D.mn, E.tb %~~% D.ct, E.tb %~~% B.fr, E.tb %~~% Cp, E.tb %~~% S.or, 
	C.sc %~~% E.tb, S.or %~~% D.mn, C.sc %~~% D.mn, S.or %~~% D.ct, C.sc %~~% D.ct, S.or %~~% B.fr, C.sc %~~% B.fr, S.or %~~% E.ln, C.sc %~~% E.ln,  E.ln %~~% D.mn, E.ln %~~% D.ct, E.ln %~~% B.fr, E.ln %~~% E.tb))

common_species_sem_green<- psem(
 glm(D.mn ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(D.ct ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(B.fr ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log),  data = zoop_species_algae1),
 glm(E.tb ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(E.ln ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(Hl ~ Bythotrephes1 + PropDaphHypo, family = gaussian(link = log), data = zoop_species_algae1),
 glm(FinalGreen ~ Bythotrephes1 + PropDaphHypo + D.mn + D.ct + B.fr + E.tb + E.ln + Hl, family = gaussian(link = log), data = zoop_species_algae1),
data = zoop_species_algae1
) 
x<-summary(update(common_species_sem_green, C.sc %~~% PropDaphHypo, 
	B.fr %~~% D.mn, D.ct %~~% D.mn, B.fr %~~% D.ct,
	E.tb %~~% D.mn, E.tb %~~% D.ct, E.tb %~~% B.fr, E.tb %~~% Cp, E.tb %~~% S.or, E.tb %~~% Hl, 
	C.sc %~~% E.tb, E.ln %~~% D.mn, E.ln %~~% D.ct, E.ln %~~% B.fr, E.ln %~~% E.tb))


#####Running GLMs to test for interactions for causal relationships with a significant effect of Bythotrephes and Daphnia Vertical Position
library(fitdistrplus)
library(MASS)


#Daphnia
glmg1<-glm(daphnia ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = mg, family =gaussian(link = "log"))
glmg2<-glm(daphnia ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = mg, family = Gamma(link = "log"))
#Gamma is a better model
descdist(daphnia1, discrete = FALSE, boot = 1000)
plot(glmg1) #decreasing resid-fit plot, increasing-scale loc plot
plot(glmg2) #looks better
#looks ok
#log likelihood ratio tests for main effects
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmg2.1<-update(glmg2, .~. -Origin)
anova(glmg2.1, glmg2, test = "Chisq") #p = 0.5852

glmg2.2<-update(glmg2, .~. -Bythotrephes:Daphnia0.3$PropDaphHypo)
anova(glmg2, glmg2.2, test = "Chisq") #p = 0.0211

#Daphnia mendotae
glmdmn1<-glm(D.mn ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = Diff2, family =gaussian(link = "log"))
glmdmn2<-glm(D.mn ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = Diff2, family = Gamma(link = "log"))
AIC(glmdct1, glmdct2)
#Gamma is a better model
descdist(Diff2$D.mn, discrete = FALSE, boot = 1000)
plot(glmdmn1) #decreasing resid-fit plot, increasing-scale loc plot
plot(glmdmn2) #looks better
#looks ok
#log likelihood ratio tests for main effects
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmdmn2.1<-update(glmdmn2, .~. -Origin)
anova(glmdmn2, glmdmn2.1) 
#Not significant p = 0.5562

glmdmn2.2<-update(glmdmn2, .~. -Bythotrephes:Daphnia0.3$PropDaphHypo)
anova(glmdmn2, glmdmn2.2, test = "Chisq")
#Interaction is not significant p = 0.1838

#Small Cladocerans
glmgds1<-glm(smallclad ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = mg, family =gaussian(link = "log"))
glmgds2<-glm(smallclad~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = mg, family = Gamma(link = "log"))
AIC(glmgds1, glmgds2) #Gamma is better
plot(glmgds1)
plot(glmgds2) #looks OK (no outliers)
#scale-location plot is increasing for log-normal, better for gamma. Gamma has better residual-fitted plot (not decreaing) as compared to log-normal
#log likelihood ratio tests for main effects
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmgds2.1<-update(glmgds2, .~. -Origin)
anova(glmgds2, glmgds2.1, test = "Chisq")
#p = 0.01073

glmgds2.2<-update(glmgds2, .~. -Bythotrephes:Daphnia0.3$PropDaphHypo)
anova(glmgds2, glmgds2.2, test = "Chisq")
#Interaction is significant p = 0.04981

#Bosmina freyi
glmdbf1<-glm(B.fr ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = Diff2, family =gaussian(link = "log"))
glmdbf2<-glm(B.fr ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = Diff2, family = Gamma(link = "log"))
AIC(glmdbf1, glmdbf2) 
descdist(Diff2$B.fr, discrete = FALSE, boot = 1000)
plot(glmdbf1) #Log normal is better than Gamma
plot(glmdbf2) 
#log likelihood ratio tests
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmdbf2.1<-update(glmdbf2, .~. -Origin)
anova(glmdbf2, glmdbf2.1, test = "Chisq") #significant effect of origin
#p <0.0001

glmdbf2.2<-update(glmdbf2, .~. -Bythotrephes:Daphnia0.3$PropDaphHypo)
anova(glmdbf2, glmdbf2.2, test = "Chisq") #significant effect of interaction despite effect of origin
#p = 0.008


#Eubosmina longispina
glmdel1<-glm(E.ln ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = Diff2, family =gaussian(link = "log"))
glmdel2<-glm(E.ln ~ Bythotrephes*Daphnia0.3$PropDaphHypo + Origin, data = Diff2, family = Gamma(link = "log"))
AIC(glmdel1, glmdel2) #Gaussian is better
descdist(Diff2$E.ln, discrete = FALSE, boot = 1000)
plot(glmdel1) #increasing scale-location (resid-fit and scale-loc plot in this one is better)
plot(glmdel2) # decreasing scale-location (no outliers)

#log likelihood ratio tests
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmdel2.1<-update(glmdel1, .~. -Origin)
anova(glmdel1, glmdel2.1, test = "Chisq") #p = 0.3633

glmdel2.2<-update(glmdel2.1, .~. -Bythotrephes:Daphnia0.3$PropDaphHypo)
anova(glmdel1, glmdel2.2, test = "Chisq") #no significant interaction p = 0.3186

glmdel2.3<-update(glmdel2.2, .~. -Bythotrephes )
anova(glmdel2.2, glmdel2.3, test = "Chisq")

#Copepods
#Cyclopoids
glmcyc1 <- glm(cyclopoids ~ Bythotrephes*Daphnia0.3$PropDaphHypo + as.numeric(as.factor(Origin)), data = mg, family = gaussian(link = "log"))
glmcyc2 <- glm(cyclopoids ~ Bythotrephes*Daphnia0.3$PropDaphHypo + as.numeric(as.factor(Origin)), data = mg, family = Gamma(link = "log"), maxit = 10)
AIC(glmcyc1, glmcyc2)
#Gamma distribution does not fit due to NaNs produced
plot(glmcyc1) #fit does not look good, outliers
plot(glmcyc2) #looks good, despite lack of convergence

#log-likelihood ratio tests
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmcyc2.1 <- update(glmcyc2, .~. -as.numeric(as.factor(Origin)))
anova(glmcyc2, glmcyc2.1, test = "Chisq") #interaction is not significant
#Origin is not significant p = 0.05076
glmcyc2.2 <- update(glmcyc2.1, .~. -Bythotrephes:Daphnia0.3$PropDaphHypo)
anova(glmcyc2.2, glmcyc2.1, test = "Chisq") #interaction is not significant (p = 0.8918)


#Cyclopds scutifer
#Cyclops scutifer with Bythotrephes and proportion of Daphnia in the hypolimnion
glmcsc1<-glm(C.sc ~ Bythotrephes*PropDaphHypo + Origin, data = Diff3, family =gaussian(link = "log"))
glmcsc2<-glm(C.sc ~ Bythotrephes*PropDaphHypo + Origin, data = Diff3, family = Gamma(link = "log"))
AIC(glmcsc1, glmcsc2)
#Gamma has the better fit
#log likelihood ratio tests for main effects
 
#testing for the effect of origin
glmcsc2.1<-update(glmcsc2, .~. -Origin)
anova(glmcsc2, glmcsc2.1) #p = 0.126 (not significant)

#testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glmcsc2.2<-update(glmcsc2.1, .~. -Bythotrephes:PropDaphHypo)
anova(glmcsc2.1, glmcsc2.2, test = "Chisq")
#Bythotrephes: Proportion of daphnia in the hypolimnion is significant (p = 0.02504)

#Total Algae - applying robust regression because of outliers
glad1<-glm(FinalTotalAlgae ~ Bythotrephes*PropDaphHypo + Origin, data = zoop_algae_1, family = Gamma(link = "log"))
glad1.1<-glm(FinalTotalAlgae ~ Bythotrephes*PropDaphHypo + Origin, data = zoop_algae_1, family = gaussian(link = "log"))
plot(glad1) #looks OK, 
plot(glad1.1) #looks OK, fit is better with Gamma 

#one outlier in gamma model, mesocosm #23. Running a robust regression with Gamma model instead
glad1.11<-glmrob(FinalTotalAlgae ~ Bythotrephes*PropDaphHypo + Origin, data = zoop_algae_1, family = Gamma(link = "log"))
glad1.12<-update(glad1.11, .~. - Origin)
anova(glad1.11, glad1.12) #p = 0.001005 significant effect of Origin

glad1.13<-update(glad1.11, .~. - Bythotrephes:PropDaphHypo)
anova(glad1.11, glad1.13) #P 0.00005 

summary(glad1.11)
covMcd(glad1.11)

fitted<-summary(glad1.11)$fitted
residuals <-summary(glad1.11)$residuals

#making diagnostic plots for robust regression
plot(residuals ~ fitted)
plot(sqrt(residuals) ~ fitted)

#getting real dispersion values for Gamma distribution
myshape <- gamma.shape(glm(fitted ~ Bythotrephes*PropDaphHypo, data = zoop_algae_1, family = Gamma(link = "log")))
predict(glad1.1, type = "response", se = T, dispersion = 1/myshape$alpha)
 dispersion = 1/myshape$alpha

#Green Algae
glgdn1<-glm(FinalGreen ~ Bythotrephes*PropDaphHypo + Origin, data = zoop_algae_1, family = Gamma(link = "log"))
plot(glgdn1) #looks OK
#log likelihood ratio tests
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
glgdn1.1<-update(glgdn1, .~. - Origin)
anova(glgdn1.1, glgdn1) #significant effect 0.0245
glgdn1.2<-update(glgdn1, .~. -Bythotrephes:PropDaphHypo)
anova(glgdn1, glgdn1.2, test = "Chisq")
#Interaction is not significant p = 0.7805
glgdn1.3<-update(glgdn1.1, .~. -Bythotrephes:PropDaphHypo)
anova(glgdn1.2, glgdn1.3, test = "Chisq")

glgdn1.4<-update(glgdn1.3, .~. -Bythotrephes)
anova(glgdn1.4, glgdn1.3, test = "Chisq")
#Diatoms
gldd1<-glm(FinalDiatoms~ Bythotrephes*PropDaphHypo + Origin, data = zoop_algae_1, family = Gamma(link = "log"))
plot(gldd1) #looks good, but outliers
#trying log normal distribution
gldd1.1<-update(gldd1, family = gaussian(link = "log")) 
plot(gldd1.1) #increasing scale-location, 

#using robust regression 
gldd1.11<-glmrob(FinalDiatoms ~ Bythotrephes*PropDaphHypo + Origin, data = zoop_algae_1, family = Gamma(link = "log"), method = "Mqle", control = glmrobMqle.control(maxit = 100)) #algorithm does not converge for Gamma, applying log-normal distribution


#running log likelihood ratio tests (Wald)
#only testing if the interactive effects of Bythotrephes and Daphnia vertical position are significant
gldd1.10<-update(gldd1.11, .~. - Origin)
anova(gldd1.11, gldd1.10, test = "Wald") #p = 0.1144
gldd1.12<-update(gldd1.10, .~. - Bythotrephes:PropDaphHypo)
anova(gldd1.11, gldd1.12, test = "Wald") #p = 0.1547
#significant effect 

####Bootstrapping analysis#############
#bootstrapping per capita change in density versus starting densities
library(boot)
#correlation between  per capita Daphnia density versus starting densities
cor(mg$daphnia, rowSums(subset(m2, Week == 0)[,c(6:32,36)]))
#-0.3404672

#bootstrapping
corr <- function(data, indices) {
  d <- data[indices,] # allows boot to select sample
  cor(d[,1], d[,2])
} 
set.seed(69)
z<-cbind(rowSums(subset(m2, Week == 0)[,c(6:32,36)]), mg$daphnia, mg$smallclad, mg$cyclopoids, mg$copepodid, mg$largeclad)
colnames(z)<-c("TotalZoop", "Daphnia", "Smallclad", "Cyclopoids", "Copepodid", "Largeclad")
z<-data.frame(z)
str(z)

#Daphnia
myBootstrap <- boot(z[,c(1:2)], corr, R=1000)

hist(myBootstrap$t, main = "Daphnia", xlab = "Correlation coefficient")
abline(v = -0.3404, lty = 2, lwd = 3, col = "red")

#Small Cladocerans
cor(mg$smallclad, rowSums(subset(m2, Week == 0)[,c(6:32,36)]))
#-0.41
myBootstrap <- boot(z[,c(1,3)], corr, R=1000)
hist(myBootstrap$t, main = "Small Cladocerans", xlab = "Correlation coefficient")
abline(v = -0.41, lty = 2, lwd = 3, col = "red")

#Large cladoceran density
cor(mg$largeclad, rowSums(subset(m2, Week == 0)[,c(6:32,36)]))
#-0.23
myBootstrap <- boot(z[,c(1,6)], corr, R=1000)
hist(myBootstrap$t, main = "Large Cladocerans", xlab = "Correlation coefficient")
abline(v = -0.23, lty = 2, lwd = 3, col = "red")

#Cyclopoid density
cor(mg$cyclopoids, rowSums(subset(m2, Week == 0)[,c(6:32,36)]))
#-0.437
myBootstrap <- boot(z[,c(1,4)], corr, R=1000)
hist(myBootstrap$t, main = "Cyclopoids", xlab = "Correlation coefficient")
abline(v = -0.437, lty = 2, lwd = 3, col = "red")
#this is very different!!!

#Copepodid density
cor(mg$copepodid, rowSums(subset(m2, Week == 0)[,c(6:32,36)]))
#-0.75
myBootstrap <- boot(z[,c(1,5)], corr, R=1000)
hist(myBootstrap$t, main = "Copepodids", xlab = "Correlation coefficient")
abline(v = -0.75, lty = 2, lwd = 3, col = "red")
dev.off()

#None of the correlations are different from random