---
title: "Environmental Data NMDS"
author: "Yaamini Venkataraman"
date: "11/30/2018"
output: html_document
---

In this script, I'll visualize differences in my environmental data between sites using a nonmetric multidimensional scaling (NMDS) and analysis of similarities (ANOSIM).

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

# Load dependencies

```{r}
#install.packages("vegan")
require(vegan)
#install.packages("cluster")
require(cluster)

source("../../biostats.R") #Load biostats file
```

```{r}
sessionInfo()
```

# Import and reformat data

## Tide data

```{r}
tideData <- read.csv("../../../../data/DNR/2017-12-13-Tidal-Data-by-Site.csv", header = TRUE, strip.white = TRUE) #Import the tide data
tideData$Date <- as.Date(tideData$Date, format = "%m/%d/%y") #Convert entries to dates
tideData$DateTime <- paste(tideData$Date, tideData$Time) #Create new DateTime column to easily merge tide and environmental data
colnames(tideData) <- c("Date", "Time", "CI-Tide", "FB-Tide", "PG-Tide", "SK-Tide", "WB-Tide", "DateTime")
head(tideData) #Confirm changes
```

## pH data

```{r}
pHDOData <- read.csv("../../../../data/DNR/2017-11-14-Environmental-Data-from-Micah.csv", header = TRUE, na.strings = "NA") #Import file with pH and DO data
colnames(pHDOData) #View column names
```

```{r}
pHData <- pHDOData[,c(2:3, 4:12)] #Subset only the bare and eelgrass outplant pH data
head(pHData) #Confirm subset
```

```{r}
colnames(pHData) <- c("Date", "Time", "WBE-pH", "WBB-pH", "SKE-pH", "SKB-pH", "PGB-pH", "CIE-pH", "CIB-pH", "FBE-pH", "FBB-pH") #Rename columns
pHData$Date <- as.Date(pHData$Date, format = "%m/%d/%y") #Convert entries to dates
pHData <- pHData[pHData$Date >= "2016-06-19", ] #The outplant only started 6/19. I need to remove all data points before this time.
pHData$DateTime <- paste(pHData$Date, pHData$Time) #Create new DateTime column to easily merge tide and environmental data
head(pHData) #Confirm changes
```

## Dissolved oxygen data

```{r}
DOData <- pHDOData[,c(2:3, 22:31)] #Subset the bare and eelgrass outplant DO data
head(DOData) #Confirm subset
```

```{r}
colnames(DOData) <- c("Date", "Time", "WBE-DO", "WBB-DO", "SKE-DO", "SKB-DO", "PGE-DO", "PGB-DO", "CIE-DO", "CIB-DO", "FBE-DO", "FBB-DO") #Rename columns
DOData$Date <- as.Date(DOData$Date, format = "%m/%d/%y") #Convert entries to dates
DOData <- DOData[DOData$Date >= "2016-06-19", ] #The outplant only started 6/19. I need to remove all data points before this time.
DOData$DateTime <- paste(DOData$Date, DOData$Time) #Create new DateTime column to easily merge tide and environmental data
head(DOData) #Confirm changes
```

## Salinity

```{r}
salinityData <- read.csv("../../../../data/DNR/2018-05-30-Fixed-Salinity-from-Micah.csv", header = TRUE, na.strings = "NA", strip.white = TRUE) #Import salinity data and remove white space from end of Date and Time columns
salinityData <- salinityData[,c(1:2, 4, 6, 8, 10, 12, 14, 16, 20)] #Subset salinity from bare and eelgrass outplants. Needed to use PGE instead of PGB since PGB has no salinity data. Also use FBE and WBE due to probe burial at bare sites.
head(salinityData) #Confirm subset
```

```{r}
colnames(salinityData) <- c("Date", "Time", "CIB-Salinity", "CIE-Salinity", "FBB-Salinity", "FBE-Salinity", "PGE-Salinity", "SKB-Salinity", "SKE-Salinity", "WBB-Salinity") #Rename columns
salinityData$Date <- as.Date(salinityData$Date, format = "%d/%m/%Y") #Convert entries to dates
salinityData <- salinityData[salinityData$Date >= "2016-06-19", ] #The outplant only started 6/19. I need to remove all data points before this time.
salinityData$DateTime <- paste(salinityData$Date, salinityData$Time) #Create new DateTime column to easily merge tide and environmental data.
head(salinityData)
```

## Temperature data

```{r}
temperatureData <- read.csv("../../../../data/DNR/2017-11-14-Environmental-Data-from-Micah.csv", header = TRUE, na.strings = "NA") #Import temperature data
colnames(temperatureData) #Get column names
```

```{r}
temperatureData <- temperatureData[,c(2:3, 32:41)] #Save temperature data from bare outplants as a new dataframe
head(temperatureData) #Confirm subset
```

```{r}
colnames(temperatureData) <- c("Date", "Time", "WBE", "WBB", "SKE", "SKB", "PGE", "PGB", "CIE", "CIB", "FBE", "FBB") #Rename columns
temperatureData$Date <- as.Date(temperatureData$Date, format = "%m/%d/%y") #Convert entries to dates
temperatureData <- temperatureData[temperatureData$Date >= "2016-06-19", ] #The outplant only started 6/19. I need to remove all data points before this time.
temperatureData$DateTime <- paste(temperatureData$Date, temperatureData$Time) #Create new DateTime column to easily merge tide and environmental data
head(temperatureData) #Confirm changes
```

```{r}
temperatureData <- temperatureData[-c(4897:4898), ] #Remove blank two rows at the end of the data
tail(temperatureData) #Confirm changes
```

# Merge tidal data with environmental data

Tide correction only needs to occur for pH, dissolved oxygen, and salinity data. Temperature data from low tide conditions is still an accurate representation of the oysters' environment.

## pH data

```{r}
pHTideData <- merge(x = pHData, y = tideData, by = "DateTime") #Merge pH and tide data
colnames(pHTideData) #Get column names
```

```{r}
pHTideData <- pHTideData[, -c(13:14)] #Remove redundant date and time columns
colnames(pHTideData) <- c("DateTime", "Date", "Time", "WBE-pH", "WBB-pH", "SKE-pH", "SKB-pH", "PGB-pH", "CIE-pH", "CIB-pH", "FBE-pH", "FBB-pH", "CI-Tide", "FB-Tide", "PG-Tide", "SK-Tide", "WB-Tide") #Change column names
head(pHTideData) #Confirm changes
```

## Dissolved oxygen 

```{r}
DOTideData <- merge(x = DOData, y = tideData, by = "DateTime") #Merge DO and tide data
colnames(DOTideData) #Get column names
```

