Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
963 lines (747 sloc) 48.2 KB
---
title: "Environmental Data NMDS"
author: "Yaamini Venkataraman"
date: "11/18/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 values for environmental data
Mean and variance values will be autocorrelated with maximum and minimum values, so I'll only calculate mean and variance.
```{r}
outplantDates <- as.Date(unique(temperatureData$Date)) #Save outplant dates as a new vector
outplantDates #Confirm vector creation
```
## pH
```{r}
head(pHTideData)
```
### Mean
```{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.
```
### Variance
```{r}
pHVariance <- data.frame("Date" = outplantDates,
"WBE-pH-var" = rep(0, length(outplantDates)),
"WBB-pH-var" = rep(0, length(outplantDates)),
"SKE-pH-var" = rep(0, length(outplantDates)),
"SKB-pH-var" = rep(0, length(outplantDates)),
"PGB-pH-var" = rep(0, length(outplantDates)),
"CIE-pH-var" = rep(0, length(outplantDates)),
"CIB-pH-var" = rep(0, length(outplantDates)),
"FBE-pH-var" = rep(0, length(outplantDates)),
"FBB-pH-var" = rep(0, length(outplantDates))) #Create an empty dataframe
head(pHVariance)
```
```{r}
for (j in 4:12) { #For each column with pH data
for(i in 1:nDates) { #For each date
pHVariance[i, j-2] <- var(pHTideData[pHTideData$Date == outplantDates[i], j], na.rm = TRUE) #Find and save the variance 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 "N/A" in the dataframe.
head(pHVariance) #Confirm dataframe filling.
```
## Dissolved oxygen
```{r}
head(DOTideData)
```
### Mean
```{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
```
### Variance
```{r}
DOVariance <- data.frame("Date" = outplantDates,
"WBE-DO-var" = rep(0, length(outplantDates)),
"WBB-DO-var" = rep(0, length(outplantDates)),
"SKE-DO-var" = rep(0, length(outplantDates)),
"SKB-DO-var" = rep(0, length(outplantDates)),
"PGE-DO-var" = rep(0, length(outplantDates)),
"PGB-DO-var" = rep(0, length(outplantDates)),
"CIE-DO-var" = rep(0, length(outplantDates)),
"CIB-DO-var" = rep(0, length(outplantDates)),
"FBE-DO-var" = rep(0, length(outplantDates)),
"FBB-DO-var" = rep(0, length(outplantDates))) #Create an empty dataframe
head(DOVariance)
```
```{r}
for (j in 4:13) { #For each column with DO data
for(i in 1:nDates) { #For each date
DOVariance[i, j-2] <- var(DOTideData[DOTideData$Date == outplantDates[i], j], na.rm = TRUE) #Find and save the variance of the designated column at each date.
}
} #The loop will cycle through each date first in one column, then move to the next column.
head(DOVariance) #Confirm dataframe filling.
```
## Salinity
```{r}
head(salinityTideData)
```
### Mean
```{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
```
### Variance
```{r}
salinityVariance <- data.frame("Date" = outplantDates,
"CIB-sal-var" = rep(0, length(outplantDates)),
"CIE-sal-var" = rep(0, length(outplantDates)),
"FBB-sal-var" = rep(0, length(outplantDates)),
"FBE-sal-var" = rep(0, length(outplantDates)),
"PGE-sal-var" = rep(0, length(outplantDates)),
"SKB-sal-var" = rep(0, length(outplantDates)),
"SKE-sal-var" = rep(0, length(outplantDates)),
"WBB-sal-var" = rep(0, length(outplantDates))) #Create an empty dataframe
head(salinityVariance)
```
```{r}
for (j in 4:11) { #For each column with salinity data
for(i in 1:nDates) { #For each date
salinityVariance[i, j-2] <- var(salinityTideData[salinityTideData$Date == outplantDates[i], j], na.rm = TRUE) #Find and save the variance 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 N/A in the dataframe.
head(salinityVariance) #Confirm dataframe filling.
```
## Temperature
```{r}
head(temperatureData)
```
```{r}
outplantDateCharacters <- as.character(outplantDates) #Convert outplant dates to characters
temperatureData$DateCharacters <- as.character(temperatureData$Date) #Convert dates to characters
```
### Mean
```{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)
```
### Variance
```{r}
temperatureVariance <- data.frame("Date" = outplantDateCharacters,
"WBE-temp-var" = rep(0, length(outplantDates)),
"WBB-temp-var" = rep(0, length(outplantDates)),
"SKE-temp-var" = rep(0, length(outplantDates)),
"SKB-temp-var" = rep(0, length(outplantDates)),
"PGE-temp-var" = rep(0, length(outplantDates)),
"PGB-temp-var" = rep(0, length(outplantDates)),
"CIE-temp-var" = rep(0, length(outplantDates)),
"CIB-temp-var" = rep(0, length(outplantDates)),
"FBE-temp-var" = rep(0, length(outplantDates)),
"FBB-temp-var" = rep(0, length(outplantDates))) #Create an empty dataframe, but be sure to use outplantDateCharacters
head(temperatureVariance)
```
```{r}
for (j in 3:12) { #For each column with temperature data
for(i in 1:nDates) { #For each date
temperatureVariance[i, j-1] <- var(temperatureData[temperatureData$DateCharacters == outplantDateCharacters[i], j], na.rm = TRUE) #Find and save the variance 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.
head(temperatureVariance) #Confirm dataframe filling. All "Inf" values are now replaced
```
```{r}
temperatureVariance$Date <- outplantDates #Replace temperatureVariance$Date with outplantDates so they are dates
head(temperatureVariance)
```
## 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
```
```{r}
dailyVarianceData <- cbind(pHVariance, DOVariance, salinityVariance, temperatureVariance) #Use cbind to merge all of the dataframes together
colnames(dailyVarianceData) #View column names
```
```{r}
dailyVarianceData <- dailyVarianceData[,-c(11, 22, 31)] #Remove all redundant "Date" columns
colnames(dailyVarianceData) #Confirm column removal
```
# Environmental data NMDS
## Transform data
```{r}
dailyEnvironmentalDataCorrected <- dailyEnvironmentalData #Duplicate dataframe
rownames(dailyEnvironmentalDataCorrected) <- dailyEnvironmentalDataCorrected$Date #Save dates as rownames
dailyEnvironmentalDataCorrected <- dailyEnvironmentalDataCorrected[, -1] #Remove date column
head(dailyEnvironmentalDataCorrected) #Confirm changes
```
```{r}
envData.log <- log(dailyEnvironmentalDataCorrected + 1) #Log transform data before analysis
envData.log <- envData.log[-c(32:34),] #Remove last three rows since the outplant ended before then
tail(envData.log)
```
```{r}
envData.log.trans <- data.frame(t(envData.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(envData.log.trans) <- rownames(envData.log) #Rename columns so there is no "X" in front of dates
head(envData.log.trans) #View transformation
```
```{r}
head(envData.log.trans) #View changes
```
## Calculate dissimilarity matrix
```{r}
envDataGowerDiss <- daisy(envData.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.
envDataGowerDiss #Confirm matrix calculation
```
## Conduct the NMDS
```{r}
nmds.scree(envDataGowerDiss, 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}
envData.log.gower.NMDS <- metaMDS(envDataGowerDiss, distance = 'euclidean', k = 2, trymax = 10000, autotransform = FALSE) #Make MDS dissimilarity matrix on log transformed data using Gower's.
envData.log.gower.NMDS$stress #Stress = 0.0344334
```
```{r}
nmds.monte(envDataGowerDiss, 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, 0.0344334, is less than the expected stress value. P-value =
```
```{r}
stressplot(envData.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 0.990 (redundant with observed stress value and p-value from the randomization test)
```
```{r}
vec.envData.log.gower <- envfit(envData.log.gower.NMDS$points, envData.log.trans, perm = 1000, na.rm = TRUE) #Calculate loadings by correlating NMDS scores with original variables
vec.envData.log.gower #Look at loadings
```
```{r}
ordiplot(envData.log.gower.NMDS, choices = c(1,2), type = "text", display = "sites", xlab = "Axis 1", ylab = "Axis 2") #Plot basic NMDS
plot(vec.envData.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(envData.log.trans$`2016-06-19`)),
"Habitat" = rep(0, times = length(envData.log.trans$`2016-06-19`)),
"Site-Habitat" = rep(0, times = length(envData.log.trans$`2016-06-19`)),
"Shape" = rep(0, times = length(envData.log.trans$`2016-06-19`)),
"Color" = rep(0, times = length(envData.log.trans$`2016-06-19`)),
"Shape2" = rep(0, times = length(envData.log.trans$`2016-06-19`)),
"Parameter" = rep(0, times = length(envData.log.trans$`2016-06-19`))) #Create a new dataframe to store plot customization information
rownames(plotCustomization) <- rownames(envData.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
```
### Parameter
```{r}
plotCustomization[grep(rownames(plotCustomization), pattern = "mean"), "Parameter"] <- "M" #Add "M" to the Parameter column for each instance of "mean" in the row names
plotCustomization[grep(rownames(plotCustomization), pattern = "var"), "Parameter"] <- "V" #Add "V" to the Parameter column for each instance of "mean" in the row names
tail(plotCustomization)
```
### Remake NMDS
#### Site only
```{r}
fig.nmds <- ordiplot(envData.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(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "CI", col = "#00A9BD") #Add confidence ellipse around the data from Case Inlet
ordiellipse(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "FB", col = "#38001C") #Add confidence ellipse around the data from Fidalgo Bay
ordiellipse(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "PG", col = "#440D82") #Add confidence ellipse around the data from Port Gamble Bay
ordiellipse(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "SK", col = "#017A74") #Add confidence ellipse around the data from Skokomish River Delta
ordiellipse(envData.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(envData.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(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "CI", col = "#00A9BD") #Add confidence ellipse around the data from Case Inlet
ordihull(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "FB", col = "#38001C") #Add confidence ellipse around the data from Fidalgo Bay
ordihull(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "PG", col = "#440D82") #Add confidence ellipse around the data from Port Gamble Bay
ordihull(envData.log.gower.NMDS, plotCustomization$Site, show.groups = "SK", col = "#017A74") #Add confidence ellipse around the data from Skokomish River Delta
ordihull(envData.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(envData.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(envData.log.gower.NMDS, plotCustomization$Habitat, show.groups = "Eelgrass", col = "grey80", lty = 1) #Add confidence ellipse around the data from eelgrass habitats
ordiellipse(envData.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(envData.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(envData.log.gower.NMDS, plotCustomization$Habitat, show.groups = "Eelgrass", col = "grey80", lty = 1) #Add confidence ellipse around the data from eelgrass habitats
ordihull(envData.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(envData.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)
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")
legend("topleft", lty = c(1, 2), legend = c("Mean", "Variance"), col = c("grey80", "grey80"), cex = 0.5, bty = "n")
ordiellipse(envData.log.gower.NMDS, plotCustomization$Parameter, show.groups = "M", col = "grey80", lty = 1) #Add confidence ellipse around the data from means
ordiellipse(envData.log.gower.NMDS, plotCustomization$Parameter, show.groups = "V", col = "grey80", lty = 2) #Add confidence ellipse around the data from means
```
```{r}
fig.nmds <- ordiplot(envData.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)
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")
legend("topleft", lty = c(1, 2), legend = c("Mean", "Variance"), col = c("grey80", "grey80"), cex = 0.5, bty = "n")
ordihull(envData.log.gower.NMDS, plotCustomization$Parameter, show.groups = "M", col = "grey80", lty = 1) #Add polygon around the data from means
ordihull(envData.log.gower.NMDS, plotCustomization$Parameter, show.groups = "V", col = "grey80", lty = 2) #Add polygon around the data from means
```
## Analysis of Similarities (ANOSIM)
### Global ANOSIM
#### Site
```{r}
siteANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$Site) #One-way ANOSIM by Site
siteANOSIM$statistic #R = -0.03428841
siteANOSIM$signif #p = 0.943
plot(siteANOSIM) #Obtain boxplots and permutation test histogram
```
#### Habitat
```{r}
habitatANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$Habitat) #One-way ANOSIM by Habitat
habitatANOSIM$statistic #R = -0.02037461
habitatANOSIM$signif #p = 0.982
plot(habitatANOSIM) #Obtain boxplots and permutation test histogram
```
#### Site and Habitat
```{r}
siteAndHabitatANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$Site.Habitat) #Two-way ANOSIM by Site and Habitat
siteAndHabitatANOSIM$statistic #R = -0.09991229
siteAndHabitatANOSIM$signif #p = 1
plot(siteAndHabitatANOSIM) #Obtain boxplots and permutation test histogram
```
#### Parameter
```{r}
parameterANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$Parameter) #One-way ANOSIM by Parameter (Means vs. Variances)
parameterANOSIM$statistic #R = 0.7378174
parameterANOSIM$signif #p = 0.001
plot(parameterANOSIM) #Obtain boxplots and permutation test histogram
```
#### Parameter and Site
```{r}
plotCustomization$Parameter.Site <- paste(plotCustomization$Parameter, "", plotCustomization$Shape) #Combine Site and Parameter
parameterSiteANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$Parameter.Site) #Two-way ANOSIM by Parameter (Means vs. Variances) and Site
parameterSiteANOSIM$statistic #R = 0.4778027
parameterSiteANOSIM$signif #p = 0.001
plot(parameterSiteANOSIM) #Obtain boxplots and permutation test histogram
```
#### Parameter and Habitat
```{r}
plotCustomization$Parameter.Habitat <- paste(plotCustomization$Parameter, "", plotCustomization$Habitat) #Combine Habitat and Parameter
parameterHabitatANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$Parameter.Habitat) #Two-way ANOSIM by Parameter (Means vs. Variances) and Habitat
parameterHabitatANOSIM$statistic #R = 0.4778027
parameterHabitatANOSIM$signif #p = 0.001
plot(parameterHabitatANOSIM) #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(envDataGowerDiss, grouping = plotCustomization$EnvironmentalVariable) #One-way ANOSIM by Environmental Variable
environmentalVariableANOSIM$statistic #R = 0.1796322
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(envDataGowerDiss, grouping = plotCustomization$EnvVariable.Site) #Two-way ANOSIM by Environmental Variable and Site
environmentalVariableSiteANOSIM$statistic #R = 0.02098092
environmentalVariableSiteANOSIM$signif #p = 0.315
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(envDataGowerDiss, grouping = plotCustomization$EnvVariable.Habitat) #Two-way ANOSIM by Environmental Variable and Habitat
environmentalVariableHabitatANOSIM$statistic #R = 0.1089352
environmentalVariableHabitatANOSIM$signif #p = 0.007
plot(environmentalVariableHabitatANOSIM) #Obtain boxplots and permutation test histogram
```
#### Environmental Variable and Parameter
```{r}
plotCustomization$EnvVariable.Parameter <- paste(plotCustomization$EnvironmentalVariable, "", plotCustomization$Parameter) #Combine Environmental Variable and Parameter
environmentalVariableParameterANOSIM <- anosim(envDataGowerDiss, grouping = plotCustomization$EnvVariable.Parameter) #Two-way ANOSIM by Environmental Variable and Habitat
environmentalVariableParameterANOSIM$statistic #R = 0.7802082
environmentalVariableParameterANOSIM$signif #p = 0.001
plot(environmentalVariableParameterANOSIM) #Obtain boxplots and permutation test histogram
```
### Pairwise ANOSIM