---
title: "*Aphrocallistes vastus* Induced Flow"
author: "Evgeni Matveev"
date: '2020-05-04'
output:
  html_document:
    code_folding: hide
    theme: lumen
    toc: true
    toc_depth: 2
    number_sections: true
    toc_float:
      
      collapsed: false
      smooth_scroll: false
  pdf_document:
    toc: true
    toc_depth: 2
    number_sections: true
params:
  BothYearsAll: "BothYearsAll-v3.csv"
  ADCP2013_1:   "ADCP2013One.csv"
  ADCP2013_2:   "ADCP2013Two.csv"
  ADCP2014_1:   "ADCP2014One.csv"
  ADCP2014_2:   "ADCP2014Two.csv"
  Calibrations: "Calibrations.csv"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning = FALSE,message=FALSE)
```


# DATA PREPARATION 
This file includes the RStudio workflow. The data was collected from thermistors paired with oxygen sensors placed in the oscula of the sponge *Aphrocallistes vastus* at the Fraser Ridge in 2013 and 2014. Preliminary data management and analysis (i.e. the alignment of different records by timestamps, calibrations, and quality flag assignment) were done in Excel. 

## Upload packages and libraries
```{r Package Upload,message=FALSE}
#Delete the '#' if packages not yet installed on computer

#install.packages("ggplot2")
#install.packages("scales")
#install.packages("ggpubr")
#install.packages("tidyverse")
#install.packages("patchwork")
#install.packages("Hmisc")
#install.packages("changepoint")
#install.packages("kableExtra")
#install.packages("knitr")
#install.packages("formattable")
#install.packages("reshape2")
#install.packages("tinytex")
#tinytex::install_tinytex() #if this installation does not work, read section below in RStudio 
  

#UPLOAD LIBRARIES

#Graphing packages
library(ggplot2)
library(scales)
library(ggpubr)
library(patchwork)

#data tidying
library(dplyr)
library(tidyverse)
library(reshape2)

#data analysis
library(Hmisc)
library(changepoint)

#Rmarkdown assembly and formatting
library(rmarkdown)
library(knitr)
library(formattable)
library(kableExtra)
library(tinytex) #if this installation does not work, read section below in RStudio 
```

```{r TinyTex Troubleshoot,include=FALSE}

#Incase TinyTeX is not working (in my case, on a mac), you could try this script for installation, but without hashtags

#if (file.access("/usr/local/bin", mode = 2) == -1) {
#        tf <- tempfile()
#        cat("#!/bin/bash\nsudo chown -R `whoami`:admin /usr/local/bin\nRscript -e #'tinytex::install_tinytex()'\n", file = tf)
#        Sys.chmod(tf, mode = "0755")
#        system(paste0('osascript -e \'tell application "Terminal" to activate\' -e \'tell #application "Terminal" to do script "', tf, '"\''))
#      } else {
#        tinytex::install_tinytex()
#      }

```

## Data Upload
All calibrated thermistor and oxygen removal data is contained within the file 'BothYearsAll-v3.csv'. The file location for the raw data of each thermistor is in the table below. 

**Sponge Name**| **Instrument** | **File Name With Raw Data**
-------| -------
Av-1 | T89_2013 | Thermoxy_8-9_2013.xlsx
Av-2 | T1415_2013  | Thermoxy_14-15_2013.xlsx
Av-3 | T1213_2013 |Thermoxy_12-13_2013.xlsx
Av-4 | T89_2014 | Thermoxy_8-9_2014.xlsx
Av-5 | T1011_2014 | Thermoxy_10-11_2014.xlsx
Av-6 | T1011B_2014 | Thermoxy_10-11B_2014.xlsx

**Upload and Format the thermistor data:**

```{r Upload and Format Data, echo=TRUE}

#Upload Data from BothYearsAll-v3.csv
AllDataFull<-read.csv(params$BothYearsAll)

#Format timestamps to POSIXct
AllDataFull$Time<-as.POSIXct(AllDataFull$Time, format="%m/%d/%Y %H:%M:%S",tz="GMT")

#Reorder instruments by year
AllDataFull$Instrument <- factor(AllDataFull$Instrument, levels=c("T89_2013",
                                                          "T1415_2013",
                                                          "T1213_2013",
                                                          "T89_2014",
                                                          "T1011_2014",
                                                          "T1011B_2014"))

#Subset Only Periods where probe was in osculum (QF=S)
AllData<-AllDataFull%>%subset(QF=="S")

#Ambient rounded to 1 cm/s intervals for binning
AllData$AmBin<-round(AllData$Am,0)

#Adding rolling median (180= 30-min median)
rollby<-180

AllData<-AllData%>%
  group_by(Instrument)%>%
  mutate(Ex_roll=rollapply(Ex,rollby,mean,align='center',fill=NA),
         Am_roll=rollapply(Am,rollby,mean,align='center',fill=NA),
         Oxy_roll=rollapply(Oxy,rollby,mean,align='center',fill=NA))

#Adding Respiration Rate: (Excurrent (cm/s) * Osculum Area (cm2) * Oxygen Uptake (µmol/s))
AllData$RespirationRate<-((AllData$Ex*AllData$OscArea)*AllData$Oxy)/1000

# ADCP data upload and format
Aqd2013One<-read.csv(params$ADCP2013_1)
Aqd2013Two<-read.csv(params$ADCP2013_2)
Aqd2014One<-read.csv(params$ADCP2014_1)
Aqd2014Two<-read.csv(params$ADCP2014_2)

Aqd2013One$Time<-as.POSIXct(Aqd2013One$Time, format="%Y-%m-%d %H:%M:%S",tz="GMT")
Aqd2013Two$Time<-as.POSIXct(Aqd2013Two$Time, format="%Y-%m-%d %H:%M:%S",tz="GMT")
Aqd2014One$Time<-as.POSIXct(Aqd2014One$Time, format="%Y-%m-%d %H:%M:%S",tz="GMT")
Aqd2014Two$Time<-as.POSIXct(Aqd2014Two$Time, format="%Y-%m-%d %H:%M:%S",tz="GMT")

```


## Graph Formatting

*(To make formatting more condensed I made common theme formats to apply to graphs later on)*

```{r Graph Formatting}