```{r}
DOTideData <- DOTideData[, -c(14:15)] #Remove redundant date and time columns
colnames(DOTideData) <- c("DateTime", "Date", "Time", "WBE-DO", "WBB-DO", "SKE-DO", "SKB-DO", "PGE-DO", "PGB-DO", "CIE-DO", "CIB-DO", "FBE-DO", "FBB-DO", "CI-Tide", "FB-Tide", "PG-Tide", "SK-Tide", "WB-Tide") #Change column names
head(DOTideData) #Confirm changes
```

## Salinity

```{r}
salinityTideData <- merge(x = salinityData, y = tideData, by = "DateTime") #Merge salinity and tide data
colnames(salinityTideData) #Get column names
```

```{r}
salinityTideData <- salinityTideData[, -c(12:13)] #Remove redundant date and time columns
colnames(salinityTideData) <- c("DateTime", "Date", "Time", "CIB-Salinity", "CIE-Salinity", "FBB-Salinity", "FBE-Salinity", "PGE-Salinity", "SKB-Salinity", "SKE-Salinity", "WBB-Salinity", "CI-Tide", "FB-Tide", "PG-Tide", "SK-Tide", "WB-Tide") #Change column names
head(salinityTideData) #Confirm changes
```

# Remove exposure times

## pH

```{r}
colnames(pHTideData)
```


```{r}
pHTideData$`WBE-pH`[pHTideData$`WB-Tide` <= 1] <- NA #Replace WBE-pH values with "NA" when tide is less than 1
pHTideData$`WBB-pH`[pHTideData$`WB-Tide` <= 1] <- NA #Replace WBB-pH values with "NA" when tide is less than 1

pHTideData$`SKE-pH`[pHTideData$`SK-Tide` <= 1] <- NA #Replace SKE-pH values with "NA" when tide is less than 1
pHTideData$`SKB-pH`[pHTideData$`SK-Tide` <= 1] <- NA #Replace SKB-pH values with "NA" when tide is less than 1

pHTideData$`PGB-pH`[pHTideData$`PG-Tide` <= 1] <- NA #Replace PGB-pH values with "NA" when tide is less than 1

pHTideData$`CIE-pH`[pHTideData$`CI-Tide` <= 1] <- NA #Replace CIE-pH values with "NA" when tide is less than 1
pHTideData$`CIB-pH`[pHTideData$`CI-Tide` <= 1] <- NA #Replace CIB-pH values with "NA" when tide is less than 1

pHTideData$`FBE-pH`[pHTideData$`FB-Tide` <= 1] <- NA #Replace FBE-pH values with "NA" when tide is less than 1
pHTideData$`FBB-pH`[pHTideData$`FB-Tide` <= 1] <- NA #Replace FBB-pH values with "NA" when tide is less than 1
```

```{r}
#Convert to numeric values

pHTideData$`WBE-pH` <- as.numeric(pHTideData$`WBE-pH`)
pHTideData$`WBB-pH` <- as.numeric(pHTideData$`WBB-pH`)

pHTideData$`SKE-pH` <- as.numeric(pHTideData$`SKE-pH`)
pHTideData$`SKB-pH` <- as.numeric(pHTideData$`SKB-pH`)

pHTideData$`PGB-pH` <- as.numeric(pHTideData$`PGB-pH`)

pHTideData$`CIE-pH` <- as.numeric(pHTideData$`CIE-pH`)
pHTideData$`CIB-pH` <- as.numeric(pHTideData$`CIB-pH`)

pHTideData$`FBE-pH` <- as.numeric(pHTideData$`FBE-pH`)
pHTideData$`FBB-pH` <- as.numeric(pHTideData$`FBB-pH`)
```

## Dissolved oxygen

```{r}
colnames(DOTideData)
```


```{r}
DOTideData$`WBE-DO`[DOTideData$`WB-Tide` <= 1] <- NA #Replace WBE-DO values with "NA" when tide is less than 1
DOTideData$`WBB-DO`[DOTideData$`WB-Tide` <= 1] <- NA #Replace WBB-DO values with "NA" when tide is less than 1

DOTideData$`SKE-DO`[DOTideData$`SK-Tide` <= 1] <- NA #Replace SKE-DO values with "NA" when tide is less than 1
DOTideData$`SKB-DO`[DOTideData$`SK-Tide` <= 1] <- NA #Replace SKB-DO values with "NA" when tide is less than 1

DOTideData$`PGE-DO`[DOTideData$`PG-Tide` <= 1] <- NA #Replace PGE-DO values with "NA" when tide is less than 1
DOTideData$`PGB-DO`[DOTideData$`PG-Tide` <= 1] <- NA #Replace PGB-DO values with "NA" when tide is less than 1

DOTideData$`CIE-DO`[DOTideData$`CI-Tide` <= 1] <- NA #Replace CIE-DO values with "NA" when tide is less than 1
DOTideData$`CIB-DO`[DOTideData$`CI-Tide` <= 1] <- NA #Replace CIB-DO values with "NA" when tide is less than 1

DOTideData$`FBE-DO`[DOTideData$`FB-Tide` <= 1] <- NA #Replace FBE-DO values with "NA" when tide is less than 1
DOTideData$`FBB-DO`[DOTideData$`FB-Tide` <= 1] <- NA #Replace FBB-DO values with "NA" when tide is less than 1
```

```{r}
#Convert to numeric values

DOTideData$`WBE-DO` <- as.numeric(DOTideData$`WBE-DO`)
DOTideData$`WBB-DO` <- as.numeric(DOTideData$`WBB-DO`)

DOTideData$`SKE-DO` <- as.numeric(DOTideData$`SKE-DO`)
DOTideData$`SKB-DO` <- as.numeric(DOTideData$`SKB-DO`)

DOTideData$`PGE-DO` <- as.numeric(DOTideData$`PGE-DO`)
DOTideData$`PGB-DO` <- as.numeric(DOTideData$`PGB-DO`)

DOTideData$`CIE-DO` <- as.numeric(DOTideData$`CIE-DO`)
DOTideData$`CIB-DO` <- as.numeric(DOTideData$`CIB-DO`)

DOTideData$`FBE-DO` <- as.numeric(DOTideData$`FBE-DO`)
DOTideData$`FBB-DO` <- as.numeric(DOTideData$`FBB-DO`)
```

## Salinity

```{r}
colnames(salinityTideData)
```

