#LUMINOXYSCAPE Analysis PART 3 # Calculate average and standard deviation/standard error (SD/SE) of maxdepths # Make some figures # Run stats # Please contact lrmccorm@ucsd.edu for questions/issues with this code. # Species abbreviations are: # Do/D. opalescens = California Market Squid, "Doryteuthis opalescens" # Ob/O. bimaculatus = Two-spot octopus, "Octopus bimaculatus" # Cg/Mg/M. gracilis = Graceful rock crab, "Metacarcinus gracilis" # formally "Cancer gracilis" # Pp/P. planipes = Tuna crab/red pelagic crab, "Pleuroncodes planipes" library("lubridate") library("ggplot2") library("dplyr") library("tidyr") library("viridis") library("viridisLite") library("respirometry") library("geosphere") library("ggmap") library("stringr") library("cmocean") setwd("~SET WD") # Add yours setwd("~/Documents/SIO/Luminoxyscape/Analysis_May21/") # Calculate Standard Error using this function: stderr <- function(x, na.rm=FALSE) { if (na.rm) x <- na.omit(x) sqrt(var(x)/length(x)) } # These load dataframes from Luminox part 2 code load("~/Documents/SIO/Luminoxyscape/Analysis_May21/SCB.sta.lims.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.SCB.lims.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/LLim_maxdepth.PFD.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_Do.0.PFD.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_Ob.0.PFD.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_Cg.0.PFD.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_Pp.0.PFD.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_V50.Do.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_V50.Ob.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_V50.Cg.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_V50.Pp.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_Met.Do.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/maxdepth_Met.Ob.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.maxdepth_V50.Do.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.maxdepth_V50.Ob.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.maxdepth_V50.Cg.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.maxdepth_V50.Pp.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.maxdepth_Met.Do.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/E.maxdepth_Met.Ob.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.SCB.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.Do.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.Ob.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.Cg.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.Pp.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.met.Do.Rdata") load("~/Documents/SIO/Luminoxyscape/Analysis_May21/Comb.met.Ob.Rdata") #Also run these: comb.Pp.sp <- comb.Pp[which(comb.Pp$Season=="Spring"),] row.names(comb.Pp.sp) <- 1:nrow(comb.Pp.sp) comb.Do.sp <- comb.Do[which(comb.Do$Season=="Spring"),] row.names(comb.Do.sp) <- 1:nrow(comb.Do.sp) comb.Ob.sp <- comb.Ob[which(comb.Ob$Season=="Spring"),] row.names(comb.Ob.sp) <- 1:nrow(comb.Ob.sp) comb.Cg.sp <- comb.Cg[which(comb.Cg$Season=="Spring"),] row.names(comb.Cg.sp) <- 1:nrow(comb.Cg.sp) # A few changes to Dist.land.... Channel island interference: RESAVED FILES SO SHOULDN'T NEED TO DO THIS AGAIN # comb.SCB$Dist.land[which(comb.SCB$Line == 86.7 & comb.SCB$Sta == 45)] <- 28 # SCB.sta$Dist.land[which(SCB.sta$Line == 86.7 & SCB.sta$Sta == 45)] <- 28 # maxdepth_L.Lim.PFD$Dist.land[which(maxdepth_L.Lim.PFD$Line == 86.7 & maxdepth_L.Lim.PFD$Sta == 45)] <- 28 # maxdepth_Do.0.PFD$Dist.land[which(maxdepth_Do.0.PFD$Line == 86.7 & maxdepth_Do.0.PFD$Sta == 45)] <- 28 # maxdepth_Ob.0.PFD$Dist.land[which(maxdepth_Ob.0.PFD$Line == 86.7 & maxdepth_Ob.0.PFD$Sta == 45)] <- 28 # maxdepth_Cg.0.PFD$Dist.land[which(maxdepth_Cg.0.PFD$Line == 86.7 & maxdepth_Cg.0.PFD$Sta == 45)] <- 28 # maxdepth_Pp.0.PFD$Dist.land[which(maxdepth_Pp.0.PFD$Line == 86.7 & maxdepth_Pp.0.PFD$Sta == 45)] <- 28 # maxdepth_V50.Do$Dist.land[which(maxdepth_V50.Do$Line == 86.7 & maxdepth_V50.Do$Sta == 45)] <- 28 # maxdepth_V50.Ob$Dist.land[which(maxdepth_V50.Ob$Line == 86.7 & maxdepth_V50.Ob$Sta == 45)] <- 28 # maxdepth_V50.Cg$Dist.land[which(maxdepth_V50.Cg$Line == 86.7 & maxdepth_V50.Cg$Sta == 45)] <- 28 # maxdepth_V50.Pp$Dist.land[which(maxdepth_V50.Pp$Line == 86.7 & maxdepth_V50.Pp$Sta == 45)] <- 28 # maxdepth_Met.Do$Dist.land[which(maxdepth_Met.Do$Line == 86.7 & maxdepth_Met.Do$Sta == 45)] <- 28 # maxdepth_Met.Ob$Dist.land[which(maxdepth_Met.Ob$Line == 86.7 & maxdepth_Met.Ob$Sta == 45)] <- 28 # E.maxdepth_V50.Do$Dist.land[which(E.maxdepth_V50.Do$St_Line == 86.7 & E.maxdepth_V50.Do$St_Station == 45)] <- 28 # E.maxdepth_V50.Ob$Dist.land[which(E.maxdepth_V50.Ob$St_Line == 86.7 & E.maxdepth_V50.Ob$St_Station == 45)] <- 28 # E.maxdepth_V50.Cg$Dist.land[which(E.maxdepth_V50.Cg$St_Line == 86.7 & E.maxdepth_V50.Cg$St_Station == 45)] <- 28 # E.maxdepth_V50.Pp$Dist.land[which(E.maxdepth_V50.Pp$St_Line == 86.7 & E.maxdepth_V50.Pp$St_Station == 45)] <- 28 # E.maxdepth_Met.Do$Dist.land[which(E.maxdepth_Met.Do$St_Line == 86.7 & E.maxdepth_Met.Do$St_Station == 45)] <- 28 # E.maxdepth_Met.Ob$Dist.land[which(E.maxdepth_Met.Ob$St_Line == 86.7 & E.maxdepth_Met.Ob$St_Station == 45)] <- 28 # comb.Do$Dist.land[which(comb.Do$Line == 86.7 & comb.Do$Sta == 45)] <- 28 # comb.Ob$Dist.land[which(comb.Ob$Line == 86.7 & comb.Ob$Sta == 45)] <- 28 # comb.Cg$Dist.land[which(comb.Cg$Line == 86.7 & comb.Cg$Sta == 45)] <- 28 # comb.Pp$Dist.land[which(comb.Pp$Line == 86.7 & comb.Pp$Sta == 45)] <- 28 # comb.met.Do$Dist.land[which(comb.met.Do$Line == 86.7 & comb.met.Do$Sta == 45)] <- 28 # comb.met.Ob$Dist.land[which(comb.met.Ob$Line == 86.7 & comb.met.Ob$Sta == 45)] <- 28 ########### SAVING AS CSV for Matlab analysis # Save as csv write.csv(SCB.sta,"SCB.sta.csv", row.names = FALSE) write.csv(md.V50.Do,"md.V50.Do.csv", row.names = FALSE) write.csv(md.V50.Ob,"md.V50.Ob.csv", row.names = FALSE) write.csv(md.V50.Cg,"md.V50.Cg.csv", row.names = FALSE) write.csv(md.V50.Pp,"md.V50.Pp.csv", row.names = FALSE) write.csv(md.L.Lim,"md.L.Lim.csv", row.names = FALSE) write.csv(md.Met.Do,"md.Met.Do.csv", row.names = FALSE) write.csv(md.Met.Ob,"md.Met.Ob.csv", row.names = FALSE) SCB.sm <- subset.data.frame(SCB.sta, select = c(Depth, OxCorr.kPa, Dist.land)) write.csv(SCB.sm,"SCB.sta.sm.csv", row.names = FALSE) ############################################################################################# ### Calculate mean max depth (md) for the dfs #### Add in SD/SE for these data max(SCB.sta$Dist.land) max(comb.SCB$Dist.land, na.rm= TRUE) # Take mean bin <- seq(0, 215, by = 2) # good for all raw.data <- maxdepth_L.Lim.PFD data <- data.frame() #blank dataframe for (d in 1:(length(bin)-1)) { i <- which(maxdepth_L.Lim.PFD$Dist.land>=bin[d] & maxdepth_L.Lim.PFD$Dist.land=bin[d] & maxdepth_V50.Do$Dist.land=bin[d] & maxdepth_V50.Ob$Dist.land=bin[d] & maxdepth_V50.Cg$Dist.land=bin[d] & maxdepth_V50.Pp$Dist.land=bin[d] & maxdepth_Met.Do$Dist.land=bin[d] & maxdepth_Met.Ob$Dist.land=tbin[y] & comb.Do$Year< tbin[y+1]) nd <- data.frame("Species" = "D. opalescens", #makes a column with just the trial name "Metric" = "V50", "Year" = tbin[y], "Depth.m" = mean(raw.data$Depth[i], na.rm= TRUE), "SD.m" = sd(raw.data$Depth[i], na.rm = TRUE), "SE.m" = stderr(raw.data$Depth[i], na.rm = TRUE)) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d y.md.Do <- data ######################################################################################################################################################## ######################################################################################################################################################## ######################################################################################################################################################## ##### D opal - Squid # All seasons ggplot(comb.Do, aes(x= comb.Do$Year, y= comb.Do$Depth, group= comb.Do$Year, fill= comb.Do$ENSO.OI))+ geom_boxplot()+ #geom_point(data= comb.Do, aes(comb.Do$Year, y= comb.Do$Depthm, color= comb.Do$ENSO.OI))+ scale_y_reverse(limits= c(300, 0), breaks=c(0, 50, 100, 150, 200, 250, 300))+ facet_wrap(facets = comb.Do$Loc)+ scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"))+ ggtitle(expression(paste(italic("D. opalescens"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Year") ggplot(comb.Do, aes(x= comb.Do$Year, y= comb.Do$Depth, group= comb.Do$Year, color= comb.Do$ENSO.OI))+ geom_point()+ #geom_point(data= comb.Do, aes(comb.Do$Year, y= comb.Do$Depthm, color= comb.Do$ENSO.OI))+ scale_y_reverse(limits= c(300, 0), breaks=c(0, 50, 100, 150, 200, 250, 300))+ facet_wrap(facets = comb.Do$Loc)+ scale_color_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "grey50"))+ ggtitle(expression(paste(italic("D. opalescens"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Year") #Just spring ggplot(comb.Do.sp, aes(x= comb.Do.sp$Year, y= comb.Do.sp$Depth, group= comb.Do.sp$Year, fill= comb.Do.sp$ENSO.Cat))+ geom_boxplot()+ #geom_point(data= comb.Do, aes(comb.Do$Year, y= comb.Do$Depthm, color= comb.Do$ENSO.OI))+ scale_y_reverse()+#limits= c(300, 0), breaks=c(0, 50, 100, 150, 200, 250, 300))+ facet_wrap(facets = comb.Do.sp$Loc)+ scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"))+ ggtitle(expression(paste(italic("D. opalescens"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Year") ggplot(comb.Do.sp, aes(x= comb.Do.sp$Year, y= comb.Do.sp$Depth, group= comb.Do.sp$Year, color= comb.Do.sp$ENSO.Cat))+ geom_point()+ #geom_point(data= comb.Do, aes(comb.Do$Year, y= comb.Do$Depthm, color= comb.Do$ENSO.OI))+ scale_y_reverse()+#limits= c(300, 0), breaks=c(0, 50, 100, 150, 200, 250, 300))+ facet_wrap(facets = comb.Do.sp$Loc)+ scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"))+ ggtitle(expression(paste(italic("D. opalescens"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Year") ################################################################################### ###### Final plots for squid # Seasonal/ENSO diff ggplot(comb.Do, aes(x=comb.Do$Season, y=comb.Do$Depth))+ geom_boxplot(fill= "grey")+ #geom_boxplot(data=comb.Do, aes(x= comb.Do$Season, y= comb.Do$Depth, fill = comb.Do$ENSO.Cat), alpha=0.7)+ #scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "fill")+ stat_summary(fun="mean", mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), geom="point", shape=20, size=5)+ stat_summary(fun.data = mean_se, fun.args = list(mult = 1), mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), width= 0.5, geom="errorbar", shape=20, size=1)+ scale_color_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "color")+ #stat_summary(fun="mean", mapping = aes(group = Season), position= position_dodge(width = .75), geom="point", shape=17, size=6, color="red", fill="red")+ scale_y_reverse(limits= c(390, 0), breaks=c(0, 50, 100, 150, 200, 250, 300, 350, 400))+ facet_wrap(facets = comb.Do$Loc)+ ggtitle(expression(paste(italic("D. opalescens"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Season") ################################################################################### ###### Stats for final figs: # D opalescens # 1. Seasonal ENSO differences in max depth #Is it normal? ggplot(comb.Do, aes(comb.Do$Depth))+ geom_histogram() #looks normal shapiro.test(comb.Do$Depth) # W = 0.98368, p-value = 6.031e-10 different than normal shapiro.test(maxdepth_V50.Do$Depth) # W = 0.96513, p-value = 8.04e-10 not normal kruskal.test(Depth~Season, data= comb.Do) dt.do <- dunn.test(x= comb.Do$Depth, g= comb.Do$Season, method = "bonferroni") # More interesting is within a season, do you see differences between ENSO yrs? comb.Do.sp <- comb.Do[which(comb.Do$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.sp) comb.Do.wi <- comb.Do[which(comb.Do$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.wi) comb.Do.su <- comb.Do[which(comb.Do$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.su) comb.Do.fa <- comb.Do[which(comb.Do$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.fa) ### Actual stats that take into account location (nearshore/offshore) comb.Do.NS <- comb.Do[which(comb.Do$Loc == "NS"),] comb.Do.OS <- comb.Do[which(comb.Do$Loc == "OS"),] shapiro.test(comb.Do$Depth) # not normal shapiro.test(comb.Do.NS$Depth) # W = 0.96417, p-value = 9.143e-08 different than normal shapiro.test(comb.Do.OS$Depth) # W = 0.98222, p-value = 4.935e-08 different than normal #NS kruskal.test(Depth~Season, data= comb.Do.NS) dt.do <- dunn.test(x= comb.Do.NS$Depth, g= comb.Do.NS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Do.NS) # More interesting is within a season, do you see differences between ENSO yrs? comb.Do.NS.sp <- comb.Do.NS[which(comb.Do.NS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.NS.sp) comb.Do.NS.wi <- comb.Do.NS[which(comb.Do.NS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.NS.wi) comb.Do.NS.su <- comb.Do.NS[which(comb.Do.NS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.NS.su) comb.Do.NS.fa <- comb.Do.NS[which(comb.Do.NS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.NS.fa) #OS kruskal.test(Depth~Season, data= comb.Do.OS) dt.do <- dunn.test(x= comb.Do.OS$Depth, g= comb.Do.OS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Do.OS) # More interesting is within a season, do you see differences between EOSO yrs? comb.Do.OS.sp <- comb.Do.OS[which(comb.Do.OS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.OS.sp) comb.Do.OS.wi <- comb.Do.OS[which(comb.Do.OS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.OS.wi) comb.Do.OS.su <- comb.Do.OS[which(comb.Do.OS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.OS.su) comb.Do.OS.fa <- comb.Do.OS[which(comb.Do.OS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Do.OS.fa) # 2. Distance from shore # Does max depth change as a function of dist.shore? kruskal.test(Depth~Dist.land, data= maxdepth_V50.Do) kruskal.test(Depth~Dist.land, data= maxdepth_L.Lim.PFD) # 3. TimeSeries #Does max depth change over time? By how much? bin <- seq(1984, 2019, by = 1) # good for all raw.data <- comb.Do.sp[which(comb.Do.sp$Loc== "NS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "D. opalescens", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Do.NS <- data # Optional data test # ggplot(mean.TS.Do.NS, aes(x= Year, y= Depth.m))+ # geom_point()+ # geom_errorbar(aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25) raw.data <- comb.Do.sp[which(comb.Do.sp$Loc== "OS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "D. opalescens", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Do.OS <- data # Optional data test # ggplot()+ # geom_point(data= mean.TS.Do.NS, aes(x= Year, y= Depth.m), color= "black")+ # geom_errorbar(data= mean.TS.Do.NS, aes(x=Year, ymin=Depth.m - Depth.se, ymax=Depth.m + Depth.se), width=0.25, color= "black")+ # geom_point(data= mean.TS.Do.OS, aes(x= Year, y= Depth.m), color= "blue")+ # geom_errorbar(data= mean.TS.Do.OS, aes(x=Year, ymin=Depth.m - Depth.se, ymax=Depth.m + Depth.se), width=0.25, color= "blue") # Stats comb.Do.sp.NS <- comb.Do.sp[which(comb.Do.sp$Loc == "NS"),] kruskal.test(Depth~Year, data= comb.Do.sp.NS) lmDo.NS = lm(Depth~Year, data = comb.Do.sp.NS) #Create the linear regression summary(lmDo.NS) #Review the results loDo.NS <- loess(data = comb.Do.sp.NS, formula = Depth~Year) summary(loDo.NS) comb.Do.sp.OS <- comb.Do.sp[which(comb.Do.sp$Loc == "OS"),] kruskal.test(Depth~Year, data= comb.Do.sp.OS) lmDo.OS = lm(Depth~Year, data = comb.Do.sp.OS) #Create the linear regression summary(lmDo.OS) #Review the results loDo.OS <- loess(data = comb.Do.sp.OS, formula = Depth~Year) summary(loDo.OS) # O bimaculatus ################################################################################### ###### Final plots for Obimac # Seasonal/ENSO diff # Final season plot ggplot(comb.Ob, aes(x=comb.Ob$Season, y=comb.Ob$Depth))+ geom_boxplot(fill= "grey")+ #geom_boxplot(data=comb.Ob, aes(x= comb.Ob$Season, y= comb.Ob$Depth, fill = comb.Ob$ENSO.Cat), alpha=0.7)+ #scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "fill")+ stat_summary(fun="mean", mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), geom="point", shape=20, size=6)+ stat_summary(fun.data = mean_se, fun.args = list(mult = 1), mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), width= 0.5, geom="errorbar", shape=20, size=1)+ scale_color_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "color")+ #stat_summary(fun="mean", mapping = aes(group = Season), position= position_dodge(width = .75), geom="point", shape=17, size=6, color="red", fill="red")+ scale_y_reverse(limits= c(390, 0), breaks=c(0, 50, 100, 150, 200, 250, 300, 350, 400))+ facet_wrap(facets = comb.Ob$Loc)+ ggtitle(expression(paste(italic("O. bimaculatus"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Season") ################################################################################### ###### Stats for final figs: # O. BIMAC # 1. Seasonal ENSO differences in max depth #Is it normal? ggplot(comb.Ob, aes(comb.Ob$Depth))+ geom_histogram() shapiro.test(comb.Ob$Depth) kruskal.test(Depth~Season, data= comb.Ob) # More interesting is within a season, do you see differences between ENSO yrs? comb.Ob.sp <- comb.Ob[which(comb.Ob$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.sp) comb.Ob.wi <- comb.Ob[which(comb.Ob$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.wi) comb.Ob.su <- comb.Ob[which(comb.Ob$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.su) comb.Ob.fa <- comb.Ob[which(comb.Ob$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.fa) ### Actual stats that take into account nearshore/offshore comb.Ob.NS <- comb.Ob[which(comb.Ob$Loc == "NS"),] comb.Ob.OS <- comb.Ob[which(comb.Ob$Loc == "OS"),] shapiro.test(comb.Ob.NS$Depth) shapiro.test(comb.Ob.OS$Depth) #NS kruskal.test(Depth~Season, data= comb.Ob.NS) dt.do <- dunn.test(x= comb.Ob.NS$Depth, g= comb.Ob.NS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Ob.NS) # More interesting is within a season, do you see differences between ENSO yrs? comb.Ob.NS.sp <- comb.Ob.NS[which(comb.Ob.NS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.NS.sp) comb.Ob.NS.wi <- comb.Ob.NS[which(comb.Ob.NS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.NS.wi) comb.Ob.NS.su <- comb.Ob.NS[which(comb.Ob.NS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.NS.su) comb.Ob.NS.fa <- comb.Ob.NS[which(comb.Ob.NS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.NS.fa) #OS kruskal.test(Depth~Season, data= comb.Ob.OS) dt.do <- dunn.test(x= comb.Ob.OS$Depth, g= comb.Ob.OS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Ob.OS) # More interesting is within a season, do you see differences between EOSO yrs? comb.Ob.OS.sp <- comb.Ob.OS[which(comb.Ob.OS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.OS.sp) comb.Ob.OS.wi <- comb.Ob.OS[which(comb.Ob.OS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.OS.wi) comb.Ob.OS.su <- comb.Ob.OS[which(comb.Ob.OS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.OS.su) comb.Ob.OS.fa <- comb.Ob.OS[which(comb.Ob.OS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Ob.OS.fa) # 2. Distance from shore # Does max depth change as a fcn of dist.shore kruskal.test(Depth~Dist.land, data= c.maxdepth_V50.Ob) kruskal.test(Depth~Dist.land, data= c.maxdepth_L.Lim.PFD) # 3. TimeSeries #Does max depth change over time? By how much? bin <- seq(1984, 2019, by = 1) # good for all raw.data <- comb.Ob.sp[which(comb.Ob.sp$Loc== "NS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "O. bimaculatus", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Ob.NS <- data # Data check # ggplot(mean.TS.Do.NS, aes(x= Year, y= Depth.m))+ # geom_point()+ # geom_errorbar(aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25) raw.data <- comb.Ob.sp[which(comb.Ob.sp$Loc== "OS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "O. bimaculatus", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Ob.OS <- data # data check # ggplot()+ # geom_point(data= mean.TS.Ob.NS, aes(x= Year, y= Depth.m), color= "black")+ # geom_errorbar(data= mean.TS.Ob.NS, aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25, color= "black")+ # geom_point(data= mean.TS.Ob.OS, aes(x= Year, y= Depth.m), color= "blue")+ # geom_errorbar(data= mean.TS.Ob.OS, aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25, color= "blue") # Stats comb.Ob.sp.NS <- comb.Ob.sp[which(comb.Ob.sp$Loc == "NS"),] kruskal.test(Depth~Year, data= comb.Ob.sp.NS) # linear model lmOb.NS = lm(Depth~Year, data = comb.Ob.sp.NS) #Create the linear regression summary(lmOb.NS) #Review the results summary(lmOb.NS)$r.squared # R sq value comb.Ob.sp.OS <- comb.Ob.sp[which(comb.Ob.sp$Loc == "OS"),] kruskal.test(Depth~Year, data= comb.Ob.sp.OS) # Kruskal-Wallis chi-squared = 116.06, df = 35, p-value = 1.279e-10 lmOb.OS = lm(Depth~Year, data = comb.Ob.sp.OS) #Create the linear regression summary(lmOb.OS) #Review the results mean(comb.Ob.sp.OS$Depth) # 192.3202 m # Determine anomaly comb.Ob.sp.OS$Anom <- (192.3 - comb.Ob.sp.OS$Depth) # ggplot(comb.Ob.sp.OS, aes(x= comb.Ob.sp.OS$Year, y= comb.Ob.sp.OS$Anom))+ # geom_point()+ # scale_y_reverse()+ # geom_hline(yintercept=0)+ # stat_smooth() ob.os <- loess(data = comb.Ob.sp.OS, formula = Depth~Year) summary(ob.os) nl.ob <- lm(Depth ~ poly(Year, 3, raw = TRUE), data = comb.Ob.sp.OS) summary(nl.ob) ggplot(comb.Ob.sp.OS, aes(x= comb.Ob.sp.OS$Year, y= comb.Ob.sp.OS$Depth))+ geom_point()+ scale_y_reverse()+ stat_smooth(method = lm, formula = y ~ poly(x, 3, raw = TRUE)) ################################################################################### ###### Final plots for M. gracilis # Seasonal/ENSO diff ggplot(comb.Cg, aes(x=comb.Cg$Season, y=comb.Cg$Depth))+ geom_boxplot(fill= "grey")+ #geom_boxplot(data=comb.Cg, aes(x= comb.Cg$Season, y= comb.Cg$Depth, fill = comb.Cg$ENSO.Cat), alpha=0.7)+ #scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "fill")+ stat_summary(fun="mean", mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), geom="point", shape=20, size=6)+ stat_summary(fun.data = mean_se, fun.args = list(mult = 1), mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), width= 0.5, geom="errorbar", shape=20, size=1)+ scale_color_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "color")+ #stat_summary(fun="mean", mapping = aes(group = Season), position= position_dodge(width = .75), geom="point", shape=17, size=6, color="red", fill="red")+ scale_y_reverse(limits= c(390, 0), breaks=c(0, 50, 100, 150, 200, 250, 300, 350, 400))+ facet_wrap(facets = comb.Cg$Loc)+ ggtitle(expression(paste(italic("M. gracilis"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Season") ###### Stats for final figs: # M gracilis # 1. Seasonal ENSO differences in max depth #Is it normal? ggplot(comb.Cg, aes(comb.Cg$Depth))+ geom_histogram() #looks normal shapiro.test(comb.Cg$Depth) shapiro.test(maxdepth_V50.Cg$Depth) kruskal.test(Depth~Season, data= comb.Cg) # More interesting is within a season, do you see differences between ENSO yrs? comb.Cg.sp <- comb.Cg[which(comb.Cg$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.sp) comb.Cg.wi <- comb.Cg[which(comb.Cg$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.wi) comb.Cg.su <- comb.Cg[which(comb.Cg$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.su) comb.Cg.fa <- comb.Cg[which(comb.Cg$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.fa) ### Actual stats that take into account nearshore/offshore comb.Cg.NS <- comb.Cg[which(comb.Cg$Loc == "NS"),] comb.Cg.OS <- comb.Cg[which(comb.Cg$Loc == "OS"),] shapiro.test(comb.Cg.NS$Depth) shapiro.test(comb.Cg.OS$Depth) #NS kruskal.test(Depth~Season, data= comb.Cg.NS) dt.do <- dunn.test(x= comb.Cg.NS$Depth, g= comb.Cg.NS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Cg.NS) # More interesting is within a season, do you see differences between ENSO yrs? comb.Cg.NS.sp <- comb.Cg.NS[which(comb.Cg.NS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.NS.sp) comb.Cg.NS.wi <- comb.Cg.NS[which(comb.Cg.NS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.NS.wi) comb.Cg.NS.su <- comb.Cg.NS[which(comb.Cg.NS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.NS.su) comb.Cg.NS.fa <- comb.Cg.NS[which(comb.Cg.NS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.NS.fa) #OS kruskal.test(Depth~Season, data= comb.Cg.OS) dt.do <- dunn.test(x= comb.Cg.OS$Depth, g= comb.Cg.OS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Cg.OS) # More interesting is within a season, do you see differences between EOSO yrs? comb.Cg.OS.sp <- comb.Cg.OS[which(comb.Cg.OS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.OS.sp) comb.Cg.OS.wi <- comb.Cg.OS[which(comb.Cg.OS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.OS.wi) comb.Cg.OS.su <- comb.Cg.OS[which(comb.Cg.OS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.OS.su) comb.Cg.OS.fa <- comb.Cg.OS[which(comb.Cg.OS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Cg.OS.fa) # 2. Distance from shore # Does max depth change as a function of dist.shore? kruskal.test(Depth~Dist.land, data= maxdepth_V50.Cg) kruskal.test(Depth~Dist.land, data= maxdepth_L.Lim.PFD) # 3. TimeSeries #Does max depth change over time? By how much? bin <- seq(1984, 2019, by = 1) # good for all raw.data <- comb.Cg.sp[which(comb.Cg.sp$Loc== "NS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "M. gracilis", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Cg.NS <- data # ggplot(mean.TS.Cg.NS, aes(x= Year, y= Depth.m))+ # geom_point()+ # geom_errorbar(aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25) raw.data <- comb.Cg.sp[which(comb.Cg.sp$Loc== "OS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "M. gracilis", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Cg.OS <- data # ggplot()+ # geom_point(data= mean.TS.Cg.NS, aes(x= Year, y= Depth.m), color= "black")+ # geom_errorbar(data= mean.TS.Cg.NS, aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25, color= "black")+ # geom_point(data= mean.TS.Cg.OS, aes(x= Year, y= Depth.m), color= "blue")+ # geom_errorbar(data= mean.TS.Cg.OS, aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25, color= "blue") # Stats comb.Cg.sp.NS <- comb.Cg.sp[which(comb.Cg.sp$Loc == "NS"),] kruskal.test(Depth~Year, data= comb.Cg.sp.NS) lmCg.NS = lm(Depth~Year, data = comb.Cg.sp.NS) #Create the linear regression summary(lmCg.NS) #Review the results comb.Cg.sp.OS <- comb.Cg.sp[which(comb.Cg.sp$Loc == "OS"),] kruskal.test(Depth~Year, data= comb.Cg.sp.OS) lmCg.OS = lm(Depth~Year, data = comb.Cg.sp.OS) #Create the linear regression summary(lmCg.OS) #Review the results ################################################################################### ###### Final plots for Plueroncodes # Seasonal/ENSO diff # Final plot ggplot(comb.Pp, aes(x=comb.Pp$Season, y=comb.Pp$Depth))+ geom_boxplot(fill= "grey")+ #geom_boxplot(data=comb.Pp, aes(x= comb.Pp$Season, y= comb.Pp$Depth, fill = comb.Pp$ENSO.Cat), alpha=0.7)+ #scale_fill_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "fill")+ stat_summary(fun="mean", mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), geom="point", shape=20, size=6)+ stat_summary(fun.data = mean_se, fun.args = list(mult = 1), mapping = aes(group = ENSO.Cat, color= ENSO.Cat), position= position_dodge(width = .75), width= 0.5, geom="errorbar", shape=20, size=1)+ scale_color_manual(name= "ENSO Index", values = c("#E69F00", "#0072B2", "white"), aesthetics = "color")+ #stat_summary(fun="mean", mapping = aes(group = Season), position= position_dodge(width = .75), geom="point", shape=17, size=6, color="red", fill="red")+ scale_y_reverse(limits= c(390, 0), breaks=c(0, 50, 100, 150, 200, 250, 300, 350, 400))+ facet_wrap(facets = comb.Pp$Loc)+ ggtitle(expression(paste(italic("P. planipes"))))+ theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ ylab("Maximum Depth of Luminoxyscape (m)")+ xlab("Season") ###### Stats for final figs: # P. planipes # 1. Seasonal ENSO differences in max depth #Is it normal? ggplot(comb.Pp, aes(comb.Pp$Depth))+ geom_histogram() #looks normal shapiro.test(comb.Pp$Depth) kruskal.test(Depth~Season, data= comb.Pp) # More interesting is within a season, do you see differences between ENSO yrs? comb.Pp.sp <- comb.Pp[which(comb.Pp$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.sp) comb.Pp.wi <- comb.Pp[which(comb.Pp$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.wi) comb.Pp.su <- comb.Pp[which(comb.Pp$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.su) comb.Pp.fa <- comb.Pp[which(comb.Pp$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.fa) ### Actual stats that take into account comb.Pp.NS <- comb.Pp[which(comb.Pp$Loc == "NS"),] comb.Pp.OS <- comb.Pp[which(comb.Pp$Loc == "OS"),] shapiro.test(comb.Pp.NS$Depth) shapiro.test(comb.Pp.OS$Depth) shapiro.test(comb.Pp$Depth) #NS kruskal.test(Depth~Season, data= comb.Pp.NS) dt.do <- dunn.test(x= comb.Pp.NS$Depth, g= comb.Pp.NS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Pp.NS) # More interesting is within a season, do you see differences between ENSO yrs? comb.Pp.NS.sp <- comb.Pp.NS[which(comb.Pp.NS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.NS.sp) comb.Pp.NS.wi <- comb.Pp.NS[which(comb.Pp.NS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.NS.wi) comb.Pp.NS.su <- comb.Pp.NS[which(comb.Pp.NS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.NS.su) comb.Pp.NS.fa <- comb.Pp.NS[which(comb.Pp.NS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.NS.fa) mean(comb.Pp.NS$Depth[which(comb.Pp.NS$Season== "Summer" & comb.Pp.NS$ENSO.Cat == "El Niño")]) sd(comb.Pp.NS$Depth[which(comb.Pp.NS$Season== "Summer" & comb.Pp.NS$ENSO.Cat == "El Niño")]) stderr(comb.Pp.NS$Depth[which(comb.Pp.NS$Season== "Summer" & comb.Pp.NS$ENSO.Cat == "El Niño")]) #OS kruskal.test(Depth~Season, data= comb.Pp.OS) dt.do <- dunn.test(x= comb.Pp.OS$Depth, g= comb.Pp.OS$Season, method = "bonferroni") kruskal.test(Depth~ENSO.Cat, data= comb.Pp.OS) # More interesting is within a season, do you see differences between EOSO yrs? comb.Pp.OS.sp <- comb.Pp.OS[which(comb.Pp.OS$Season== "Spring"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.OS.sp) comb.Pp.OS.wi <- comb.Pp.OS[which(comb.Pp.OS$Season== "Winter"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.OS.wi) comb.Pp.OS.su <- comb.Pp.OS[which(comb.Pp.OS$Season== "Summer"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.OS.su) comb.Pp.OS.fa <- comb.Pp.OS[which(comb.Pp.OS$Season== "Fall"),] kruskal.test(Depth~ENSO.Cat, data= comb.Pp.OS.fa) # 2. Distance from shore # Does max depth change as a fcn of dist.shore? kruskal.test(Depth~Dist.land, data= maxdepth_V50.Pp) kruskal.test(Depth~Dist.land, data= maxdepth_L.Lim.PFD) kruskal.test(Depth~Loc, data= maxdepth_V50.Pp) mean(maxdepth_V50.Pp$Depth[which(maxdepth_V50.Pp$Loc == "NS")]) mean(maxdepth_V50.Pp$Depth[which(maxdepth_V50.Pp$Loc == "OS")]) kruskal.test(Depth~Loc, data= maxdepth_V50.Do) mean(maxdepth_V50.Do$Depth[which(maxdepth_V50.Do$Loc == "NS")]) mean(maxdepth_V50.Do$Depth[which(maxdepth_V50.Do$Loc == "OS")]) kruskal.test(Depth~Loc, data= maxdepth_V50.Ob) mean(maxdepth_V50.Ob$Depth[which(maxdepth_V50.Ob$Loc == "NS")]) mean(maxdepth_V50.Ob$Depth[which(maxdepth_V50.Ob$Loc == "OS")]) kruskal.test(Depth~Loc, data= maxdepth_V50.Cg) mean(maxdepth_V50.Cg$Depth[which(maxdepth_V50.Cg$Loc == "NS")]) mean(maxdepth_V50.Cg$Depth[which(maxdepth_V50.Cg$Loc == "OS")]) # Data check # ggplot(maxdepth_V50.Pp, aes(x= maxdepth_V50.Pp$Loc, y= maxdepth_V50.Pp$Depth))+ # geom_boxplot() # # ggplot(maxdepth_V50.Do, aes(x= maxdepth_V50.Do$Loc, y= maxdepth_V50.Do$Depth))+ # geom_boxplot() # 3. TimeSeries #Does max depth change over time? By how much? bin <- seq(1984, 2019, by = 1) # good for all raw.data <- comb.Pp.sp[which(comb.Pp.sp$Loc== "NS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "P. planipes", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Pp.NS <- data # ggplot(mean.TS.Pp.NS, aes(x= Year, y= Depth.m))+ # geom_point()+ # geom_errorbar(aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25) raw.data <- comb.Pp.sp[which(comb.Pp.sp$Loc== "OS"),] data <- data.frame() #blank dataframe for (d in 1:(length(bin))) { i <- which(raw.data$Year == bin[d]) nd <- data.frame("Species" = "P. planipes", #makes a column with just the trial name "Season" = "Spring", "Year" = bin[d], "Depth.m" = mean(raw.data$Depth[i]), "Depth.sd" = sd(raw.data$Depth[i]), "Depth.se" = stderr(raw.data$Depth[i])) data <- rbind(data, nd) #replace data dataframe with new data sticking on } #end d mean.TS.Pp.OS <- data # ggplot()+ # geom_point(data= mean.TS.Pp.NS, aes(x= Year, y= Depth.m), color= "black")+ # geom_errorbar(data= mean.TS.Pp.NS, aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25, color= "black")+ # geom_point(data= mean.TS.Pp.OS, aes(x= Year, y= Depth.m), color= "blue")+ # geom_errorbar(data= mean.TS.Pp.OS, aes(x=Year, ymin=Depth.m - Depth.sd, ymax=Depth.m + Depth.sd), width=0.25, color= "blue") # Stats comb.Pp.sp.NS <- comb.Pp.sp[which(comb.Pp.sp$Loc == "NS"),] kruskal.test(Depth~Year, data= comb.Pp.sp.NS) lmPp.NS = lm(Depth~Year, data = comb.Pp.sp.NS) #Create the linear regression summary(lmPp.NS) #Review the results comb.Pp.sp.OS <- comb.Pp.sp[which(comb.Pp.sp$Loc == "OS"),] kruskal.test(Depth~Year, data= comb.Pp.sp.OS) lmPp.OS = lm(Depth~Year, data = comb.Pp.sp.OS) #Create the linear regression summary(lmPp.OS) #Review the results kruskal.test(Depth~ENSO.Cat, data= comb.Pp.sp.NS) kruskal.test(Depth~ENSO.Cat, data= comb.Pp.sp.OS) # ggplot(comb.Pp.sp, aes(x= comb.Pp.sp$ENSO.Cat, y= comb.Pp.sp$Depth))+ # geom_point()+ # stat_summary(fun="mean", mapping = aes(group = ENSO.Cat), position= position_dodge(width = .75), geom="point", shape=20, size=3, color="red", fill="red")+ # scale_y_reverse()+ # facet_wrap(facets = comb.Pp.sp$Loc) ############################################################ ########### Light Lim vs. o2 lim vs. both vs. no lim (Fig 1 McCormick et al.) ############################################# #### FINAL ANALYSIS FOR HISTOGRAMS # D opalescens SCB.sta$Lim.Do <- NA SCB.sta$Lim.Do[which(SCB.sta$Light.lim == FALSE & SCB.sta$Do.50 == TRUE)] <- "Irradiance Limited" SCB.sta$Lim.Do[which(SCB.sta$Light.lim == FALSE & SCB.sta$Do.50 == FALSE)] <- "Both Limited" SCB.sta$Lim.Do[which(SCB.sta$Light.lim == TRUE & SCB.sta$Do.50 == FALSE)] <- "O2 Limited" SCB.sta$Lim.Do[which(SCB.sta$Light.lim == TRUE & SCB.sta$Do.50 == TRUE)] <- "Not Limited" SCB.hist.Do <- SCB.sta[-c(which(is.na(SCB.sta$Lim.Do))),] SCB.hist.Do$Lim.Do <- as.factor(SCB.hist.Do$Lim.Do) SCB.hist.Do$Lim.Do <- ordered(SCB.hist.Do$Lim.Do, levels= c("Not Limited", "O2 Limited", "Irradiance Limited", "Both Limited")) # O. bimaculatus SCB.sta$Lim.Ob <- NA SCB.sta$Lim.Ob[which(SCB.sta$Light.lim == FALSE & SCB.sta$Ob.50 == TRUE)] <- "Irradiance Limited" SCB.sta$Lim.Ob[which(SCB.sta$Light.lim == FALSE & SCB.sta$Ob.50 == FALSE)] <- "Both Limited" SCB.sta$Lim.Ob[which(SCB.sta$Light.lim == TRUE & SCB.sta$Ob.50 == FALSE)] <- "O2 Limited" SCB.sta$Lim.Ob[which(SCB.sta$Light.lim == TRUE & SCB.sta$Ob.50 == TRUE)] <- "Not Limited" SCB.hist.Ob <- SCB.sta[-c(which(is.na(SCB.sta$Lim.Ob))),] SCB.hist.Ob$Lim.Ob <- as.factor(SCB.hist.Ob$Lim.Ob) SCB.hist.Ob$Lim.Ob <- ordered(SCB.hist.Ob$Lim.Ob, levels= c("Not Limited", "O2 Limited", "Irradiance Limited", "Both Limited")) # M. gracilis SCB.sta$Lim.Cg <- NA SCB.sta$Lim.Cg[which(SCB.sta$Light.lim == FALSE & SCB.sta$Cg.50 == TRUE)] <- "Irradiance Limited" SCB.sta$Lim.Cg[which(SCB.sta$Light.lim == FALSE & SCB.sta$Cg.50 == FALSE)] <- "Both Limited" SCB.sta$Lim.Cg[which(SCB.sta$Light.lim == TRUE & SCB.sta$Cg.50 == FALSE)] <- "O2 Limited" SCB.sta$Lim.Cg[which(SCB.sta$Light.lim == TRUE & SCB.sta$Cg.50 == TRUE)] <- "Not Limited" SCB.hist.Cg <- SCB.sta[-c(which(is.na(SCB.sta$Lim.Cg))),] SCB.hist.Cg$Lim.Cg <- as.factor(SCB.hist.Cg$Lim.Cg) SCB.hist.Cg$Lim.Cg <- ordered(SCB.hist.Cg$Lim.Cg, levels= c("Not Limited", "O2 Limited", "Irradiance Limited", "Both Limited")) # P. planipes SCB.sta$Lim.Pp <- NA SCB.sta$Lim.Pp[which(SCB.sta$Light.lim == FALSE & SCB.sta$Pp.50 == TRUE)] <- "Irradiance Limited" SCB.sta$Lim.Pp[which(SCB.sta$Light.lim == FALSE & SCB.sta$Pp.50 == FALSE)] <- "Both Limited" SCB.sta$Lim.Pp[which(SCB.sta$Light.lim == TRUE & SCB.sta$Pp.50 == FALSE)] <- "O2 Limited" SCB.sta$Lim.Pp[which(SCB.sta$Light.lim == TRUE & SCB.sta$Pp.50 == TRUE)] <- "Not Limited" SCB.hist.Pp <- SCB.sta[-c(which(is.na(SCB.sta$Lim.Pp))),] SCB.hist.Pp$Lim.Pp <- as.factor(SCB.hist.Pp$Lim.Pp) SCB.hist.Pp$Lim.Pp <- ordered(SCB.hist.Pp$Lim.Pp, levels= c("Not Limited", "O2 Limited", "Irradiance Limited", "Both Limited")) # Get counts count(SCB.hist.Do[which(SCB.hist.Do$Lim.Do == "Irradiance Limited"),]) count(SCB.hist.Do[which(SCB.hist.Do$Lim.Do == "O2 Limited"),]) count(SCB.hist.Do[which(SCB.hist.Do$Lim.Do == "Not Limited"),]) count(SCB.hist.Do[which(SCB.hist.Do$Lim.Do == "Both Limited"),]) count(SCB.hist.Ob[which(SCB.hist.Ob$Lim.Ob == "Irradiance Limited"),]) count(SCB.hist.Ob[which(SCB.hist.Ob$Lim.Ob == "O2 Limited"),]) count(SCB.hist.Ob[which(SCB.hist.Ob$Lim.Ob == "Not Limited"),]) count(SCB.hist.Ob[which(SCB.hist.Ob$Lim.Ob == "Both Limited"),]) count(SCB.hist.Cg[which(SCB.hist.Cg$Lim.Cg == "Irradiance Limited"),]) count(SCB.hist.Cg[which(SCB.hist.Cg$Lim.Cg == "O2 Limited"),]) count(SCB.hist.Cg[which(SCB.hist.Cg$Lim.Cg == "Not Limited"),]) count(SCB.hist.Cg[which(SCB.hist.Cg$Lim.Cg == "Both Limited"),]) count(SCB.hist.Pp[which(SCB.hist.Pp$Lim.Pp == "Irradiance Limited"),]) count(SCB.hist.Pp[which(SCB.hist.Pp$Lim.Pp == "O2 Limited"),]) count(SCB.hist.Pp[which(SCB.hist.Pp$Lim.Pp == "Not Limited"),]) count(SCB.hist.Pp[which(SCB.hist.Pp$Lim.Pp == "Both Limited"),]) # NEW DB for species % limitation: do.lim <- c(0.06, 50.63, 3.76, 45.55) ob.lim <- c(1.9, 17.94, 1.92, 78.24) cg.lim <- c(0.38, 37.67, 3.44, 58.51) pp.lim <- c(2.19, 15.14, 1.63, 81.04) spp.lim.sum <- rbind.data.frame(do.lim, ob.lim, cg.lim, pp.lim) colnames(spp.lim.sum) <- c("Irradiance Limited", "O2 Limited", "Both Limited", "Not Limited") spp.lim.sum$Spp <- NA spp.lim.sum$Spp <- c("D. opalescens", "O. bimaculatus", "M. gracilis", "P. planipes") do.lim <- rbind.data.frame(0.06, 50.63, 3.76, 45.55) colnames(do.lim) <- c("Percent") do.lim$Cat <- c("Irradiance Limited", "O2 Limited", "Both Limited", "Not Limited") do.lim$Spp <- "D. opalescens" ob.lim <- rbind.data.frame(1.9, 17.94, 1.92, 78.24) colnames(ob.lim) <- c("Percent") ob.lim$Cat <- c("Irradiance Limited", "O2 Limited", "Both Limited", "Not Limited") ob.lim$Spp <- "O. bimaculatus" cg.lim <- rbind.data.frame(0.38, 37.67, 3.44, 58.51) colnames(cg.lim) <- c("Percent") cg.lim$Cat <- c("Irradiance Limited", "O2 Limited", "Both Limited", "Not Limited") cg.lim$Spp <- "M. gracilis" pp.lim <- rbind.data.frame(2.19, 15.14, 1.63, 81.04) colnames(pp.lim) <- c("Percent") pp.lim$Cat <- c("Irradiance Limited", "O2 Limited", "Both Limited", "Not Limited") pp.lim$Spp <- "P. planipes" spp.lim.sum <- rbind.data.frame(do.lim, ob.lim, cg.lim, pp.lim) spp.lim.sum$Cat <- as.factor(spp.lim.sum$Cat) spp.lim.sum$Cat <- ordered(spp.lim.sum$Cat, levels= c("Not Limited", "O2 Limited", "Irradiance Limited", "Both Limited")) spp.lim.sum$Spp <- as.factor(spp.lim.sum$Spp) spp.lim.sum$Spp <- ordered(spp.lim.sum$Spp, levels= c("D. opalescens", "O. bimaculatus", "M. gracilis", "P. planipes")) ##### OKAY. NOW CHANGE EACH DEPTH BIN TO PERCENT OF EACH LIMITATION TYPE. # depth bins of 10m # for each depth bin, count each type of limitation / total number of points in that depth bin # then save, species, depth bin, fraction each limitation, total n in bin # D. opalescens test.df <- SCB.hist.Do bin <- seq(0, 770, by = 10) # good for all raw.data <- test.df data <- data.frame() #blank dataframe for (d in 1:(length(bin)-1)) { i <- which(raw.data$Depth>=bin[d] & raw.data$Depth=bin[d] & raw.data$Depth=bin[d] & raw.data$Depth=bin[d] & raw.data$Depth