#Flow Themes
FlowTheme<-list(
  geom_line(aes(y=Am),colour="#c8dcf0ff"),
  geom_line(aes(y=Ex),colour="#fad3b6"),
  geom_line(aes(y=Am_roll),colour="#2476C4"),
  geom_line(aes(y=Ex_roll),colour="#F1944B"),
  theme_classic(),
  scale_x_datetime(breaks = "2 hour",labels = date_format("%H:%M",tz="GMT")),
  labs(x="Time (hh:mm)",y=expression('Flow Rate (cm s'^-1*')')),
  theme(axis.text.x = element_text(angle=45,vjust=0.5,hjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent")))

LongFlowTheme<-list(
  geom_line(aes(y=Am),colour="#c8dcf0ff"),
  geom_line(aes(y=Ex),colour="#fad3b6"),
  geom_line(aes(y=Am_roll),colour="#2476C4"),
  geom_line(aes(y=Ex_roll),colour="#F1944B"),
  theme_classic(),
  scale_x_datetime(breaks = "5 hour",labels = date_format("%H:%M",tz="GMT")),
  labs(x="Time (hh:mm)",y=expression('Flow Rate (cm s'^-1*')')),
  theme(axis.text.x = element_text(colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent")))

ShortFlowTheme<-list(
  geom_line(aes(y=Am),colour="#c8dcf0ff"),
  geom_line(aes(y=Ex),colour="#fad3b6"),
  geom_line(aes(y=Am_roll),colour="#2476C4"),
  geom_line(aes(y=Ex_roll),colour="#F1944B"),
  theme_classic(),
  scale_x_datetime(breaks = "1 hour",labels = date_format("%H:%M",tz="GMT")),
  labs(x="Time (hh:mm)",y=expression('Flow Rate (cm s'^-1*')')),
  theme(axis.text.x = element_text(colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent")))

ArrestFlowTheme<-list(
  geom_line(aes(y=Am),colour="#2476C4"), 
  geom_line(aes(y=Ex),colour="#F1944B"),
  theme_classic(),
  scale_x_datetime(breaks = "15 mins",labels = date_format("%H:%M",tz="GMT")),
  labs(x="Time (dd/mm hh:mm)",y=expression('Flow Rate (cm s'^-1*')')),
  theme(axis.text.x = element_text(angle=45,vjust=0,hjust=0,colour="black",size=8,face="bold"),
        text = element_text(size=9,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))
)

#Oxygen Theme
OxyTheme <- list(
  geom_hline(yintercept=0,linetype=2,alpha=0.25),
  geom_line(aes(y=Oxy),colour="#c5eae2ff"),
  geom_line(aes(y=Oxy_roll),colour="#19AD8E"),
  labs(y=expression('Oxygen Removal \n (µmol O'[2]*'L'[p]*''^-1*')',sep="")),
  scale_x_datetime(breaks = "2 hour",labels = date_format("%H:%M",tz="GMT")),
  theme_classic(),
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.line.x=element_blank(),
        text=element_text(face="bold",size=7.5,colour="black"),
        axis.text.y=element_text(size=10,colour="black"),
        axis.line.y=element_line(size=1),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))
)

BackScatterTheme<-list(
    geom_line(aes(y=Backscatter),colour="black"),
    labs(y="Backscatter \n (counts)"),
    scale_x_datetime(breaks = "2 hour",labels = date_format("%H:%M",tz="GMT")),
    theme_classic(),
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.line.x=element_blank(),
          text=element_text(size=7.5,colour="black"),
          axis.text.y=element_text(size=7,colour="black"),
          axis.line.y=element_line(size=1),
          panel.background = element_rect(fill = "transparent"),
          plot.background = element_rect(fill = "transparent"))
)

#ADCP DATA FORMATTING

#Heatmap Format 
HeatmapFlowFormat<-list(scale_fill_gradientn(colours = rev(rainbow(5)),limits=c(0,70)),
                      theme_classic(),
                      theme(axis.text.x = element_text(colour="black",size=10,angle=90),
                            text = element_text(size=12,colour = "black",face="bold"),
                            line=element_line(size=1),
                            legend.title = element_text(size=9,vjust=1),
                            axis.text.y=element_text(size=10,colour="black"),
                            legend.position=c(0.8,0.9),
                            legend.direction = "horizontal",
                            panel.background = element_rect(fill = "transparent"),
                            plot.background = element_rect(fill = "transparent")))

HeatmapBackscatterFormat<-list(scale_fill_gradientn(colours = rev(rainbow(5))),
                      theme_classic(),
                      theme(axis.text.x = element_text(colour="black",size=10,angle=90),
                            text = element_text(size=12,colour = "black",face="bold"),
                            line=element_line(size=1),
                            legend.title = element_text(size=9,vjust=1),
                            axis.text.y=element_text(size=10,colour="black"),
                            legend.position=c(0.8,0.9),
                            legend.direction = "horizontal",
                            panel.background = element_rect(fill = "transparent"),
                            plot.background = element_rect(fill = "transparent")))

#Flow Distribution Format 
DistributionFormat<-list(geom_density(aes(colour=Profile,group=Profile),size=1),
  geom_density(size=1.3,colour="black"),
  theme_classic(),
  scale_colour_distiller(palette = "RdPu"),
  theme(axis.text.x = element_text(colour="black",size=10,angle=90),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        legend.position=c(0.8,0.8),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent")))

```

#**FIGURE 2 Ambient: AMBIENT FLOW CONDITIONS**

##Figure 2 Ambient A: *ADCP 2014-1 Flow Heatmap*

```{r Figure Ambient A}

SpongeHeights<-AllData%>%
                  group_by(Instrument)%>%
                  summarise(HeightAB=mean(HeightAB,na.rm = T)/100)

MeanSpongeHeight<-mean(SpongeHeights$HeightAB,na.rm=T)

HeatmapFlowAqd2014One<-ggplot(Aqd2014One,aes(x=Time,y=as.numeric(Profile)))+
  geom_tile(aes(fill=Flow))+
  geom_hline(yintercept=MeanSpongeHeight,linetype=2,col="darkgrey")+
  labs(x="Time",y="Height from Bottom (ft)",fill="Flow Speed (cm/s)",tag="A.")+
  scale_x_datetime(breaks = "5 hour",labels = date_format("%d-%m-%Y %H:%M",tz="GMT"))+
  HeatmapFlowFormat

#View Plot
HeatmapFlowAqd2014One


#Save Plot
ggsave("Figure Ambient A.pdf",
       path="./PDF Outputs",
       plot=HeatmapFlowAqd2014One,
       device="pdf",
       width=9,
       height=6,
       useDingbats=FALSE)


```

##Figure 2 Ambient B: *ADCP 2014-1 Backscatter Heatmap*

```{r Figure 2 Ambient B}

HeatmapBackscatterAqd2014One<-ggplot(Aqd2014One,aes(x=Time,y=as.numeric(Profile)))+
  geom_tile(aes(fill=Backscatter))+
  geom_hline(yintercept=MeanSpongeHeight,linetype=2,col="darkgrey")+
  labs(x="Time",y="Height from Bottom (ft)",fill="Backscatter (counts)",tag="B.")+
  scale_x_datetime(breaks = "5 hour",labels = date_format("%d-%m-%Y %H:%M",tz="GMT"))+
  HeatmapBackscatterFormat

#View Plot
HeatmapBackscatterAqd2014One

ggsave("Figure Ambient B.pdf",
       path="./PDF Outputs",
       plot=HeatmapBackscatterAqd2014One,
       width=9,
       height=6,
       useDingbats=FALSE)

ggplot(subset(Aqd2014One,Profile==0.5),aes(x=Time))+
  geom_line(aes(y=Flow))+
  geom_line(aes(y=Backscatter/10),col="red")


#When is the perfect correlation between backscatter and ambient? 
#NOTE: dayshift does not represent day, just that I took code I wrote for something else and was too lazy to replace. Thanks for reading. 
dayshift0<-Aqd2014One%>%
  subset(Profile==2)%>%
  select(Flow,Backscatter)

maxlength<-length(dayshift0$Flow)

lag<-c(0:200)
nn<-1
corcoeff<-0
dayshift1cor<- for(i in lag){
  corcoeff[nn]<-cor(dayshift0$Flow[1:(maxlength-i)],dayshift0$Backscatter[(1+i):maxlength])
  nn<-nn+1
  }

corbyshift<-data.frame(lag,corcoeff)

plot(corbyshift)

lagwithhighcor<-corbyshift[which.max(corbyshift$corcoeff),]
perflag<-lagwithhighcor$lag

paste("The best correlation between ambient and backscatter is",perflag*0.5,"minutes apart.")

```

##Figure 2 Ambient C: 2014 ADCP-1 flow by profile

```{r flow by profile, message=FALSE}
Aqd2014One_FlowByProfile<-Aqd2014One%>%
  summarise("0-1 mab median"     =median(Flow[Profile>=0 & Profile<=1]),
            "0-1 mab mean"     =mean(Flow[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Flow[Profile>=0 & Profile<=1]),
            "1-2 mab median"     =median(Flow[Profile>1 &  Profile<=2]),
            "1-2 mab mean"     =mean(Flow[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Flow[Profile>1 &  Profile<=2]),
            "2-3 mab median"     =median(Flow[Profile>2 &  Profile<=3]),
            "2-3 mab mean"     =mean(Flow[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Flow[Profile>2 &  Profile<=3]),
            "3-4 mab median"     =median(Flow[Profile>3 &  Profile<=4]),
            "3-4 mab mean"     =mean(Flow[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Flow[Profile>3 &  Profile<=4]),
            "4-5 mab median"     =median(Flow[Profile>4 &  Profile<=5]),
            "4-5 mab mean"     =mean(Flow[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Flow[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014One Flow",.before="0-1 mab median")

Aqd2014One_FlowByProfileB<-Aqd2014One%>%
  group_by(Profile)%>%
  summarise("medianflow"=median(Flow),"meanflow"=mean(Flow),FlowSD=sd(Flow))%>%
  add_column(Instrument="Aqd2014One Flow")


Aqd2014One_FlowByProfilePlot<-ggplot(Aqd2014One_FlowByProfileB,aes(x=medianflow,y=Profile))+
  geom_line(colour="black")+
  geom_point(colour="black")+
  geom_ribbon(aes(xmax=medianflow+FlowSD,xmin=medianflow-FlowSD),alpha=0.25,
              colour="NA",linetype=2,
              fill="grey")+
  theme_classic()+
  labs(x="Median Flow (cm/s)",y="Profile (mab)",tag="C.")+
  theme(legend.position = "none")+
  theme(axis.text.x = element_text(vjust=0.5,colour="black",size=10),
        text = element_text(size=15,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))

Aqd2014One_FlowByProfilePlot

ggsave("Figure Ambient C.pdf",
       path="./PDF Outputs",
       plot=Aqd2014One_FlowByProfilePlot,
       width=5,
       height=5,
       useDingbats=FALSE)



```

##Figure 2 Ambient D: 2014 ADCP-1 backscatter by profile
```{r Backscatter by profile, message=FALSE}
Aqd2014One_BackscatterByProfile<-Aqd2014One%>%
   summarise("0-1 mab mean"     =mean(Flow[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Flow[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Flow[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Flow[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Flow[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Flow[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Flow[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Flow[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Flow[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Flow[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014One Backscatter",.before="0-1 mab mean")

Aqd2014One_BackscatterByProfileB<-Aqd2014One%>%
  group_by(Profile)%>%
  summarise("meanbackscatter"=mean(Backscatter),BackscatterSD=sd(Backscatter))%>%
  add_column(Instrument="Aqd2014One Backscatter")

Aqd2014One_BackscatterByProfilePlot<-ggplot(Aqd2014One_BackscatterByProfileB,aes(x=meanbackscatter,y=Profile))+
  geom_line(colour="black")+
  geom_point(colour="black")+
  geom_ribbon(aes(xmax=meanbackscatter+BackscatterSD,xmin=meanbackscatter-BackscatterSD),alpha=0.25,
              colour="NA",linetype=2,
              fill="grey")+
  theme_classic()+
  labs(x="mean Backscatter (counts)",y="Profile (mab)",tag="D.")+
  theme(legend.position = c(0.8,0.8))+
  theme(axis.text.x = element_text(vjust=0.5,colour="black",size=10),
        text = element_text(size=15,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))

Aqd2014One_BackscatterByProfilePlot

ggsave("Figure Ambient D.pdf",
       path="./PDF Outputs",
       plot=Aqd2014One_BackscatterByProfilePlot,
       width=5,
       height=5,
       useDingbats=FALSE)

```


##Figure 2 Ambient E: *Ambient Thermistor Flow Distribution*

```{r Figure 2 Ambient E, warning=FALSE}

AmbientDistribution<-ggplot(AllData,aes(x=Am))+
  geom_density(col="grey",aes(group=Instrument),size=1)+
  geom_density(size=1)+
  theme_classic()+
  scale_color_brewer(palette = "Set2")+
  scale_x_continuous(breaks=seq(0,40,5))+
  labs(x="Ambient Flow Speed (cm/s)",y="Frequency")+
  theme(axis.text.x = element_text(vjust=0.5,colour="black",size=10),
        text = element_text(size=15,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))

AmbientDistribution

ggsave("Figure Ambient E.pdf",
       path="./PDF Outputs",
       plot=AmbientDistribution,
       width=6,
       height=6,
       useDingbats=FALSE)

#The code below shows that the most frequent ambient flow is at around 1 cm/s
AllAmDensity<-density(AllData$Am,na.rm=T)
AllAmDensityDf<-cbind(Ams=AllAmDensity$x,Freq=AllAmDensity$y)%>%as.data.frame()
MaxAmDensity<-subset(AllAmDensityDf,Freq==max(Freq))

paste("The most frequent ambient flow throughout all thermistors, with a frequency of",round(MaxAmDensity$Freq,2), ", was", round(MaxAmDensity$Ams,2),"cm/s")

```

##Aquadopp Temperatures 
```{r Temperatures, message= FALSE}

paste("the mean temperature from 2013 aquadopp 1 was ",round(mean(Aqd2013One$Temp,na.rm=T),2),"°C +/-",round(sd(Aqd2013One$Temp,na.rm=T),2),"SD")

paste("the mean temperature from 2013 aquadopp 2 was ",round(mean(Aqd2013Two$Temp,na.rm=T),2),"°C +/-",round(sd(Aqd2013Two$Temp,na.rm=T),2),"SD")

paste("the mean temperature from 2014 aquadopp 1 was ",round(mean(Aqd2014One$Temp,na.rm=T),2),"°C +/-",round(sd(Aqd2014One$Temp,na.rm=T),2),"SD")

paste("the mean temperature from 2014 aquadopp 2 was ",round(mean(Aqd2014Two$Temp,na.rm=T),2),"°C +/-",round(sd(Aqd2014Two$Temp,na.rm=T),2),"SD")

# Temperature Profiles, just in case 

Aqd2013OneTemp<-ggplot(subset(Aqd2013One,Temp!="NA"),aes(x=Time,y=Temp))+
  geom_line()

Aqd2013TwoTemp<-ggplot(subset(Aqd2013Two,Temp!="NA"),aes(x=Time,y=Temp))+
  geom_line()

Aqd2014OneTemp<-ggplot(subset(Aqd2014One,Temp!="NA"),aes(x=Time,y=Temp))+
  geom_line()

Aqd2014TwoTemp<-ggplot(subset(Aqd2014Two,Temp!="NA"),aes(x=Time,y=Temp))+
  geom_line()

```

## ADCP Ambient vs. Backscatter
```{r Binned ADCP vs. Backscatter}

BackscatterByAmBin<-Aqd2014One%>%
  subset(Profile<0.5)%>%
  group_by(Flow=round(Flow,0))%>%
  summarise(MedianBackscatter=median(Backscatter),
            BackscatterLower=quantile(na.omit(Backscatter),probs=0.25),
            BackscatterUpper=quantile(na.omit(Backscatter),probs=0.75),
            n=n())

BackscatterByAmBinPlot<-ggplot(subset(BackscatterByAmBin,n>50),aes(x=Flow,y=MedianBackscatter))+
  geom_ribbon(aes(ymin=BackscatterLower,ymax=BackscatterUpper),alpha=0.15,colour=NA)+
  geom_point(col="black")+
  geom_line( col="black")+
  theme_classic()+
  labs(x=expression('Ambient ADCP Flow Rate (cm s'^-1*')'),
       y="Backscatter (counts)")+
  theme(legend.position="none",
        strip.background = element_blank())

ggsave("ADCP AmBin vs. Backscatter.pdf",
       path="./PDF Outputs/Supplementary",
       plot=BackscatterByAmBinPlot,
       width=6,
       height=6,
       useDingbats=FALSE)

BackscatterByAmBinPlot

```

#**FIGURE 3 BinOxyEx: EXCURRENT AND OXYGEN REMOVAL BY AMBIENT BINS**

Figure 3 BinOxyEx summarizes the data from all the sponge excurrent and oxygen removal records into means plots binned by 1 cm/s ambient flow bins. 

###Figure 3 BinOxyEx *Sorting the data*

This groups the data by Instrument, and for each Ambient Bin gives the Median, Lower Quartile (25%), Upper Quartile (75%), Secondary Lower Quartile (15%), Secondary Upper Quartile (85%), and sample number (n).


```{r Sorting Ex and Oxy data By Abient Bins,warning=FALSE}

#Calculating upper and lower (0.75, 0.25) quantiles, and secondary upper and lower quantiles (0.85,0.15)
AllDataMedsErrs<-AllData%>%
  group_by(Instrument,AmBin)%>%
  summarise(MedianOxy=median(Oxy),
            OxyLower=quantile(na.omit(Oxy),probs=0.25),
            OxyUpper=quantile(na.omit(Oxy),probs=0.75),
            OxyLower2=quantile(na.omit(Oxy),probs=0.15),
            OxyUpper2=quantile(na.omit(Oxy),probs=0.85),
            MedianEx=median(Ex),
            ExLower=quantile(na.omit(Ex),probs=0.25),
            ExUpper=quantile(na.omit(Ex),probs=0.75),
            ExLower2=quantile(na.omit(Ex),probs=0.15),
            ExUpper2=quantile(na.omit(Ex),probs=0.85),
            MedianResp=median(RespirationRate),
            RespLower=quantile(na.omit(RespirationRate),probs=0.25),
            RespUpper=quantile(na.omit(RespirationRate),probs=0.75),
            RespLower2=quantile(na.omit(RespirationRate),probs=0.15),
            RespUpper2=quantile(na.omit(RespirationRate),probs=0.85),
            n=n())    

```

###Figure 3 BinOxyEx A: *Sorting excurrent speed by ambient bins*

Excurrent speed by ambient bins (note that only bins with more than 50 samples were used): 

```{r Excurrent By Am Bin Plot}

ExByBins<-ggplot(subset(AllDataMedsErrs,n>50),aes(x=AmBin,y=MedianEx))+
  geom_ribbon(aes(ymin=ExLower,ymax=ExUpper),alpha=0.15,colour=NA)+
  geom_line( col="black")+
  geom_point(col="black")+
  theme_classic()+
  scale_colour_brewer(palette = "Dark2")+
  scale_fill_brewer(  palette = "Dark2")+
  facet_wrap(vars(Instrument),scales = "free")+
  labs(x=expression('Ambient Flow Speed (cm s'^-1*')'),
       y=expression('Excurrent Flow Speed (cm s'^-1*')'),tag="A.")+
  theme(legend.position="none",
        strip.background = element_blank(),
        axis.title.x = element_blank())+
  guides(col = guide_legend(ncol = 8))

```

###Figure 3 BinOxyEx B: *Sorting oxygen removal by ambient bins*

Oxygen by ambient bins (note that only bins with more than 50 samples were used):

```{r Oxygen By Ambient Bin}

#Plotting oxygen by ambient bins
OxyRemByBins<-ggplot(subset(AllDataMedsErrs,n>50),aes(x=AmBin,y=MedianOxy))+
  geom_ribbon(aes(ymin=OxyLower,ymax=OxyUpper),alpha=0.15,colour=NA)+
  geom_line(colour="white")+
  geom_point(pch=21,fill="white")+
  theme_classic()+
  scale_colour_brewer(palette = "Dark2")+
  scale_fill_brewer(  palette = "Dark2")+
  facet_wrap(vars(Instrument),scales = "free")+
  labs(x=expression('Ambient Flow Speed (cm s'^-1*')'),
       y=expression('Oxygen Removal (µmol O'[2]*' L'[p]*''^-1*')'),tag="B.")+
  theme(legend.position="bottom",
        strip.background = element_blank())+
  guides(col = guide_legend(ncol = 8))

```

###Figure 3 BinOxyEx C?: *Sorting respiration rate by ambient bins*

Respiration rate by ambient bins (note that only bins with more than 50 samples were used):

```{r Oxygen Respiration By Ambient Bin}
RespByBins<-ggplot(subset(AllDataMedsErrs,n>50),aes(x=AmBin,y=MedianResp))+
  geom_ribbon(aes(ymin=RespLower,ymax=RespUpper),alpha=0.15,colour=NA)+
  geom_line(colour="darkgrey")+
  geom_point(pch=21,fill="darkgrey")+
  theme_classic()+
  scale_colour_brewer(palette = "Dark2")+
  scale_fill_brewer(  palette = "Dark2")+
  facet_wrap(vars(Instrument),scales = "free")+
  labs(x=expression('Ambient Flow Speed (cm s'^-1*')'),
       y=expression('Respiration Rate (µmol O'[2]*' s'^-1*')'))+
  theme(legend.position="bottom",
        strip.background = element_blank())+
  guides(col = guide_legend(ncol = 8))

ggsave("Respiration By Bins.pdf",
       path="./PDF Outputs/Supplementary",
       plot=RespByBins,
       width=8,
       height=6,
       useDingbats=FALSE)

RespByBins

```

##Figure 3 BinOxyEx: *Assembly*

Combining both A and B together using 'patchwork'
```{r Fig 3 BinOxyEx A and B, fig.width=8, fig.height=10}
#Put Plot Together
BinOxyEx<-(ExByBins/OxyRemByBins)+ plot_layout(widths = 8, heights = 10)

#Save Plot
ggsave("Figure 3 BinOxyEx.pdf",
       path="./PDF Outputs",
       plot=BinOxyEx,
       width=8,
       height=10,
       useDingbats=FALSE)

#View Plot
BinOxyEx
```

##Figure 3 BinOxyEx: *Description*

**Figure 3 BinOxyEx. Excurrent flow (cm/s) (A) and oxygen removal (µmol) (B) sorted into 1 cm/s Ambient flow bins**. Only bins where the number of observations was >50 were used. The points represent the median, the dark shading represents 25-75 percentiles, and the light shading represents 15-85 percentiles. **A)** In most cases, excurrent flow was not related to ambient flow. In T89_2014 the excurrent flow increases with ambient flow, but based on Fig 4 T89Sample this is more likely a behavioural response rather than passive flow. **B)** The oxygen removal (µmol) generally decreases with increase in ambient flow. This is true for all sponges except T1213_2013. 


#**FIGURE 4 T89Sample: T89_2014 SAMPLE RECORD**

## Figure 4 T89Sample A: *Full flow and oxygen record*
This is a sample plot of data collected from thermistors and oxygen sensors. This record was taken during 2014 from the sponge termed T8-9. 

### Figure 4 T89Sample A *Flow Plot*

```{r Flow Plot, warning=FALSE}

#Subsetting sponge T89_2014 for Figure 4 T89Sample
T89_2014<-subset(AllData,Instrument=="T89_2014")

#Graph of T8-9 Ambient (blue) and Excurrent (orange) Flow 
#(For FlowTheme details see Graph Formatting above)

T89Flow<-ggplot(T89_2014,aes(x=Time))+LongFlowTheme

#View Flow Plot
T89Flow

```

### Figure 4 T89Sample A *Oxygen Removal*

```{r Oxygen Plot, warning=FALSE}

#Oxygen plot. 
#(For OxyTheme details see Graph Formatting above)
T89OxyRem<-ggplot(T89_2014, aes(x=Time))+OxyTheme

#View oxygen plot 
T89OxyRem
```

### Figure 4 T89Sample A *Flow Categories*
Making  bars that represent the flow categories. Flow categories were determined based on changepoint analysis (see Supplementary for Details). Flow category cutoffs were based on the criteria in the table below, which were chosen based on a visual inspection of the flow graph.

**Flow Category**| **Ambient**  | **Excurrent**
-------| ------- | -------
Low    | < 2 cm/s | < 2 cm/s
High   | > 5 cm/s | > 4 cm/s

```{r Changepoint Code, include=FALSE}
#THIS CODE IS INCLUDED BELOW, IN THE SUPPLEMENARY, BUT HAS TO BE HERE SO THE FOLLOWING SECTION FUNCTIONS 

#Subset T8-9
T89_2014<-subset(AllData,Instrument=="T89_2014")

#Excurrent Changepoints with minimum segement length of 180 (30 minutes). 
T89_2014_ExChp<-cpt.meanvar(T89_2014$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T89_2014_MVEx<-as.data.frame(T89_2014_ExChp@param.est)  #mean variance
T89_2014_ChpLocEx<-as.data.frame(T89_2014_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T89_2014_AmChp<-cpt.meanvar(T89_2014$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T89_2014_MVAm<-as.data.frame(T89_2014_AmChp@param.est) #mean variance
T89_2014_ChpLocAm<-as.data.frame(T89_2014_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T89_2014_ChangepointAmData<-as.data.frame(c(T89_2014_ChpLocAm,T89_2014_MVAm))
T89_2014_ChangepointExData<-as.data.frame(c(T89_2014_ChpLocEx,T89_2014_MVEx))

#Rename columns so they're ex and am specific 
colnames(T89_2014_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T89_2014_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T89_2014_ChangepointExData$ChangeStart<-c(1,head(T89_2014_ChangepointExData$ExChp.cpts,-1))
T89_2014_ChangepointAmData$ChangeStart<-c(1,head(T89_2014_ChangepointAmData$AmChp.cpts,-1))
T89_2014$count<-seq.int(nrow(T89_2014))


#Aligning changepoint means/variances to original data
T89_2014_ChpCatEx<-merge(T89_2014,T89_2014_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T89_2014_ChpCat<-merge(T89_2014_ChpCatEx,T89_2014_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T89_2014_ChpCat<-fill_(T89_2014_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

#Calculating Percentile cutoffs (not used, but a method that was attempted)
#T89_2014ExHighCut<-T89_2014%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T89_2014ExLowCut <-T89_2014%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T89_2014AmHighCut<-T89_2014%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T89_2014AmLowCut <-T89_2014%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()

#CHANGING HIGH AND LOW CUT-OFFS CAN BE DONE HERE

T89_2014ExHighCut<-4
T89_2014ExLowCut <-2
T89_2014AmHighCut<-5
T89_2014AmLowCut <-2

#Adding Changepoint CategoriesT89_2014
T89_2014_ChpCat$ExCat <-  ifelse(T89_2014_ChpCat$ExMean <=  T89_2014ExLowCut, 'Low', 
                          ifelse(T89_2014_ChpCat$ExMean <=  T89_2014ExHighCut, 'Med', 'High'))
T89_2014_ChpCat$AmCat <-  ifelse(T89_2014_ChpCat$AmMean <=  T89_2014AmLowCut, 'Low', 
                          ifelse(T89_2014_ChpCat$AmMean <=  T89_2014AmHighCut, 'Med', 'High'))

#Remove areas with high variance (more than half of the mean)
T89_2014_ChpCat$ExCat<-ifelse(T89_2014_ChpCat$ExVar<(T89_2014_ChpCat$ExMean/2),
                              T89_2014_ChpCat$ExCat,'NA')
T89_2014_ChpCat$AmCat<-ifelse(T89_2014_ChpCat$AmVar<(T89_2014_ChpCat$AmMean/2),
                              T89_2014_ChpCat$AmCat,'NA')

#Make plots visualizing the changepoints and categories 

T89_2014_AmChangepoints<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T89_2014AmHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T89_2014AmLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="A.")

T89_2014_ExChangepoints<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T89_2014ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T89_2014ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category",tag="B.")

#Visualize Plots
T89_2014_ChangepointsPlot<-T89_2014_AmChangepoints/T89_2014_ExChangepoints


T89_2014_ChpCat%>%subset(AmCat!="NA")%>%ggplot(.,aes(y=Oxy,x=Ex,col=ExCat))+geom_point()+geom_smooth()+facet_wrap(~AmCat)

```

Adding Flow Category Bars: 

```{r Flow Category Bars}
ExcurrentHighCategory<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_rug(size=0.05,sides="b",length = unit(5, "npc"),aes(colour=ExCat))+
  scale_colour_manual(values=c("#F1944B","NA","NA"))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.line.x=element_blank(),
        text=element_blank(),
        axis.text.y=element_text(size=10,colour="black"),
        axis.line.y=element_line(size=1),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        legend.position = "none",
        plot.margin = unit(c(0,0,0,0), "cm"))

ExcurrentLowCategory<-ggplot(data=subset(T89_2014_ChpCat),aes(x=Time))+
  geom_rug(size=0.05,sides="b",length = unit(5, "npc"),aes(colour=ExCat))+
  scale_colour_manual(values=c("NA","#B65509","NA"))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.line.x=element_blank(),
        text=element_blank(),
        axis.text.y=element_text(size=10,colour="black"),
        axis.line.y=element_line(size=1),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        legend.position = "none",
        plot.margin = unit(c(0,0,0,0), "cm"))

AmbientHighCategory<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_rug(size=0.05,aes(colour=AmCat),length = unit(5, "npc"))+
  scale_colour_manual(values=c("#2476C4","NA","NA","NA"))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.line.x=element_blank(),
        text=element_blank(),
        axis.text.y=element_text(size=10,colour="black"),
        axis.line.y=element_line(size=1),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        legend.position = "none",
        plot.margin = unit(c(0,0,0,0), "cm"))

AmbientLowCategory<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_rug(size=0.05,aes(colour=AmCat),length = unit(5, "npc"))+
  scale_colour_manual(values=c("NA","#137393","NA","NA"))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.line.x=element_blank(),
        text=element_blank(),
        axis.text.y=element_text(size=10,colour="black"),
        axis.line.y=element_line(size=1),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        legend.position = "none",
        plot.margin = unit(c(0,0,0,0), "cm"))

#Put the bars together
Categories<-ggarrange(ExcurrentHighCategory,
                      AmbientHighCategory,
                      ExcurrentLowCategory,
                      AmbientLowCategory,
                      align="hv",
                      nrow=4,
                      heights=c(0.15,0.15,0.15,0.15))

```

###Figure 4 T89Sample A *Assembly* 

```{r Put It Together, message=FALSE, warning=FALSE}
#Combining all elements of Fig 4 T89Sample A 
T89Plot<-ggarrange(T89OxyRem,Categories,T89Flow,align="v",nrow=3,heights=c(1,0.4,2))

#Save Plot
ggsave("Figure 4 T89SampleA.pdf",
       path="./PDF Outputs",
       plot=T89Plot,
       width=30,
       height=20,
       units="cm",
       device="pdf")

#View full Fig 4 T89Sample A
T89Plot
```

##Figure 4 T89Sample B: *Probe placement photo*

Photo of thermistor placemet in the osculum of *A. vastus*. The thermistor T89_2014 is on the top of the photo (unused thermistor on the bottom). 


![Original photo of T89 (top instrument) in osculum](./R1808_00166.jpg)

**Placement Details**

* Yellow: Ambient probe (flow and oxygen)
* Red: Excurrent probe (flow and oxygen)
* Depth: 169m
* Coordinates: N49° 9.3905, W123° 22.9528
* Dive #: 1808
* Observation #: 00166


##Figure 4 T89Sample C: *Low ambient flow, high excurrent flow*
This is the period when ambient flow was low, but excurrent was high. Based on the changepoint bands in figure 4 T89Sample A we see that there is a period like this between 3 AM and 5 AM on Nov 1st. 

```{r Low Am High Ex, message=FALSE, warning=FALSE}
#Subseting period where the excurrent was consistently high, 
#and the ambient was consistently low (between 3 and 5 on the 1st of Nov)

T89ExHighAmLow<-T89_2014 %>% 
  filter(Time>= as.POSIXct('11/01/2014 3:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/01/2014 5:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

#Making the same graphs as Fig 4 T89SampleA, but only for the subsetted time
T89ExHighAmLowFlow<-ggplot(T89ExHighAmLow ,aes(x=Time))+ShortFlowTheme+lims(y=c(0,10))

#Same thing but for oxygen removal
T89ExHighAmLowOxyRem<-ggplot(T89ExHighAmLow, aes(x=Time))+OxyTheme+lims(y=c(0,5))

#Combining flow and oxygen removal graphs 
T89HighExLowAmPlot<-ggarrange(T89ExHighAmLowOxyRem,
                              T89ExHighAmLowFlow,
                              align="v",
                              nrow=2,
                              heights=c(1,2))

#Save Plot
ggsave("Figure T89_Sample C.pdf",
       path="./PDF Outputs", 
       plot=T89HighExLowAmPlot,
       width=10,
       height=20,
       units="cm",
       device="pdf")

#View plot
T89HighExLowAmPlot
```

##Figure 4 T89Sample D: *High ambient and excurrent*

Period when ambient flow and excurrent flow were consistently high, between 9AM and 11AM on Nov 2nd. 

```{r High Ambient High Excurrent, message=FALSE, warning=FALSE}
#Subestting a period where ambient was high AND excurrent was high 
# (between 9 and 11 on the 2nd of Nov, 2014)
T89ExHighAmHigh<-T89_2014 %>% 
  filter(Time>= as.POSIXct('11/02/2014 09:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/02/2014 11:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

#Same thing as above 
T89ExHighAmHighFlow<-ggplot(T89ExHighAmHigh ,aes(x=Time))+ShortFlowTheme+lims(y=c(0,10))

T89ExHighAmHighOxyRem<-ggplot(T89ExHighAmHigh, aes(x=Time))+OxyTheme+lims(y=c(0,5))

#Put flow and oxygen together
T89HighExHighAmPlot<-ggarrange(T89ExHighAmHighOxyRem,
                               T89ExHighAmHighFlow,
                               align="v",
                               nrow=2,
                               heights=c(1,2))

#Save Plot
ggsave("Figure 4 T89SampleD.pdf",path="./PDF Outputs",plot=T89HighExHighAmPlot,width=10,height=20,units="cm",device="pdf")

#View Plot 
T89HighExHighAmPlot
```

##Figure 4 T89Sample E: *Excurrent is independent of ambient*

This plot shows that the excurrent flow is not "passive":

  1. There are prolonged periods when the ambient flow decreases, while the excurrent flow increases. 
  2. There are periods prolonged where the ambient flow increases without affecting the excurrent flow. 

```{r Excurrent is Not Ambient, message=FALSE, warning=FALSE}
#Subset time between 22:00 on the 2nd, to 1:00 on the 3rd of Nov, 2014
T89ExIsNotAm<-T89_2014 %>% 
  filter(Time>= as.POSIXct('11/02/2014 22:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/03/2014 01:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

#Flow graph for that time period
T89ExIsNotAmFlow<-ggplot(T89ExIsNotAm ,aes(x=Time))+ShortFlowTheme+lims(y=c(0,10))

#Oxygen Graph for that time period
T89ExIsNotAmOxyRem<-ggplot(T89ExIsNotAm, aes(x=Time))+OxyTheme+lims(y=c(0,5))

#putting the graphs together
T89ExIsNotAmPlot<-ggarrange(T89ExIsNotAmOxyRem,
                            T89ExIsNotAmFlow,
                            align="v",
                            nrow=2,
                            heights=c(1,2))

#Save Plot
ggsave("Figure 4 T89SampleE.pdf",path="./PDF Outputs",
       plot=T89ExIsNotAmPlot,
       width=10,
       height=20,
       units="cm",
       device="pdf")

#View Plot 
T89ExIsNotAmPlot

```

##Figure 4 T89Sample F: *Boxplots of oxygen removal and excurrent at low/high ambient*

This boxplot plot shows that at the same excurrent flow, oxygen removal was on average lower when the ambient flow was high than when the ambient flow was low.

```{r Boxplots at same excurrent but different ambients, message=FALSE, warning=FALSE}
# 1.5 hours in each category 
HighExDiffAm<-T89_2014_ChpCat %>% 
  filter(Time>= as.POSIXct('11/01/2014 3:15:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/01/2014 4:45:00', format="%m/%d/%Y %H:%M:%S",tz="GMT")|
         Time>= as.POSIXct('11/02/2014 9:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/02/2014 10:30:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

HighExDiffAm$AmCat<-factor(HighExDiffAm$AmCat,levels=c("Low","High"))

# Adding N, mean, max, and min above categorical x-axis categories [not done, but included just in case]
# Troubleshoot: (make sure to manually select the full code for it to work)

stat_box_data <- function(y, upper_limit = median(y)) {
  return( 
    data.frame(
      y = 0.95 * upper_limit,
      label = paste('N =', length(y), ',',
                    'M =', round(mean(y), 1), ',', "\n",
                    'max=',round(max(y),1),',',
                    'min=',round(min(y),1))))}

#Oxygen Removal Boxplot plot 
#to include ALL changepoint segements, replace HighExDiffAm with T89_2014_ChpCat
T89OxyBoxplot<-ggplot(subset(T89_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                           aes(x=AmCat,y=Oxy))+
  geom_boxplot(fill="#19AD8E",colour="black")+
  stat_summary(fun=mean,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Oxygen Uptake (uMol/s)",tag="iii.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))


#Excurrent flow is not different
T89ExBoxplot<-ggplot(subset(T89_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                    aes(x=AmCat,y=Ex))+
  geom_boxplot(fill="#F1944B",colour="black",position="dodge2")+
  stat_summary(fun=mean,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Excurrent (cm/s)",tag="i.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+lims(y=c(-0.5,10))

#Ambient flow is different
T89AmBoxplot<-ggplot(subset(T89_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                    aes(x=AmCat,y=Am))+
  geom_boxplot(aes(y=Am),fill="#2476C4",colour="black")+
  stat_summary(fun=mean,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",tag="ii.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.line.y=element_blank(),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+lims(y=c(-0.5,10))

#Add to plot to add mean and sd. 
#stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+

#Add to plot to add statistics summary
#    stat_summary(
#    fun.data = stat_box_data, 
#    geom = "text", 
#    hjust = 0.25,
#    vjust = -1,angle=90)


#Assemble Plot
HighExBoxPlot<-ggarrange(T89ExBoxplot,T89AmBoxplot,T89OxyBoxplot,align="v",ncol=3,heights=c(1,1,1))

#Save Plot
ggsave("Figure 4 T89SampleF.pdf",
       path="./PDF Outputs",
       plot=HighExBoxPlot,
       width=15,
       height=10,
       units="cm",
       device="pdf")

#View Plot
HighExBoxPlot
```

###Table S1: T89 Box Plot Summary Statistics Table

```{r Table S1 Box Plot Table}
#Troubleshoot: might need to detatch(package:plyr)

#Time subset to 1.5 h
BoxplotSumStat<-HighExDiffAm%>%
  filter(row_number() %% 20 == 1)%>%
            group_by("Ambient Category"=AmCat)%>%
            summarise("Mean Ambient"=mean(Am),
                      "Ambient SD"=sd(Am),
                      "Mean Excurrent"=mean(Ex),
                      "Excurrent SD"=sd(Ex),
                      "Mean Oxygen"=mean(Oxy),
                      "Oxygen SD"= sd(Oxy),
                      n=n())

#JUST IN CASE: Only categories where excurrent is high and ambient is either low or high within the selected timeframes
BoxplotSumStatCI<-BoxplotSumStat%>%
  add_column("Oxygen CI"=round(qt(0.975,df=BoxplotSumStat$n-1)*
  BoxplotSumStat$`Oxygen SD`/sqrt(BoxplotSumStat$n),3),.after="Oxygen SD")%>%
  add_column("Excurrent CI"=round(qt(0.975,df=BoxplotSumStat$n-1)*
  BoxplotSumStat$`Excurrent SD`/sqrt(BoxplotSumStat$n),3),.after="Excurrent SD")%>%
  add_column("Ambient CI"=round(qt(0.975,df=BoxplotSumStat$n-1)*
  BoxplotSumStat$`Ambient SD`/sqrt(BoxplotSumStat$n),4),.after="Mean Ambient")


#The table below has all the periods where excurrent is high but ambient is high or low, and was the one used
BoxplotSumStatAll<-T89_2014_ChpCat%>%
  filter(row_number() %% 20 == 1)%>%
            subset((AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High") %>%
            group_by("Ambient Category"=AmCat)%>%
            summarise("Mean Ambient"=mean(Am),
                      "Ambient SD"=sd(Am),
                      "Mean Excurrent"=mean(Ex),
                      "Excurrent SD"=sd(Ex),
                      "Mean Oxygen"=mean(Oxy),
                      "Median Oxygen"=median(Oxy),
                      "Oxygen SD"= sd(Oxy),
                      n=n())

#Adding confidence intervals
BoxplotSumStatAllCI<-BoxplotSumStatAll%>%
  add_column("Oxygen CI"=round(qt(0.975,df=BoxplotSumStatAll$n-1)*
  BoxplotSumStatAll$`Oxygen SD`/sqrt(BoxplotSumStatAll$n),3),.after="Oxygen SD")%>%
  add_column("Excurrent CI"=round(qt(0.975,df=BoxplotSumStatAll$n-1)*
  BoxplotSumStatAll$`Excurrent SD`/sqrt(BoxplotSumStatAll$n),3),.after="Excurrent SD")%>%
  add_column("Ambient CI"=round(qt(0.975,df=BoxplotSumStatAll$n-1)*
  BoxplotSumStatAll$`Ambient SD`/sqrt(BoxplotSumStatAll$n),4),.after="Mean Ambient")

BoxplotTable<-kable_styling(kable(BoxplotSumStatAllCI),bootstrap_options = c("striped", "hover"))

BoxplotTable

write.csv(BoxplotSumStatAllCI,"./Tables/Table S1-T89-Boxplot.csv")

paste("Mean oxygen uptake under high ambient flow is",
      round((BoxplotSumStatAllCI$`Mean Oxygen`[1]/BoxplotSumStatAllCI$`Mean Oxygen`[2])*100,2),
      "% that of low ambient flow.")

#VISUALISATION OF CONFIDENCE INTERVALS 

#ggplot(BoxplotSumStatAllCI,aes(x=`Ambient Category`))+
#  geom_pointrange(aes(y=`Mean Excurrent`,ymax=`Mean Excurrent`+`Excurrent CI`,ymin=`Mean Excurrent`-`Excurrent #CI`))
#
#ggplot(BoxplotSumStatAllCI,aes(x=`Ambient Category`))+
#  geom_pointrange(aes(y=`Mean Ambient`,ymax=`Mean Ambient`+`Ambient CI`,ymin=`Mean Ambient`-`Ambient CI`))
#
#ggplot(BoxplotSumStatAllCI,aes(x=`Ambient Category`))+
#  geom_pointrange(aes(y=`Mean Oxygen`,ymax=`Mean Oxygen`+`Oxygen CI`,ymin=`Mean Oxygen`-`Oxygen CI`))
#
```

##Figure 4 T89Sample *Assembly*
Figure 4 T89Sample was assembled in Inkscape (v 1.0). 


##Figure 4 T89Sample *Caption*

**Figure 4 T89Sample. Sample record from a thermistor probe paired with oxygen sensors placed in the osculum of *Aphrocallistes vastus*. A.** Full 2 day 16 hr flow and oxygen removal record. Teal: oxygen removal (uMol), measured as the difference between ambient oxygen levels and oxygen levels in the sponge osculum. Orange: flow velocity (cm s^-1^) coming out of the sponge osculum. Blue: flow velocity (cm s^-1^) around the sponge. Bars: periods of consistent flow separated into low (Ex<2 cm s^-1^, Am<2 cm s^-1^) and high (Ex>4 cm s^-1^, Am>5 cm s^-1^) determined by a mean-variance changepoint analysis. The record contains several tidal cycles denoted as T1-4. **B.** A photograph showing the positioning of ambient (yellow) and excurrent (red) probes along with sponge location information. **C.** Period where ambient flow was consistently low, while excurrent flow was consistently high. **D.** Period where ambient and excurrent flow were both consistently high. Note the difference in oxygen removal between B and C (See also F) **E.** Although excurrent flow seems to follow ambient flow (A, Supp. Fig _), this is likely a behavioural response. There are periods where ambient flow increases while excurrent flow decreases (D-I), and periods where ambient flow decreases without affecting excurrent flow (D-II). **F.** Violin plot comparing oxygen removal (teal) under low and high ambient flows during periods of consistently high excurrent flow (orange). Black points represent mean, black bars represent standard deviation, and white points represent median. 



#**FIGURE 3: ARRESTS**

These two examples of arrests come from sponge T1011B_2014. There are further examples in the Supplementary section below. 

##Figure 3 *Arrest plots*

This code is selecting the 1011B arrests. See full flow records under Supplementary.

```{r Arrests 1011B}

#T10-11B 2014-Arrest1
T1011B_2014Arrest1<-subset(AllData,Instrument=="T1011B_2014") %>% 
  filter(Time>= as.POSIXct('11/06/2014 18:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/06/2014 20:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

#Arrest 1 plots
T1011B_2014Arrest1Flow<-ggplot(data=subset(T1011B_2014Arrest1) ,aes(x=Time))+
  ArrestFlowTheme+theme(text = element_text(size=15))

T1011B_2014Arrest1OxyRem<-ggplot(data=subset(T1011B_2014Arrest1) , aes(x=Time))+
  OxyTheme

T1011B_2014Arrest1BackScatter<-ggplot(data=subset(T1011B_2014Arrest1) , aes(x=Time))+
  BackScatterTheme+labs(tag="A.",title=element_blank())

#Arrest 1 plot assembly
T1011B_2014Arrest1Plot<-ggarrange(T1011B_2014Arrest1BackScatter,
                                  T1011B_2014Arrest1OxyRem,
                                  T1011B_2014Arrest1Flow,
                                  align="v",
                                  nrow=3,
                                  heights=c(0.6,0.6,2))

#Arrest2 subsetting time
T1011B_2014Arrest2<-subset(AllData,Instrument=="T1011B_2014") %>% 
  filter(Time>= as.POSIXct('11/06/2014 11:05:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & 
         Time<= as.POSIXct('11/06/2014 11:30:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

#Arrest2 plots 
T1011B_2014Arrest2Flow<-ggplot(data=subset(T1011B_2014Arrest2) ,aes(x=Time))+
  ArrestFlowTheme

T1011B_2014Arrest2OxyRem<-ggplot(data=subset(T1011B_2014Arrest2) , aes(x=Time))+
  OxyTheme

T1011B_2014Arrest2BackScatter<-ggplot(data=subset(T1011B_2014Arrest2) , aes(x=Time))+
  BackScatterTheme+labs(tag="B.",title=element_blank())

#Arrest 2 plot assembly
T1011B_2014Arrest2Plot<-ggarrange(T1011B_2014Arrest2BackScatter,
                                  T1011B_2014Arrest2OxyRem,
                                  T1011B_2014Arrest2Flow,
                                  align="v",
                                  nrow=3,
                                  heights=c(0.6,0.6,2))

#Visualize both plots
Figure3<-(T1011B_2014Arrest1Plot|T1011B_2014Arrest2Plot)+plot_layout(widths=3,heights=6)

Figure3



ggsave("Figure3.pdf",
       path="./PDF Outputs",
       plot=Figure3,
       width=6,
       height=6,
       useDingbats=FALSE)

ggsave("T1011B_2014Arrest1Flow.pdf",
       path="./PDF Outputs",
       plot=T1011B_2014Arrest1Flow,
       width=3,
       height=3,
       useDingbats=FALSE)


```

##Figure 3 *Caption*
**Figure 3. Examples of sponge arrests** Sponges show distinct behaviours not related to ambient flow. Here we show a clogging event (left) and a single arrest (right) as described in Grant et al. (2018) recorded from sponge T1011B_2014. It is suggested by Grant et al. (2018) that these are responses to increased sediment levels, which is possible in these records based on the backscatter record (black). More examples of arrrests can be seen in the supplementary. 




#**TABLE 1: SUMMARY OF EXCURRENT AND OXYGEN BY SPONGE**

This is a summary table of the length of recording, average excurrent (+/- CI), and average oxygen removal (+/- CI). A subsample rate of 1 in 20 samples (every 3.3 min) was taken based on a Partial Autocorrelation Function (PACF) analysis (see Supplementary).

Calulating a summary table with confidence intervals
```{r Summary Table with CI, warning=FALSE}

#One in every 20 data points is taken based on a partial autocorrelation analysis (supplementary)
SummaryTable1meanSD<-AllData%>%
  filter(row_number() %% 60 == 1)%>%
  subset(Instrument!="NA"&Ex>0)%>%
  group_by(Instrument)%>%
  summarise("Time Span (h)"= round(difftime(max(Time), min(Time), unit = "hours"),2),
            "Excurrent Mean (cm/s)"=round(mean(Ex),2),
            "Excurrent SD"=round(sd(Ex),2),
            "Oxygen Mean (µmol)"=round(mean(Oxy),2),
            "Oxygen SD"=round(sd(Oxy),2),
            "Respiration Mean"=round(mean(RespirationRate),2),
            "Respiration SD"=round(sd(RespirationRate),2),
             n=n())

#Adding Sponge Dimensions
SummaryTable1meanSDb<-SummaryTable1meanSD%>%
  add_column("Diameter(cm)"=c(3.6,2.6,3.5,4.8,8.3,7.6),.after="Instrument")%>%
  add_column("Cross Sectional Area (cm2)"=c(10.18,5.3,9.58,17.8,53.65,45.1),.after="Diameter(cm)")%>%
  add_column("Height Above Bottom (cm)"=c(33.1,38.1,24,39.5,54.1,33.7),
           .after="Cross Sectional Area (cm2)")

#Adding Oscular Flow Rate
SummaryTable1meanSDc<-SummaryTable1meanSDb%>%
  add_column("Oscular Flow Rate (L/s)"=round((SummaryTable1meanSDb$`Excurrent Mean (cm/s)`
                                             *SummaryTable1meanSDb$`Cross Sectional Area (cm2)`)
                                             /1000,3),.after="Excurrent SD")%>%
  add_column("Oscular Flow Rate SD (L/s)"=round((SummaryTable1meanSDb$`Excurrent SD`
                                             *SummaryTable1meanSDb$`Cross Sectional Area (cm2)`)/1000,3),.after="Oscular Flow Rate (L/s)")

#Adding 95% confidence intervals
SummaryTable1meanCI<-SummaryTable1meanSDc%>%
  group_by(Instrument)%>%
  add_column("Oxygen CI"=round(qt(0.975,df=SummaryTable1meanSDc$n-1)*
  SummaryTable1meanSDc$`Oxygen SD`/sqrt(SummaryTable1meanSDc$n),3),.after="Oxygen SD")%>%
  add_column("Excurrent CI"=round(qt(0.975,df=SummaryTable1meanSDc$n-1)*
  SummaryTable1meanSDc$`Excurrent SD`/sqrt(SummaryTable1meanSDc$n),3),.after="Excurrent SD")%>%
  add_column("Oscular Flow Rate CI"=round(qt(0.975,df=SummaryTable1meanSDc$n-1)*
  SummaryTable1meanSDc$`Oscular Flow Rate SD (L/s)`/sqrt(SummaryTable1meanSDc$n),4),.after="Oscular Flow Rate SD (L/s)")%>%
    add_column("Respiration Rate CI"=round(qt(0.975,df=SummaryTable1meanSDc$n-1)*
  SummaryTable1meanSDc$`Respiration SD`/sqrt(SummaryTable1meanSDc$n),4),.after="Respiration SD")

kable_styling(kable(SummaryTable1meanCI,align = "c"),bootstrap_options = c("striped", "hover"))

write.csv(SummaryTable1meanCI,"./Tables/Table 1-Summary.csv")

MeanExcurrent<-mean(SummaryTable1meanCI$`Excurrent Mean (cm/s)`)
ExcurrentCI<-round(qt(0.975,df=6-1)*
  sd(SummaryTable1meanCI$`Excurrent Mean (cm/s)`)/sqrt(6),3)

MeanOscFlowRate<-mean(SummaryTable1meanCI$`Oscular Flow Rate (L/s)`)
OscFlowRateCI<-round(qt(0.975,df=6-1)*
  sd(SummaryTable1meanCI$`Oscular Flow Rate (L/s)`)/sqrt(6),3)

MeanOxygen<-mean(SummaryTable1meanCI$`Oxygen Mean (µmol)`)
OxygenCI<-round(qt(0.975,df=6-1)*
  sd(SummaryTable1meanCI$`Oxygen Mean (µmol)`)/sqrt(6),3)

MeanRespirationRate<-mean(SummaryTable1meanCI$`Respiration Mean`)
RespirationCI<-round(qt(0.975,df=6-1)*
  sd(SummaryTable1meanCI$`Respiration Mean`)/sqrt(6),3)

paste("On average, A. vastus pumps at a rate of",round(MeanExcurrent,2)," cm/s ±",ExcurrentCI,"CI, processing on average",round(MeanOscFlowRate,3),"L/s ±",OscFlowRateCI,"CI, and removes on average",round(MeanOxygen,2),"µmol O2 /L ±",OxygenCI,"CI of Oxygen, with a respiration rate of ",MeanRespirationRate,"±",RespirationCI,"µmol Lp-1 s-1")

```
#**SUPPLEMENTARY**

##Figure 1 Supplementary: Binned Data Table

This is the data used to generate Figure 1. Specific code used to make the data is by Fig 1.

Upper: upper 25 percentile 
Lower: lower 25 percentile

The figure also includes upper and lower 15 percentiles. 

Note: any bins with an n<50 were not used for the figure, but are shown here. 

```{r Binned Data,warning=FALSE,message=FALSE}
AllBinOxyExData<-kable_styling(kable(AllDataMedsErrs),bootstrap_options = c("striped", "hover"))

AllDataMedsErrsB<-AllDataMedsErrs[c(1,2:5,8:10,13)]

BinOxyExDataSummarized<-kable_styling(kable(AllDataMedsErrsB),bootstrap_options = c("striped", "hover"))

BinOxyExDataSummarized

```

##Figure Supp. 1: All Flow and Oxygen Graphs

Thermistor 8-9 from 2014 was used in Figure 4 T89Sample as an example of the flow and oxygen data that suggests induced flow (excurrent increases with ambient, oxygen decreases). 

The flow and oxygen data from the other 5 instruments is shown below. 

###*Code*: 

```{r Supp. Flow and Oxygen Code, warning=FALSE, message=FALSE}

#Thermistor 8-9, 2013
T89_2013Flow<-ggplot(subset(AllData,Instrument=="T89_2013"),aes(x=Time))+
  FlowTheme+
  theme(axis.title.y = element_text(size=8),axis.text.y = element_text(size=6))

T89_2013Oxygen<-ggplot(subset(AllData,Instrument=="T89_2013"),aes(x=Time))+
  OxyTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T89_2013Backscatter<-ggplot(subset(AllData,Instrument=="T89_2013"),aes(x=Time))+
  BackScatterTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T89_2013Full<-T89_2013Backscatter/T89_2013Oxygen/T89_2013Flow+
  plot_layout(heights = c(8,8,24))+
  plot_annotation(title="A. T89_2013")

ggsave("T89_2013Full.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T89_2013Full,
       width=8,
       height=6,
       useDingbats=FALSE)

#Thermistor 14-15, 2013
T1415_2013Flow<-ggplot(subset(AllData,Instrument=="T1415_2013"),aes(x=Time))+
  FlowTheme+
  theme(axis.title.y = element_text(size=8),axis.text.y = element_text(size=6))

T1415_2013Oxygen<-ggplot(subset(AllData,Instrument=="T1415_2013"),aes(x=Time))+
  OxyTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1415_2013Backscatter<-ggplot(subset(AllData,Instrument=="T1415_2013"),aes(x=Time))+
  BackScatterTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1415_2013Full<-T1415_2013Backscatter/T1415_2013Oxygen/T1415_2013Flow+
  plot_layout(heights = c(8,8,24))+
  plot_annotation(title="B. T1415_2013")

ggsave("T1415_2013Full.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T1415_2013Full,
       width=8,
       height=6,
       useDingbats=FALSE)


#Thermistor 12-13, 2013
T1213_2013Flow<-ggplot(subset(AllData,Instrument=="T1213_2013"),aes(x=Time))+
  FlowTheme+
  theme(axis.title.y = element_text(size=8),axis.text.y = element_text(size=6))

T1213_2013Oxygen<-ggplot(subset(AllData,Instrument=="T1213_2013"),aes(x=Time))+
  OxyTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1213_2013Backscatter<-ggplot(subset(AllData,Instrument=="T1213_2013"),aes(x=Time))+
  BackScatterTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1213_2013Full<-T1213_2013Backscatter/T1213_2013Oxygen/T1213_2013Flow+
  plot_layout(heights = c(8,8,24))+
  plot_annotation(title="C. T1213_2013")

ggsave("T1213_2013Full.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T1213_2013Full,
       width=8,
       height=6,
       useDingbats=FALSE)


#Thermistor 10-11, 2014
T1011_2014Flow<-ggplot(subset(AllData,Instrument=="T1011_2014"),aes(x=Time))+
  FlowTheme+
  theme(axis.title.y = element_text(size=8),axis.text.y = element_text(size=6))

T1011_2014Oxygen<-ggplot(subset(AllData,Instrument=="T1011_2014"),aes(x=Time))+
  OxyTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1011_2014Backscatter<-ggplot(subset(AllData,Instrument=="T1011_2014"),aes(x=Time))+
  BackScatterTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1011_2014Full<-T1011_2014Backscatter/T1011_2014Oxygen/T1011_2014Flow+
  plot_layout(heights = c(8,8,24))+
  plot_annotation(title="D. T1011_2014")

ggsave("T1011_2014Full.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T1011_2014Full,
       width=8,
       height=6,
       useDingbats=FALSE)

#Thermistor 10-11, 2014
T1011B_2014Flow<-ggplot(subset(AllData,Instrument=="T1011B_2014"),aes(x=Time))+
  FlowTheme+
  theme(axis.title.y = element_text(size=8),axis.text.y = element_text(size=6))

T1011B_2014Oxygen<-ggplot(subset(AllData,Instrument=="T1011B_2014"),aes(x=Time))+
  OxyTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1011B_2014Backscatter<-ggplot(subset(AllData,Instrument=="T1011B_2014"),aes(x=Time))+
  BackScatterTheme+
  theme(axis.title.y = element_text(size=6),axis.text.y = element_text(size=6))

T1011B_2014Full<-T1011B_2014Backscatter/T1011B_2014Oxygen/T1011B_2014Flow+
  plot_layout(heights = c(8,8,24))+
  plot_annotation(title="E. T1011B_2014")

ggsave("T1011B_2014Full.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T1011B_2014Full,
       width=8,
       height=6,
       useDingbats=FALSE)

```

###*Summary of Graphs*

```{r Supp. Flow and Oxygen Summary Graphs,  warning=FALSE, message=FALSE}
AllFloxy<-list(T89_2013Full,T1415_2013Full,T1213_2013Full,T1011_2014Full,T1011B_2014Full)

#Save all to a single PDF
pdf("./PDF Outputs/Supplementary/All Flow and Oxygen.pdf",useDingbats = F)
AllFloxy
dev.off()


```

##Figure Supp. 4: Arrest Plots

Sponges behave by arresting their filtration (see Grant et al, 2018). A few arrests are shown as examples below.

```{r Arrest Plots, fig.width=3,  warning=FALSE, message=FALSE}
#Arrests

#T89_2013: Arrest 1
T89_2013Arrest1<-subset(AllData,Instrument=="T89_2013") %>% filter(Time>= as.POSIXct('10/11/2013 1:45:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('10/11/2013 2:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T89_2013Arrest1Flow<-       ggplot(data=subset(T89_2013Arrest1) ,aes(x=Time))+ArrestFlowTheme

T89_2013Arrest1OxyRem<-     ggplot(data=subset(T89_2013Arrest1) , aes(x=Time))+OxyTheme

T89_2013Arrest1BackScatter<-ggplot(data=subset(T89_2013Arrest1) , aes(x=Time))+
  BackScatterTheme+labs(title="T89_2013")

T89_2013Arrest1Plot<-ggarrange(T89_2013Arrest1BackScatter,T89_2013Arrest1OxyRem,T89_2013Arrest1Flow,align="v",nrow=3,heights=c(0.5,0.5,2))


#T1415_2013: Arrest 1

T1415_2013Arrest1<-subset(AllData,Instrument=="T1415_2013") %>% filter(Time>= as.POSIXct('10/11/2013 19:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('10/11/2013 20:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T1415_2013Arrest1Flow<-       ggplot(data=subset(T1415_2013Arrest1) ,aes(x=Time))+ArrestFlowTheme

T1415_2013Arrest1OxyRem<-     ggplot(data=subset(T1415_2013Arrest1) , aes(x=Time))+OxyTheme

T1415_2013Arrest1BackScatter<-ggplot(data=subset(T1415_2013Arrest1) , aes(x=Time))+
  BackScatterTheme+labs(title="T1415_2013")

T1415_2013Arrest1Plot<-ggarrange(T1415_2013Arrest1BackScatter,T1415_2013Arrest1OxyRem,T1415_2013Arrest1Flow,align="v",nrow=3,heights=c(0.5,0.5,2))


#T1213_2013: Arrest1

T1213_2013Arrest1<-subset(AllData,Instrument=="T1213_2013") %>% filter(Time>= as.POSIXct('10/12/2013 5:20:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('10/12/2013 6:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T1213_2013Arrest1Flow<-       ggplot(data=subset(T1213_2013Arrest1) ,aes(x=Time))+ArrestFlowTheme

T1213_2013Arrest1OxyRem<-     ggplot(data=subset(T1213_2013Arrest1) , aes(x=Time))+OxyTheme

T1213_2013Arrest1BackScatter<-ggplot(data=subset(T1213_2013Arrest1) , aes(x=Time))+
  BackScatterTheme+labs(title="T1213_2013")

T1213_2013Arrest1Plot<-ggarrange(T1213_2013Arrest1BackScatter,
                                 T1213_2013Arrest1OxyRem,
                                 T1213_2013Arrest1Flow,align="v",nrow=3,heights=c(0.5,0.5,2))


#T10-11 2014: Arrest 1

T1011_2014Arrest1<-subset(AllData,Instrument=="T1011_2014") %>% filter(Time>= as.POSIXct('11/03/2014 19:45:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('11/03/2014 21:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T1011_2014Arrest1Flow<-ggplot(data=subset(T1011_2014Arrest1) ,aes(x=Time))+ArrestFlowTheme

T1011_2014Arrest1OxyRem<-ggplot(data=subset(T1011_2014Arrest1) , aes(x=Time))+OxyTheme

T1011_2014Arrest1BackScatter<-ggplot(data=subset(T1011_2014Arrest1) , aes(x=Time))+
  BackScatterTheme+labs(title="T1011_2014")

T1011_2014Arrest1Plot<-ggarrange(T1011_2014Arrest1BackScatter,T1011_2014Arrest1OxyRem,T1011_2014Arrest1Flow,align="v",nrow=3,heights=c(0.5,0.5,2))


#T10-11 2014: Arrest 2
T1011_2014Arrest2<-subset(AllData,Instrument=="T1011_2014") %>% filter(Time>= as.POSIXct('11/04/2014 4:30:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('11/04/2014 5:15:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T1011_2014Arrest2Flow<-ggplot(data=subset(T1011_2014Arrest2) ,aes(x=Time))+ArrestFlowTheme

T1011_2014Arrest2OxyRem<-ggplot(data=subset(T1011_2014Arrest2) , aes(x=Time))+OxyTheme

T1011_2014Arrest2BackScatter<-ggplot(data=subset(T1011_2014Arrest2) , aes(x=Time))+
  BackScatterTheme+labs(title="T1011_2014")

T1011_2014Arrest2Plot<-ggarrange(T1011_2014Arrest2BackScatter,T1011_2014Arrest2OxyRem,T1011_2014Arrest2Flow,align="v",nrow=3,heights=c(0.5,0.5,2))


#T10-11B 2014-Arrest1
T1011B_2014Arrest1<-subset(AllData,Instrument=="T1011B_2014") %>% filter(Time>= as.POSIXct('11/06/2014 18:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('11/06/2014 20:00:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T1011B_2014Arrest1Flow<-ggplot(data=subset(T1011B_2014Arrest1) ,aes(x=Time))+ArrestFlowTheme

T1011B_2014Arrest1OxyRem<-ggplot(data=subset(T1011B_2014Arrest1) , aes(x=Time))+OxyTheme

T1011B_2014Arrest1BackScatter<-ggplot(data=subset(T1011B_2014Arrest1) , aes(x=Time))+
  BackScatterTheme+labs(title="T1011B_2014")

T1011B_2014Arrest1Plot<-ggarrange(T1011B_2014Arrest1BackScatter,T1011B_2014Arrest1OxyRem,T1011B_2014Arrest1Flow,align="v",nrow=3,heights=c(0.5,0.5,2))


#Arrest2
T1011B_2014Arrest2<-subset(AllData,Instrument=="T1011B_2014") %>% filter(Time>= as.POSIXct('11/06/2014 11:05:00', format="%m/%d/%Y %H:%M:%S",tz="GMT") & Time<= as.POSIXct('11/06/2014 11:30:00', format="%m/%d/%Y %H:%M:%S",tz="GMT"))

T1011B_2014Arrest2Flow<-ggplot(data=subset(T1011B_2014Arrest2) ,aes(x=Time))+ArrestFlowTheme

T1011B_2014Arrest2OxyRem<-ggplot(data=subset(T1011B_2014Arrest2) , aes(x=Time))+OxyTheme

T1011B_2014Arrest2BackScatter<-ggplot(data=subset(T1011B_2014Arrest2) , aes(x=Time))+
BackScatterTheme+labs(title="T1011B_2014")

T1011B_2014Arrest2Plot<-ggarrange(T1011B_2014Arrest2BackScatter,T1011B_2014Arrest2OxyRem,T1011B_2014Arrest2Flow,align="v",nrow=3,heights=c(0.5,0.5,2))

#Arrest Summary
T89_2013Arrest1Plot

T1415_2013Arrest1Plot

T1213_2013Arrest1Plot

T1011_2014Arrest1Plot
T1011_2014Arrest2Plot

T1011B_2014Arrest1Plot
T1011B_2014Arrest2Plot

AllArrests<-c(T89_2013Arrest1Plot,
                 T1415_2013Arrest1Plot,
                 T1213_2013Arrest1Plot,
                 T1011_2014Arrest1Plot,
                 T1011_2014Arrest2Plot,
                 T1011B_2014Arrest1Plot,
                 T1011B_2014Arrest2Plot)

#Save all to a single PDF

ggsave(paste("T1011B_2014Arrest2Plot",".pdf"),
       path="./PDF Outputs/Supplementary",
       plot=T1011B_2014Arrest2Plot,
       useDingbats=FALSE,
       width=3,
       height=7)



```

##Figure Supp. 2: Changepoint Analysis of T89_2014

###*Generating Changepoints*

Changepoints were generated to compare oxygen removal during periods where excurrent flow was consistently high, but ambient flow was either low or high. 

Using the code below you: 

1. Determine mean-variance changepoints using the "changepoints" package. The minimum segment length is set to a half hour (180 points at 10 s intervals). You can change the segment length by changing `minseglen=`. The changepoints are placed where the mean and/or variance changes. 

2. Align the mean and variance to the T8-9_2014 data and fill down so you get continuous lines representing segments. 

3. Sort the ambient and excurrent segments into "Low" and "High" categories based on the table below. The segments are assigned categories based on the mean, and segments with high variance were taken out (decided that if variance was more than half of the mean we should not use them). 

**Flow Category**| **Ambient**  | **Excurrent**
-------| ------- | -------
Low    | < 2 cm/s | < 2 cm/s
High   | > 5 cm/s | > 4 cm/s


4. Visualize the changepoints. 

```{r Changepoints in 89,  warning=FALSE, message=FALSE}

#Subset T8-9
T89_2014<-subset(AllData,Instrument=="T89_2014")

#Excurrent Changepoints with minimum segement length of 180 (30 minutes)
T89_2014_ExChp<-cpt.meanvar(T89_2014$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T89_2014_MVEx<-as.data.frame(T89_2014_ExChp@param.est)  #mean variance
T89_2014_ChpLocEx<-as.data.frame(T89_2014_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T89_2014_AmChp<-cpt.meanvar(T89_2014$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T89_2014_MVAm<-as.data.frame(T89_2014_AmChp@param.est) #mean variance
T89_2014_ChpLocAm<-as.data.frame(T89_2014_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T89_2014_ChangepointAmData<-as.data.frame(c(T89_2014_ChpLocAm,T89_2014_MVAm))
T89_2014_ChangepointExData<-as.data.frame(c(T89_2014_ChpLocEx,T89_2014_MVEx))

#Rename columns so they're ex and am specific 
colnames(T89_2014_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T89_2014_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T89_2014_ChangepointExData$ChangeStart<-c(1,head(T89_2014_ChangepointExData$ExChp.cpts,-1))
T89_2014_ChangepointAmData$ChangeStart<-c(1,head(T89_2014_ChangepointAmData$AmChp.cpts,-1))
T89_2014$count<-seq.int(nrow(T89_2014))


#Aligning changepoint means/variances to original data
T89_2014_ChpCatEx<-merge(T89_2014,T89_2014_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T89_2014_ChpCat<-merge(T89_2014_ChpCatEx,T89_2014_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T89_2014_ChpCat<-fill_(T89_2014_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

#Adding Changepoint Categories
T89_2014_ChpCat$ExCat <- ifelse(T89_2014_ChpCat$ExMean <= 2, 'Low', 
                          ifelse(T89_2014_ChpCat$ExMean <= 4, 'Med', 'High'))
T89_2014_ChpCat$AmCat <- ifelse(T89_2014_ChpCat$AmMean <= 2, 'Low', 
                          ifelse(T89_2014_ChpCat$AmMean <= 5, 'Med', 'High'))

#Calculating Percentile cutoffs
#T89_2014ExHighCut<-T89_2014%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T89_2014ExLowCut <-T89_2014%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T89_2014AmHighCut<-T89_2014%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T89_2014AmLowCut <-T89_2014%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()

T89_2014ExHighCut<-4
T89_2014ExLowCut <-2
T89_2014AmHighCut<-5
T89_2014AmLowCut <-2

#Adding Changepoint CategoriesT89_2014
T89_2014_ChpCat$ExCat <-  ifelse(T89_2014_ChpCat$ExMean <=  T89_2014ExLowCut, 'Low', 
                          ifelse(T89_2014_ChpCat$ExMean <=  T89_2014ExHighCut, 'Med', 'High'))
T89_2014_ChpCat$AmCat <-  ifelse(T89_2014_ChpCat$AmMean <=  T89_2014AmLowCut, 'Low', 
                          ifelse(T89_2014_ChpCat$AmMean <=  T89_2014AmHighCut, 'Med', 'High'))

#Remove areas with high variance (more than half of the mean)
T89_2014_ChpCat$ExCat<-ifelse(T89_2014_ChpCat$ExVar<(T89_2014_ChpCat$ExMean/2),
                              T89_2014_ChpCat$ExCat,'NA')
T89_2014_ChpCat$AmCat<-ifelse(T89_2014_ChpCat$AmVar<(T89_2014_ChpCat$AmMean/2),
                              T89_2014_ChpCat$AmCat,'NA')

#Make plots visualizing the changepoints and categories 

T89_2014_AmChangepoints<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T89_2014AmHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T89_2014AmLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="A.")

T89_2014_ExChangepoints<-ggplot(T89_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T89_2014ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T89_2014ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category",tag="B.")

#Visualize Plots
T89_2014_ChangepointsPlot<-T89_2014_AmChangepoints/T89_2014_ExChangepoints

T89_2014_ChangepointsPlot

ggsave("T89 Changepoints.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T89_2014_ChangepointsPlot,
       width=12,
       height=12,
       useDingbats=FALSE)
```

##Other Changepoints

Changepoints for other thermistors were generated to see if the pattern seen above holds up. In just two of the below we see a similar pattern (T1415_2013 and T1011B_2014), and it looks considerable only in 14-15_2013. 

T89_2013
```{r T89_2013 Changepoints}
#Subset T14-15 2013
T89_2013<-subset(AllData,Instrument=="T89_2013")

#Excurrent Changepoints with minimum segement length of 180 (30 minutes)
T89_2013_ExChp<-cpt.meanvar(T89_2013$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T89_2013_MVEx<-as.data.frame(T89_2013_ExChp@param.est)  #mean variance
T89_2013_ChpLocEx<-as.data.frame(T89_2013_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T89_2013_AmChp<-cpt.meanvar(T89_2013$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T89_2013_MVAm<-as.data.frame(T89_2013_AmChp@param.est) #mean variance
T89_2013_ChpLocAm<-as.data.frame(T89_2013_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T89_2013_ChangepointAmData<-as.data.frame(c(T89_2013_ChpLocAm,T89_2013_MVAm))
T89_2013_ChangepointExData<-as.data.frame(c(T89_2013_ChpLocEx,T89_2013_MVEx))

#Rename columns so they're ex and am specific 
colnames(T89_2013_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T89_2013_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T89_2013_ChangepointExData$ChangeStart<-c(1,head(T89_2013_ChangepointExData$ExChp.cpts,-1))
T89_2013_ChangepointAmData$ChangeStart<-c(1,head(T89_2013_ChangepointAmData$AmChp.cpts,-1))
T89_2013$count<-seq.int(nrow(T89_2013))


#Aligning changepoint means/variances to original data
T89_2013_ChpCatEx<-merge(T89_2013,T89_2013_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T89_2013_ChpCat<-merge(T89_2013_ChpCatEx,T89_2013_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T89_2013_ChpCat<-fill_(T89_2013_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

#Adding Changepoint Categories
T89_2013_ChpCat$ExCat <- ifelse(T89_2013_ChpCat$ExMean <= 2, 'Low', 
                          ifelse(T89_2013_ChpCat$ExMean <= 4, 'Med', 'High'))
T89_2013_ChpCat$AmCat <- ifelse(T89_2013_ChpCat$AmMean <= 2, 'Low', 
                          ifelse(T89_2013_ChpCat$AmMean <= 5, 'Med', 'High'))

#Calculating Percentile cutoffs
#T89_2013ExHighCut<-T89_2013%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T89_2013ExLowCut <-T89_2013%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T89_2013AmHighCut<-T89_2013%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T89_2013AmLowCut <-T89_2013%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()

#Manual Cutoffs
T89_2013ExHighCut<-7
T89_2013ExLowCut <-5
T89_2013AmHighCut<-7
T89_2013AmLowCut <-6

#Adding Changepoint Categories
T89_2013_ChpCat$ExCat <- ifelse( T89_2013_ChpCat$ExMean <= T89_2013ExLowCut, 'Low', 
                          ifelse(T89_2013_ChpCat$ExMean <= T89_2013ExHighCut, 'Med', 'High'))
T89_2013_ChpCat$AmCat <- ifelse( T89_2013_ChpCat$AmMean <= T89_2013AmLowCut, 'Low', 
                          ifelse(T89_2013_ChpCat$AmMean <= T89_2013AmHighCut, 'Med', 'High'))

#Make plots visualizing the changepoints and categories 

T89_2013_AmChangepoints<-ggplot(T89_2013_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T89_2013ExHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T89_2013ExLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="A.")

T89_2013_ExChangepoints<-ggplot(T89_2013_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T89_2013ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T89_2013ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category",tag="B.")

#Visualize Plots
T89_2013_AmChangepoints/T89_2013_ExChangepoints

T89_2013_OxyViolin<-ggplot(subset(T89_2013_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                           aes(x=AmCat,y=Oxy))+
  geom_violin(fill="#19AD8E",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Oxygen Removal \n (uMol)",tag="a.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
  stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

#Excurrent flow is not different
T89_2013_ExViolin<-ggplot(subset(T89_2013_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                    aes(x=AmCat,y=Ex))+
  geom_violin(fill="red",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="High Excurrent \n Flow Category (cm/s)",tag="b.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
    stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

T89_2013_OxyViolin|T89_2013_ExViolin

```

T1415_2013
```{r T1415_2013 Changepoints}
#Subset T14-15 2013
T1415_2013<-subset(AllData,Instrument=="T1415_2013")

#Excurrent Changepoints with minimum segement length of 180 (30 minutes)
T1415_2013_ExChp<-cpt.meanvar(T1415_2013$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T1415_2013_MVEx<-as.data.frame(T1415_2013_ExChp@param.est)  #mean variance
T1415_2013_ChpLocEx<-as.data.frame(T1415_2013_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T1415_2013_AmChp<-cpt.meanvar(T1415_2013$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T1415_2013_MVAm<-as.data.frame(T1415_2013_AmChp@param.est) #mean variance
T1415_2013_ChpLocAm<-as.data.frame(T1415_2013_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T1415_2013_ChangepointAmData<-as.data.frame(c(T1415_2013_ChpLocAm,T1415_2013_MVAm))
T1415_2013_ChangepointExData<-as.data.frame(c(T1415_2013_ChpLocEx,T1415_2013_MVEx))

#Rename columns so they're ex and am specific 
colnames(T1415_2013_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T1415_2013_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T1415_2013_ChangepointExData$ChangeStart<-c(1,head(T1415_2013_ChangepointExData$ExChp.cpts,-1))
T1415_2013_ChangepointAmData$ChangeStart<-c(1,head(T1415_2013_ChangepointAmData$AmChp.cpts,-1))
T1415_2013$count<-seq.int(nrow(T1415_2013))


#Aligning changepoint means/variances to original data
T1415_2013_ChpCatEx<-merge(T1415_2013,T1415_2013_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T1415_2013_ChpCat<-merge(T1415_2013_ChpCatEx,T1415_2013_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T1415_2013_ChpCat<-fill_(T1415_2013_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

#Adding Changepoint Categories
T1415_2013_ChpCat$ExCat <- ifelse(T1415_2013_ChpCat$ExMean <= 2, 'Low', 
                          ifelse(T1415_2013_ChpCat$ExMean <= 4, 'Med', 'High'))
T1415_2013_ChpCat$AmCat <- ifelse(T1415_2013_ChpCat$AmMean <= 2, 'Low', 
                          ifelse(T1415_2013_ChpCat$AmMean <= 5, 'Med', 'High'))

#Changepoint cutoffs (top and bottom 30%)
#T1415_2013ExHighCut<-T1415_2013%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T1415_2013ExLowCut <-T1415_2013%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T1415_2013AmHighCut<-T1415_2013%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T1415_2013AmLowCut <-T1415_2013%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()


T1415_2013ExHighCut<-2.5
T1415_2013ExLowCut <-2.25
T1415_2013AmHighCut<-5
T1415_2013AmLowCut <-4
  
#Adding Changepoint Categories
T1415_2013_ChpCat$ExCat <- ifelse( T1415_2013_ChpCat$ExMean <= T1415_2013ExLowCut, 'Low', 
                            ifelse(T1415_2013_ChpCat$ExMean <= T1415_2013ExHighCut, 'Med', 'High'))
T1415_2013_ChpCat$AmCat <- ifelse( T1415_2013_ChpCat$AmMean <= T1415_2013AmLowCut, 'Low', 
                          ifelse(  T1415_2013_ChpCat$AmMean <= T1415_2013AmHighCut, 'Med', 'High'))


#Make plots visualizing the changepoints and categories 

T1415_2013_AmChangepoints<-ggplot(T1415_2013_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T1415_2013AmHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T1415_2013AmLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="A.")

T1415_2013_ExChangepoints<-ggplot(T1415_2013_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T1415_2013ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T1415_2013ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category")

#Visualize Plots
T1415_2013_AmChangepoints/T1415_2013_ExChangepoints

T1415_2013_OxyViolin<-ggplot(subset(T1415_2013_ChpCat,
                           (AmCat=="Low"|
                            AmCat=="High")&
                            ExCat=="High"),
                           aes(x=AmCat,y=Oxy))+
  geom_violin(fill="#19AD8E",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Oxygen Removal \n (uMol)",tag="B.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
  stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

#Excurrent flow is not different
T1415_2013_ExViolin<-ggplot(subset(T1415_2013_ChpCat,
                           (AmCat=="Low"|
                            AmCat=="High")&
                            ExCat=="High"),
                    aes(x=AmCat,y=Ex))+
  geom_violin(fill="red",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="High Excurrent \n Flow Category (cm/s)")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
    stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

T1415_2013_OxyViolin|T1415_2013_ExViolin

###Violin Plot Table for T1415_2013  

#Troubleshoot: might need to detatch(package:plyr)

T1415_Violin_Table<-T1415_2013_ChpCat%>%
 filter(row_number() %% 20 == 1)%>%
            subset((AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High") %>%
            group_by("Ambient Category"=AmCat)%>%
            summarise("Mean Ambient"=mean(Am),
                      "Ambient SD"=sd(Am),
                      "Mean Excurrent"=mean(Ex),
                      "Excurrent SD"=sd(Ex),
                      "Mean Oxygen"=mean(Oxy),
                      "Oxygen SD"= sd(Oxy),
                      n=n())

T1415_Violin_TableCI<-T1415_Violin_Table%>%
  add_column("Oxygen CI"=round(qt(0.975,df=T1415_Violin_Table$n-1)*
  T1415_Violin_Table$`Oxygen SD`/sqrt(T1415_Violin_Table$n),3),.after="Mean Oxygen")%>%
  add_column("Excurrent CI"=round(qt(0.975,df=T1415_Violin_Table$n-1)*
  T1415_Violin_Table$`Excurrent SD`/sqrt(T1415_Violin_Table$n),3),.after="Mean Excurrent")%>%
  add_column("Ambient CI"=round(qt(0.975,df=T1415_Violin_Table$n-1)*
  T1415_Violin_Table$`Ambient SD`/sqrt(T1415_Violin_Table$n),4),.after="Mean Ambient")

T1415_Violin_TableB<-kable_styling(kable(T1415_Violin_TableCI),bootstrap_options = c("striped", "hover"))

T1415_Violin_TableB

paste("For T1415_2013, mean oxygen uptake under high ambient flow is",
      round((T1415_Violin_Table$'Mean Oxygen'[1]/T1415_Violin_Table$'Mean Oxygen'[2])*100,2),
      "% that of low ambient flow.")

```


T1213_2013
```{r T1213_2013 Changepoints}
#Subset T1213_2013 
T1213_2013<-subset(AllData,Instrument=="T1213_2013")

#Excurrent Changepoints with minimum segement length of 180 (30 minutes)
T1213_2013_ExChp<-cpt.meanvar(T1213_2013$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T1213_2013_MVEx<-as.data.frame(T1213_2013_ExChp@param.est)  #mean variance
T1213_2013_ChpLocEx<-as.data.frame(T1213_2013_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T1213_2013_AmChp<-cpt.meanvar(T1213_2013$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T1213_2013_MVAm<-as.data.frame(T1213_2013_AmChp@param.est) #mean variance
T1213_2013_ChpLocAm<-as.data.frame(T1213_2013_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T1213_2013_ChangepointAmData<-as.data.frame(c(T1213_2013_ChpLocAm,T1213_2013_MVAm))
T1213_2013_ChangepointExData<-as.data.frame(c(T1213_2013_ChpLocEx,T1213_2013_MVEx))

#Rename columns so they're ex and am specific 
colnames(T1213_2013_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T1213_2013_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T1213_2013_ChangepointExData$ChangeStart<-c(1,head(T1213_2013_ChangepointExData$ExChp.cpts,-1))
T1213_2013_ChangepointAmData$ChangeStart<-c(1,head(T1213_2013_ChangepointAmData$AmChp.cpts,-1))
T1213_2013$count<-seq.int(nrow(T1213_2013))


#Aligning changepoint means/variances to original data
T1213_2013_ChpCatEx<-merge(T1213_2013,T1213_2013_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T1213_2013_ChpCat<-merge(T1213_2013_ChpCatEx,T1213_2013_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T1213_2013_ChpCat<-fill_(T1213_2013_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

#Adding Changepoint Categories
T1213_2013_ChpCat$ExCat <- ifelse(T1213_2013_ChpCat$ExMean <= 2, 'Low', 
                          ifelse(T1213_2013_ChpCat$ExMean <= 4, 'Med', 'High'))
T1213_2013_ChpCat$AmCat <- ifelse(T1213_2013_ChpCat$AmMean <= 2, 'Low', 
                          ifelse(T1213_2013_ChpCat$AmMean <= 5, 'Med', 'High'))

##Calculating Percentile cutoffs
#T1213_2013ExHighCut<-T1213_2013%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T1213_2013ExLowCut <-T1213_2013%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T1213_2013AmHighCut<-T1213_2013%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T1213_2013AmLowCut <-T1213_2013%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()

T1213_2013ExHighCut<-7
T1213_2013ExLowCut <-5
T1213_2013AmHighCut<-9
T1213_2013AmLowCut <-5

#Adding Changepoint Categories
T1213_2013_ChpCat$ExCat <- ifelse(T1213_2013_ChpCat$ExMean <= T1213_2013ExLowCut, 'Low', 
                          ifelse(T1213_2013_ChpCat$ExMean <= T1213_2013ExHighCut, 'Med', 'High'))
T1213_2013_ChpCat$AmCat <- ifelse(T1213_2013_ChpCat$AmMean <= T1213_2013AmLowCut, 'Low', 
                          ifelse(T1213_2013_ChpCat$AmMean <= T1213_2013AmHighCut, 'Med', 'High'))

#Make plots visualizing the changepoints and categories 

T1213_2013_AmChangepoints<-ggplot(T1213_2013_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T1213_2013AmHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T1213_2013AmLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="A.")

T1213_2013_ExChangepoints<-ggplot(T1213_2013_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T1213_2013ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T1213_2013ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category",tag="B.")

#Visualize Plots
T1213_2013_AmChangepoints/T1213_2013_ExChangepoints

T1213_2013_OxyViolin<-ggplot(subset(T1213_2013_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                           aes(x=AmCat,y=Oxy))+
  geom_violin(fill="#19AD8E",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Oxygen Removal \n (uMol)",tag="a.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
  stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

#Excurrent flow is not different

T1213_2013_ExViolin<-ggplot(subset(T1213_2013_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                    aes(x=AmCat,y=Ex))+
  geom_violin(fill="red",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="High Excurrent \n Flow Category (cm/s)",tag="b.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
    stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

T1213_2013_OxyViolin|T1213_2013_ExViolin
```


T1011_2014
```{r T1011_2014 Changepoints}
#Subset T1011_2014 
T1011_2014<-subset(AllData,Instrument=="T1011_2014")


#Excurrent Changepoints with minimum segement length of 180 (30 minutes)
T1011_2014_ExChp<-cpt.meanvar(T1011_2014$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T1011_2014_MVEx<-as.data.frame(T1011_2014_ExChp@param.est)  #mean variance
T1011_2014_ChpLocEx<-as.data.frame(T1011_2014_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T1011_2014_AmChp<-cpt.meanvar(T1011_2014$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T1011_2014_MVAm<-as.data.frame(T1011_2014_AmChp@param.est) #mean variance
T1011_2014_ChpLocAm<-as.data.frame(T1011_2014_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T1011_2014_ChangepointAmData<-as.data.frame(c(T1011_2014_ChpLocAm,T1011_2014_MVAm))
T1011_2014_ChangepointExData<-as.data.frame(c(T1011_2014_ChpLocEx,T1011_2014_MVEx))

#Rename columns so they're ex and am specific 
colnames(T1011_2014_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T1011_2014_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T1011_2014_ChangepointExData$ChangeStart<-c(1,head(T1011_2014_ChangepointExData$ExChp.cpts,-1))
T1011_2014_ChangepointAmData$ChangeStart<-c(1,head(T1011_2014_ChangepointAmData$AmChp.cpts,-1))
T1011_2014$count<-seq.int(nrow(T1011_2014))


#Aligning changepoint means/variances to original data
T1011_2014_ChpCatEx<-merge(T1011_2014,T1011_2014_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T1011_2014_ChpCat<-merge(T1011_2014_ChpCatEx,T1011_2014_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T1011_2014_ChpCat<-fill_(T1011_2014_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

##Calculating Percentile cutoffs
#T1011_2014ExHighCut<-T1011_2014%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T1011_2014ExLowCut <-T1011_2014%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T1011_2014AmHighCut<-T1011_2014%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T1011_2014AmLowCut <-T1011_2014%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()

T1011_2014ExHighCut<-5.5
T1011_2014ExLowCut <-5
T1011_2014AmHighCut<-6
T1011_2014AmLowCut <-4

#Adding Changepoint Categories
T1011_2014_ChpCat$ExCat <- ifelse(T1011_2014_ChpCat$ExMean <= T1011_2014ExLowCut, 'Low', 
                           ifelse(T1011_2014_ChpCat$ExMean <= T1011_2014ExHighCut, 'Med', 'High'))
T1011_2014_ChpCat$AmCat <- ifelse(T1011_2014_ChpCat$AmMean <= T1011_2014AmLowCut, 'Low', 
                           ifelse(T1011_2014_ChpCat$AmMean <= T1011_2014AmHighCut, 'Med', 'High'))


#Make plots visualizing the changepoints and categories 

T1011_2014_AmChangepoints<-ggplot(T1011_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T1011_2014AmHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T1011_2014AmLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="A.")

T1011_2014_ExChangepoints<-ggplot(T1011_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T1011_2014ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T1011_2014ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category",tag="B.")

#Visualize Plots
T1011_2014_AmChangepoints/T1011_2014_ExChangepoints

T1011_2014_OxyViolin<-ggplot(subset(T1011_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                           aes(x=AmCat,y=Oxy))+
  geom_violin(fill="#19AD8E",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Oxygen Removal \n (uMol)",tag="a.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
  stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

#Excurrent flow is not different

T1011_2014_ExViolin<-ggplot(subset(T1011_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                    aes(x=AmCat,y=Ex))+
  geom_violin(fill="red",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="High Excurrent \n Flow Category (cm/s)",tag="b.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black",face="bold"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
    stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

T1011_2014_OxyViolin|T1011_2014_ExViolin



```

T1011B_2014
```{r T1011B_2014 Changepoints}
#Subset T14-15 2013
T1011B_2014<-subset(AllData,Instrument=="T1011B_2014")

#Excurrent Changepoints with minimum segement length of 180 (30 minutes)
T1011B_2014_ExChp<-cpt.meanvar(T1011B_2014$Ex,method="PELT",minseglen = 180)

#Save excurrent changepoints as dataframe
T1011B_2014_MVEx<-as.data.frame(T1011B_2014_ExChp@param.est)  #mean variance
T1011B_2014_ChpLocEx<-as.data.frame(T1011B_2014_ExChp@cpts)   #changepoint location

#Ambient Changepoints with minimum segement length of 180 (30 minutes)
T1011B_2014_AmChp<-cpt.meanvar(T1011B_2014$Am,method="PELT",minseglen = 180)

#save ambient changepoints as dataframe
T1011B_2014_MVAm<-as.data.frame(T1011B_2014_AmChp@param.est) #mean variance
T1011B_2014_ChpLocAm<-as.data.frame(T1011B_2014_AmChp@cpts)  #changepoint location 

#Merging ambient and excurrent changepoint means/variances to changepoint locations
require(data.table)
T1011B_2014_ChangepointAmData<-as.data.frame(c(T1011B_2014_ChpLocAm,T1011B_2014_MVAm))
T1011B_2014_ChangepointExData<-as.data.frame(c(T1011B_2014_ChpLocEx,T1011B_2014_MVEx))

#Rename columns so they're ex and am specific 
colnames(T1011B_2014_ChangepointAmData)<-c("AmChp.cpts","AmMean","AmVar")
colnames(T1011B_2014_ChangepointExData)<-c("ExChp.cpts","ExMean","ExVar")

#Align changepoint location of mean/variance to where segments start, rather than end
T1011B_2014_ChangepointExData$ChangeStart<-c(1,head(T1011B_2014_ChangepointExData$ExChp.cpts,-1))
T1011B_2014_ChangepointAmData$ChangeStart<-c(1,head(T1011B_2014_ChangepointAmData$AmChp.cpts,-1))
T1011B_2014$count<-seq.int(nrow(T1011B_2014))


#Aligning changepoint means/variances to original data
T1011B_2014_ChpCatEx<-merge(T1011B_2014,T1011B_2014_ChangepointExData[2:4],by.x="count",
                   by.y="ChangeStart",all.x=T)

T1011B_2014_ChpCat<-merge(T1011B_2014_ChpCatEx,T1011B_2014_ChangepointAmData[2:4],by.x="count",
                 by.y="ChangeStart",all.x=T)

#The means/variances at changepoints have been added, now we fill down
T1011B_2014_ChpCat<-fill_(T1011B_2014_ChpCat,c("ExMean","ExVar","AmMean","AmVar"))

#Changepoint cutoffs (top and bottom 30%)
#T1011B_2014ExHighCut<-T1011B_2014%>%top_frac(0.30,wt=Ex)%>%summarise(min(Ex))%>%as.numeric()
#T1011B_2014ExLowCut <-T1011B_2014%>%top_frac(-0.30,wt=Ex)%>%summarise(max(Ex))%>%as.numeric()
#
#T1011B_2014AmHighCut<-T1011B_2014%>%top_frac(0.30,wt=Am)%>%summarise(min(Am))%>%as.numeric()
#T1011B_2014AmLowCut <-T1011B_2014%>%top_frac(-0.30,wt=Am)%>%summarise(max(Am))%>%as.numeric()

#Manual cutoffs based on visual inspection of the graph
T1011B_2014ExHighCut<-4
T1011B_2014ExLowCut <-3
T1011B_2014AmHighCut<-5.5
T1011B_2014AmLowCut <-3
  
#Adding Changepoint CategoriesT1011B_2014
T1011B_2014_ChpCat$ExCat <- ifelse( T1011B_2014_ChpCat$ExMean <= T1011B_2014ExLowCut, 'Low', 
                            ifelse( T1011B_2014_ChpCat$ExMean <= T1011B_2014ExHighCut, 'Med', 'High'))
T1011B_2014_ChpCat$AmCat <- ifelse( T1011B_2014_ChpCat$AmMean <= T1011B_2014AmLowCut, 'Low', 
                          ifelse(   T1011B_2014_ChpCat$AmMean <= T1011B_2014AmHighCut, 'Med', 'High'))

#Make plots visualizing the changepoints and categories 

T1011B_2014_AmChangepoints<-ggplot(T1011B_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Am),alpha=0.25)+
  geom_ribbon(aes(ymin=AmMean-AmVar,ymax=AmMean+AmVar),fill="blue",alpha=0.3)+
  geom_line(aes(y=AmMean),col="blue")+
  geom_rug(aes(col=AmCat),sides="t")+
  scale_color_manual(values=c("lightblue","darkblue","NA","NA"))+
  geom_hline(yintercept=T1011B_2014AmHighCut,col="lightblue",alpha=0.8)+
  geom_hline(yintercept=T1011B_2014AmLowCut,col="darkblue",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Ambient Flow \n cm/s",col="Flow Category",tag="C.")

T1011B_2014_ExChangepoints<-ggplot(T1011B_2014_ChpCat,aes(x=Time))+
  geom_line(aes(y=Ex),alpha=0.25)+
  geom_ribbon(aes(ymin=ExMean-ExVar,ymax=ExMean+ExVar),fill="red",alpha=0.3)+
  geom_line(aes(y=ExMean),col="red")+
  geom_rug(aes(col=ExCat),sides="t")+
  scale_color_manual(values=c("red","darkred","NA","NA"))+
  geom_hline(yintercept=T1011B_2014ExLowCut,col="darkred",alpha=0.8)+
  geom_hline(yintercept=T1011B_2014ExHighCut,col="red",alpha=0.8)+
  theme_classic()+
  labs(x="Time",y="Excurrent Flow \n cm/s",col="Flow Category")

#Visualize Plots
T1011B_2014_AmChangepoints/T1011B_2014_ExChangepoints

T1011B_2014_OxyViolin<-ggplot(subset(T1011B_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                           aes(x=AmCat,y=Oxy))+
  geom_violin(fill="#19AD8E",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="Oxygen Removal \n (uMol)",tag="D.")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
  stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

#Excurrent flow is not different
T1011B_2014_ExViolin<-ggplot(subset(T1011B_2014_ChpCat,
                           (AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High"),
                    aes(x=AmCat,y=Ex))+
  geom_violin(fill="red",colour="NA",alpha=0.5)+
  stat_summary(fun.data=mean_sdl, geom="pointrange", color="black",size=1)+
  stat_summary(fun=median,geom="point", color="white",size=2)+
  theme_light()+
  labs(x="Ambient \n Flow Category",y="High Excurrent \n Flow Category (cm/s)")+
  theme_classic()+
  theme(axis.text.x = element_text(angle=0,vjust=0.5,colour="black",size=10),
        text = element_text(size=12,colour = "black"),
        line=element_line(size=1),
        axis.text.y=element_text(size=10,colour="black"),
        panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(0,0,0,0), "cm"))+
    stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.25,
    vjust = -1,angle=90)

T1011B_2014_OxyViolin|T1011B_2014_ExViolin

###Violin Plot Table for T1011B_2014  

#Troubleshoot: might need to detatch(package:plyr)

T1011B_2014_Violin_Table<-T1011B_2014_ChpCat%>%
  filter(row_number() %% 20 == 1)%>%
            subset((AmCat=="Low"|
                           AmCat=="High")&
                           ExCat=="High") %>%
            group_by("Ambient Category"=AmCat)%>%
            summarise("Mean Ambient"=mean(Am),
                      "Ambient SD"=sd(Am),
                      "Mean Excurrent"=mean(Ex),
                      "Excurrent SD"=sd(Ex),
                      "Mean Oxygen"=mean(Oxy),
                      "Oxygen SD"= sd(Oxy),
                      n=n())

T1011B_2014_Violin_TableCI<-T1011B_2014_Violin_Table%>%
  add_column("Oxygen CI"=round(qt(0.975,df=T1011B_2014_Violin_Table$n-1)*
  T1011B_2014_Violin_Table$`Oxygen SD`/sqrt(T1011B_2014_Violin_Table$n),3),.after="Mean Oxygen")%>%
  add_column("Excurrent CI"=round(qt(0.975,df=T1011B_2014_Violin_Table$n-1)*
  T1011B_2014_Violin_Table$`Excurrent SD`/sqrt(T1011B_2014_Violin_Table$n),3),.after="Mean Excurrent")%>%
  add_column("Ambient CI"=round(qt(0.975,df=T1011B_2014_Violin_Table$n-1)*
  T1011B_2014_Violin_Table$`Ambient SD`/sqrt(T1011B_2014_Violin_Table$n),4),.after="Mean Ambient")

T1011B_2014_Violin_TableCI_B<-kable_styling(kable(T1011B_2014_Violin_TableCI),bootstrap_options = c("striped", "hover"))

T1011B_2014_Violin_TableCI_B

paste("For T1011B_2014, mean oxygen uptake under high ambient flow is",
      round((T1011B_2014_Violin_Table$'Mean Oxygen'[1]/T1011B_2014_Violin_Table$'Mean Oxygen'[2])*100,2),
      "% that of low ambient flow.")

```

###Figure Supp. 3/Table S2: Changepoints for T14-15 and T10-11B Combined

Figure S3
```{r T1415 and 1011B changepoint figure}

#Violin Plots
T1415_2013_ViolinPlot <-(T1415_2013_OxyViolin|T1415_2013_ExViolin)+plot_layout(widths=5,heights=5)
T1011B_2014_ViolinPlot<-(T1011B_2014_OxyViolin|T1011B_2014_ExViolin)+plot_layout(widths=5,heights=5)

#Changepoints
T1415_2013_ChpPlot <-(T1415_2013_AmChangepoints/T1415_2013_ExChangepoints)+plot_layout(widths=5,heights=5)
T1011B_2014_ChpPlot <-(T1011B_2014_AmChangepoints/T1011B_2014_ExChangepoints)+plot_layout(widths=5,heights=5)

T1415_2013_T1011B_2014_ChpAndViolin<-(T1415_2013_ChpPlot|T1415_2013_ViolinPlot)/(T1011B_2014_ChpPlot|T1011B_2014_ViolinPlot)

ggsave("FigureS3.pdf",
       path="./PDF Outputs/Supplementary",
       plot=T1415_2013_T1011B_2014_ChpAndViolin,
       width=15,
       height=12,
       useDingbats=FALSE)

```

Table S2:  T14-15 and T10-11B violin plot summary
```{r Violin table 1415 and 1011B}

T1415_Violin_TableC<-T1415_Violin_TableCI%>%add_column(Instrument="T1415_2013",.before="Ambient Category")%>%
  mutate_at(3:11,round,2)

T1011B_2014_Violin_TableC<-T1011B_2014_Violin_TableCI%>%add_column(Instrument="T1011B_2014",.before="Ambient Category")%>%
  mutate_at(3:11,round,2)

ViolinTable_1415_and_1011B<-rbind(T1415_Violin_TableC,T1011B_2014_Violin_TableC)

write.csv(ViolinTable_1415_and_1011B,"./Tables/Table S2-ViolinTable_1415_and_1011B.csv")

kable_styling(kable(ViolinTable_1415_and_1011B),bootstrap_options = c("striped", "hover"),position="center")

paste("For T1415_2013, mean oxygen uptake under high ambient flow is",
      round((T1415_Violin_Table$'Mean Oxygen'[1]/T1415_Violin_Table$'Mean Oxygen'[2])*100,2),
      "% that of low ambient flow.")

paste("For T1011B_2014, mean oxygen uptake under high ambient flow is",
      round((T1011B_2014_Violin_Table$'Mean Oxygen'[1]/T1011B_2014_Violin_Table$'Mean Oxygen'[2])*100,2),
      "% that of low ambient flow.")
```

##ADCP vs Ambient Thermistors
Note: ADCP-1 in 2013 had a faulty flow record and should not be considered. 
```{r ADCP vs Ambient,warning=FALSE,message=FALSE}

AmVsADCP1<-ggplot(subset(AllData,Instrument!="NA"),aes(x=ADCP1,y=Am))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="red")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Ambient Thermistor \n cm/s",x="ADCP1 (cm/s)",
       title="Ambient Thermistors vs. ADCP-1")

AmVsADCP2<-ggplot(subset(AllData,Instrument!="NA"),aes(x=ADCP2,y=Am))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="red")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Ambient Thermistor \n cm/s",x="ADCP2 (cm/s)",
       title="Ambient Thermistors vs. ADCP-2")

AmVsBackscatter<-ggplot(subset(AllData,Instrument!="NA"),aes(x=Backscatter,y=Am))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="blue")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Ambient Thermistor \n cm/s",x="Backscatter (counts)",
       title="Ambient Thermistors vs. Backscatter")
  
AmVsADCP1
AmVsADCP2
AmVsBackscatter

```

##ADCP Flow vs. Excurrent Thermistor 

Note: ADCP-1 in 2013 had a faulty flow record and should not be considered. 

```{r ADCP vs Excurrent,warning=FALSE,message=FALSE}

ExVsADCP1<-ggplot(subset(AllData,Instrument!="NA"),aes(x=ADCP1,y=Ex))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="red")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Excurrent Thermistor \n cm/s",x="ADCP1 (cm/s)",
       title="Excurrent Thermistors vs. ADCP-1")

ExVsADCP2<-ggplot(subset(AllData,Instrument!="NA"),aes(x=ADCP2,y=Ex))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="red")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Excurrent Thermistor \n cm/s",x="ADCP2 (cm/s)",
       title="Excurrent Thermistors vs. ADCP-2")

ExVsBackscatter<-ggplot(subset(AllData,Instrument!="NA"),aes(x=Backscatter,y=Ex))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="blue")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Excurrent Thermistor \n cm/s",x="Backscatter (counts)",
       title="Excurrent Thermistors vs. Backscatter")

ExVsADCP1
ExVsADCP2
ExVsBackscatter
```

##Ambient Thermistor vs. Excurrent Thermistor
```{r Ambient vs Excurrent,warning=FALSE,message=FALSE}

ExVsAm<-ggplot(subset(AllData,Instrument!="NA"),aes(x=Am,y=Ex))+
  geom_point(alpha=0.25,size=0.5)+
  geom_smooth(method="lm",se=F,col="red")+
  facet_wrap(vars(Instrument),scales="free")+
  theme_classic()+
  labs(y="Excurrent Thermistor \n cm/s",x="Ambient Thermistor (cm/s)",
       title="Excurrent Thermistors vs. Ambient Thermistors")


ExVsAm

```

##Tables: Correlation Coefficient and Slope 

Pearsons Correlation Coefficient (R) and slope for linear relationships between: 

- Excurrent~Ambient
- Oxygen Uptake~Ambient
- Oxygen Uptake~Excurrent
- Ambient Thermistor~ADCP1

(others might need to be added)

```{r Correlation and Slope,warning=FALSE,message=FALSE}

#Excurrent is positively related to ambient in only one of the instruments, 
# and negatively related in another
ExAm_SlopesAndCorrels<-AllData%>%
  subset(Instrument!="NA")%>%
  group_by(Instrument)%>%
  summarise(Slope = lm(Ex ~ Am,na.action=na.omit)$coeff[2],
            R2 = cor(Ex, Am, method="pearson")^2)%>%
  mutate_at(2:3, round, 3)%>%
  add_column("Correlation Pair"="Excurrent~Ambient",.after="Instrument")

#oxygen is sightly negatively correlated to ambient flow in each case. 
OxyAm_SlopesAndCorrels<-AllData%>%
  subset(Instrument!="NA")%>%
  group_by(Instrument)%>%
  summarise(Slope = lm(Oxy ~ Am,na.action=na.omit)$coeff[2],
            R2 = cor(Oxy, Am, method="pearson")^2)%>%
  mutate_at(2:3, round, 3)%>%
  add_column("Correlation Pair"="Oxygen Uptake~Ambient",.after="Instrument")

#oxygen is negatively correlated to excurrent flow in each case 
# (this is why we only looked) at periods of consistently high excurrent flow
OxyEx_SlopesAndCorrels<-AllData%>%
  subset(Instrument!="NA")%>%
  group_by(Instrument)%>%
  summarise(Slope = lm(Oxy ~ Ex,na.action=na.omit)$coeff[2],
            R2 = cor(Oxy, Ex, method="pearson")^2)%>%
  mutate_at(2:3, round, 3)%>%
  add_column("Correlation Pair"="Oxygen Uptake~Excurrent",.after="Instrument")

#ADCP2 was used because ADCP1 was bad in 2013. 
#ADCP and Ambient, especially in 2014, are positively correlated to each other. 
AmADCP2_SlopesAndCorrels<-AllData%>%
  subset(Instrument!="NA")%>%
  group_by(Instrument)%>%
  summarise(Slope = lm(Am ~ ADCP2,na.action=na.omit)$coeff[2],
            R2 = cor(Am, ADCP2, method="pearson",use="complete.obs")^2)%>%
  mutate_at(2:3, round, 3)%>%
  add_column("Correlation Pair"="Ambient Thermistor~ADCP1",.after="Instrument")

SlopeCorrelationTable<-rbind(ExAm_SlopesAndCorrels,
                             OxyAm_SlopesAndCorrels,
                             OxyEx_SlopesAndCorrels,
                             AmADCP2_SlopesAndCorrels)

write.csv(SlopeCorrelationTable,"./Tables/Table S5-Slope Correlation.csv")


SlopeCorrelationTable%>%
  group_by(`Correlation Pair`)%>%
  mutate(Slope=cell_spec(Slope,color=ifelse(Slope > 0, "black", "grey")),
         R2=color_tile("white","red")(R2))%>%
  kable("html",escape=F)%>%
  kable_styling(position = "center",bootstrap_options = c("striped", "hover"))

```


##Table: Proportion of Flows above 5, 10, and 15 cm/s

This table shows that flows rarely exceed 15 cm/s. Thermistor flows are generally lower than ADCP flows, potentially due to probe placement. 

```{r Flows Above 5 10 and 15,warning=FALSE,message=FALSE}
#Ambient Thermistors
 AmbientThermistor_Proportions<-AllData%>%
  group_by(Instrument)%>%
  summarise("% Below 5"      =length(Am[Am>=0 & Am<=5])/length(Am[Am>=0])*100,
            "% Between 5-10" =length(Am[Am>=5 & Am<=10])/length(Am[Am>=0])*100,
            "% Between 10-15"=length(Am[Am>=10 & Am<=15])/length(Am[Am>=0])*100,
            "% Above 15"     =length(Am[Am>=15 & Am<=Inf])/length(Am[Am>=0])*100,
            "Median"         =median(Am),
            "SD"             =sd(Am))

AmbientThermistor_Proportions[,-1]<-round(AmbientThermistor_Proportions[,-1],1)

AmbientThermistor_Proportions<-AmbientThermistor_Proportions[1:6,]

#ADCP 2013 1 !! This record is corrupted and was not used !!
Aqd2013One_Proportions<- Aqd2013One%>%
  filter(Profile==0.5)%>%
  summarise("% Below 5"      =length(Flow[Flow>=0 & Flow<=5])/length(Flow[Flow>=0])*100,
            "% Between 5-10" =length(Flow[Flow>=5 & Flow<=10])/length(Flow[Flow>=0])*100,
            "% Between 10-15"=length(Flow[Flow>=10 & Flow<=15])/length(Flow[Flow>=0])*100,
            "% Above 15"     =length(Flow[Flow>=15 & Flow<=Inf])/length(Flow[Flow>=0])*100,
            "Median"         =median(Flow),
            "SD"             =sd(Flow))%>%
  round(2)%>%
  add_column(Instrument="Aqd2013One",.before="% Below 5")


#ADCP 2013 2
Aqd2013Two_Proportions<- Aqd2013Two%>%
  filter(Profile==0.5)%>%
  summarise("% Below 5"      =length(Flow[Flow>=0 & Flow<=5])/length(Flow[Flow>=0])*100,
            "% Between 5-10" =length(Flow[Flow>=5 & Flow<=10])/length(Flow[Flow>=0])*100,
            "% Between 10-15"=length(Flow[Flow>=10 & Flow<=15])/length(Flow[Flow>=0])*100,
            "% Above 15"     =length(Flow[Flow>=15 & Flow<=Inf])/length(Flow[Flow>=0])*100,
            "Median"         =median(Flow),
            "SD"             =sd(Flow))%>%
  round(2)%>%
  add_column(Instrument="Aqd2013Two",.before="% Below 5")
 
 #ADCP 2014 1
Aqd2014One_Proportions<- Aqd2014One%>%
  filter(Profile==0.5)%>%
  summarise("% Below 5"      =length(Flow[Flow>=0 & Flow<=5])/length(Flow[Flow>=0])*100,
            "% Between 5-10" =length(Flow[Flow>=5 & Flow<=10])/length(Flow[Flow>=0])*100,
            "% Between 10-15"=length(Flow[Flow>=10 & Flow<=15])/length(Flow[Flow>=0])*100,
            "% Above 15"     =length(Flow[Flow>=15 & Flow<=Inf])/length(Flow[Flow>=0])*100,
            "Median"         =median(Flow),
            "SD"             =sd(Flow))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014One",.before="% Below 5")

  #ADCP 2014 2
Aqd2014Two_Proportions<- Aqd2014Two%>%
  filter(Profile==0.5)%>%
  summarise("% Below 5"      =length(Flow[Flow>=0 & Flow<=5])/length(Flow[Flow>=0])*100,
            "% Between 5-10" =length(Flow[Flow>=5 & Flow<=10])/length(Flow[Flow>=0])*100,
            "% Between 10-15"=length(Flow[Flow>=10 & Flow<=15])/length(Flow[Flow>=0])*100,
            "% Above 15"     =length(Flow[Flow>=15 & Flow<=Inf])/length(Flow[Flow>=0])*100,
            "Median"         =median(Flow),
            "SD"             =sd(Flow))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014Two",.before="% Below 5")


AmbitentFlowConditions<-rbind(AmbientThermistor_Proportions,Aqd2013Two_Proportions,Aqd2014One_Proportions,Aqd2014Two_Proportions)

write.csv(AmbitentFlowConditions,"./Tables/Table S3-Ambient Flow Conditions.csv")

AmbientFlowConditionsTable<-kable_styling(kable(AmbitentFlowConditions),position = "center",bootstrap_options = c("striped", "hover"))

AmbientFlowConditionsTable

 


```

##Ambient Oxygen Removal
```{r Ambient Oxygen Removal,warning=FALSE,message=FALSE}

AllData%>%subset(Instrument!="T1011_2014")%>%group_by(Instrument)%>%summarise("Mean Oxygen"=mean(OxyAm,na.rm=T),"SD"=sd(OxyAm))



```

##Table: ADCP flow and backscatter by profile

Sorted by 1 m intervals to make it easier to read, however ADCP records at 0.25 m profiles. 

This table and plot shows that as the distance from the seafloor increases, the ambient flow increases slightly and the backscatter decreases dramatically.

```{r ADCP profile table }

Aqd2013Two_FlowByProfile<-Aqd2013Two%>%
  summarise("0-1 mab mean"     =mean(Flow[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Flow[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Flow[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Flow[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Flow[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Flow[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Flow[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Flow[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Flow[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Flow[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2013Two Flow",.before="0-1 mab mean")

Aqd2013Two_FlowByProfileB<-Aqd2013Two%>%
  group_by(Profile)%>%
  summarise("meanflow"=mean(Flow),FlowSD=sd(Flow))%>%
  add_column(Instrument="Aqd2013Two Flow")

Aqd2013Two_BackscatterByProfile<-Aqd2013Two%>%
  summarise("0-1 mab mean"     =mean(Backscatter[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Backscatter[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Backscatter[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Backscatter[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Backscatter[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Backscatter[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Backscatter[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Backscatter[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Backscatter[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Backscatter[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2013Two Backscatter",.before="0-1 mab mean")

Aqd2013Two_BackscatterByProfileB<-Aqd2013Two%>%
  group_by(Profile)%>%
  summarise("meanbackscatter"=mean(Backscatter),BackscatterSD=sd(Backscatter))%>%
  add_column(Instrument="Aqd2013Two Backscatter")

Aqd2014One_FlowByProfile<-Aqd2014One%>%
  summarise("0-1 mab mean"     =mean(Flow[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Flow[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Flow[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Flow[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Flow[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Flow[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Flow[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Flow[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Flow[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Flow[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014One Flow",.before="0-1 mab mean")

Aqd2014One_FlowByProfileB<-Aqd2014One%>%
  group_by(Profile)%>%
  summarise("meanflow"=mean(Flow),FlowSD=sd(Flow))%>%
  add_column(Instrument="Aqd2014One Flow")

Aqd2014One_BackscatterByProfile<-Aqd2014One%>%
  summarise("0-1 mab mean"     =mean(Backscatter[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Backscatter[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Backscatter[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Backscatter[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Backscatter[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Backscatter[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Backscatter[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Backscatter[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Backscatter[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Backscatter[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014One Backscatter",.before="0-1 mab mean")

Aqd2014One_BackscatterByProfileB<-Aqd2014One%>%
  group_by(Profile)%>%
  summarise("meanbackscatter"=mean(Backscatter),BackscatterSD=sd(Backscatter))%>%
  add_column(Instrument="Aqd2014One Backscatter")

Aqd2014Two_FlowByProfile<-Aqd2014Two%>%
  summarise("0-1 mab mean"     =mean(Flow[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Flow[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Flow[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Flow[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Flow[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Flow[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Flow[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Flow[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Flow[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Flow[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014Two Flow",.before="0-1 mab mean")

Aqd2014Two_FlowByProfileB<-Aqd2014Two%>%
  group_by(Profile)%>%
  summarise("meanflow"=mean(Flow),FlowSD=sd(Flow))%>%
  add_column(Instrument="Aqd2014Two Flow")

Aqd2014Two_BackscatterByProfile<-Aqd2014Two%>%
  summarise("0-1 mab mean"     =mean(Backscatter[Profile>=0 & Profile<=1]),
            "0-1 mab SD"     =sd(Backscatter[Profile>=0 & Profile<=1]),
            "1-2 mab mean"     =mean(Backscatter[Profile>1 &  Profile<=2]),
            "1-2 mab SD"     =sd(Backscatter[Profile>1 &  Profile<=2]),
            "2-3 mab mean"     =mean(Backscatter[Profile>2 &  Profile<=3]),
            "2-3 mab SD"     =sd(Backscatter[Profile>2 &  Profile<=3]),
            "3-4 mab mean"     =mean(Backscatter[Profile>3 &  Profile<=4]),
            "3-4 mab SD"     =sd(Backscatter[Profile>3 &  Profile<=4]),
            "4-5 mab mean"     =mean(Backscatter[Profile>4 &  Profile<=5]),
            "4-5 mab SD"     =sd(Backscatter[Profile>4 &  Profile<=5]))%>%
  round(2)%>%
  add_column(Instrument="Aqd2014Two Backscatter",.before="0-1 mab mean")

Aqd2014Two_BackscatterByProfileB<-Aqd2014Two%>%
  group_by(Profile)%>%
  summarise("meanbackscatter"=mean(Backscatter),BackscatterSD=sd(Backscatter))%>%
  add_column(Instrument="Aqd2014Two Backscatter")

ADCPby1mProfile<-rbind(Aqd2013Two_FlowByProfile,
                     Aqd2014One_FlowByProfile,
                     Aqd2014Two_FlowByProfile,
                     Aqd2013Two_BackscatterByProfile,
                     Aqd2014One_BackscatterByProfile,
                     Aqd2014Two_BackscatterByProfile)

ADCPFlowbyProfileFull<-rbind(Aqd2013Two_FlowByProfileB,
                             Aqd2014One_FlowByProfileB,
                             Aqd2014Two_FlowByProfileB)

ADCPBackscatterbyProfileFull<-rbind(Aqd2013Two_BackscatterByProfileB,
                                    Aqd2014One_BackscatterByProfileB,
                                    Aqd2014Two_BackscatterByProfileB)
  
write.csv(ADCPby1mProfile,"./Tables/Table S4-ADCP Flow and Backscatter by Profile.csv")

kable_styling(kable(ADCPby1mProfile,caption="mean flow and backscatter by profile"),position = "center",bootstrap_options = c("striped", "hover"))

paste("At 0-1 mab the mean flow is on average",round((mean(ADCPby1mProfile$'0-1 mab'[1:3])
                                            /mean(ADCPby1mProfile$'4-5 mab'[1:3])),2)*100,
      "% that of mean flow rate at 4-5 mab")

paste("At 0-1 mab the mean backscatter is on average",round((mean(ADCPby1mProfile$'0-1 mab'[4:6])
                                            /mean(ADCPby1mProfile$'4-5 mab'[4:6])),2)*100,
      "% that of mean backscatter at 4-5 mab")   
      
FlowByProfile<-ggplot(ADCPFlowbyProfileFull,aes(y=meanflow,x=Profile))+
  geom_line(aes(colour=Instrument))+
  geom_ribbon(aes(ymax=meanflow+FlowSD,ymin=meanflow-FlowSD,
                  fill=Instrument),alpha=0.25)+
  theme_classic()+
  labs(y="mean Flow (cm/s)",x="Profile (mab)",tag="A.")+
  theme(legend.position = "none")

BackscatterByProfile<-ggplot(ADCPBackscatterbyProfileFull,aes(y=meanbackscatter,x=Profile))+
  geom_line(aes(colour=Instrument))+
  geom_ribbon(aes(ymax=meanbackscatter+BackscatterSD,ymin=meanbackscatter-BackscatterSD,
                  fill=Instrument),alpha=0.25)+
  theme_classic()+
  labs(y="mean Backscatter (counts)",x="Profile (mab)",tag="B.")+
  theme(legend.position = c(0.8,0.8))

ADCPbyProfilePlot<-(FlowByProfile|BackscatterByProfile)+plot_annotation(title="Flow and Backscatter by profile (line is mean, shading is ± SD)")

ggsave("Flow and Backscatter by Profile (all instruments).pdf",
       path="./PDF Outputs/Supplementary",
       plot=ADCPbyProfilePlot,
       width=12,
       height=6,
       useDingbats=FALSE)

ADCPbyProfilePlot

```

##Figure 2 Ambient Supp: *ADCP 2013-2 Flow Distribution*

```{r, Figure 2 Ambient Supp 2013-2, warning=FALSE}

DistributionAqd2013Two<-ggplot(Aqd2013Two,aes(x=Flow))+
  lims(x=c(0,60))+
  labs(x="2013 ADCP 2 Flow (cm/s)",y="Frequency")+
  DistributionFormat

DistributionAqd2013Two

ggsave("Figure 2[Ambient] Supp: 2013 flow distribution.pdf",
       path="./PDF Outputs/Supplementary",
       plot=DistributionAqd2013Two,
       width=4,
       height=4,
       useDingbats=FALSE)

```

##Figure 2 Ambient Supp: *ADCP 2014-2 Flow Distribution*

```{r, Figure 2 Ambient Supp 2014-2, warning=FALSE}

DistributionAqd2014Two<-ggplot(Aqd2014Two,aes(x=Flow))+
  lims(x=c(0,60))+
  labs(x="2014 ADCP 2 Flow (cm/s)",y="Frequency")+
  DistributionFormat

DistributionAqd2014Two

ggsave("Figure Ambient Supp: 2014 flow distribution.pdf",
       path="./PDF Outputs/Supplementary",
       plot=DistributionAqd2014Two,
       width=4,
       height=4,
       useDingbats=FALSE)
```

##Calibration Curves
```{r Calibrations}

Calibrations<-read.csv(params$Calibrations,fileEncoding="UTF-8-BOM")

Calibrations$Instrument<-factor(Calibrations$Instrument, levels=c("T89_2013",
                                                          "T1415_2013",
                                                          "T1213_2013",
                                                          "T89_2014",
                                                          "T1011_2014",
                                                          "T1011B_2014"))

Calibrations%>%group_by(Instrument,ExAm)%>%summarise(length(Vector)/6)

ExCalibrationPlots<-ggplot(subset(Calibrations,ExAm=="Ex"),aes(x=LnVector,y=Thermistor))+
  geom_point(col="#F1944B")+
  geom_smooth(method="lm",se=F,col="black")+
  labs(x="ln(Vector) cm/s",y="Excurrent Termistor (V)")+
  theme_classic()+
  facet_wrap(vars(Instrument),scales="free")

AmCalibrationPlots<-ggplot(subset(Calibrations,ExAm=="Am"),aes(x=LnVector,y=Thermistor))+
  geom_point(col="#2476C4")+
  geom_smooth(method="lm",se=F,col="black")+
  labs(x="ln(Vector) cm/s",y="Ambient Termistor (V)")+
  theme_classic()+
  facet_wrap(vars(Instrument),scales="free")



ExCalibrationPlots
AmCalibrationPlots

```

##Calculating Partial Autocorrelation Factor
```{r Partial Autocorrelation Factor, include=FALSE}

# 1. AMBIENT

#Select entire section below to make pacf graphs for each instrument  
Instruments<-c("T89_2013","T1415_2013","T1213_2013","T89_2014","T1011_2014","T1011B_2014")
for (i in Instruments) {
Inst<-AllData%>%
  subset(Instrument==i)

pacf(Inst$Am,lag=100,main="",plot=T)
title(paste(i," Ambient"),line=-2)
}

#At a lag of about 20 the pacf decreases below 0.05 in for all instruments for Ambient


# 2. EXCURRENT


#Select entire section below to make pacf graphs for each instrument  
Instruments<-c("T89_2013","T1415_2013","T1213_2013","T89_2014","T1011_2014","T1011B_2014")
for (i in Instruments) {
Inst<-subset(AllData,Instrument==i)

pacf(Inst$Ex,lag=100,main="",plot=T)
title(paste(i," Excurrent"),line=-2)
}

#At a lag of about 20 the pacf decreases below 0.05 in for all instruments for Excurrent

#. OXYGEN REMOVAL 

#Select entire section below to make pacf graphs for each instrument  
Instruments<-c("T89_2013","T1415_2013","T1213_2013","T89_2014","T1011_2014","T1011B_2014")
for (i in Instruments) {
Inst<-subset(AllData,Instrument==i)

pacf(Inst$Oxy,lag=100,main="",plot=T)
title(paste(i," Oxy"),line=-2)
}

#At a lag of less than 20 the pacf decreases below 0.05 in for all instruments for Oxygen

# a subset of every 20 points was chosen for independent samples. 
```

Confidence interval calculation for each instrument by ambient bins
```{r Confidence Intervals,include=FALSE}

ACFcorrectedMeans<-AllData%>%
  filter(row_number() %% 20 == 1)%>%
  group_by(Instrument,AmBin)%>%
  summarise(MeanOxy=mean(Oxy),
            SdOxy=sd(Oxy),
            MeanEx=mean(Ex),
            SdEx=sd(Ex),
            MeanResp=mean(RespirationRate),
            SdResp=sd(RespirationRate),
            n=n())

#Calculating Confidence Intervals for a t distribution
ACFcorrectedMeans$OxyCI<-qt(0.975,df=ACFcorrectedMeans$n-1)*
  ACFcorrectedMeans$SdOxy/sqrt(ACFcorrectedMeans$n)

ACFcorrectedMeans$ExCI<-qt(0.975,df=ACFcorrectedMeans$n-1)*
  ACFcorrectedMeans$SdEx/sqrt(ACFcorrectedMeans$n)

ACFcorrectedMeans$RespCI<-qt(0.975,df=ACFcorrectedMeans$n-1)*
  ACFcorrectedMeans$SdResp/sqrt(ACFcorrectedMeans$n)

ExByBinsACFcorrencted<-ggplot(subset(ACFcorrectedMeans,n>20),aes(x=AmBin,y=MeanEx))+
  geom_ribbon(aes(ymax=MeanEx+ExCI,ymin=MeanEx-ExCI),fill="lightgrey")+
  geom_line( col="black")+
  geom_point(col="black")+
  theme_classic()+
  scale_colour_brewer(palette = "Dark2")+
  scale_fill_brewer(  palette = "Dark2")+
  facet_wrap(vars(Instrument),scales = "free")+
  labs(x=expression('Ambient Flow Rate (cm s'^-1*')'),
       y=expression('Osculum Jet Speed (cm s'^-1*')'),tag="A.")+
  theme(legend.position="none",
        strip.background = element_blank(),
        axis.title.x = element_blank())+
  guides(col = guide_legend(ncol = 8))

OxyByBinsACFcorrencted<-ggplot(subset(ACFcorrectedMeans,n>20),aes(x=AmBin,y=MeanOxy))+
  geom_ribbon(aes(ymax=MeanOxy+OxyCI,ymin=MeanOxy-OxyCI),fill="lightgrey")+
  geom_line( col="black")+
  geom_point(col="black")+
  theme_classic()+
  scale_colour_brewer(palette = "Dark2")+
  scale_fill_brewer(  palette = "Dark2")+
  facet_wrap(vars(Instrument),scales = "free")+
  labs(x=expression('Ambient Flow Rate (cm s'^-1*')'),
       y=expression('Oxygen Removal (µmol O'[2]*' L'[p]*''^-1*')'),tag="B.")+
  theme(legend.position="none",
        strip.background = element_blank(),
        axis.title.x = element_blank())+
  guides(col = guide_legend(ncol = 8))

ExByBinsACFcorrencted/OxyByBinsACFcorrencted

```

#Extra tid bits (probably to be deleted)

```{r Summary Table with SD, warning=FALSE}

SummaryTable1<-AllData%>%
  subset(Instrument!="NA"&Ex>0)%>%
  group_by(Instrument)%>%
  summarise("Time Span (h)"= round(difftime(max(Time), min(Time), unit = "hours"),2),
            "Excurrent Median (cm/s)"=round(median(Ex),2),
            "Excurrent SD"=round(sd(Ex),2),
            "Oxygen Median (umol)"=round(median(Oxy),2),
            "Oxygen SD"=round(sd(Oxy),2))

#Adding Sponge Dimensions
SummaryTable1b<-SummaryTable1%>%
  add_column("Diameter(cm)"=c(3.6,2.6,3.5,4.8,8.3,7.6),.after="Instrument")%>%
  add_column("Cross Sectional Area (cm2)"=c(10.18,5.3,9.58,17.8,53.65,45.1),.after="Diameter(cm)")%>%
  add_column("Height Above Bottom (cm)"=c(33.1,38.1,24,39.5,54.1,33.7),
           .after="Cross Sectional Area (cm2)")

#Adding Oscular Flow Rate
SummaryTable1c<-SummaryTable1b%>%
  add_column("Oscular Flow Rate (L/s)"=round((SummaryTable1b$`Excurrent Median (cm/s)`
                                             *SummaryTable1b$`Cross Sectional Area (cm2)`)
                                             /1000,3),.after="Excurrent SD")%>%
  add_column("Oscular Flow Rate SD (L/s)"=round((SummaryTable1b$`Excurrent SD`
                                             *SummaryTable1b$`Cross Sectional Area (cm2)`)/1000,3),.after="Oscular Flow Rate (L/s)")
  

kable_styling(kable(SummaryTable1c,align = "c"),bootstrap_options = c("striped", "hover"))

write.csv(SummaryTable1c,"./Tables/Table 1-Summary.csv")

```

###*Calculating propagated error*

Equation used was based on http://ipl.physics.harvard.edu/wp-uploads/2013/03/PS3_Error_Propagation_sp13.pdf

Equation used is: 

$\sqrt{\sum{SD^2}}$

```{r Propagated Error, warning=FALSE,message=FALSE}

#Exurrent mean 
ExMean<-mean(SummaryTable1c$`Excurrent Median (cm/s)`)

#Exurrent error
ExPropErr<-SummaryTable1c$`Excurrent Median (cm/s)` %>% 
  sapply(function(x) x^2)%>%
  sum() %>%
  sqrt()

#Oxygen removal mean 
OxyMean<-mean(SummaryTable1c$`Oxygen Median (umol)`)

#Oxygen Removal propagated error
OxyPropErr<-SummaryTable1c$`Oxygen Median (umol)`%>% 
  sapply(function(x) x^2)%>%
  sum() %>%
  sqrt()

#Oscular Flow Rate Mean
OscFlowRateMean<-mean(SummaryTable1c$`Oscular Flow Rate (L/s)`)

#Oscular flow rate propagated error
OscFlowRatePropErr<-SummaryTable1c$`Oscular Flow Rate SD (L/s)`%>% 
  sapply(function(x) x^2)%>%
  sum() %>%
  sqrt()

cat("Excurrent ± PropErr =",round(ExMean,2),"±",round(ExPropErr,2))

cat("Oxygen Removal ± PropErr =",round(OxyMean,2),"±",round(OxyPropErr,2))

cat("Oscular Flow Rate ± PropErr =",round(OscFlowRateMean,3),"±",round(OscFlowRatePropErr,2))
```