```{r}
salinityTideData$`CIB-Salinity`[salinityTideData$`CIB-Salinity` <= 1] <- NA #Replace CIB-Salinity values with "NA" when tide is less than 1
salinityTideData$`CIE-Salinity`[salinityTideData$`CIB-Salinity` <= 1] <- NA #Replace CIB-Salinity values with "NA" when tide is less than 1

salinityTideData$`FBB-Salinity`[salinityTideData$`FB-Salinity` <= 1] <- NA #Replace FBB-Salinity values with "NA" when tide is less than 1
salinityTideData$`FBE-Salinity`[salinityTideData$`FB-Salinity` <= 1] <- NA #Replace FBE-Salinity values with "NA" when tide is less than 1

salinityTideData$`PGE-Salinity`[salinityTideData$`PG-Salinity` <= 1] <- NA #Replace PGE-Salinity values with "NA" when tide is less than 1

salinityTideData$`SKB-Salinity`[salinityTideData$`SK-Salinity` <= 1] <- NA #Replace SKB-Salinity values with "NA" when tide is less than 1
salinityTideData$`SKE-Salinity`[salinityTideData$`SK-Salinity` <= 1] <- NA #Replace SKE-Salinity values with "NA" when tide is less than 1

salinityTideData$`WBB-Salinity`[salinityTideData$`WB-Salinity` <= 1] <- NA #Replace WBB-Salinity values with "NA" when tide is less than 1
```

```{r}
#Convert to numeric values
salinityTideData$`CIB-Salinity` <- as.numeric(salinityTideData$`CIB-Salinity`)
salinityTideData$`CIE-Salinity` <- as.numeric(salinityTideData$`CIE-Salinity`)

salinityTideData$`FBB-Salinity` <- as.numeric(salinityTideData$`FBB-Salinity`)
salinityTideData$`FBE-Salinity` <- as.numeric(salinityTideData$`FBE-Salinity`)

salinityTideData$`PGE-Salinity` <- as.numeric(salinityTideData$`PGE-Salinity`)

salinityTideData$`SKB-Salinity` <- as.numeric(salinityTideData$`SKB-Salinity`)
salinityTideData$`SKE-Salinity` <- as.numeric(salinityTideData$`SKE-Salinity`)

salinityTideData$`WBB-Salinity` <- as.numeric(salinityTideData$`WBB-Salinity`)
```

# Remove outliers

## pH

```{r}
nSites <- 17 #Sites are from columns 4 to 17
for(i in 4:nSites) { #For individual site-habitat data
  upperBound <- as.numeric((quantile(pHTideData[, i], na.rm = TRUE)[4]) + (1.5*(quantile(pHTideData[, i], na.rm = TRUE)[4] - quantile(pHTideData[, i], na.rm = TRUE)[2]))) #Calculate upper bound
  lowerBound <- as.numeric((quantile(pHTideData[, i], na.rm = TRUE)[2]) - (1.5*(quantile(pHTideData[, i], na.rm = TRUE)[4] - quantile(pHTideData[, i], na.rm = TRUE)[2]))) #Calculate lower bound
  pHTideData[, i][pHTideData[, i] > upperBound] <- NA #Replace any values higher than upper bound with NA
  pHTideData[, i][pHTideData[, i] < lowerBound] <- NA #Replace any values lower than upper bound with NA
} #Replace outliers with NA values
```

## Dissolved oxygen

```{r}
nSites <- 18 #Sites are from columns 4 to 18
for(i in 4:nSites) { #For individual site data
  upperBound <- as.numeric((quantile(DOTideData[, i], na.rm = TRUE)[4]) + (1.5*(quantile(DOTideData[, i], na.rm = TRUE)[4] - quantile(DOTideData[, i], na.rm = TRUE)[2]))) #Calculate upper bound
  lowerBound <- 0 #Dissolved oxygen content cannot be less than zero
  DOTideData[, i][DOTideData[, i] > upperBound] <- NA #Replace any values higher than upper bound with NA
  DOTideData[, i][DOTideData[, i] < lowerBound] <- NA #Replace any values lower than upper bound with NA
} #Replace outliers with NA values
```

## Salinity

```{r}
nSites <- 16 #Sites are from columns 4 to 16
for(i in 4:nSites) { #For individual site data
  upperBound <- as.numeric((quantile(salinityTideData[, i], na.rm = TRUE)[4]) + (1.5*(quantile(salinityTideData[, i], na.rm = TRUE)[4] - quantile(salinityTideData[, i], na.rm = TRUE)[2]))) #Calculate upper bound
  lowerBound <- as.numeric((quantile(salinityTideData[, i], na.rm = TRUE)[2]) - (1.5*(quantile(salinityTideData[, i], na.rm = TRUE)[4] - quantile(salinityTideData[, i], na.rm = TRUE)[2]))) #Calculate lower bound
  salinityTideData[, i][salinityTideData[, i] > upperBound] <- NA #Replace any values higher than upper bound with NA
  salinityTideData[, i][salinityTideData[, i] < lowerBound] <- NA #Replace any values lower than upper bound with NA
} #Replace outliers with NA values
```

#Export corrected data

```{r}
#write.csv(pHTideData, "2018-11-18-pH-Tide-Data-Corrected.csv") #Export pH data
#write.csv(DOTideData, "2018-11-18-DO-Tide-Data-Corrected.csv") #Export DO data
#write.csv(salinityTideData, "2018-11-18-Salinity-Tide-Data-Corrected.csv") #Export salinity data
#write.csv(temperatureData, "2018-11-18-Temperature-Corrected.csv") #Export temperature data
```

# Obtain daily mean values for environmental data

```{r}
outplantDates <- as.Date(unique(temperatureData$Date)) #Save outplant dates as a new vector
outplantDates #Confirm vector creation
```

## pH

```{r}
head(pHTideData)
```

```{r}
pHMean <- data.frame("Date" = outplantDates,
                     "WBE-pH-mean" = rep(0, length(outplantDates)),
                     "WBB-pH-mean" = rep(0, length(outplantDates)),
                     "SKE-pH-mean" = rep(0, length(outplantDates)),
                     "SKB-pH-mean" = rep(0, length(outplantDates)),
                     "PGB-pH-mean" = rep(0, length(outplantDates)),
                     "CIE-pH-mean" = rep(0, length(outplantDates)),
                     "CIB-pH-mean" = rep(0, length(outplantDates)),
                     "FBE-pH-mean" = rep(0, length(outplantDates)),
                     "FBB-pH-mean"  = rep(0, length(outplantDates))) #Create an empty dataframe
head(pHMean)
```

```{r}
nDates <- length(outplantDates)
for (j in 4:12) { #For each column with pH data
  for(i in 1:nDates) { #For each date
    pHMean[i, j-2] <- mean(pHTideData[pHTideData$Date == outplantDates[i], j], na.rm = TRUE) #Find and save the mean value of the designated column at each date.
  }
} #The loop will cycle through each date first in one column, then move to the next column. Any days with no non-missing data will save as "NaN" in the dataframe.
is.na(pHMean) <- sapply(pHMean, is.nan) #Replace all "NaN" with N/A
head(pHMean) #Confirm dataframe filling. All "NaN" values are now replaced with N/A.
```

