#### Import environmental data, plot all data to inspect for abnormalities # ==> Temperature (T) # ==> Dissolved Oxygen (DO) # ==> pH (temp-adjusted) # ==> Salinity # ==> Tide height (estim. from http://tbone.biol.sc.edu/). library(reshape) library(plotly) #open plotly program package setwd("~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics") Env.Data <- data.frame(read.csv(file="~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/data/Environmental/EnvData-Master.csv", stringsAsFactors=F, header=T, na.strings = "")) # pH pH.Data <- Env.Data[,c(1,grep("pH\\.", colnames(Env.Data)))] names(pH.Data) <- c("DateTime", "WBE", "WBB", "SKE", "SKB", "PGE", "PGB", "CIE", "CIB", "FBE", "FBB") pH.Data.melted <- melt(pH.Data, id="DateTime") pH.Data.melted.noNA <- pH.Data.melted[which(!is.na(pH.Data.melted$value)),] pH.series <- plot_ly(data = pH.Data.melted.noNA, x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% #generate plotly plot layout(title="pH across sites, 2016 DNR outplant", yaxis = list(title = 'pH (total scale)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(pH.series), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-pH-series.html") pH.box <- plot_ly(data = pH.Data.melted.noNA, x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="pH across sites, 2016 DNR outplant", yaxis = list(title = 'pH (total scale)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(pH.box), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-pH-box.html") # Dissolved Oxygen DO.Data <- Env.Data[,c(1,grep("do\\.", colnames(Env.Data)))] DO.Data.noOuts <- DO.Data[which(DO.Data$FBE < 50),] #Remove FBE values over 50 names(DO.Data) <- c("DateTime", "WBE", "WBB", "SKE", "SKB", "PGE", "PGB", "CIE", "CIB", "FBE", "FBB") DO.Data.melted <- melt(DO.Data, id="DateTime") DO.Data.melted.noNA <- DO.Data.melted[which(!is.na(DO.Data.melted$value) & DO.Data.melted$value > 0),] DO.series <- plot_ly(data = DO.Data.melted.noNA, x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% layout(title="Dissolved Oxygen across sites, 2016 DNR outplant)", yaxis = list(title = 'DO (mg/L)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(DO.series), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-DO-series.html") DO.box <- plot_ly(data = DO.Data.melted.noNA, x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="Dissolved Oxygen across sites, 2016 DNR outplant)", yaxis = list(title = 'DO (mg/L)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(DO.box), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-DO-box.html") # Temperature T.Data <- Env.Data[,c(1,grep("doT", colnames(Env.Data)))] names(T.Data) <- c("DateTime", "WBE", "WBB", "SKE", "SKB", "PGE", "PGB", "CIE", "CIB", "FBE", "FBB") T.Data.melted <- melt(T.Data, id="DateTime") T.Data.melted.noNA <- T.Data.melted[which(!is.na(T.Data.melted$value)),] T.series <- plot_ly(data = T.Data.melted.noNA, x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% layout(title="Temperature (from DO sensor) across sites, 2016 DNR outplant", yaxis = list(title = 'Temperature (C)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(T.series), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Temp-series.html") T.box <- plot_ly(data = T.Data.melted.noNA, x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="Temperature (from DO sensor) across sites, 2016 DNR outplant", yaxis = list(title = 'Temperature (C)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(T.box), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Temp-box.html") # Salinity S.Data <- Env.Data[,c(1,grep("ctS", colnames(Env.Data)))] names(S.Data) <- c("DateTime", "CIB", "CIE", "FBB", "FBE", "PGE", "SKE", "SKB", "WBB", "WBE") S.Data.melted <- melt(S.Data, id="DateTime") S.Data.melted$value <- as.numeric(levels(S.Data.melted$value))[S.Data.melted$value] S.Data.melted.noNA <- S.Data.melted[which(!is.na(S.Data.melted$value)),] S.series <- plot_ly(data = S.Data.melted.noNA, x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% layout(title="Salinity across sites, 2016 DNR outplant", yaxis = list(title = 'Salinity'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(S.series), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Salinity-series.html") S.box <- plot_ly(data = S.Data.melted.noNA, x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="Salinity across sites, 2016 DNR outplant", yaxis = list(title = 'Salinity'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(S.box), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Salinity-box.html") # Tidal height Tide.Data <- Env.Data[,c(1,grep("Tide", colnames(Env.Data)))] names(Tide.Data) <- c("DateTime", "FB", "PG", "CI", "WB") Tide.Data.melted <- melt(Tide.Data, id="DateTime") # Tide.Data.melted$value <- as.numeric(levels(Tide.Data.melted$value))[Tide.Data.melted$value] Tide.Data.melted.noNA <- Tide.Data.melted[which(!is.na(Tide.Data.melted$value)),] Tide.series <- plot_ly(data = Tide.Data.melted.noNA, x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% layout(title="Tidal height across sites, 2016 DNR outplant", yaxis = list(title = 'Tidal Height'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(Tide.series), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Tide-series.html") Tide.box <- plot_ly(data = Tide.Data.melted.noNA, x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="Tidal height across sites, 2016 DNR outplant", yaxis = list(title = 'Tidal Height'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(Tide.box), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Tide-box.html") #### Remove erroneous data points # Sampling Dates: # ==> Case Inlet: July 19th; cut off env. data at -1.0 (estimate)- "never fully exposed" # ==> Willapa Bay: July 20th-22nd (?); cut off env. data at 1.5 (estimate)- "fully exposed/dry" # ==> Port Gamble: July 20th-22nd; cut off env. data at 0 (estimate)- "fully exposed/dry" # ==> Fidalgo Bay: July 20th-22nd; cut off env. data at -1.25 (estimate)- "never fully exposed" # ==> Skokomish: July 22nd # Create master melted dataframe with all environmental data # First, add another column "metric" to each env. data frame pH.Data.melted.noNA$metric <- c("pH") T.Data.melted.noNA$metric <- c("Temperature") DO.Data.melted.noNA$metric <- c("DO") S.Data.melted.noNA$metric <- c("Salinity") Tide.Data.melted.noNA$metric <- c("Tide") # Identify any time points in pH, DO & S that where probes were exposed. NOTE: likely need to update the threshold based on Micah/Alex's input. Tide.sites <- as.factor(c("FB", "PG", "CI", "WB")) Tide.exposed <- c(-1.25, 0, -1.0, 1.5) #tide depth at which probes are exposed (estimated) submerged.times = list() for (i in 1:length(Tide.sites)) { submerged.times[[i]] <- Tide.Data.melted.noNA[which(Tide.Data.melted.noNA$value>=Tide.exposed[i] & Tide.Data.melted.noNA$variable %in% Tide.sites[[i]]), ] } # Filter pH, DO & S at time points where probes were exposed Tide.location <- as.factor(c("FBB", "FBE", "PGB", "PGE", "CIB", "CIE", "WBB", "WBE")) J <- c(1,1,2,2,3,3,4,4) submerged.data = list() for (i in 1:length(Tide.location)) { submerged.data[[i]] <- rbind(pH.Data.melted.noNA[which(pH.Data.melted.noNA$DateTime %in% submerged.times[[J[i]]]$DateTime & pH.Data.melted.noNA$variable %in% Tide.location[[i]]), ], DO.Data.melted.noNA[which(DO.Data.melted.noNA$DateTime %in% submerged.times[[J[i]]]$DateTime & DO.Data.melted.noNA$variable %in% Tide.location[[i]]), ], S.Data.melted.noNA[which(S.Data.melted.noNA$DateTime %in% submerged.times[[J[i]]]$DateTime & S.Data.melted.noNA$variable %in% Tide.location[[i]]), ]) } # Combine results from the previous for loop for pH, DO & Salinity, and the Temp and Tide data (unedited) into one master env. data frame (long form) Env.Data.Master <- do.call(rbind, submerged.data) Env.Data.Master <- rbind(Env.Data.Master, T.Data.melted.noNA, Tide.Data.melted.noNA) #this is a dataframe with env. data, screened for times when pH, DO & S probes were exposed Env.Data.Master$metric <- as.factor(Env.Data.Master$metric) # Remove DO data from FBE after 6/24 @ 08:40:00, as the probe clearly malfunctioned after that time. Env.Data.Master <- subset(Env.Data.Master, !(variable=="FBE" & metric=="DO" & DateTime > "06/24/16 08:40:00")) # Remove Salinity data where probes clearly malfunctioned (identified via plots) Env.Data.Master <- subset(Env.Data.Master, !((variable=="CIE" & metric=="Salinity") | (variable=="FBB" & metric=="Salinity" & DateTime > "07/03/16 09:50:00") | (variable=="WBB" & metric=="Salinity" & DateTime > "06/25/16 05:30:00"))) # Let's plot the low-tide scrubbed salinity data. We should see that extreme low readings (~0) are gone Salinity.series.noLows <- plot_ly(data = subset(Env.Data.Master, metric=="Salinity"), x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% #generate plotly plot layout(title="Salinity across sites (low tides removed), 2016 DNR outplant", yaxis = list(title = 'Salinity (ppt)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(Salinity.series.noLows), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Salinity-series-noLowTides.html") Salinity.box.noLows <- plot_ly(data = subset(Env.Data.Master, metric=="Salinity"), x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="Salinity across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'Salinity (ppt)'), legend = list(x=.95, y=.2)) htmlwidgets::saveWidget(as_widget(Salinity.box.noLows), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Salinity-box-noLowTides.html") # Identify and remove outliers from pH, DO & Salinity data. Apply Tukey's method of removing outlying values, where values outside the inner fence removed: Env.Data.Master.noOuts <- Env.Data.Master # pH Data for(i in 1:length(Tide.location)) { #For individual site data IQR <- quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "pH" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[4] - quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "pH" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[2] upperBound <- as.numeric(quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "pH" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[4] + 1.5*IQR) #calculate upper bound lowerBound <- as.numeric(quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "pH" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[2] - 1.5*IQR) #calculate lower bound Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% Tide.location[i] & Env.Data.Master.noOuts$value > upperBound), "value"] <- NA #replace values higher than upper bound with NA Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% Tide.location[i] & Env.Data.Master.noOuts$value < lowerBound), "value"] <- NA #replace values lower than lower bound with NA } # DO Data for(i in 1:length(Tide.location)) { #For individual site data IQR <- quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "DO" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[4] - quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "DO" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[2] upperBound <- as.numeric(quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "DO" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[4] + 1.5*IQR) #replace values higher than upper bound with NA lowerBound <- 0 #DO cannot be less than 0 Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$variable %in% Tide.location[i] & Env.Data.Master.noOuts$value > upperBound), "value"] <- NA Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$variable %in% Tide.location[i] & Env.Data.Master.noOuts$value < lowerBound), "value"] <- NA } # Salinity Data for(i in 1:length(Tide.location)) { #For individual site data IQR <- quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "Salinity" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[4] - quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "Salinity" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[2] upperBound <- as.numeric(quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "Salinity" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[4] + 1.5*IQR) #calculate upper bound lowerBound <- as.numeric(quantile(Env.Data.Master[which(Env.Data.Master$metric %in% "Salinity" & Env.Data.Master$variable %in% Tide.location[i]),"value"], na.rm=TRUE)[2] - 1.5*IQR) #calculate lower bound Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$variable %in% Tide.location[i] & Env.Data.Master.noOuts$value > upperBound), "value"] <- NA #replace values higher than upper bound with NA Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$variable %in% Tide.location[i] & Env.Data.Master.noOuts$value < lowerBound), "value"] <- NA #replace values lower than lower bound with NA } # Remove the NA entries Env.Data.Master.noOuts <- Env.Data.Master.noOuts[which(!is.na(Env.Data.Master.noOuts$value)),] # Remove SK entries Env.Data.Master.noOuts <- subset(Env.Data.Master.noOuts, variable!="SKE") Env.Data.Master.noOuts <- subset(Env.Data.Master.noOuts, variable!="SKB") # Save this outlier-scrubbed dataset as .csv write.csv(file="~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/EnvData-Melted-NoOutliers.csv",Env.Data.Master.noOuts, col.names = T, row.names=F) # Let's plot the outlier-scrubbed data , to check it out after outliers have been removed. # ==> pH pH.series.noOuts <- plot_ly(data = subset(Env.Data.Master.noOuts, metric=="pH"), x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% #generate plotly plot layout(title="pH across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'pH (total scale)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(pH.series.noOuts), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-pH-series-noOutliers.html") pH.box.noOuts <- plot_ly(data = subset(Env.Data.Master.noOuts, metric=="pH"), x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="pH across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'pH (total scale)'), legend = list(x=.97, y=.98)) htmlwidgets::saveWidget(as_widget(pH.box.noOuts), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-pH-box-noOutliers.html") # Generate time-series of daily means, variance pH # convert to date/time and retain as a new field Env.Data.Master.noOuts$DateTime <- as.POSIXct(strptime(Env.Data.Master.noOuts$DateTime, format="%m/%d/%Y %H:%M:%S", tz=Sys.timezone())) # date in the format: YearMonthDay Hour:Minute library(dplyr) # use dplyr and mutate to add a day column to your data Env.Data.Master.noOuts_daily <- Env.Data.Master.noOuts %>% mutate(Day = as.Date(DateTime, format = "%Y-%m-%d")) #Env.Data.Master.noOuts-daily Env.Data.Master.noOuts_daily <- subset(Env.Data.Master.noOuts, !((metric=="Tide"))) %>% mutate(Day = as.Date(DateTime, format = "%Y-%m-%d")) %>% group_by(Day, variable, metric, Region, Site, Habitat) %>% # group by the day column summarise(daily.mean=mean(value), daily.sd=sd(value), daily.var=var(value)) %>% na.omit() group.colors <- c(WB = "sienna1", CI = "goldenrod1", PG ="steelblue2", FB = "royalblue3") ### JUST NEED TO SWAP EELGRASS & BARE LINES # Plot daily pH means ggplot(data=subset(Env.Data.Master.noOuts_daily, (metric=="pH" & (variable=="FBE" | variable=="FBB"))), aes(x=Day,y=daily.mean,colour=variable,group=variable)) + geom_line(size=2, aes(linetype=variable), color="royalblue3") + scale_linetype_manual(values=c("dashed", "solid")) + geom_ribbon(aes(ymax=daily.mean + daily.sd, ymin=daily.mean-daily.sd, alpha=0.5), colour=NA, fill = "grey70") + theme_light() + theme(plot.title = element_text(size=24, face="bold"), axis.text.y=element_text(size=18, angle=45, face="bold"), axis.text.x=element_text(size=18, face="bold"), axis.title=element_blank(), legend.position = "none", panel.background = element_blank()) + ggtitle("Fidalgo Bay pH") + ylim(6.84, 8.24) + xlim("0016-06-01"," 0016-07-19") ggplot(data=subset(Env.Data.Master.noOuts_daily, (metric=="pH" & (variable=="PGE" | variable=="PGB"))), aes(x=Day,y=daily.mean,colour=variable,group=variable)) + geom_line(size=2, aes(linetype=variable), color="steelblue2") + scale_linetype_manual(values=c("solid")) + geom_ribbon(aes(ymax=daily.mean + daily.sd, ymin=daily.mean-daily.sd, alpha=0.5), colour=NA, fill = "grey70") + theme_light() + theme(plot.title = element_text(size=24, face="bold"), axis.text.y=element_text(size=18, angle=45, face="bold"), axis.text.x=element_text(size=18, face="bold"), axis.title=element_blank(), legend.position = "none", panel.background = element_blank()) + ggtitle("Port Gamble Bay pH") + ylim(6.84, 8.24) + xlim("0016-06-01"," 0016-07-19") ggplot(data=subset(Env.Data.Master.noOuts_daily, (metric=="pH" & (variable=="CIE" | variable=="CIB"))), aes(x=Day,y=daily.mean,colour=variable,group=variable)) + geom_line(size=2, aes(linetype=variable), color="goldenrod1") + scale_linetype_manual(values=c("dashed", "solid")) + geom_ribbon(aes(ymax=daily.mean + daily.sd, ymin=daily.mean-daily.sd, alpha=0.5), colour=NA, fill = "grey70") + theme_light() + theme(plot.title = element_text(size=24, face="bold"), axis.text.y=element_text(size=18, angle=45, face="bold"), axis.text.x=element_text(size=18, face="bold"), axis.title=element_blank(), legend.position = "none", panel.background = element_blank()) + ggtitle("Case Inlet pH") + ylim(6.84, 8.24) + xlim("0016-06-01"," 0016-07-19") ggplot(data=subset(Env.Data.Master.noOuts_daily, (metric=="pH" & (variable=="WBE" | variable=="WBB"))), aes(x=Day,y=daily.mean,colour=variable,group=variable)) + geom_line(size=2, aes(linetype=variable), color="sienna1") + scale_linetype_manual(values=c("dashed", "solid")) + geom_ribbon(aes(ymax=daily.mean + daily.sd, ymin=daily.mean-daily.sd, alpha=0.5), colour=NA, fill = "grey70") + theme_light() + theme(plot.title = element_text(size=24, face="bold"), axis.text.y=element_text(size=18, angle=45, face="bold"), axis.text.x=element_text(size=18, face="bold"), axis.title=element_blank(), legend.position = "none", panel.background = element_blank()) + ggtitle("Willapa Bay pH") + ylim(6.84, 8.24) + xlim("0016-06-01"," 0016-07-19") # Try nested lm - anova results summary(lm(daily.mean ~ 1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) #grand pH mean anova(lm(daily.mean ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) summary(lm(daily.mean ~ Habitat-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) #mean daily mean pH in each habitat summary(lm(daily.mean ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) # mean daily mean pH at each location anova(lm(daily.mean ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "DO"),])) summary(lm(daily.mean ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "DO"),]))# mean daily mean DO anova(lm(daily.mean ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Temperature"),])) summary(lm(daily.mean ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Temperature"),])) # mean daliy mean temperature anova(lm(daily.mean ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Salinity"),])) summary(lm(daily.mean ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Salinity"),])) # mean daily mean salinity # daily means for bays (not separated by habitat) summary(lm(value ~ Site-1, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO"),]))# mean daily mean DO summary(lm(value ~ Site-1, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature"),])) # mean daliy mean temperature summary(lm(value ~ Site-1, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity"),])) # mean daily mean salinity anova(lm(daily.sd ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) summary(lm(daily.var ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) #mean daily sd pH anova(lm(daily.sd ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "DO"),])) summary(lm(daily.sd ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "DO"),])) #mean daily sd DO anova(lm(daily.sd ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Temperature"),])) summary(lm(daily.sd ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Temperature"),])) #mean daily sd Temp anova(lm(daily.sd ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Salinity"),])) summary(lm(daily.sd ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Salinity"),])) #mean daily sd Salinity anova(lm(daily.var ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) summary(lm(daily.var ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "pH"),])) #mean daily variance pH anova(lm(daily.var ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "DO"),])) summary(lm(daily.var ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "DO"),])) #mean daily variance DO anova(lm(daily.var ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Temperature"),])) summary(lm(daily.var ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Temperature"),])) #mean daily variance Temp anova(lm(daily.var ~ Site + Site/Habitat, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Salinity"),])) summary(lm(daily.var ~ variable-1, data=Env.Data.Master.noOuts_daily[which(Env.Data.Master.noOuts_daily$metric %in% "Salinity"),])) #mean daily variance Salinity # ==> Dissolved Oxygen # Let's plot the data again, to check it out after outliers have been removed. DO.series.noOuts <- plot_ly(data = subset(Env.Data.Master.noOuts, metric=="DO"), x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% #generate plotly plot layout(title="DO across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'DO (mg/L)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(DO.series.noOuts), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-DO-series-noOutliers.html") DO.box.noOuts <- plot_ly(data = subset(Env.Data.Master.noOuts, metric=="DO"), x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="DO across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'DO (total scale)'), legend = list(x=.95, y=.97)) htmlwidgets::saveWidget(as_widget(DO.box.noOuts), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-DO-box-noOutliers.html") # ==> Salinity # Let's plot the data again, to check it out after outliers have been removed. Salinity.series.noOuts <- plot_ly(data = subset(Env.Data.Master.noOuts, metric=="Salinity"), x = ~DateTime, y = ~value, type="scatter", mode="lines", color=~variable, hovertext=~value) %>% #generate plotly plot layout(title="Salinity across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'Salinity (ppt)'), legend = list(x=.95, y=.95)) htmlwidgets::saveWidget(as_widget(Salinity.series.noOuts), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Salinity-series-noOutliers.html") Salinity.box.noOuts <- plot_ly(data = subset(Env.Data.Master.noOuts, metric=="Salinity"), x = ~variable, y = ~value, type="box", color=~variable) %>% layout(title="Salinity across sites (outliers removed), 2016 DNR outplant", yaxis = list(title = 'Salinity (ppt)'), legend = list(x=.95, y=.2)) htmlwidgets::saveWidget(as_widget(Salinity.box.noOuts), "~/Documents/Roberts Lab/Paper-DNR-Geoduck-Proteomics/analyses/Environmental/June2016-Outplant-Salinity-box-noOutliers.html") # Run QQplots to assess normality on raw data par(mfrow = c(2, 3)) for (i in 1:length(Env.parameters)) { qqnorm(Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% Env.parameters[i]),"value"], main = Env.parameters[[i]], xlab = "Theoretical Quantiles", ylab = "Parameter Value", plot.it = TRUE) qqline(Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% Env.parameters[i]),"value"]) } # Add Habitat identifier to dataframe, since some environmental parameters did not differ between habitat within sites Env.Data.Master.noOuts$Habitat <- Env.Data.Master.noOuts$variable Env.Data.Master.noOuts$Habitat <- gsub("FBB|PGB|WBB|CIB", "Bare", Env.Data.Master.noOuts$Habitat) Env.Data.Master.noOuts$Habitat <- gsub("FBE|PGE|WBE|CIE", "Eelgrass", Env.Data.Master.noOuts$Habitat) Env.Data.Master.noOuts$Habitat <- gsub("FB|PG|WB|CI", "Nope", Env.Data.Master.noOuts$Habitat) Env.Data.Master.noOuts[Env.Data.Master.noOuts == "Nope"] <- NA Env.Data.Master.noOuts$Habitat <- as.factor(Env.Data.Master.noOuts$Habitat) # Add Site identifier to dataframe, since some environmental parameters did not differ between habitat within sites Env.Data.Master.noOuts$Site <- Env.Data.Master.noOuts$variable Env.Data.Master.noOuts$Site <- gsub("FBB|FBE|FB", "FB", Env.Data.Master.noOuts$Site) Env.Data.Master.noOuts$Site <- gsub("PGB|PGE|PG", "PG", Env.Data.Master.noOuts$Site) Env.Data.Master.noOuts$Site <- gsub("WBB|WBE|WB", "WB", Env.Data.Master.noOuts$Site) Env.Data.Master.noOuts$Site <- gsub("CIB|CIE|CI", "CI", Env.Data.Master.noOuts$Site) Env.Data.Master.noOuts$Site <- as.factor(Env.Data.Master.noOuts$Site) # Add North & South regions identifier to dataframe, since that's how the proteins differed Env.Data.Master.noOuts$Region <- Env.Data.Master.noOuts$variable Env.Data.Master.noOuts$Region <- gsub("FBB|FBE|PGB|PGE|FB|PG", "North", Env.Data.Master.noOuts$Region) Env.Data.Master.noOuts$Region <- gsub("WBB|WBE|CIB|CIE|WB|CI", "South", Env.Data.Master.noOuts$Region) Env.Data.Master.noOuts$Region <- as.factor(Env.Data.Master.noOuts$Region) # ENVIRONMENTAL DATA STATISTICS # 0. Assume unbalanced data sets across the board due to outlier scrubbing # 1. Assess equal variances via Bartlett Test of Homogeneity of Variances for each environmental parameter # PH bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "FB"),]) #unequal variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "CI"),]) #unequal variance (but not as sig.) bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "WB"),]) #unequal variance # DISSOLVED OXYGEN bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "FB"),]) #unequal variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "PG"),]) #unequal variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "CI"),]) #unequal variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "WB"),]) #unequal variance # TEMPERATURE bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "FB"),]) #EQUAL variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "PG"),]) #unequal variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "CI"),]) #EQUAL variance bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "WB"),]) #unequal variance # SALINITY bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$Site %in% "FB"),]) #EQUAL variance (marginal...) bartlett.test(value~variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$Site %in% "WB"),]) #unequal variance install.packages("FSA") library(FSA) # 2. RUN STATS TO COMPARE DATA BETWEEN HABITATS, WITHIN BAYS # STAT SELECTED BASED ON UNEQUAL VARIANCES, PARAMETRIC OR NOT # PH = normal distribution, unequal variance t.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "FB"),]$value, g=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "FB"),]$Habitat,var.equal=FALSE) t.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "CI"),]$value, g=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "FB"),]$Habitat,var.equal=FALSE) t.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "WB"),]$value, g=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$Site %in% "FB"),]$Habitat,var.equal=FALSE) # DISSOLVED OXYGEN: nonparametric, unequal variance = Kolmogorov-Smirnov Test ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "FB" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "FB" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "PG" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "PG" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "CI" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "CI" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "WB" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO" & Env.Data.Master.noOuts$Site %in% "WB" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) # TEMPERATURE: nonparametric, equal variance = krusgal-wallis; unequal=kolmogorov-smirnov kruskal.test(value ~ Habitat, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "FB"),]) ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "PG" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "PG" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) kruskal.test(value ~ Habitat, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "CI"),]) ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "WB" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature" & Env.Data.Master.noOuts$Site %in% "WB" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) # SALINITY kruskal.test(value ~ Habitat, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$Site %in% "FB"),]) ks.test(x=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$Site %in% "WB" & Env.Data.Master.noOuts$Habitat %in% "Bare"),]$value, y=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity" & Env.Data.Master.noOuts$Site %in% "WB" & Env.Data.Master.noOuts$Habitat %in% "Eelgrass"),]$value) #### If it's OK to assume equal variance ... # Run Stats on Environmental Data - if we can assume equal variance ... summary(aov(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH"),])) TukeyHSD(aov(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH"),])) kruskal.test(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO"),]) dunnTest(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO"),], method="bonferroni") kruskal.test(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature"),]) dunnTest(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature"),], method="bonferroni") kruskal.test(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity"),]) dunnTest(value ~ variable, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity"),], method="bonferroni") # Run stats to compare environmental data between bays, TukeyHSD(aov(value ~ Site, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH"),])) dunnTest(value ~ Site, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO"),], method="bonferroni") dunnTest(value ~ Site, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature"),], method="bonferroni") dunnTest(value ~ Site, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity"),], method="bonferroni") # Run stats to compare environmental data between north vs. south, TukeyHSD(aov(value ~ Region, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH"),])) dunnTest(value ~ Region, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO"),], method="bonferroni") dunnTest(value ~ Region, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature"),], method="bonferroni") dunnTest(value ~ Region, data=Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity"),], method="bonferroni") # PLAY AROUND with time series data # PH test <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "FBB"),]$value test1 <- ts(test, frequency=24*60/10) test2 <- decompose(test1) plot(test2) test3 <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "FBE"),]$value test4 <- ts(test3, frequency=24*60/10) test5 <- decompose(test4) plot(test5) test6 <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "PGB"),]$value test7 <- ts(test6, frequency=24*60/10) test8 <- decompose(test7) plot(test8) test9 <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "CIB"),]$value test10 <- ts(test9, frequency=24*60/10) test11 <- decompose(test10) plot(test11) test12 <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "CIE"),]$value test13<- ts(test12, frequency=24*60/10) test14 <- decompose(test13) plot(test14) test15 <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "WBB"),]$value test16 <- ts(test15, frequency=24*60/10) test17 <- decompose(test16) plot(test17) test18 <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH" & Env.Data.Master.noOuts$variable %in% "WBE"),]$value test19<- ts(test18, frequency=24*60/10) test20 <- decompose(test19) plot(test20) # CORRELATION PLOTS ! # Isolate dataframes per environmental variable pH <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "pH"),] DO <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "DO"),] Temperature <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Temperature"),] Salinity <- Env.Data.Master.noOuts[which(Env.Data.Master.noOuts$metric %in% "Salinity"),] # Make dataframe wide format pH.wide <- dcast(pH, DateTime ~ variable) DO.wide <- dcast(DO, DateTime ~ variable) Temperature.wide <- dcast(Temperature, DateTime ~ variable) Salinity.wide <- dcast(Salinity, DateTime ~ variable) # Correlation plots, between habitats within bays # pH Correlation Plots ggscatter(data=pH.wide, x="FBE", y="FBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="FBE pH", ylab="FBB pH", main="Fidalgo Bay pH Correlation Plot") ggscatter(data=pH.wide, x="WBE", y="WBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="WBE pH", ylab="WBB pH", main="Willapa Bay pH Correlation Plot") ggscatter(data=pH.wide, x="CIE", y="CIB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="CIE pH", ylab="CIB pH", main="Case Inlet pH Correlation Plot") # DO Correlation Plots ggscatter(data=DO.wide, x="FBE", y="FBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="FBE DO", ylab="FBB DO", main="Fidalgo Bay DO Correlation Plot") ggscatter(data=DO.wide, x="PGE", y="PGB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="PGE DO", ylab="PGB DO", main="Port Gamble Bay DO Correlation Plot") ggscatter(data=DO.wide, x="WBE", y="WBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="WBE DO", ylab="WBB DO", main="Willapa Bay DO Correlation Plot") ggscatter(data=DO.wide, x="CIE", y="CIB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="CIE DO", ylab="CIB DO", main="Case Inlet DO Correlation Plot") # Temperature Correlation Plots ggscatter(data=Temperature.wide, x="FBE", y="FBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="FBE Temperature", ylab="FBB Temperature", main="Fidalgo Bay Temperature Correlation Plot") ggscatter(data=Temperature.wide, x="PGE", y="PGB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="PGE Temperature", ylab="PGB Temperature", main="Port Gamble Bay Temperature Correlation Plot") ggscatter(data=Temperature.wide, x="WBE", y="WBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="WBE Temperature", ylab="WBB Temperature", main="Willapa Bay Temperature Correlation Plot") ggscatter(data=Temperature.wide, x="CIE", y="CIB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="CIE Temperature", ylab="CIB Temperature", main="Case Inlet Temperature Correlation Plot") # Salinity Correlation Plots ggscatter(data=Salinity.wide, x="FBE", y="FBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="FBE Salinity", ylab="FBB Salinity", main="Fidalgo Bay Salinity Correlation Plot") ggscatter(data=Salinity.wide, x="WBE", y="WBB", add="reg.line", conf.int = TRUE, cor.coef=TRUE, cor.method="pearson", xlab="WBE Salinity", ylab="WBB Salinity", main="Willapa Bay Salinity Correlation Plot") # Random stats # Survival statistics do.variance <- data.frame(rbind(c(22.5, 11.7, 7.2, 2.0), c(31.4, 13.3, 4.7, 2.3)), row.names = c("EELGRASS", "BARE"), stringsAsFactors = F) names(do.variance) <- c("FB","PG","CI","WB") chisq.test(t(do.variance), simulate.p.value = TRUE) #not significant ?!