## Dissolved oxygen

```{r}
head(DOTideData)
```

```{r}
DOMean <- data.frame("Date" = outplantDates,
                        "WBE-DO-mean" = rep(0, length(outplantDates)),
                        "WBB-DO-mean" = rep(0, length(outplantDates)),
                        "SKE-DP-mean" = rep(0, length(outplantDates)),
                        "SKB-DO-mean" = rep(0, length(outplantDates)),
                        "PGE-DO-mean" = rep(0, length(outplantDates)),
                        "PGB-DO-mean" = rep(0, length(outplantDates)),
                        "CIE-DO-mean" = rep(0, length(outplantDates)),
                        "CIB-DO-mean" = rep(0, length(outplantDates)),
                        "FBE-DO-mean" = rep(0, length(outplantDates)),
                        "FBB-DO-mean" = rep(0, length(outplantDates))) #Create an empty dataframe
head(DOMean)
```

```{r}
for (j in 4:13) { #For each column with DO data
  for(i in 1:nDates) { #For each date
    DOMean[i, j-2] <- mean(DOTideData[DOTideData$Date == outplantDates[i], j], na.rm = TRUE) #Find and save the mean value of the designated column at each date.
  }
} #The loop will cycle through each date first in one column, then move to the next column. Any days with no non-missing data will save as "NaN" in the dataframe.
is.na(DOMean) <- sapply(DOMean, is.nan) #Replace all "NaN" with N/A
head(DOMean) #Confirm dataframe filling. All "Inf" values are now replaced
```

## Salinity

```{r}
head(salinityTideData)
```

```{r}
salinityMean <- data.frame("Date" = outplantDates,
                           "CIB-sal-mean" = rep(0, length(outplantDates)),
                           "CIE-sal-mean" = rep(0, length(outplantDates)),
                           "FBB-sal-mean" = rep(0, length(outplantDates)),
                           "FBE-sal-mean" = rep(0, length(outplantDates)),
                           "PGE-sal-mean" = rep(0, length(outplantDates)),
                           "SKB-sal-mean" = rep(0, length(outplantDates)),
                           "SKE-sal-mean" = rep(0, length(outplantDates)),
                           "WBB-sal-mean" = rep(0, length(outplantDates))) #Create an empty dataframe
head(salinityMean)
```

```{r}
for (j in 4:11) { #For each column with salinity data
  for(i in 1:nDates) { #For each date
    salinityMean[i, j-2] <- mean(salinityTideData[salinityTideData$Date == outplantDates[i], j], na.rm = TRUE) #Find and save the mean value of the designated column at each date.
  }
} #The loop will cycle through each date first in one column, then move to the next column. Any days with no non-missing data will save as "NaN" in the dataframe.
is.na(salinityMean) <- sapply(salinityMean, is.nan) #Replace all "NaN" with N/A
head(salinityMean) #Confirm dataframe filling. All "NaN" values are now replaced
```

## Temperature

```{r}
head(temperatureData)
```

```{r}
outplantDateCharacters <- as.character(outplantDates) #Convert outplant dates to characters
temperatureData$DateCharacters <- as.character(temperatureData$Date) #Convert dates to characters
```

```{r}
temperatureMean <- data.frame("Date" = outplantDateCharacters,
                              "WBE-temp-mean" = rep(0, length(outplantDates)),
                              "WBB-temp-mean" = rep(0, length(outplantDates)),
                              "SKE-temp-mean" = rep(0, length(outplantDates)),
                              "SKB-temp-mean" = rep(0, length(outplantDates)),
                              "PGE-temp-mean" = rep(0, length(outplantDates)),
                              "PGB-temp-mean" = rep(0, length(outplantDates)),
                              "CIE-temp-mean" = rep(0, length(outplantDates)),
                              "CIB-temp-mean" = rep(0, length(outplantDates)),
                              "FBE-temp-mean" = rep(0, length(outplantDates)),
                              "FBB-temp-mean" = rep(0, length(outplantDates))) #Create an empty dataframe, but be sure to use outplantDateCharacters
head(temperatureMean)
```

```{r}
for (j in 3:12) { #For each column with temperature data
  for(i in 1:nDates) { #For each date
    temperatureMean[i, j-1] <- mean(temperatureData[temperatureData$DateCharacters == outplantDateCharacters[i], j], na.rm = TRUE) #Find and save the mean value of the designated column at each date.
  }
} #The loop will cycle through each date first in one column, then move to the next column. Any days with no non-missing data will save as "Inf" in the dataframe.
is.na(temperatureMean) <- sapply(temperatureMean, is.nan) #Replace all "NaN" with N/A
head(temperatureMean) #Confirm dataframe filling. All "NaN" values are now replaced
```

```{r}
temperatureMean$Date <- outplantDates #Replace temperatureMean$Date with outplantDates so they are dates 
head(temperatureMean)
```

## Merge all dataframes together
I'm going to ordinate means and variances separately, so they need to be in seaprate dataframes.

```{r}
dailyMeanData <- cbind(pHMean, DOMean, salinityMean, temperatureMean) #Use cbind to merge all of the dataframes together
colnames(dailyMeanData) #View column names
```

```{r}
dailyMeanData <- dailyMeanData[,-c(11, 22, 31)] #Remove all redundant "Date" columns
colnames(dailyMeanData) #Confirm column removal
```

# Environmental data NMDS

## Transform data

```{r}
dailyMeanDataCorrected <- dailyMeanData #Duplicate dataframe
rownames(dailyMeanDataCorrected) <- dailyMeanDataCorrected$Date #Save dates as rownames
dailyMeanDataCorrected <- dailyMeanDataCorrected[, -1] #Remove date column
head(dailyMeanDataCorrected) #Confirm changes
```

```{r}
meanData.log <- log(dailyMeanDataCorrected + 1) #Log transform data before analysis
meanData.log <- meanData.log[-c(32:34),] #Remove last three rows since the outplant ended before then
tail(meanData.log)
```

```{r}
meanData.log.trans <- data.frame(t(meanData.log)) #Transpose dataframe so objects are environmental conditions at site-habitat combinations, and descriptors are dates. This is the data I will use for the NMDS.
colnames(meanData.log.trans) <- rownames(meanData.log) #Rename columns so there is no "X" in front of dates
head(meanData.log.trans) #View transformation
```

```{r}
#write.csv(meanData.log.trans, "2018-12-01-Mean-Environmental-Data-Log-Transformed.csv", row.names = TRUE) #Write out transformed data for future analyses
```


## Calculate dissimilarity matrix

```{r}
meanDataGowerDiss <- daisy(meanData.log.trans, metric = "gower") #Calculate a dissimilarity (distance) matrix based on Gower's coefficient. Since vegan cannot handle missing values, use the daisy function from the cluster library to calculate the dissimilarity matrix to use in all vegan functions.
meanDataGowerDiss #Confirm matrix calculation
```

## Conduct the NMDS

```{r}
nmds.scree(meanDataGowerDiss, distance = "euclidean", k = 10, autotransform = FALSE, trymax = 20) #Create a screeplot to compare the stress for solutions across different k values from 2 to 10. Use 20 different random start configurations. As the number of ordination axes increases, stress is minimized because the NMDS algorithm is trying to represent p dimensional data in k dimensions. There is an elbow at 2 ordination axes, which tells me that 2 axes should suffice for the analysis.
```

```{r}
meanData.log.gower.NMDS <- metaMDS(meanDataGowerDiss, distance = 'euclidean', k = 2, trymax = 10000, autotransform = FALSE) #Make MDS dissimilarity matrix on log transformed data using Gower's.
meanData.log.gower.NMDS$stress #Stress = 0.01740668
```

```{r}
#Use code from Jennifer Gardner to define a new function, nmds.montejg, to calculate a p-value with a permutation test for missing data.

nmds.montejg <-
  function(dissim,rawdata,k,nmds.distance='bray',daisy.distance='gower', trymax=50,autotransform=FALSE,
           trace=0,zerodist='add',perm=99,col.hist='blue',col.line='red',
           lty=2,las=1,lab=c(5,5,4),...){
    
    library(vegan)
    library(MASS)
    
    z<-metaMDS(comm=dissim,k=k,distance=nmds.distance,trymax=trymax,
               autotransform=autotransform,trace=trace,...) #nmds analysis on dissimilarity matrix
    z.stress<-z$stress #get stress
    y.stress<-rep(0,perm)
    
    for(i in 1:perm){
      y<-apply(rawdata,2,sample) #permute raw data matrix
      y<- daisy(y, metric = daisy.distance) # calculate gower's for permuted data
      y<-metaMDS(comm=y,k=k,distance=nmds.distance,trymax=trymax,
                 autotransform=autotransform,trace=trace,...) #nmds analysis on dissimilarity matrix
      y.stress[i]<-y$stress #get stress
    }
    n<-sum(y.stress<=z.stress) #compute number of random runs with stress < observed
    p.value<-(1+n)/(1+perm) #compute p-value
    
    xmin<-min(z.stress,min(y.stress))
    xmax<-max(z.stress,max(y.stress))
    hist(y.stress,col=col.hist,las=las,lab=lab,
         xaxs='i',yaxs='i',xlim=c(xmin,xmax),xlab='Stress',
         main=paste('Random Permutation Distribution of Stress for',k,'Dimensions',sep=' '),...)
    abline(v=z.stress,col=col.line,lty=lty,lwd=2,...)
    
    cat('Randomization Test of Stress:\n')
    cat('Permutation stress values:\n')
    print(y.stress)
    z<-rbind('Observed stress'=z.stress,'P-value'=p.value,'Times Permute Less Than Observed'=n)
    return(z)
  }
```


```{r}
nmds.montejg(dissim = meanDataGowerDiss, rawdata = meanData.log.trans, nmds.distance = 'euclidean', k = 2, autotransform = FALSE, trymax = 20) #Perform a randomization test to determine if the solution for 2 dimensions is significant. The observed stress value is less than the expected stress value. P-value = 0.01000000
```

```{r}
stressplot(meanData.log.gower.NMDS) #Make Shepard plot to visualize the relationship between original dissimilarities (distance matrix) and distnaces in ordination space. The non-metric R-squared value is 1 (redundant with observed stress value and p-value from the randomization test)
```

```{r}
vec.meanData.log.gower <- envfit(meanData.log.gower.NMDS$points, meanData.log.trans, perm = 1000, na.rm = TRUE) #Calculate loadings by correlating NMDS scores with original variables
vec.meanData.log.gower #Look at loadings. All are equally significant.
```

```{r}
ordiplot(meanData.log.gower.NMDS, choices = c(1,2), type = "text", display = "sites", xlab = "Axis 1", ylab = "Axis 2") #Plot basic NMDS
plot(vec.meanData.log.gower, p.max = 0.001, col = 'blue') #Plot loadings that are significant at the 0.001 level
```

## Customize NMDS

```{r}
plotCustomization <- data.frame("Site" = rep(0, times = length(meanData.log.trans$`2016-06-19`)),
                                "Habitat" = rep(0, times = length(meanData.log.trans$`2016-06-19`)),
                                "Site-Habitat" = rep(0, times = length(meanData.log.trans$`2016-06-19`)),
                                "Shape" = rep(0, times = length(meanData.log.trans$`2016-06-19`)),
                                "Color" = rep(0, times = length(meanData.log.trans$`2016-06-19`)),
                                "Shape2" = rep(0, times = length(meanData.log.trans$`2016-06-19`))) #Create a new dataframe to store plot customization information
rownames(plotCustomization) <- rownames(meanData.log.trans) #Use rownames from NMDS data as rownames for the customization information
head(plotCustomization) #Confirm creation
```

### Site information

```{r}
plotCustomization[grep(rownames(plotCustomization), pattern = "CI"), "Site"] <- "CI" #For each rowname with "CI", add "CI" to the site column
plotCustomization[grep(rownames(plotCustomization), pattern = "FB"), "Site"] <- "FB" #For each rowname with "FB", add "FB" to the site column
plotCustomization[grep(rownames(plotCustomization), pattern = "PG"), "Site"] <- "PG" #For each rowname with "PG", add "PG" to the site column
plotCustomization[grep(rownames(plotCustomization), pattern = "SK"), "Site"] <- "SK" #For each rowname with "SK", add "SK" to the site column
plotCustomization[grep(rownames(plotCustomization), pattern = "WB"), "Site"] <- "WB" #For each rowname with "WB", add "WB" to the site column
head(plotCustomization) #Confirm changes
```

### Habitat information

```{r}
plotCustomization[grep(rownames(plotCustomization), pattern = "CIB"), "Habitat"] <- "Bare" #For each rowname with "CIB", add "Bare" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "FBB"), "Habitat"] <- "Bare" #For each rowname with "FBB", add "Bare" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "PGB"), "Habitat"] <- "Bare" #For each rowname with "PGB", add "Bare" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "SKB"), "Habitat"] <- "Bare" #For each rowname with "SKB", add "Bare" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "WBB"), "Habitat"] <- "Bare" #For each rowname with "WBB", add "Bare" to the habitat column
head(plotCustomization) #Confirm changes
```

```{r}
plotCustomization[grep(rownames(plotCustomization), pattern = "CIE"), "Habitat"] <- "Eelgrass" #For each rowname with "CIE", add "Eelgrass" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "FBE"), "Habitat"] <- "Eelgrass" #For each rowname with "FBE", add "Eelgrass" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "PGE"), "Habitat"] <- "Eelgrass" #For each rowname with "PGE", add "Eelgrass" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "SKE"), "Habitat"] <- "Eelgrass" #For each rowname with "SKE", add "Eelgrass" to the habitat column
plotCustomization[grep(rownames(plotCustomization), pattern = "WBE"), "Habitat"] <- "Eelgrass" #For each rowname with "WBE", add "Eelgrass" to the habitat column
head(plotCustomization) #Confirm changes
```

### Site and habitat information

```{r}
plotCustomization$Site.Habitat <- paste(plotCustomization$Site, "", plotCustomization$Habitat) #Combine "Site" and "Habitat" columns for Site.Habitat
head(plotCustomization) #Confirm changes
```

### Shape assignment

#### Bare vs. Eelgrass

```{r}
plotCustomization[grep(plotCustomization$Habitat, pattern = "Eelgrass"), "Shape"] <- 16 #Add 16 (filled circles) to the Shape column for each instance of "Eelgrass" in the Habitat column
plotCustomization[grep(plotCustomization$Habitat, pattern = "Bare"), "Shape"] <- 1 #Add 1 (open circles) to the Shape column for each instance of "Bare" in the Habitat column
head(plotCustomization) #Confirm changes
```

#### Habitat and environmental variable

```{r}
#pH
plotCustomization[grep(rownames(plotCustomization), pattern = "B.pH"), "Shape2"] <- 0 #Add 0 "open square" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names
plotCustomization[grep(rownames(plotCustomization), pattern = "E.pH"), "Shape2"] <- 22 #Add 22 "closed square" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names

#Dissolved oxygen
plotCustomization[grep(rownames(plotCustomization), pattern = "B.DO"), "Shape2"] <- 1 #Add 1 "open circle" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names
plotCustomization[grep(rownames(plotCustomization), pattern = "E.DO"), "Shape2"] <- 21 #Add 21 "closed circle" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names

#Salinity
plotCustomization[grep(rownames(plotCustomization), pattern = "B.sal"), "Shape2"] <- 5 #Add 5 "open diamond" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names
plotCustomization[grep(rownames(plotCustomization), pattern = "E.sal"), "Shape2"] <- 23 #Add 23 "closed diamond" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names

#Temperature
plotCustomization[grep(rownames(plotCustomization), pattern = "B.temp"), "Shape2"] <- 6 #Add 6 "open triangle" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names
plotCustomization[grep(rownames(plotCustomization), pattern = "E.temp"), "Shape2"] <- 25 #Add 25 "closed triangle" to the Shape2 column for each instance pH measurement from a bare habitat identified in the row names

tail(plotCustomization) #Confirm changes
```

### Color assignment

```{r}
plotCustomization[grep(plotCustomization$Site, pattern = "CI"), "Color"] <- "#00A9BD" #Add "#00A9BD" to the Color column for each instance of "CI" in the Site column
plotCustomization[grep(plotCustomization$Site, pattern = "FB"), "Color"] <- "#38001C" #Add "#38001C" to the Color column for each instance of "FB" in the Site column
plotCustomization[grep(plotCustomization$Site, pattern = "PG"), "Color"] <- "#440D82" #Add "#440D82" to the Color column for each instance of "PG" in the Site column
plotCustomization[grep(plotCustomization$Site, pattern = "SK"), "Color"] <- "#017A74" #Add "#017A74" to the Color column for each instance of "SK" in the Site column
plotCustomization[grep(plotCustomization$Site, pattern = "WB"), "Color"] <- "#EB8B0C" #Add "#EB8B0C" to the Color column for each instance of "WB" in the Site column
head(plotCustomization) #Confirm changes
```

### Remake NMDS

#### Site only

```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "Axis 1", ylab = "Axis 2", cex = 0.5) #Save NMDS as a new object
points(fig.nmds, "sites", col = plotCustomization$Color, pch = 16)

ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "CI", col = "#00A9BD") #Add confidence ellipse around the data from Case Inlet
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "FB", col = "#38001C") #Add confidence ellipse around the data  from Fidalgo Bay
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "PG", col = "#440D82") #Add confidence ellipse around the data from Port Gamble Bay
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "SK", col = "#017A74") #Add confidence ellipse around the data from Skokomish River Delta
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "WB", col = "#EB8B0C") #Add confidence ellipse around the data from Willapa Bay

legend("topright", pch = c(rep(x = 16, times = 5)), legend=c('Case Inlet', "Fidalgo Bay", "Port Gamble Bay", "Skokomish River Delta", "Willapa Bay"), col=c('#00A9BD', '#38001C', '#440D82', '#017A74', '#EB8B0C'), cex = 0.5, bty = "n")
```

```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "Axis 1", ylab = "Axis 2", cex = 0.5) #Save NMDS as a new object
points(fig.nmds, "sites", col = plotCustomization$Color, pch = 16)

ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "CI", col = "#00A9BD") #Add confidence ellipse around the data from Case Inlet
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "FB", col = "#38001C") #Add confidence ellipse around the data  from Fidalgo Bay
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "PG", col = "#440D82") #Add confidence ellipse around the data from Port Gamble Bay
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "SK", col = "#017A74") #Add confidence ellipse around the data from Skokomish River Delta
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "WB", col = "#EB8B0C") #Add confidence ellipse around the data from Willapa Bay

legend("topright", pch = c(rep(x = 16, times = 5)), legend=c('Case Inlet', "Fidalgo Bay", "Port Gamble Bay", "Skokomish River Delta", "Willapa Bay"), col=c('#00A9BD', '#38001C', '#440D82', '#017A74', '#EB8B0C'), cex = 0.5, bty = "n")
```

#### Site and Habitat

```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "Axis 1", ylab = "Axis 2", cex = 0.5) #Save NMDS as a new object
points(fig.nmds, "sites", col = plotCustomization$Color, pch = plotCustomization$Shape)

ordiellipse(meanData.log.gower.NMDS, plotCustomization$Habitat, show.groups = "Eelgrass", col = "grey80", lty = 1) #Add confidence ellipse around the data from eelgrass habitats
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Habitat, show.groups = "Bare", col = "grey80", lty = 2) #Add confidence ellipse around the data  from bare habitats

legend("topright", pch = c(rep(x = 1, times = 6), 16), legend=c('Case Inlet', "Fidalgo Bay", "Port Gamble Bay", "Skokomish River Delta", "Willapa Bay", "Bare", "Eelgrass"), col=c('#00A9BD', '#38001C', '#440D82', '#017A74', '#EB8B0C', "grey80", "grey80"), cex = 0.5, bty = "n")
```

```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "Axis 1", ylab = "Axis 2", cex = 0.5) #Save NMDS as a new object
points(fig.nmds, "sites", col = plotCustomization$Color, pch = plotCustomization$Shape)

ordihull(meanData.log.gower.NMDS, plotCustomization$Habitat, show.groups = "Eelgrass", col = "grey80", lty = 1) #Add confidence ellipse around the data from eelgrass habitats
ordihull(meanData.log.gower.NMDS, plotCustomization$Habitat, show.groups = "Bare", col = "grey80", lty = 2) #Add confidence ellipse around the data  from bare habitats

legend("topright", pch = c(rep(x = 1, times = 6), 16), legend=c('Case Inlet', "Fidalgo Bay", "Port Gamble Bay", "Skokomish River Delta", "Willapa Bay", "Bare", "Eelgrass"), col=c('#00A9BD', '#38001C', '#440D82', '#017A74', '#EB8B0C', "grey80", "grey80"), cex = 0.5, bty = "n")
```

#### Site, habitat, and environmental variables

```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "", ylab = "", cex = 0.5, xaxt = "n", yaxt = "n") #Save NMDS as a new object
points(fig.nmds, "sites", col = plotCustomization$Color, pch = plotCustomization$Shape2, bg = plotCustomization$Color)
axis(side = 1, labels = TRUE, col = "grey80", cex.axis = 0.75)
mtext(side = 1, text = "NMDS1", line = 2)
axis(side = 2, labels = TRUE, col = "grey80", cex.axis = 0.75)
mtext(side = 2, text = "NMDS2", line = 2)
box(col = "grey80")
legend("topleft", pch = c(rep(x = 1, times = 6), 16, 0, 1, 5, 6), legend=c('Case Inlet', "Fidalgo Bay", "Port Gamble Bay", "Skokomish River Delta", "Willapa Bay", "Bare", "Eelgrass", "pH", "Dissolved Oxygen", "Salinity", "Temperature"), col=c('#00A9BD', '#38001C', '#440D82', '#017A74', '#EB8B0C', "black", "black", "black", "black", "black", "black"), cex = 0.4, bty = "n")

ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "CI", col = "#00A9BD88") #Add confidence ellipse around the data from Case Inlet
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "FB", col = "#38001C88") #Add confidence ellipse around the data  from Fidalgo Bay
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "PG", col = "#440D8288") #Add confidence ellipse around the data from Port Gamble Bay
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "SK", col = "#017A7488") #Add confidence ellipse around the data from Skokomish River Delta
ordiellipse(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "WB", col = "#EB8B0C88") #Add confidence ellipse around the data from Willapa Bay
```


```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "Axis 1", ylab = "Axis 2", cex = 0.5) #Save NMDS as a new object
points(fig.nmds, "sites", col = plotCustomization$Color, pch = plotCustomization$Shape2, bg = plotCustomization$Color)
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "CI", col = "#00A9BD") #Add confidence ellipse around the data from Case Inlet
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "FB", col = "#38001C") #Add confidence ellipse around the data  from Fidalgo Bay
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "PG", col = "#440D82") #Add confidence ellipse around the data from Port Gamble Bay
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "SK", col = "#017A74") #Add confidence ellipse around the data from Skokomish River Delta
ordihull(meanData.log.gower.NMDS, plotCustomization$Site, show.groups = "WB", col = "#EB8B0C") #Add confidence ellipse around the data from Willapa Bay
legend("topright", pch = c(rep(x = 1, times = 6), 16, 0, 1, 5, 6), legend=c('Case Inlet', "Fidalgo Bay", "Port Gamble Bay", "Skokomish River Delta", "Willapa Bay", "Bare", "Eelgrass", "pH", "Dissolved Oxygen", "Salinity", "Temperature"), col=c('#00A9BD', '#38001C', '#440D82', '#017A74', '#EB8B0C', "black", "black", "black", "black", "black", "black"), cex = 0.5, bty = "n")
```

#### Loadings only

```{r}
sigLoadings <- envfit(meanData.log.gower.NMDS$points, meanData.log.trans[,c(1:3, 5, 18:19)], perm = 1000, na.rm = TRUE) #Only calculate loadings simper identified as driving differences between FB-WB and SK-WB. See ANOSIM section for simper results.
sigLoadings #View loadings
```


```{r}
fig.nmds <- ordiplot(meanData.log.gower.NMDS, choices=c(1,2), type = "none", display = "sites", xlab = "", ylab = "", cex = 0.5, xaxt = "n", yaxt = "n") #Save NMDS as a new object
plot(sigLoadings, col = 'grey20') #Plot loadings that simper determined were significant
axis(side = 1, labels = TRUE, col = "grey80", cex.axis = 0.75)
mtext(side = 1, text = "NMDS1", line = 2)
axis(side = 2, labels = TRUE, col = "grey80", cex.axis = 0.75)
mtext(side = 2, text = "NMDS2", line = 2)
box(col = "grey80")
```

## Analysis of Similarities (ANOSIM)

### Global ANOSIM

#### Site

```{r}
siteANOSIM <- anosim(meanData.log.trans, grouping = plotCustomization$Site) #One-way ANOSIM by Site
siteANOSIM$statistic #R = -0.07799145
siteANOSIM$signif #p = 0.968
plot(siteANOSIM) #Obtain boxplots and permutation test histogram
```

#### Habitat

```{r}
habitatANOSIM <- anosim(meanData.log.trans, grouping = plotCustomization$Habitat) #One-way ANOSIM by Habitat
habitatANOSIM$statistic #R = -0.03833658
habitatANOSIM$signif #p = 0.925
plot(habitatANOSIM) #Obtain boxplots and permutation test histogram
```

#### Site and Habitat

```{r}
siteAndHabitatANOSIM <- anosim(meanData.log.trans, grouping = plotCustomization$Site.Habitat) #Two-way ANOSIM by Site and Habitat
siteAndHabitatANOSIM$statistic #R = -0.2203411
siteAndHabitatANOSIM$signif #p = 1
plot(siteAndHabitatANOSIM) #Obtain boxplots and permutation test histogram
```

#### Environmental Variable

```{r}
plotCustomization$EnvironmentalVariable <- rep(0, times = length(plotCustomization$Site))
plotCustomization[grep(rownames(plotCustomization), pattern = "pH"), "EnvironmentalVariable"] <- "pH"
plotCustomization[grep(rownames(plotCustomization), pattern = "DO"), "EnvironmentalVariable"] <- "DO"
plotCustomization[grep(rownames(plotCustomization), pattern = "sal"), "EnvironmentalVariable"] <- "S"
plotCustomization[grep(rownames(plotCustomization), pattern = "temp"), "EnvironmentalVariable"] <- "T"
tail(plotCustomization)
```

```{r}
environmentalVariableANOSIM <- anosim(meanData.log.trans, grouping = plotCustomization$EnvironmentalVariable) #One-way ANOSIM by Environmental Variable
environmentalVariableANOSIM$statistic #R = 0.8339003
environmentalVariableANOSIM$signif #p = 0.001
plot(environmentalVariableANOSIM) #Obtain boxplots and permutation test histogram
```

#### Environmental Variable and Site

```{r}
plotCustomization$EnvVariable.Site <- paste(plotCustomization$EnvironmentalVariable, "", plotCustomization$Site) #Combine Environmental Variable and Site
environmentalVariableSiteANOSIM <- anosim(meanData.log.trans, grouping = plotCustomization$EnvVariable.Site) #Two-way ANOSIM by Environmental Variable and Site
environmentalVariableSiteANOSIM$statistic #R = 0.7675
environmentalVariableSiteANOSIM$signif #p = 0.001
plot(environmentalVariableSiteANOSIM) #Obtain boxplots and permutation test histogram
```

#### Environmental Variable and Habitat

```{r}
plotCustomization$EnvVariable.Habitat <- paste(plotCustomization$EnvironmentalVariable, "", plotCustomization$Habitat) #Combine Environmental Variable and Habitat
environmentalVariableHabitatANOSIM <- anosim(meanData.log.trans, grouping = plotCustomization$EnvVariable.Habitat) #Two-way ANOSIM by Environmental Variable and Habitat
environmentalVariableHabitatANOSIM$statistic #R = 0.698505
environmentalVariableHabitatANOSIM$signif #p = 0.001
plot(environmentalVariableHabitatANOSIM) #Obtain boxplots and permutation test histogram
```

### Pairwise ANOSIM

```{r}
pairwiseData <- merge(x = meanData.log.trans, y = plotCustomization, by = "row.names") #Merge dataframes
rownames(pairwiseData) <- pairwiseData$Row.names #Make the column the rownames
pairwiseData <- pairwiseData[, -1] #Remove Sample.Number column
head(pairwiseData) #Confirm merge
```

#### Environmental Variable and Habitat

##### pH: Bare vs. Eelgrass

```{r}
datapH <- pairwiseData[grep(rownames(pairwiseData), pattern = "pH"),] #Subset data
head(datapH) #Confirm subset
```

```{r}
pHHabitatANOSIM <- anosim(datapH[,1:31], grouping = datapH$Habitat) #Conduct ANOSIM
pHHabitatANOSIM$statistic #R = 0.2125
pHHabitatANOSIM$signif #p = 0.135
plot(pHHabitatANOSIM) #Obtain boxplots and permutation test histogram
```

##### Dissolved oxygen: Bare vs. Eelgrass

```{r}
dataDO <- pairwiseData[grep(rownames(pairwiseData), pattern = "DO"),] #Subset data
head(dataDO) #Confirm subset
```

```{r}
DOHabitatANOSIM <- anosim(dataDO[,1:31], grouping = dataDO$Habitat) #Conduct ANOSIM
DOHabitatANOSIM$statistic #R = -0.14375
DOHabitatANOSIM$signif #p = 0.922
plot(DOHabitatANOSIM) #Obtain boxplots and permutation test histogram
```

##### Salinity: Bare vs. Eelgrass

```{r}
dataSalinity <- pairwiseData[grep(rownames(pairwiseData), pattern = "sal"),] #Subset data
head(dataSalinity) #Confirm subset
```

```{r}
salinityHabitatANOSIM <- anosim(dataSalinity[,1:31], grouping = dataSalinity$Habitat) #Conduct ANOSIM
salinityHabitatANOSIM$statistic #R = -0.1666667
salinityHabitatANOSIM$signif #p = 0.971
plot(salinityHabitatANOSIM) #Obtain boxplots and permutation test histogram
```

##### Temperature: Bare vs. Eelgrass

```{r}
dataTemperature <- pairwiseData[grep(rownames(pairwiseData), pattern = "temp"),] #Subset data
head(dataTemperature) #Confirm subset
```

```{r}
temperatureHabitatANOSIM <- anosim(dataTemperature[,1:31], grouping = dataTemperature$Habitat) #Conduct ANOSIM
temperatureHabitatANOSIM$statistic #R = -0.192
temperatureHabitatANOSIM$signif #p = 0.905
plot(temperatureHabitatANOSIM) #Obtain boxplots and permutation test histogram
```

#### Environmental Variable and Site

##### pH: All sites

```{r}
pHSiteANOSIM <- anosim(datapH[,1:31], grouping = datapH$Site) #Conduct ANOSIM
pHSiteANOSIM$statistic #R = 0.328125
pHSiteANOSIM$signif #p = 0.151
plot(pHSiteANOSIM) #Obtain boxplots and permutation test histogram
```

##### Dissolved oxygen: All sites

```{r}
DOSiteANOSIM <- anosim(dataDO[,1:31], grouping = dataDO$Site) #Conduct ANOSIM
DOSiteANOSIM$statistic #R = 0.40625
DOSiteANOSIM$signif #p = 0.053
plot(DOSiteANOSIM) #Obtain boxplots and permutation test histogram
```

```{r}
simperDOSite <- simper(dataDO[,1:31], group = dataDO$Site) #Identify what could be driving significant differences
summary(simperDOSite)
```

##### Salinity: All sites

```{r}
salinitySiteANOSIM <- anosim(dataSalinity[,1:31], grouping = dataSalinity$Site) #Conduct ANOSIM
salinitySiteANOSIM$statistic #R = -0.01333333
salinitySiteANOSIM$signif #p = 0.482
plot(salinitySiteANOSIM) #Obtain boxplots and permutation test histogram
```

##### Temperature: All sites

```{r}
temperatureSiteANOSIM <- anosim(dataTemperature[,1:31], grouping = dataTemperature$Site) #Conduct ANOSIM
temperatureSiteANOSIM$statistic #R = 1
temperatureSiteANOSIM$signif #p = 0.002
plot(temperatureSiteANOSIM) #Obtain boxplots and permutation test histogram
```

```{r}
simperTemperatureSite <- simper(dataTemperature[,1:31], group = dataTemperature$Site) #Identify what could be driving significant differences
summary(simperTemperatureSite)
```