diff --git a/phyto-final/CSVs/EMP_Stations_All.csv b/phyto-final/CSVs/EMP_Stations_All.csv deleted file mode 100644 index 8616247..0000000 --- a/phyto-final/CSVs/EMP_Stations_All.csv +++ /dev/null @@ -1,94 +0,0 @@ -Station ID,Location,Latitude,Longitude,StationType -ANH,San Joaquin River at Antioch,38.017854,-121.802939,Continuous -C10A,San Joaquin River near Vernalis @ SJR Club,37.67934,-121.2647,Continuous -C10A,San Joaquin River near Vernalis @ SJR Club,37.67934,-121.2647,Discrete -C10A,San Joaquin River near Vernalis @ SJR Club,37.67934,-121.2647,Phytoplankton -C3A,Sacramento River @ Hood,38.36771,-121.5205,Continuous -C3A,Sacramento River @ Hood,38.36771,-121.5205,Discrete -C3A,Sacramento River @ Hood,38.36771,-121.5205,Phytoplankton -C9,West Canal @ Clifton Court Intake,37.83095,-121.554,Discrete -C9,West Canal @ Clifton Court Intake,37.83095,-121.554,Phytoplankton -D10,Sacramento River @ Chipps Island,38.04631,-121.9183,Discrete -D10,Sacramento River @ Chipps Island,38.04631,-121.9183,Phytoplankton -D10,Sacramento River @ Chipps Island,38.04631,-121.9183,Zooplankton -D12,San Joaquin River @ Antioch Ship Channel,38.02161,-121.8063,Benthic -D12,San Joaquin River @ Antioch Ship Channel,38.02161,-121.8063,Discrete -D12,San Joaquin River @ Antioch Ship Channel,38.02161,-121.8063,Phytoplankton -D12,San Joaquin River @ Antioch Ship Channel,38.02161,-121.8063,Zooplankton -D16,San Joaquin River @ Twitchell Island,38.0969,-121.6691,Discrete -D16,San Joaquin River @ Twitchell Island,38.0969,-121.6691,Phytoplankton -D16,San Joaquin River @ Twitchell Island,38.0969,-121.6691,Zooplankton -D19,Frank's Tract near Russo's Landing,38.04376,-121.6148,Discrete -D19,Frank's Tract near Russo's Landing,38.04376,-121.6148,Phytoplankton -D22,Sacramento River @ Emmaton,38.08453,-121.7391,Benthic -D22,Sacramento River @ Emmaton,38.08453,-121.7391,Discrete -D22,Sacramento River @ Emmaton,38.08453,-121.7391,Phytoplankton -D22,Sacramento River @ Emmaton,38.08453,-121.7391,Zooplankton -D24,Sacramento River below Rio Vista Bridge,38.15778,-121.6847,Continuous -D26,San Joaquin River @ Potato Slough,38.07664,-121.5669,Benthic -D26,San Joaquin River @ Potato Slough,38.07664,-121.5669,Discrete -D26,San Joaquin River @ Potato Slough,38.07664,-121.5669,Phytoplankton -D26,San Joaquin River @ Potato Slough,38.07664,-121.5669,Zooplankton -D28A,Old River @ Rancho Del Rio,37.97048,-121.573,Benthic -D28A,Old River @ Rancho Del Rio,37.97048,-121.573,Discrete -D28A,Old River @ Rancho Del Rio,37.97048,-121.573,Phytoplankton -D28A,Old River @ Rancho Del Rio,37.97048,-121.573,Zooplankton -D4,Sacramento River above Point Sacramento,38.06248,-121.8205,Benthic -D4,Sacramento River above Point Sacramento,38.06248,-121.8205,Discrete -D4,Sacramento River above Point Sacramento,38.06248,-121.8205,Phytoplankton -D4,Sacramento River above Point Sacramento,38.06248,-121.8205,Zooplankton -D41 ,San Pablo Bay near Pinole Point,38.03022,-122.3729,Benthic -D41 ,San Pablo Bay near Pinole Point,38.03022,-122.3729,Discrete -D41 ,San Pablo Bay near Pinole Point,38.03022,-122.3729,Phytoplankton -D41 ,San Pablo Bay near Pinole Point,38.03022,-122.3729,Zooplankton -D41A,San Pablo Bay near Mouth of Petaluma River,38.08472,-122.3907,Benthic -D41A,San Pablo Bay near Mouth of Petaluma River,38.08472,-122.3907,Discrete -D41A,San Pablo Bay near Mouth of Petaluma River,38.08472,-122.3907,Phytoplankton -D41A,San Pablo Bay near Mouth of Petaluma River,38.08472,-122.3907,Zooplankton -D6,Suisun Bay @ Bulls Head nr. Martinez,38.04436,-122.1177,Benthic -D6,Suisun Bay @ Bulls Head nr. Martinez,38.04436,-122.1177,Discrete -D6,Suisun Bay @ Bulls Head nr. Martinez,38.04436,-122.1177,Phytoplankton -D6,Suisun Bay @ Bulls Head nr. Martinez,38.04436,-122.1177,Zooplankton -D7,Grizzly Bay @ Dolphin nr. Suisun Slough,38.11714,-122.0397,Discrete -D7,Grizzly Bay @ Dolphin nr. Suisun Slough,38.11714,-122.0397,Phytoplankton -D7,Grizzly Bay @ Dolphin nr. Suisun Slough,38.11714,-122.0397,Zooplankton -D8,Suisun Bay off Middle Point nr. Nichols,38.05992,-121.99,Discrete -D8,Suisun Bay off Middle Point nr. Nichols,38.05992,-121.99,Phytoplankton -D8,Suisun Bay off Middle Point nr. Nichols,38.05992,-121.99,Zooplankton -FRK,Frank's Tract Mid Tract,38.046499,-121.598063,Continuous -GZL,Grizzly Bay,38.124281,-122.037962,Continuous -HON,Honker Bay,38.072229,-121.93718,Continuous -MAL,Sacramento River at Mallard Island,38.042799,-121.92013,Continuous -MD10A,Disappointment Slough @ Bishop Cut,38.04226,-121.4199,Discrete -MD10A,Disappointment Slough @ Bishop Cut,38.04226,-121.4199,Phytoplankton -MD10A,Disappointment Slough @ Bishop Cut,38.04226,-121.4199,Zooplankton -MRZ,Suisun Bay - Martinez,38.02762,-122.14052,Continuous -MSD,San Joaquin River at Mossdale Bridge,37.786144,-121.306474,Continuous -NZ002,Carquinez Strait near Glencove Harbor,38.06529,-122.2152,Discrete -NZ002,Carquinez Strait near Glencove Harbor,38.06529,-122.2152,Phytoplankton -NZ002,Carquinez Strait near Glencove Harbor,38.06529,-122.2152,Zooplankton -NZ004,Carquinez Strait near Ozol Pier,38.03576,-122.1616,Discrete -NZ004,Carquinez Strait near Ozol Pier,38.03576,-122.1616,Phytoplankton -NZ004,Carquinez Strait near Ozol Pier,38.03576,-122.1616,Zooplankton -NZ032,"Montezuma Slough, 2nd bend from mouth",38.16991,-122.0211,Discrete -NZ032,"Montezuma Slough, 2nd bend from mouth",38.16991,-122.0211,Phytoplankton -NZ032,"Montezuma Slough, 2nd bend from mouth",38.16991,-122.0211,Zooplankton -NZ068,Sacramento River below Rio Vista Bridge,38.14272,-121.6895,Benthic -NZ068,Sacramento River below Rio Vista Bridge,38.14272,-121.6895,Discrete -NZ068,Sacramento River below Rio Vista Bridge,38.14272,-121.6895,Phytoplankton -NZ068,Sacramento River below Rio Vista Bridge,38.14272,-121.6895,Zooplankton -NZ325,San Pablo Bay near Rock Wall and Light 15,38.05798,-122.2919,Discrete -NZ325,San Pablo Bay near Rock Wall and Light 15,38.05798,-122.2919,Phytoplankton -NZ325,San Pablo Bay near Rock Wall and Light 15,38.05798,-122.2919,Zooplankton -NZS42,Suisun Slough @ Volanti Slough,38.18045,-122.0476,Benthic -NZS42,Suisun Slough @ Volanti Slough,38.18045,-122.0476,Discrete -NZS42,Suisun Slough @ Volanti Slough,38.18045,-122.0476,Phytoplankton -NZS42,Suisun Slough @ Volanti Slough,38.18045,-122.0476,Zooplankton -P8,San Joaquin River @ Buckley Cove,37.97817,-121.3823,Discrete -P8,San Joaquin River @ Buckley Cove,37.97817,-121.3823,Phytoplankton -P8,San Joaquin River @ Buckley Cove,37.97817,-121.3823,Zooplankton -PPT,SJR - Prisoner's Point,38.056296,-121.549973,Continuous -RRI,SJR - Prisoner-Rough & Ready,37.962726,-121.365603,Continuous -RYC,Suisun Bay - Cutoff Near Ryer,38.083961,-121.995641,Continuous -SSI,Sacramento R Nr Sherman Island,38.074153,-121.761764,Continuous -TWI,San Joaquin River at Twitchell Island,38.097455,-121.668718,Continuous diff --git a/phyto-final/CSVs/station_regions.csv b/phyto-final/CSVs/station_regions.csv new file mode 100644 index 0000000..31d0b81 --- /dev/null +++ b/phyto-final/CSVs/station_regions.csv @@ -0,0 +1,40 @@ +StationCode,Region +C3A,Northern Interior Delta +D24,Northern Interior Delta +NZ068,Northern Interior Delta +C9,Southern Interior Delta +C10A,Southern Interior Delta +MD10A,Southern Interior Delta +P8,Southern Interior Delta +D16,Central Delta +D19,Central Delta +D26,Central Delta +D28A,Central Delta +D4,Confluence +D10,Confluence +D12,Confluence +D22,Confluence +D7,Grizzly and Suisun Bay +D8,Grizzly and Suisun Bay +NZ032,Grizzly and Suisun Bay +NZS42,Grizzly and Suisun Bay +D6,San Pablo Bay +D41,San Pablo Bay +D41A,San Pablo Bay +NZ002,San Pablo Bay +NZ004,San Pablo Bay +NZ325,San Pablo Bay +EZ2,Entrapment Zone +EZ2-SJR,Entrapment Zone +EZ6,Entrapment Zone +EZ6-SJR,Entrapment Zone +RMB,Upstream +RCS,Upstream +RD22,Upstream +I80,Upstream +LIS,Upstream +STTD,Upstream +BL5,Downstream +LIB,Downstream +RYI,Downstream +RVB,Downstream diff --git a/phyto-final/CSVs/station_regions_EMP.csv b/phyto-final/CSVs/station_regions_EMP.csv deleted file mode 100644 index 9644c5e..0000000 --- a/phyto-final/CSVs/station_regions_EMP.csv +++ /dev/null @@ -1,30 +0,0 @@ -StationCode,Region -C3A,Northern.Interior.Delta -D24,Northern.Interior.Delta -NZ068,Northern.Interior.Delta -C9,Southern.Interior.Delta -C10A,Southern.Interior.Delta -MD10A,Southern.Interior.Delta -P8,Southern.Interior.Delta -D16,Central.Delta -D19,Central.Delta -D26,Central.Delta -D28A,Central.Delta -D4,Confluence -D10,Confluence -D12,Confluence -D22,Confluence -D7,Grizzly.and.Suisun.Bay -D8,Grizzly.and.Suisun.Bay -NZ032,Grizzly.and.Suisun.Bay -NZS42,Grizzly.and.Suisun.Bay -D6,San.Pablo.Bay -D41,San.Pablo.Bay -D41A,San.Pablo.Bay -NZ002,San.Pablo.Bay -NZ004,San.Pablo.Bay -NZ325,San.Pablo.Bay -EZ2,Entrapment.Zone -EZ2-SJR,Entrapment.Zone -EZ6,Entrapment.Zone -EZ6-SJR,Entrapment.Zone diff --git a/phyto-final/RData/df_phyto.RData b/phyto-final/RData/df_phyto.RData index 0da9403..f2ec4c7 100644 Binary files a/phyto-final/RData/df_phyto.RData and b/phyto-final/RData/df_phyto.RData differ diff --git a/phyto-final/RData/df_phyto_FASTR_gen.RData b/phyto-final/RData/df_phyto_FASTR_gen.RData new file mode 100644 index 0000000..3cfb30d Binary files /dev/null and b/phyto-final/RData/df_phyto_FASTR_gen.RData differ diff --git a/phyto-final/RData/df_phyto_FASTR_gen_BV.RData b/phyto-final/RData/df_phyto_FASTR_gen_BV.RData new file mode 100644 index 0000000..46166f4 Binary files /dev/null and b/phyto-final/RData/df_phyto_FASTR_gen_BV.RData differ diff --git a/phyto-final/RData/df_phyto_FASTR_grp.RData b/phyto-final/RData/df_phyto_FASTR_grp.RData new file mode 100644 index 0000000..a1119eb Binary files /dev/null and b/phyto-final/RData/df_phyto_FASTR_grp.RData differ diff --git a/phyto-final/RData/df_phyto_FASTR_grp_BV.RData b/phyto-final/RData/df_phyto_FASTR_grp_BV.RData new file mode 100644 index 0000000..fb561e4 Binary files /dev/null and b/phyto-final/RData/df_phyto_FASTR_grp_BV.RData differ diff --git a/phyto-final/RData/df_phyto_FASTR_sum.RData b/phyto-final/RData/df_phyto_FASTR_sum.RData new file mode 100644 index 0000000..2346261 Binary files /dev/null and b/phyto-final/RData/df_phyto_FASTR_sum.RData differ diff --git a/phyto-final/RData/df_phyto_NMDS.RData b/phyto-final/RData/df_phyto_NMDS.RData new file mode 100644 index 0000000..f23b2ad Binary files /dev/null and b/phyto-final/RData/df_phyto_NMDS.RData differ diff --git a/phyto-final/analyses/phyto_NMDS_stress.csv b/phyto-final/analyses/phyto_NMDS_stress.csv new file mode 100644 index 0000000..48079a2 --- /dev/null +++ b/phyto-final/analyses/phyto_NMDS_stress.csv @@ -0,0 +1,7 @@ +Year,stress +2014,0.17163426035862414 +2015,0.21308563651147788 +2016,0.1944142356455662 +2017,0.18958791263284083 +2018,0.1938840696502086 +2019,0.17471661483123263 diff --git a/phyto-final/phyto-analyses-FASTR.Rmd b/phyto-final/phyto-analyses-FASTR.Rmd index ccbc399..8ad06d0 100644 --- a/phyto-final/phyto-analyses-FASTR.Rmd +++ b/phyto-final/phyto-analyses-FASTR.Rmd @@ -1,5 +1,5 @@ --- -title: "Phytoplankton Data Processing for FASTR" +title: "FASTR - Phytoplankton Analyses" author: "Ted Flynn" date: "2024-03-18" output: html_document @@ -16,581 +16,27 @@ suppressWarnings(suppressMessages(library(vegan))) suppressWarnings(suppressMessages(library(janitor))) suppressWarnings(suppressMessages(library(here))) -# Set visual theme in ggplot -theme_set(theme_bw()) - -``` - -## Import EMP Data - -Import phytoplankton data collected by the Environmental Monitoring Program. - -```{r import EMP, echo=FALSE} - -# Import EMP data files -phyto_files_EMP <- list.files(path = here("phyto-final","data","EMP","csv"), - pattern = "\\.csv", - full.names = T) - -df_phyto_EMP <- map(phyto_files_EMP, ~read_csv(.x, show_col_types = FALSE)) %>% - list_rbind() - -# Read in files with non-standard headers individually -df_Dec2021 <- read_csv(here("phyto-final","data","EMP","oddballs","December 2021.csv")) -df_Nov2021 <- read_csv(here("phyto-final","data","EMP","oddballs","November 2021.csv")) -df_Sep2013 <- read_csv(here("phyto-final","data","EMP","oddballs","September 2013.csv")) -df_Nov2013 <- read_csv(here("phyto-final","data","EMP","oddballs","November 2013.csv")) - -# Combine like oddball dfs -df_phyto2013 <- bind_rows(df_Sep2013, df_Nov2013) -df_phyto2021 <- bind_rows(df_Dec2021, df_Nov2021) - -# Remove individual dfs -rm(df_Dec2021) -rm(df_Nov2021) -rm(df_Nov2013) -rm(df_Sep2013) - -# Rename headers to match standard BSA headers Oddballs actually have the -# "correct" name of Total Cells rather than the incorrect "Number of cells per -# unit" - -df_phyto2013 <- df_phyto2013 %>% - rename("Number of cells per unit" = "Total Cells Counted") - -df_phyto2021 <- df_phyto2021 %>% - rename("Number of cells per unit" = "Total Number of Cells") %>% - rename("Unit Abundance" = "Unit Abundance (# of Natural Units)") - -# Combine oddball files with others -df_phyto_EMP <- bind_rows(df_phyto_EMP, df_phyto2013) -df_phyto_EMP <- bind_rows(df_phyto_EMP, df_phyto2021) - -# Remove unneeded dfs -rm(df_phyto2013) -rm(df_phyto2021) - -# Remove empty rows -df_phyto_EMP <- df_phyto_EMP %>% filter_all(any_vars(!is.na(.))) - -# Correct GALD, which is imported into two separate columns -# Test to see if NAs are 'either/or' and that there aren't some rows with a -# value in both GALD and GALD 1 - -sum(is.na(df_phyto_EMP$GALD)) # Total is 5880 -sum(is.na(df_phyto_EMP$`GALD 1`)) # Total is 8072 - -# Sum of NAs is 13952 which is the same as the number of rows in the df. -# This shows that there aren't any rows with two values so that we can -# Combine them without any issues. - -# Move GALD header -df_phyto_EMP <- df_phyto_EMP %>% relocate(`GALD 1`, .after = GALD) - -# Combine both GALD columns -df_phyto_EMP <- df_phyto_EMP %>% - rowwise() %>% - mutate(GALD.Tot = sum(c_across(GALD:`GALD 1`), na.rm = TRUE)) - -# Remove old GALD columns and rename GALD.Tot -df_phyto_EMP <- df_phyto_EMP %>% - select(!(GALD:`GALD 1`)) %>% - rename("GALD" = "GALD.Tot") - -# Clean up column names -df_phyto_EMP <- df_phyto_EMP %>% clean_names(case = "big_camel") - -df_phyto_EMP <- df_phyto_EMP %>% rename("GALD" = "Gald") - -# Remove blank columns -df_phyto_EMP <- df_phyto_EMP %>% select_if(~ !all(is.na(.))) - -# Remove columns that just have the method code "Phyto" as well as -# pre-calculated organisms per mL -df_phyto_EMP <- df_phyto_EMP %>% select(SampleDate:Biovolume10,GALD) - -# Add column to indicate what study it came from -df_phyto_EMP <- df_phyto_EMP %>% mutate(Study = "EMP", .after = StationCode) - -``` - -## Import phyto data collected for FASTR project - -```{r import FASTR} -# Import AEU data files (comment out when finished) - -phyto_files_AEU <- list.files(path = here("phyto-final","data","csv"), - pattern = "\\.csv", - full.names = T) - -df_phyto_FASTR <- map(phyto_files_AEU, ~read_csv(.x, show_col_types = FALSE)) %>% - list_rbind() - -# Remove empty rows -df_phyto_FASTR <- df_phyto_FASTR %>% filter_all(any_vars(!is.na(.))) - -# Remove weird row with only a single zero in biovolume -df_phyto_FASTR <- df_phyto_FASTR %>% drop_na(MethodCode) - -# Clean up column names -df_phyto_FASTR <- df_phyto_FASTR %>% clean_names(case = "big_camel") - -# Filter out samples from other projects. Read in station names with flags and -# merge together. -df_keepers <- read_csv(here("phyto-final","CSVs", "station_names_flagged.csv")) -df_phyto_FASTR <- left_join(df_phyto_FASTR, df_keepers) -rm(df_keepers) - -# Remove data unrelated to project -df_phyto_FASTR <- df_phyto_FASTR %>% filter(Flag == "Keep") - -# Remove flag column -df_phyto_FASTR <- df_phyto_FASTR %>% select(!(Flag)) - -# Confirm station IDs -unique(df_phyto_FASTR$StationCode) # Shows only 10 FASTR stations -table(df_phyto_FASTR$StationCode) - -# Remove blank columns -df_phyto_FASTR <- df_phyto_FASTR %>% select_if(~ !all(is.na(.))) - -# Remove individual measurement columns as they are only present in a small -# number of samples. -df_phyto_FASTR <- df_phyto_FASTR %>% select(!(Dimension:DepthM)) - -# Remove secondary and tertiary GALD measurements as they don't vary much and -# aren't present in EMP data -df_phyto_FASTR <- df_phyto_FASTR %>% select(!(Gald2:Gald3)) - -# Correct GALD, which is imported into two separate columns. Test to see if NAs -# are 'either/or' and that there aren't some rows with a value in both GALD and -# GALD 1 - -sum(is.na(df_phyto_FASTR$Gald)) # Total is 1791 -sum(is.na(df_phyto_FASTR$Gald1)) # Total is 4535 - -# Sum of NAs is 6326 which is the same as the number of rows in the df. This -# shows that there aren't any rows with two values so that we can Combine them -# without any issues. - -# Move Gald1 column -df_phyto_FASTR <- df_phyto_FASTR %>% relocate(Gald1, .after = Gald) - -# Combine both GALD columns -df_phyto_FASTR <- df_phyto_FASTR %>% - rowwise() %>% - mutate(GALD.Tot = sum(c_across(Gald:Gald1), na.rm = TRUE)) - -# Check if rows have two NAs or no NAs in the GALD columns -# test <- df_phyto_FASTR %>% -# mutate(GALD.Value = case_when(is.na(Gald) & is.na(Gald1) ~ "Fix", -# TRUE ~ "Okay")) - -# Remove old GALD columns and rename GALD.Tot -df_phyto_FASTR <- df_phyto_FASTR %>% - select(!(Gald:Gald1)) %>% - rename("GALD" = "GALD.Tot") - -# Remove MethodCode column that just says "Phyto" -df_phyto_FASTR <- df_phyto_FASTR %>% select(!(MethodCode)) - -# Add column to indicate what study it came from -df_phyto_FASTR <- df_phyto_FASTR %>% mutate(Study = "FASTR", .after = StationCode) -``` - -## Combine EMP and AEU data - -```{r combine data} -df_phyto <- bind_rows(df_phyto_EMP, df_phyto_FASTR) - -# Average all 10 biovolume measurements for each taxon -df_phyto <- df_phyto %>% - rowwise() %>% - mutate(BV.Avg = mean(c_across(Biovolume1:Biovolume10), na.rm = T)) %>% - select(!(Biovolume1:Biovolume10)) # Remove Individual Biovolume Columns - -# Remove unneeded columns -df_phyto <- df_phyto %>% select(!c("BsaTin","DiatomSoftBody")) -df_phyto <- df_phyto %>% select(!(ColonyFilamentIndividualGroupCode:Shape)) -df_phyto <- df_phyto %>% select(!(VolumeReceivedML:NumberOfFieldsCounted)) - -# Fix samples with missing times -# Just replace with 12:00:00 b/c aren't doing any time-based analyses -df_phyto <- df_phyto %>% replace_na(list(SampleTime = "12:00:00")) - -# Get dates in the right format. Some are 1/1/14 and others 1/1/2014. -df_phyto$SampleDate <- mdy(df_phyto$SampleDate) - -# Combine date and time column -df_phyto <- df_phyto %>% unite(DateTime, c("SampleDate","SampleTime"), sep = " ") #, remove = FALSE, na.rm = FALSE) - -df_phyto$DateTime <- as_datetime(df_phyto$DateTime, - tz = "US/Pacific", - format = c("%Y-%m-%d %H:%M:%OS")) - -# Check for missing dates -df_phyto %>% filter(is.na(DateTime)) # No missing dates - -# Correct BSA header -df_phyto <- df_phyto %>% rename("TotalCells" = "NumberOfCellsPerUnit") - -# Calculate Unit Density & Biovolume Density -df_phyto <- df_phyto %>% - mutate(Units.per.mL = UnitAbundance * Factor) %>% - mutate(BV.um3.per.mL= TotalCells * BV.Avg * Factor) %>% - mutate(Cells.per.mL = TotalCells * Factor) - -# Remove columns used for density calcs -# Remove Biovolume Columns -df_phyto <- df_phyto %>% select(!(c(Factor,UnitAbundance:TotalCells))) - -# Add column for year and month for highlighting data -df_phyto <- df_phyto %>% - mutate(Year = year(df_phyto$DateTime)) %>% - mutate(Month = month(df_phyto$DateTime, label = T)) - -# Order month in calendar order rather than (default) alphabetical -df_phyto$Month = factor(df_phyto$Month, levels = month.abb) - -# Reorder date/time columns -df_phyto <- df_phyto %>% - relocate(Year, .after = DateTime) %>% - relocate(Month, .after = DateTime) - -# Convert years to factors for plotting -df_phyto$Year <- as.factor(df_phyto$Year) - -``` - -## Clean Combined Data - -```{r data cleaning, echo = FALSE} -# Fix EMP site names -df_phyto$StationCode <- gsub("EZ6 SAC","EZ6",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ6SAC","EZ6",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ6SJR","EZ6-SJR",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ2SAC","EZ2",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ2 SAC","EZ2",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ2SJR","EZ2-SJR",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ2 SJR","EZ2-SJR",df_phyto$StationCode) -df_phyto$StationCode <- gsub("EZ6 SJR","EZ6-SJR",df_phyto$StationCode) -df_phyto$StationCode <- gsub("D16-Twitchell","D16",df_phyto$StationCode) -df_phyto$StationCode <- gsub("D16-Twitchel","D16",df_phyto$StationCode) -df_phyto$StationCode <- gsub("D16 - Twitchell","D16",df_phyto$StationCode) -df_phyto$StationCode <- gsub("D16 Twitchell","D16",df_phyto$StationCode) -df_phyto$StationCode <- gsub("NZ328","NZ325",df_phyto$StationCode) # Typo in August 2019 -df_phyto$StationCode <- gsub("C3A-HOOD","C3A",df_phyto$StationCode) -df_phyto$StationCode <- gsub("C3A Hood","C3A",df_phyto$StationCode) -df_phyto$StationCode <- gsub("C3A- Hood","C3A",df_phyto$StationCode) -df_phyto$StationCode <- gsub("C3A-Hood","C3A",df_phyto$StationCode) -df_phyto$StationCode <- gsub("NZ542","NZS42",df_phyto$StationCode) -df_phyto$StationCode <- gsub("E26","EZ6",df_phyto$StationCode) -df_phyto$StationCode <- gsub("E22","EZ2",df_phyto$StationCode) # Typo in May 2018 - -# Fix AEU site names -df_phyto$StationCode <- gsub("SHER","SHR",df_phyto$StationCode) -df_phyto$StationCode <- gsub("BL-5","BL5",df_phyto$StationCode) - -# Remove extra Microcystis tows at D19 -df_phyto <- df_phyto %>% filter(StationCode != "D19 MC Tow") - -# In Fall 2016, taxonomists began classifying the species -# Chroococcus microscopicus as Eucapsis microscopica. This is one of the most -# dominant species in this samples, so all taxa previously classified as -# C. microscopicus will be re-named E. microscopica - -df_phyto <- df_phyto %>% - mutate(Taxon = case_when(Taxon == 'Chroococcus microscopicus' ~ 'Eucapsis microscopica', - TRUE ~ Taxon)) %>% - mutate(Genus = case_when(Taxon == 'Eucapsis microscopica' ~ 'Eucapsis', - TRUE ~ Genus)) - -# The taxon Plagioselmis lacustris is inconsistently named, appearing sometimes as -# Rhodomonas lacustris. Change to Rhodomonas lacustris to avoid confusion. - -df_phyto <- df_phyto %>% - mutate(Taxon = case_when(Taxon == 'Plagioselmis lacustris' ~ 'Rhodomonas lacustris', - TRUE ~ Taxon)) %>% - mutate(Genus = case_when(Taxon == 'Rhodomonas lacustris' ~ 'Rhodomonas', - TRUE ~ Genus)) - -# Correct the genus label for a Chlorella entry -df_phyto$Genus <- gsub("cf Chlorella","Chlorella",df_phyto$Genus) - -# Add Taxonomy Data------------------------------------------------------------- - -# Read in CSV with manually-added WoRMS classification -df_taxa <- read_csv(here("phyto-final","CSVs","phyto_group_classification.csv"), - show_col_types = FALSE) - -df_phyto <- left_join(df_phyto, df_taxa) - -# Check if any groups are missing -sum(is.na(df_phyto$Group)) # none missing -rm(df_taxa) - -# Reorder Group classification column -df_phyto <- df_phyto %>% - relocate(Group, .before = Taxon) - -``` - -## Analyze FASTR Biovolume Data for Outliers - -```{r outlier analysis} - -# Subset combined dataset to FASTR-only -df_phyto_FASTR <- df_phyto %>% - filter(Study == "FASTR") - -# Identify if outliers are present -boxplot(df_phyto_FASTR$BV.um3.per.mL) - -# Identify 0.1% and 99.9% percentiles of the data -quartiles <- quantile(df_phyto_FASTR$BV.um3.per.mL, - probs=c(0.001, 0.999), - na.rm = TRUE) - -# List upper cutoff (99.9%) -cutoff <- quartiles[2] - -# Filter FASTR dataset to show all taxa above this 99.9% cutoff -df_outliers <- df_phyto_FASTR %>% - filter(BV.um3.per.mL > cutoff) - -df_outliers %>% - select(DateTime,Taxon,BV.um3.per.mL) %>% - kable(caption = "Al")## 4 of top 6 are Spirogyra. - -df_phyto_FASTR %>% filter(Genus == "Spirogyra")# Only 5 total samples w/ Spirogyra - -# Decision: Remove Spirogyra from datasets -# 4/5 are above 99.9% cutoff, 5th is also highly abundant -# These are benthic algae that clump, more likely to be "bullseye" samples -# that got a big blob or clump - -# Remove Spirogyra from FASTR and combined dataset -df_phyto_FASTR <- df_phyto_FASTR %>% - filter(Genus != "Spirogyra") - -df_phyto <- df_phyto %>% - filter(Genus != "Spirogyra") - -## Remove 2013 data from FASTR dataset (very limited) -df_phyto_FASTR <- df_phyto_FASTR %>% filter(Year != 2013) - -``` - -## Review Outliers in FASTR Data - -```{r spirogyra, echo = FALSE} - - - -``` - - -## Estimate Phytoplankton Biomass - -```{r biomass calcs, echo = FALSE} - -# Convert biovolume to biomass using relationships from Menden-Deuer and Lussard -# (2000) doi: 10.4319/lo.2000.45.3.0569 -# Only relevant to group-level data -# Units of BV.Density are um^3 per mL - -# Calculate Biomass (pg-C per cell) from Biovolume (um^3 per cell) -df_phyto_FASTR <- df_phyto_FASTR %>% - mutate(Biomass.pg.C = case_when(Group == "Diatoms" ~ 0.288 * (BV.Avg^0.811), - Group != "Diatoms" ~ 0.216 * (BV.Avg^0.939))) - -# Convert pg to ug (ug-C per L) -df_phyto_FASTR <- df_phyto_FASTR %>% - mutate(Biomass.ug.C = Biomass.pg.C / 10^6, .keep = "unused") - -# Calculate Biomass Density (pg-C per mL) -df_phyto_FASTR <- df_phyto_FASTR %>% - mutate(BM.ug.per.mL = Biomass.ug.C * Cells.per.mL) - - - -``` - -## Add FASTR-Specific Metadata - -```{r FASTR metadata} - -# Import Stations listed as Upstream and Downstream -df_region <- read_csv(here("phyto-final","CSVs","upstream_downstream_stations.csv"), - show_col_types = FALSE) - -df_region <- df_region %>% - rename("Region" = "UpDown") %>% - rename("StationCode" = "Site") - -df_region$Region <- factor(df_region$Region, - levels = c("Upstream","Downstream")) - -## Add Region to data frames -df_phyto_FASTR <- left_join(df_phyto_FASTR, df_region) - -df_phyto_FASTR <- df_phyto_FASTR %>% - relocate(Region, .after = StationCode) - -# Cat Pien's Code for Merging Flow Designations ##### -# -# FlowDesignation <- read_csv("FlowDatesDesignations.csv") -# Update with 45 days limit on either end -# Update 5/27/22 TF added >= to two of the mutate commands to avoid removing -# samples falling on the PreFlowEnd and PostFlowStart dates. - -# Upload flow dates and categories -df_flow_designation <- read_csv(here("phyto-final","CSVs","FlowDatesDesignations_45days.csv"), - show_col_types = FALSE) - -df_flow_designation$Year <- as.factor(df_flow_designation$Year) - -# Format dates as dates -df_flow_designation$PreFlowStart <- mdy(df_flow_designation$PreFlowStart) -df_flow_designation$PreFlowEnd <- mdy(df_flow_designation$PreFlowEnd) -df_flow_designation$PostFlowStart <- mdy(df_flow_designation$PostFlowStart) -df_flow_designation$PostFlowEnd <- mdy(df_flow_designation$PostFlowEnd) - -# Remove column with net flow days -df_flow_designation <- df_flow_designation %>% - select(!(NetFlowDays)) - -# Save file for use in making plots -# save(FlowDesignation, file = "RData/FlowDesignation.RData") - -# Merge data from FlowDesignation Table (Water Year Type, Flow days and type) -df_phyto_FASTR <- inner_join(df_phyto_FASTR, df_flow_designation, - by = "Year") - -# Filter only Pre-During-Post Flow Action Data. -df_phyto_FASTR <- df_phyto_FASTR %>% - mutate(ActionPhase = ifelse(DateTime > PreFlowStart & DateTime < PreFlowEnd, "Before", NA)) %>% - mutate(ActionPhase = replace(ActionPhase, DateTime >= PreFlowEnd & DateTime < PostFlowStart, "During")) %>% - mutate(ActionPhase = replace(ActionPhase, DateTime >= PostFlowStart & DateTime < PostFlowEnd, "After")) %>% - filter(!is.na(ActionPhase)) %>% - select(-c(PreFlowStart:PostFlowEnd)) - -# Order the new ActionPhase so that it plots in the order Pre < During < Post -phase.order <- c("Before","During","After") -df_phyto_FASTR$ActionPhase <- factor(as.character(df_phyto_FASTR$ActionPhase), - levels = phase.order) - -## Add in Flow Pulse Category Data -df_pulse_category <- read_csv(here("phyto-final","CSVs","FlowPulseType.csv"), - show_col_types = FALSE) - -df_pulse_category <- df_pulse_category %>% - rename("FlowPulseCategory" = "FlowPulseType") - -df_pulse_category$Year <- as.factor(df_pulse_category$Year) - -# Add Flow Pulse Category to data frames to be exported -df_phyto_FASTR <- left_join(df_phyto_FASTR, df_pulse_category, - by = "Year") - -# List of stations used for FASTR -stations <- c("RMB","RCS","RD22","I80","LIS","STTD","BL5","LIB","RYI","RVB") - -# Re-order stations from North to South -df_phyto_FASTR$StationCode <- factor(as.character(df_phyto_FASTR$StationCode), - levels = stations) - -## Relocate metadata columns -df_phyto_FASTR <- df_phyto_FASTR %>% - relocate(FlowPulseCategory, .before = ActionPhase) %>% - relocate(WYType:ActionPhase, .before = GALD) - -# Save df to use for making plots -# save(df_phyto, file = "RData/df_phyto.RData") - -``` - -## Create New Data Frames summarizing by Genus and Algal Group - -```{r summarize phyto data} -# Summarize by algal group -df_phyto_FASTR_grp <- df_phyto_FASTR %>% - group_by(Year, Month, Region, ActionPhase, DateTime, StationCode, Group) %>% - summarize(across(Units.per.mL:BM.ug.per.mL, ~sum(.x, na.rm = TRUE))) %>% - ungroup - -## Summarize by genus -df_phyto_FASTR_gen <- df_phyto_FASTR %>% - group_by(Year, Month, Region, ActionPhase, DateTime, StationCode, Genus) %>% - summarize(across(Units.per.mL:BM.ug.per.mL, ~sum(.x, na.rm = TRUE))) %>% - ungroup - -## Set group display order -group.order <- c("Diatoms", - "Cyanobacteria", - "Green Algae", - "Cryptophytes", - "Ciliates", - "Dinoflagellates", - "Golden Algae", - "Other") - -df_phyto_FASTR_grp$Group <- factor(as.character(df_phyto_FASTR_grp$Group), - levels = group.order) - -## Summarize totals -df_phyto_FASTR_sum <- df_phyto_FASTR %>% - group_by(Year, Month, DateTime, StationCode, Region, WYType, FlowPulseType, FlowPulseCategory, ActionPhase) %>% - summarize(across(Units.per.mL:BM.ug.per.mL, ~sum(.x, na.rm = TRUE))) %>% - rename("Total.Units.per.mL" = "Units.per.mL") %>% - rename("Total.Cells.per.mL" = "Cells.per.mL") %>% - rename("Total.BV.per.mL" = "BV.um3.per.mL") %>% - rename("Total.BM.per.mL" = "BM.ug.per.mL") %>% - ungroup - -# Separate out data frames by density type and add zeros -# Biovolume Density + Group -df_phyto_FASTR_grp_BV <- df_phyto_FASTR_grp %>% - select(Year:Group,BV.um3.per.mL) - -# Add zeros for taxa that weren't detected -temp <- pivot_wider(df_phyto_FASTR_grp_BV, - names_from = Group, - values_from = BV.um3.per.mL, - values_fill = 0) - -df_phyto_FASTR_grp_BV <- pivot_longer(temp, - cols = Cyanobacteria:last_col(), - names_to = "Group", - values_to = "BV.um3.per.mL") - - - +# Set location to save plots +plots <- here("phyto-final","plots") ``` -## Calculate LCEFA Abundance in FASTR Data - -```{r LCEFA calcs} -## Calculate LCEFA composition based on method in Galloway & Winder (2015) -## doi: 10.1371/journal.pone.0130053 - -df_phyto_FASTR_grp <- df_phyto_FASTR_grp %>% - mutate(LCEFA.per.mL = case_when(Group == 'Diatoms' ~ BM.ug.per.mL*2.77/100, - Group == 'Cyanobacteria' ~ BM.ug.per.mL*0.02/100, - Group == 'Green Algae' ~ BM.ug.per.mL*0.52/100, - Group == 'Cryptophytes' ~ BM.ug.per.mL*2.13/100, - Group == 'Dinoflagellates' ~ BM.ug.per.mL*3.82/100, - TRUE ~ 0)) +## Import data from Processing Scripts +```{r import data} +# Read in RData files processed in the data processing Rmd file +load(file = here("phyto-final","RData","df_phyto.RData")) +load(file = here("phyto-final","RData","df_phyto_FASTR_sum.RData")) +load(file = here("phyto-final","RData","df_phyto_FASTR_gen.RData")) +load(file = here("phyto-final","RData","df_phyto_FASTR_gen_BV.RData")) +load(file = here("phyto-final","RData","df_phyto_FASTR_grp.RData")) +load(file = here("phyto-final","RData","df_phyto_FASTR_grp_BV.RData")) +load(file = here("phyto-final","RData","df_phyto_NMDS.RData")) ``` ## Create Plots for Paper -```{r Total Biovolume Plots} - -# Set plot output directory -plots <- here("phyto-final","plots") +```{r Total Biovolume Plots, echo=FALSE, message=FALSE} # Plot Upstream and Downstream Total BV by ActionPhase for each year fig1 <- ggplot(data = df_phyto_FASTR_sum, @@ -736,70 +182,9 @@ for (year in years) { ``` -## Multivariate Calculations - -```{r multivariate calcs} - -## Calculate NMDS axes -## Create Biovolume-only data frame at genus level -df_phyto_FASTR_gen_BV <- df_phyto_FASTR_gen %>% select(Year:Genus,BV.um3.per.mL) - -## Generate NMDS data with metaMDS by each year separately - -#df_phyto_FASTR_gen_BV$Year <- as.integer(df_phyto_FASTR_gen_BV$Year) +## NMDS Plots -years <- unique(df_phyto_FASTR_gen_BV$Year) -years <- sort(years, decreasing = F, na.last = T) - -## Create blank data frame to fill stresses in -stresses <- data.frame(Year = years, stress = NA) -ls_dfs <- list() - -for (i in 1:length(years)) { - genw <- pivot_wider(df_phyto_FASTR_gen_BV, - names_from = "Genus", - values_from = "BV.um3.per.mL", - values_fill = 0) - - genw <- genw %>% filter(Year == years[i]) - - #look at number of observations per station - #table(genw$StationCode) - - # Calculate the nMDS using vegan - # A good rule of thumb: stress < 0.05 provides an excellent representation in reduced dimensions, - # < 0.1 is great, < 0.2 is good/ok, and stress < 0.3 provides a poor representation. - phyto.NMDS <- metaMDS( - comm = genw[c(7:139)], - distance = "bray", - k = 3, - trymax = 500 - #trace = F, - #autotransform = F - ) - - stresses$stress[which(stresses$Year == years[i])] <- phyto.NMDS$stress - - #look at Shepard plot which shows scatter around the regression between the interpoint distances - #in the final configuration (i.e., the distances between each pair of communities) against their - #original dissimilarities. - stressplot(phyto.NMDS) - - # Using the scores function from vegan to extract the site scores and convert to a data.frame - data.scores <- as_tibble(scores(phyto.NMDS, display = "sites")) - - # Combine metadata with NMDS data scores to plot in ggplot - meta <- genw %>% select(1:6) - meta <- cbind(meta, data.scores) - - # Read in years as a character otherwise it shows up as a number and gets displayed as a gradient - meta$Year <- as.character(meta$Year) - - ls_dfs[[i]] <- meta - -} - -df_phyto_NMDS <- do.call(rbind, ls_dfs) +```{r NMDS, echo=FALSE, message=FALSE} # Plot NMDS Results @@ -961,5 +346,17 @@ bd <- betadisper(dm_phyto_gen, genw$Region) anova(bd) permutest(bd) +``` +## Create map of EMP stations for Supplemental Material + +```{r, maps, echo=FALSE, message=FALSE} + +suppressWarnings(suppressMessages(library(deltamapr))) +suppressWarnings(suppressMessages(library(sf))) +suppressWarnings(suppressMessages(library(ggrepel))) +suppressWarnings(suppressMessages(library(maps))) + + + ``` diff --git a/phyto-final/phyto-data-processing-FASTR.Rmd b/phyto-final/phyto-data-processing-FASTR.Rmd new file mode 100644 index 0000000..bb5155c --- /dev/null +++ b/phyto-final/phyto-data-processing-FASTR.Rmd @@ -0,0 +1,659 @@ +--- +title: "FASTR - Phytoplankton Data Processing" +author: "Ted Flynn" +date: "2024-03-21" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) + +suppressWarnings(suppressMessages(library(tidyverse))) +suppressWarnings(suppressMessages(library(knitr))) +suppressWarnings(suppressMessages(library(lubridate))) +suppressWarnings(suppressMessages(library(RColorBrewer))) +suppressWarnings(suppressMessages(library(vegan))) +suppressWarnings(suppressMessages(library(janitor))) +suppressWarnings(suppressMessages(library(here))) + +``` + +## Import phyto data collected for FASTR project + +These steps import raw phytoplankton data collected by DWR's Aquatic Ecology group as part of the FASTR project and combines them into a single data frame. + +```{r import FASTR, echo=FALSE, message=FALSE} +# Import AEU data files (comment out when finished) + +files_AEU <- list.files(path = here("phyto-final","data","csv"), + pattern = "\\.csv", + full.names = T) + +df_phyto_AEU <- map(files_AEU, ~read_csv(.x, show_col_types = FALSE)) %>% + list_rbind() + +# Remove empty rows +df_phyto_AEU <- df_phyto_AEU %>% filter_all(any_vars(!is.na(.))) + +# Remove weird row with only a single zero in biovolume +df_phyto_AEU <- df_phyto_AEU %>% drop_na(MethodCode) + +# Clean up column names +df_phyto_AEU <- df_phyto_AEU %>% clean_names(case = "big_camel") + +# Filter out samples from other projects. Read in station names with flags and +# merge together. +df_keepers <- read_csv(here("phyto-final","CSVs", "station_names_flagged.csv"), + show_col_types = FALSE) + +df_phyto_AEU <- left_join(df_phyto_AEU, df_keepers) + +rm(df_keepers) + +# Remove data unrelated to project +df_phyto_AEU <- df_phyto_AEU %>% filter(Flag == "Keep") + +# Remove flag column +df_phyto_AEU <- df_phyto_AEU %>% select(!(Flag)) + +# Confirm station IDs +unique(df_phyto_AEU$StationCode) # Shows only 10 FASTR stations +table(df_phyto_AEU$StationCode) + +# Remove blank columns +df_phyto_AEU <- df_phyto_AEU %>% select_if(~ !all(is.na(.))) + +# Remove individual measurement columns as they are only present in a small +# number of samples. +df_phyto_AEU <- df_phyto_AEU %>% select(!(Dimension:DepthM)) + +# Remove secondary and tertiary GALD measurements as they don't vary much and +# aren't present in EMP data +df_phyto_AEU <- df_phyto_AEU %>% select(!(Gald2:Gald3)) + +# Correct GALD, which is imported into two separate columns. Test to see if NAs +# are 'either/or' and that there aren't some rows with a value in both GALD and +# GALD 1 + +sum(is.na(df_phyto_AEU$Gald)) # Total is 1791 +sum(is.na(df_phyto_AEU$Gald1)) # Total is 4535 + +# Sum of NAs is 6326 which is the same as the number of rows in the df. This +# shows that there aren't any rows with two values so that we can Combine them +# without any issues. + +# Move Gald1 column +df_phyto_AEU <- df_phyto_AEU %>% relocate(Gald1, .after = Gald) + +# Combine both GALD columns +df_phyto_AEU <- df_phyto_AEU %>% + rowwise() %>% + mutate(GALD.Tot = sum(c_across(Gald:Gald1), na.rm = TRUE)) + +# Check if rows have two NAs or no NAs in the GALD columns +# test <- df_phyto_AEU %>% +# mutate(GALD.Value = case_when(is.na(Gald) & is.na(Gald1) ~ "Fix", +# TRUE ~ "Okay")) + +# Remove old GALD columns and rename GALD.Tot +df_phyto_AEU <- df_phyto_AEU %>% + select(!(Gald:Gald1)) %>% + rename("GALD" = "GALD.Tot") + +# Remove MethodCode column that just says "Phyto" +df_phyto_AEU <- df_phyto_AEU %>% select(!(MethodCode)) + +# Add column to indicate what study it came from +df_phyto_AEU <- df_phyto_AEU %>% mutate(Study = "FASTR", .after = StationCode) + +``` + +## Import EMP Data + +Import phytoplankton data collected by the Environmental Monitoring Program. + +```{r import EMP, echo=FALSE, message=FALSE} + +# Import EMP data files +phyto_files_EMP <- list.files(path = here("phyto-final","data","EMP","csv"), + pattern = "\\.csv", + full.names = T) + +df_phyto_EMP <- map(phyto_files_EMP, ~read_csv(.x, show_col_types = FALSE)) %>% + list_rbind() + +# Read in files with non-standard headers individually +df_Dec2021 <- read_csv(here("phyto-final","data","EMP","oddballs","December 2021.csv"), + show_col_types = FALSE) +df_Nov2021 <- read_csv(here("phyto-final","data","EMP","oddballs","November 2021.csv"), + show_col_types = FALSE) +df_Sep2013 <- read_csv(here("phyto-final","data","EMP","oddballs","September 2013.csv"), + show_col_types = FALSE) +df_Nov2013 <- read_csv(here("phyto-final","data","EMP","oddballs","November 2013.csv"), + show_col_types = FALSE) + +# Combine like oddball dfs +df_phyto2013 <- bind_rows(df_Sep2013, df_Nov2013) +df_phyto2021 <- bind_rows(df_Dec2021, df_Nov2021) + +# Remove individual dfs +rm(df_Dec2021) +rm(df_Nov2021) +rm(df_Nov2013) +rm(df_Sep2013) + +# Rename headers to match standard BSA headers Oddballs actually have the +# "correct" name of Total Cells rather than the incorrect "Number of cells per +# unit" + +df_phyto2013 <- df_phyto2013 %>% + rename("Number of cells per unit" = "Total Cells Counted") + +df_phyto2021 <- df_phyto2021 %>% + rename("Number of cells per unit" = "Total Number of Cells") %>% + rename("Unit Abundance" = "Unit Abundance (# of Natural Units)") + +# Combine oddball files with others +df_phyto_EMP <- bind_rows(df_phyto_EMP, df_phyto2013) +df_phyto_EMP <- bind_rows(df_phyto_EMP, df_phyto2021) + +# Remove unneeded dfs +rm(df_phyto2013) +rm(df_phyto2021) + +# Remove empty rows +df_phyto_EMP <- df_phyto_EMP %>% filter_all(any_vars(!is.na(.))) + +# Correct GALD, which is imported into two separate columns +# Test to see if NAs are 'either/or' and that there aren't some rows with a +# value in both GALD and GALD 1 + +sum(is.na(df_phyto_EMP$GALD)) # Total is 5880 +sum(is.na(df_phyto_EMP$`GALD 1`)) # Total is 8072 + +# Sum of NAs is 13952 which is the same as the number of rows in the df. +# This shows that there aren't any rows with two values so that we can +# Combine them without any issues. + +# Move GALD header +df_phyto_EMP <- df_phyto_EMP %>% relocate(`GALD 1`, .after = GALD) + +# Combine both GALD columns +df_phyto_EMP <- df_phyto_EMP %>% + rowwise() %>% + mutate(GALD.Tot = sum(c_across(GALD:`GALD 1`), na.rm = TRUE)) + +# Remove old GALD columns and rename GALD.Tot +df_phyto_EMP <- df_phyto_EMP %>% + select(!(GALD:`GALD 1`)) %>% + rename("GALD" = "GALD.Tot") + +# Clean up column names +df_phyto_EMP <- df_phyto_EMP %>% clean_names(case = "big_camel") + +df_phyto_EMP <- df_phyto_EMP %>% rename("GALD" = "Gald") + +# Remove blank columns +df_phyto_EMP <- df_phyto_EMP %>% select_if(~ !all(is.na(.))) + +# Remove columns that just have the method code "Phyto" as well as +# pre-calculated organisms per mL +df_phyto_EMP <- df_phyto_EMP %>% select(SampleDate:Biovolume10,GALD) + +# Add column to indicate what study it came from +df_phyto_EMP <- df_phyto_EMP %>% mutate(Study = "EMP", .after = StationCode) + +``` + +## Combine EMP and AEU data + +```{r combine data, echo=FALSE, message=FALSE} + +df_phyto <- bind_rows(df_phyto_EMP, df_phyto_AEU) + +# Average all 10 biovolume measurements for each taxon +df_phyto <- df_phyto %>% + rowwise() %>% + mutate(BV.Avg = mean(c_across(Biovolume1:Biovolume10), na.rm = T)) %>% + select(!(Biovolume1:Biovolume10)) # Remove Individual Biovolume Columns + +# Remove unneeded columns +df_phyto <- df_phyto %>% select(!c("BsaTin","DiatomSoftBody")) +df_phyto <- df_phyto %>% select(!(ColonyFilamentIndividualGroupCode:Shape)) +df_phyto <- df_phyto %>% select(!(VolumeReceivedML:NumberOfFieldsCounted)) + +# Fix samples with missing times +# Just replace with 12:00:00 b/c aren't doing any time-based analyses +df_phyto <- df_phyto %>% replace_na(list(SampleTime = "12:00:00")) + +# Get dates in the right format. Some are 1/1/14 and others 1/1/2014. +df_phyto$SampleDate <- mdy(df_phyto$SampleDate) + +# Combine date and time column +df_phyto <- df_phyto %>% unite(DateTime, c("SampleDate","SampleTime"), sep = " ") #, remove = FALSE, na.rm = FALSE) + +df_phyto$DateTime <- as_datetime(df_phyto$DateTime, + tz = "US/Pacific", + format = c("%Y-%m-%d %H:%M:%OS")) + +# Check for missing dates +df_phyto %>% filter(is.na(DateTime)) # No missing dates + +# Correct BSA header +df_phyto <- df_phyto %>% rename("TotalCells" = "NumberOfCellsPerUnit") + +# Calculate Unit Density & Biovolume Density +df_phyto <- df_phyto %>% + mutate(Units.per.mL = UnitAbundance * Factor) %>% + mutate(BV.um3.per.mL= TotalCells * BV.Avg * Factor) %>% + mutate(Cells.per.mL = TotalCells * Factor) + +# Remove columns used for density calcs +# Remove Biovolume Columns +df_phyto <- df_phyto %>% select(!(c(Factor,UnitAbundance:TotalCells))) + +# Add column for year and month for highlighting data +df_phyto <- df_phyto %>% + mutate(Year = year(df_phyto$DateTime)) %>% + mutate(Month = month(df_phyto$DateTime, label = T)) + +# Order month in calendar order rather than (default) alphabetical +df_phyto$Month = factor(df_phyto$Month, levels = month.abb) + +# Reorder date/time columns +df_phyto <- df_phyto %>% + relocate(Year, .after = DateTime) %>% + relocate(Month, .after = DateTime) + +# Remove 2013 data (FASTR data is very limited) +df_phyto <- df_phyto %>% filter(Year != 2013) + +``` + +## Fix Station Names + +```{r data cleaning, echo=FALSE, message=FALSE} + +# Fix EMP site names +df_phyto$StationCode <- gsub("EZ6 SAC","EZ6",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ6SAC","EZ6",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ6SJR","EZ6-SJR",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ2SAC","EZ2",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ2 SAC","EZ2",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ2SJR","EZ2-SJR",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ2 SJR","EZ2-SJR",df_phyto$StationCode) +df_phyto$StationCode <- gsub("EZ6 SJR","EZ6-SJR",df_phyto$StationCode) +df_phyto$StationCode <- gsub("D16-Twitchell","D16",df_phyto$StationCode) +df_phyto$StationCode <- gsub("D16-Twitchel","D16",df_phyto$StationCode) +df_phyto$StationCode <- gsub("D16 - Twitchell","D16",df_phyto$StationCode) +df_phyto$StationCode <- gsub("D16 Twitchell","D16",df_phyto$StationCode) +df_phyto$StationCode <- gsub("NZ328","NZ325",df_phyto$StationCode) # Typo in August 2019 +df_phyto$StationCode <- gsub("C3A-HOOD","C3A",df_phyto$StationCode) +df_phyto$StationCode <- gsub("C3A Hood","C3A",df_phyto$StationCode) +df_phyto$StationCode <- gsub("C3A- Hood","C3A",df_phyto$StationCode) +df_phyto$StationCode <- gsub("C3A-Hood","C3A",df_phyto$StationCode) +df_phyto$StationCode <- gsub("NZ542","NZS42",df_phyto$StationCode) +df_phyto$StationCode <- gsub("E26","EZ6",df_phyto$StationCode) +df_phyto$StationCode <- gsub("E22","EZ2",df_phyto$StationCode) # Typo in May 2018 + +# Fix AEU site names +df_phyto$StationCode <- gsub("SHER","SHR",df_phyto$StationCode) +df_phyto$StationCode <- gsub("BL-5","BL5",df_phyto$StationCode) + +# Remove extra Microcystis tows at D19 +df_phyto <- df_phyto %>% filter(StationCode != "D19 MC Tow") + +``` + +## Add Taxonomy Data and Fix Errors + +In Fall 2016, taxonomists began classifying the species Chroococcus microscopicus as Eucapsis microscopica. This is one of the most dominant species in this samples, so all taxa previously classified asC. microscopicus will be re-named E. microscopica. Also, the taxon Plagioselmis lacustris is inconsistently named, appearing sometimes as Rhodomonas lacustris. Change to Rhodomonas lacustris to avoid having these show up as separate genera. + +We also use a CSV of group classification names matched up with the genera occurring here to add higher-level taxonomic resolution. + +```{r taxonomy, echo=FALSE, message=FALSE} + +# Rename Chroococcus microscopicus as Eucapsis microscopica +df_phyto <- df_phyto %>% + mutate(Taxon = case_when(Taxon == 'Chroococcus microscopicus' ~ 'Eucapsis microscopica', + TRUE ~ Taxon)) %>% + mutate(Genus = case_when(Taxon == 'Eucapsis microscopica' ~ 'Eucapsis', + TRUE ~ Genus)) + +# Rename Plagioselmis lacustris as Rhodomonas lacustris +df_phyto <- df_phyto %>% + mutate(Taxon = case_when(Taxon == 'Plagioselmis lacustris' ~ 'Rhodomonas lacustris', + TRUE ~ Taxon)) %>% + mutate(Genus = case_when(Taxon == 'Rhodomonas lacustris' ~ 'Rhodomonas', + TRUE ~ Genus)) + +# Correct the genus label for a Chlorella entry +df_phyto$Genus <- gsub("cf Chlorella","Chlorella",df_phyto$Genus) + +# Read in CSV with manually-added WoRMS classification +df_taxa <- read_csv(here("phyto-final","CSVs","phyto_group_classification.csv"), + show_col_types = FALSE) + +df_phyto <- left_join(df_phyto, df_taxa) + +# Check if any groups are missing +sum(is.na(df_phyto$Group)) # none missing +rm(df_taxa) + +# Reorder Group classification column +df_phyto <- df_phyto %>% + relocate(Group, .before = Taxon) + +# Read in Region data for EMP stations +df_region <- read_csv(here("phyto-final","CSVs","station_regions.csv"), + show_col_types = FALSE) + +# Combine region data and phyto data +df_phyto <- left_join(df_phyto, df_region) + +# Reorder column +df_phyto <- df_phyto %>% relocate(Region, .after = StationCode) + +# check if there are any NAs in Region after the join +table(is.na(df_phyto$Region)) # no NAs + +# Remove data frames with old data +rm(df_region) +rm(df_phyto_EMP) +rm(df_phyto_AEU) + +``` +## Estimate Phytoplankton Biomass and LCEFA Abundance + +Convert biovolume to biomass using relationships from Menden-Deuer and Lussard (2000, doi: 10.4319/lo.2000.45.3.0569). This is only relevant to group-level data. Use biovolume measurements to calculate the abundance of long-chain essential fatty acids for each group based on the method in Galloway & Winder (2015, doi: 10.1371/journal.pone.0130053). + + +```{r biomass calcs, echo=FALSE, message=FALSE} + +# Calculate Biomass (pg-C per cell) from Biovolume (um^3 per cell) +df_phyto <- df_phyto %>% + mutate(Biomass.pg.C = case_when(Group == "Diatoms" ~ 0.288 * (BV.Avg^0.811), + Group != "Diatoms" ~ 0.216 * (BV.Avg^0.939))) + +# Convert pg to ug (ug-C per mL) +df_phyto <- df_phyto %>% + mutate(Biomass.ug.C = Biomass.pg.C / 10^6, .keep = "unused") + +# Calculate Biomass Density (pg-C per mL) +df_phyto <- df_phyto %>% + mutate(BM.ug.per.mL = Biomass.ug.C * Cells.per.mL) + +df_phyto <- df_phyto %>% + mutate(LCEFA.per.mL = case_when(Group == 'Diatoms' ~ BM.ug.per.mL*2.77/100, + Group == 'Cyanobacteria' ~ BM.ug.per.mL*0.02/100, + Group == 'Green Algae' ~ BM.ug.per.mL*0.52/100, + Group == 'Cryptophytes' ~ BM.ug.per.mL*2.13/100, + Group == 'Dinoflagellates' ~ BM.ug.per.mL*3.82/100, + TRUE ~ 0)) + +``` + +## Analyze FASTR Data for Outliers + +```{r outlier analysis, echo=FALSE, message=FALSE} + +# Subset combined dataset to FASTR-only +df_phyto_FASTR <- df_phyto %>% filter(Study == "FASTR") + +# Identify if outliers are present +boxplot(df_phyto_FASTR$BV.um3.per.mL) + +# Identify 0.1% and 99.9% percentiles of the data +quartiles <- quantile(df_phyto_FASTR$BV.um3.per.mL, + probs=c(0.001, 0.999), + na.rm = TRUE) + +# List upper cutoff (99.9%) +cutoff <- quartiles[2] + +# Filter FASTR dataset to show all taxa above this 99.9% cutoff +df_outliers <- df_phyto_FASTR %>% + filter(BV.um3.per.mL > cutoff) + +df_outliers %>% + select(DateTime,Taxon,BV.um3.per.mL) %>% + kable(caption = "Al")## 4 of top 6 are Spirogyra. + +df_phyto_FASTR %>% filter(Genus == "Spirogyra")# Only 5 total samples w/ Spirogyra + +# Decision: Remove Spirogyra from datasets +# 4/5 are above 99.9% cutoff, 5th is also highly abundant +# These are benthic algae that clump, more likely to be "bullseye" samples +# that got a big blob or clump + +# Remove Spirogyra from FASTR and combined dataset +df_phyto_FASTR <- df_phyto_FASTR %>% + filter(Genus != "Spirogyra") + +df_phyto <- df_phyto %>% + filter(Genus != "Spirogyra") + +``` + +## Add FASTR-Specific Metadata + +This section adds metadata specific to the FASTR project to a subset of the data. EMP data is not included in this data frame as the categories do not apply to it. + +Code for Merging Flow Designations was provided by Cat Pien and includes 45 days limit on either end of the flow pulse. Code was tweaked by TF in May 2022 to include greater-than or equal-to signs avoid inadvertently removing samples falling on the exact dates of pre- and post-flow phases. + +```{r FASTR metadata, echo=FALSE, message=FALSE} + +# Set levels so that Upstream stations appear before Downstream in graphs +df_phyto_FASTR$Region <- factor(df_phyto_FASTR$Region, + levels = c("Upstream","Downstream")) + +# Import flow pulse dates and categories +df_flow_designation <- read_csv(here("phyto-final", + "CSVs", + "FlowDatesDesignations_45days.csv"), + show_col_types = FALSE) + +# Remove column with net flow days. Won't be used in any analyses. +df_flow_designation <- df_flow_designation %>% select(!(NetFlowDays)) + +# Format dates as dates in each of the 4 date columns +df_flow_designation <- df_flow_designation %>% + mutate(across(PreFlowStart:PostFlowEnd, mdy)) + +# Merge data from FlowDesignation Table (Water Year Type, Flow days and type) +df_phyto_FASTR <- inner_join(df_phyto_FASTR, + df_flow_designation, + by = "Year") + +# Label samples with the appropriate phase (Pre-During-Post) and filter +# FASTR phyto data to include only samples collected within the flow pulse window +df_phyto_FASTR <- df_phyto_FASTR %>% + mutate(ActionPhase = ifelse(DateTime > PreFlowStart & DateTime < PreFlowEnd, "Before", NA)) %>% + mutate(ActionPhase = replace(ActionPhase, DateTime >= PreFlowEnd & DateTime < PostFlowStart, "During")) %>% + mutate(ActionPhase = replace(ActionPhase, DateTime >= PostFlowStart & DateTime < PostFlowEnd, "After")) %>% + filter(!is.na(ActionPhase)) %>% + select(-c(PreFlowStart:PostFlowEnd)) + +# Order the new ActionPhase so that it plots in the order Pre < During < Post +phase.order <- c("Before","During","After") +df_phyto_FASTR$ActionPhase <- factor(as.character(df_phyto_FASTR$ActionPhase), + levels = phase.order) + +# Add in Flow Pulse Category Data +df_pulse_category <- read_csv(here("phyto-final","CSVs","FlowPulseType.csv"), + show_col_types = FALSE) + +df_pulse_category <- df_pulse_category %>% + rename("FlowPulseCategory" = "FlowPulseType") + +# Add Flow Pulse Category to phytoplankton data frame +df_phyto_FASTR <- left_join(df_phyto_FASTR, df_pulse_category, by = "Year") + +# Create a list of stations used for FASTR +stations <- c("RMB","RCS","RD22","I80","LIS","STTD","BL5","LIB","RYI","RVB") + +# Re-order stations from North to South +df_phyto_FASTR$StationCode <- factor(as.character(df_phyto_FASTR$StationCode), + levels = stations) +## Relocate metadata columns +df_phyto_FASTR <- df_phyto_FASTR %>% + relocate(FlowPulseCategory, .before = ActionPhase) %>% + relocate(WYType:ActionPhase, .before = GALD) %>% + relocate(Biomass.ug.C, .after = BV.Avg) + + +``` + +## Create New Data Frames summarizing by Genus and Algal Group + +```{r summarize phyto data, echo=FALSE, message=FALSE} + +# Summarize by algal group +df_phyto_FASTR_grp <- df_phyto_FASTR %>% + group_by(Year, Month, Region, ActionPhase, DateTime, StationCode, Group) %>% + summarize(across(Units.per.mL:LCEFA.per.mL, ~sum(.x, na.rm = TRUE))) %>% + ungroup + +# Summarize by genus +df_phyto_FASTR_gen <- df_phyto_FASTR %>% + group_by(Year, Month, Region, ActionPhase, DateTime, StationCode, Genus) %>% + summarize(across(Units.per.mL:LCEFA.per.mL, ~sum(.x, na.rm = TRUE))) %>% + ungroup + +# Set group display order +group.order <- c("Diatoms", + "Cyanobacteria", + "Green Algae", + "Cryptophytes", + "Ciliates", + "Dinoflagellates", + "Golden Algae", + "Other") + +df_phyto_FASTR_grp$Group <- factor(as.character(df_phyto_FASTR_grp$Group), + levels = group.order) + +# Summarize totals +df_phyto_FASTR_sum <- df_phyto_FASTR %>% + group_by(Year, Month, DateTime, StationCode, Region, WYType, FlowPulseType, FlowPulseCategory, ActionPhase) %>% + summarize(across(Units.per.mL:LCEFA.per.mL, ~sum(.x, na.rm = TRUE))) %>% + rename("Total.Units.per.mL" = "Units.per.mL") %>% + rename("Total.BV.per.mL" = "BV.um3.per.mL") %>% + rename("Total.Cells.per.mL" = "Cells.per.mL") %>% + rename("Total.BM.per.mL" = "BM.ug.per.mL") %>% + rename("Total.LCEFA.per.mL" = "LCEFA.per.mL") %>% + ungroup + +# Separate out data frames by density type and add zeros +# Biovolume Density + Group +df_phyto_FASTR_grp_BV <- df_phyto_FASTR_grp %>% + select(Year:Group,BV.um3.per.mL) + +# Add zeros for taxa that weren't detected +temp <- pivot_wider(df_phyto_FASTR_grp_BV, + names_from = Group, + values_from = BV.um3.per.mL, + values_fill = 0) + +df_phyto_FASTR_grp_BV <- pivot_longer(temp, + cols = Cyanobacteria:last_col(), + names_to = "Group", + values_to = "BV.um3.per.mL") + +# Create Biovolume-only data frame at genus level +df_phyto_FASTR_gen_BV <- df_phyto_FASTR_gen %>% + select(Year:Genus,BV.um3.per.mL) + +temp <- pivot_wider(df_phyto_FASTR_gen_BV, + names_from = Genus, + values_from = BV.um3.per.mL, + values_fill = 0) + +df_phyto_FASTR_gen_BV <- pivot_longer(temp, + cols = Chlorella:last_col(), + names_to = "Genus", + values_to = "BV.um3.per.mL") + +rm(temp) + +``` + +## Multivariate Calculations +Use the vegan package in R to calculate the Bray-Curtis similarity coefficients for phytoplankton community comparisons in each year (2014 - 2019). A good rule of thumb: stress < 0.05 provides an excellent representation in reduced dimensions, < 0.1 is great, < 0.2 is good/ok, and stress < 0.3 provides a poor representation. + +```{r multivariate, echo=FALSE, message=FALSE} + +set.seed(2390) + +# Generate NMDS data with metaMDS by each year separately +years <- unique(df_phyto_FASTR_gen_BV$Year) +years <- sort(years, decreasing = F, na.last = T) + +# Create blank data frame to fill stresses in +stresses <- data.frame(Year = years, stress = NA) +ls_dfs <- list() + +for (i in 1:length(years)) { + genw <- pivot_wider(df_phyto_FASTR_gen_BV, + names_from = "Genus", + values_from = "BV.um3.per.mL", + values_fill = 0) + + genw <- genw %>% + filter(Year == years[i]) + + # Calculate the nMDS using vegan + + phyto.NMDS <- metaMDS( + comm = genw[c(7:139)], + distance = "bray", + k = 3, + trymax = 500 + #trace = F, + #autotransform = F + ) + + stresses$stress[which(stresses$Year == years[i])] <- phyto.NMDS$stress + + #look at Shepard plot which shows scatter around the regression between the interpoint distances + #in the final configuration (i.e., the distances between each pair of communities) against their + #original dissimilarities. + stressplot(phyto.NMDS) + + # Using the scores function from vegan to extract the site scores and convert to a data.frame + data.scores <- as_tibble(scores(phyto.NMDS, display = "sites")) + + # Combine metadata with NMDS data scores to plot in ggplot + meta <- genw %>% select(1:6) + meta <- cbind(meta, data.scores) + + # Read in years as a character otherwise it shows up as a number and gets displayed as a gradient + meta$Year <- as.character(meta$Year) + + ls_dfs[[i]] <- meta + +} + + +df_phyto_NMDS <- do.call(rbind, ls_dfs) + +# Output stresses for us in plotting NMDS +write_csv(stresses, file = here("phyto-final","analyses","phyto_NMDS_stress.csv")) + +``` + +## Save RData files for use in analysis + +```{r save RData, echo=FALSE, message=FALSE} + +## Save data files +save(df_phyto, file = here("phyto-final","RData","df_phyto.RData")) +save(df_phyto_FASTR_sum, file = here("phyto-final","RData","df_phyto_FASTR_sum.RData")) +save(df_phyto_FASTR_gen, file = here("phyto-final","RData","df_phyto_FASTR_gen.RData")) +save(df_phyto_FASTR_gen_BV, file = here("phyto-final","RData","df_phyto_FASTR_gen_BV.RData")) +save(df_phyto_FASTR_grp, file = here("phyto-final","RData","df_phyto_FASTR_grp.RData")) +save(df_phyto_FASTR_grp_BV, file = here("phyto-final","RData","df_phyto_FASTR_grp_BV.RData")) +save(df_phyto_NMDS, file = here("phyto-final","RData","df_phyto_NMDS.RData")) + +``` + diff --git a/phyto-final/phyto-data-processing-FASTR.html b/phyto-final/phyto-data-processing-FASTR.html new file mode 100644 index 0000000..d61381c --- /dev/null +++ b/phyto-final/phyto-data-processing-FASTR.html @@ -0,0 +1,1802 @@ + + + + +
+ + + + + + + + + + +These steps import raw phytoplankton data collected by DWR’s Aquatic +Ecology group as part of the FASTR project and combines them into a +single data frame.
+## [1] "LIS" "STTD" "RVB" "RCS" "LIB" "RD22" "RYI" "BL5" "I80" "RMB"
+##
+## BL5 I80 LIB LIS RCS RD22 RMB RVB RYI STTD
+## 325 587 500 1413 440 564 27 544 656 1270
+## [1] 1791
+## [1] 4535
+Import phytoplankton data collected by the Environmental Monitoring +Program.
+## [1] 5880
+## [1] 8072
+## # A tibble: 0 × 11
+## # ℹ 11 variables: DateTime <dttm>, StationCode <chr>, Study <chr>,
+## # Factor <dbl>, Taxon <chr>, Genus <chr>, Species <chr>, UnitAbundance <dbl>,
+## # NumberOfCellsPerUnit <dbl>, GALD <dbl>, BV.Avg <dbl>
+In Fall 2016, taxonomists began classifying the species Chroococcus +microscopicus as Eucapsis microscopica. This is one of the most dominant +species in this samples, so all taxa previously classified asC. +microscopicus will be re-named E. microscopica. Also, the taxon +Plagioselmis lacustris is inconsistently named, appearing sometimes as +Rhodomonas lacustris. Change to Rhodomonas lacustris to avoid having +these show up as separate genera.
+We also use a CSV of group classification names matched up with the +genera occurring here to add higher-level taxonomic resolution.
+## [1] 0
+##
+## FALSE
+## 19895
+Convert biovolume to biomass using relationships from Menden-Deuer +and Lussard (2000, doi: 10.4319/lo.2000.45.3.0569). This is only +relevant to group-level data. Use biovolume measurements to calculate +the abundance of long-chain essential fatty acids for each group based +on the method in Galloway & Winder (2015, doi: +10.1371/journal.pone.0130053).
+# Subset combined dataset to FASTR-only
+df_phyto_FASTR <- df_phyto %>% filter(Study == "FASTR")
+
+# Identify if outliers are present
+boxplot(df_phyto_FASTR$BV.um3.per.mL)
+
+# Identify 0.1% and 99.9% percentiles of the data
+quartiles <- quantile(df_phyto_FASTR$BV.um3.per.mL,
+ probs=c(0.001, 0.999),
+ na.rm = TRUE)
+
+# List upper cutoff (99.9%)
+cutoff <- quartiles[2]
+
+# Filter FASTR dataset to show all taxa above this 99.9% cutoff
+df_outliers <- df_phyto_FASTR %>%
+ filter(BV.um3.per.mL > cutoff)
+
+df_outliers %>%
+ select(DateTime,Taxon,BV.um3.per.mL) %>%
+ kable(caption = "Al")## 4 of top 6 are Spirogyra.
+DateTime | +Taxon | +BV.um3.per.mL | +
---|---|---|
2015-08-27 10:04:00 | +Spirogyra sp. | +53284341 | +
2015-09-01 11:55:00 | +Spirogyra sp. | +57390494 | +
2016-03-24 08:26:00 | +cf. Cyclotella sp. | +45565715 | +
2016-07-07 09:40:00 | +Cyclotella sp. | +31747497 | +
2016-09-13 12:17:00 | +Spirogyra sp. | +256335927 | +
2017-09-13 12:10:00 | +Spirogyra sp. | +76015935 | +
2018-08-02 10:15:00 | +Ulnaria ulna | +65089153 | +
df_phyto_FASTR %>% filter(Genus == "Spirogyra")# Only 5 total samples w/ Spirogyra
+## # A tibble: 5 × 18
+## DateTime Month Year StationCode Region Study Group Taxon Genus
+## <dttm> <ord> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
+## 1 2015-08-27 10:04:00 Aug 2015 LIB Downstream FASTR Gree… Spir… Spir…
+## 2 2015-09-01 11:55:00 Sep 2015 BL5 Downstream FASTR Gree… Spir… Spir…
+## 3 2016-08-02 12:53:00 Aug 2016 RVB Downstream FASTR Gree… Spir… Spir…
+## 4 2016-09-13 12:17:00 Sep 2016 RYI Downstream FASTR Gree… Spir… Spir…
+## 5 2017-09-13 12:10:00 Sep 2017 BL5 Downstream FASTR Gree… Spir… Spir…
+## # ℹ 9 more variables: Species <chr>, GALD <dbl>, BV.Avg <dbl>,
+## # Units.per.mL <dbl>, BV.um3.per.mL <dbl>, Cells.per.mL <dbl>,
+## # Biomass.ug.C <dbl>, BM.ug.per.mL <dbl>, LCEFA.per.mL <dbl>
+# Decision: Remove Spirogyra from datasets
+# 4/5 are above 99.9% cutoff, 5th is also highly abundant
+# These are benthic algae that clump, more likely to be "bullseye" samples
+# that got a big blob or clump
+
+# Remove Spirogyra from FASTR and combined dataset
+df_phyto_FASTR <- df_phyto_FASTR %>%
+ filter(Genus != "Spirogyra")
+
+df_phyto <- df_phyto %>%
+ filter(Genus != "Spirogyra")
+This section adds metadata specific to the FASTR project to a subset +of the data. EMP data is not included in this data frame as the +categories do not apply to it.
+Code for Merging Flow Designations was provided by Cat Pien and +includes 45 days limit on either end of the flow pulse. Code was tweaked +by TF in May 2022 to include greater-than or equal-to signs avoid +inadvertently removing samples falling on the exact dates of pre- and +post-flow phases.
+# Set levels so that Upstream stations appear before Downstream in graphs
+df_phyto_FASTR$Region <- factor(df_phyto_FASTR$Region,
+ levels = c("Upstream","Downstream"))
+
+# Import flow pulse dates and categories
+df_flow_designation <- read_csv(here("phyto-final",
+ "CSVs",
+ "FlowDatesDesignations_45days.csv"),
+ show_col_types = FALSE)
+
+# Remove column with net flow days. Won't be used in any analyses.
+df_flow_designation <- df_flow_designation %>% select(!(NetFlowDays))
+
+# Format dates as dates in each of the 4 date columns
+df_flow_designation <- df_flow_designation %>%
+ mutate(across(PreFlowStart:PostFlowEnd, mdy))
+
+# Merge data from FlowDesignation Table (Water Year Type, Flow days and type)
+df_phyto_FASTR <- inner_join(df_phyto_FASTR,
+ df_flow_designation,
+ by = "Year")
+
+# Label samples with the appropriate phase (Pre-During-Post) and filter
+# FASTR phyto data to include only samples collected within the flow pulse window
+df_phyto_FASTR <- df_phyto_FASTR %>%
+ mutate(ActionPhase = ifelse(DateTime > PreFlowStart & DateTime < PreFlowEnd, "Before", NA)) %>%
+ mutate(ActionPhase = replace(ActionPhase, DateTime >= PreFlowEnd & DateTime < PostFlowStart, "During")) %>%
+ mutate(ActionPhase = replace(ActionPhase, DateTime >= PostFlowStart & DateTime < PostFlowEnd, "After")) %>%
+ filter(!is.na(ActionPhase)) %>%
+ select(-c(PreFlowStart:PostFlowEnd))
+
+# Order the new ActionPhase so that it plots in the order Pre < During < Post
+phase.order <- c("Before","During","After")
+df_phyto_FASTR$ActionPhase <- factor(as.character(df_phyto_FASTR$ActionPhase),
+ levels = phase.order)
+
+# Add in Flow Pulse Category Data
+df_pulse_category <- read_csv(here("phyto-final","CSVs","FlowPulseType.csv"),
+ show_col_types = FALSE)
+
+df_pulse_category <- df_pulse_category %>%
+ rename("FlowPulseCategory" = "FlowPulseType")
+
+# Add Flow Pulse Category to phytoplankton data frame
+df_phyto_FASTR <- left_join(df_phyto_FASTR, df_pulse_category, by = "Year")
+
+# Create a list of stations used for FASTR
+stations <- c("RMB","RCS","RD22","I80","LIS","STTD","BL5","LIB","RYI","RVB")
+
+# Re-order stations from North to South
+df_phyto_FASTR$StationCode <- factor(as.character(df_phyto_FASTR$StationCode),
+ levels = stations)
+## Relocate metadata columns
+df_phyto_FASTR <- df_phyto_FASTR %>%
+ relocate(FlowPulseCategory, .before = ActionPhase) %>%
+ relocate(WYType:ActionPhase, .before = GALD) %>%
+ relocate(Biomass.ug.C, .after = BV.Avg)
+# Summarize by algal group
+df_phyto_FASTR_grp <- df_phyto_FASTR %>%
+ group_by(Year, Month, Region, ActionPhase, DateTime, StationCode, Group) %>%
+ summarize(across(Units.per.mL:LCEFA.per.mL, ~sum(.x, na.rm = TRUE))) %>%
+ ungroup
+## `summarise()` has grouped output by 'Year', 'Month', 'Region', 'ActionPhase',
+## 'DateTime', 'StationCode'. You can override using the `.groups` argument.
+# Summarize by genus
+df_phyto_FASTR_gen <- df_phyto_FASTR %>%
+ group_by(Year, Month, Region, ActionPhase, DateTime, StationCode, Genus) %>%
+ summarize(across(Units.per.mL:LCEFA.per.mL, ~sum(.x, na.rm = TRUE))) %>%
+ ungroup
+## `summarise()` has grouped output by 'Year', 'Month', 'Region', 'ActionPhase',
+## 'DateTime', 'StationCode'. You can override using the `.groups` argument.
+# Set group display order
+group.order <- c("Diatoms",
+ "Cyanobacteria",
+ "Green Algae",
+ "Cryptophytes",
+ "Ciliates",
+ "Dinoflagellates",
+ "Golden Algae",
+ "Other")
+
+df_phyto_FASTR_grp$Group <- factor(as.character(df_phyto_FASTR_grp$Group),
+ levels = group.order)
+
+# Summarize totals
+df_phyto_FASTR_sum <- df_phyto_FASTR %>%
+ group_by(Year, Month, DateTime, StationCode, Region, WYType, FlowPulseType, FlowPulseCategory, ActionPhase) %>%
+ summarize(across(Units.per.mL:LCEFA.per.mL, ~sum(.x, na.rm = TRUE))) %>%
+ rename("Total.Units.per.mL" = "Units.per.mL") %>%
+ rename("Total.BV.per.mL" = "BV.um3.per.mL") %>%
+ rename("Total.Cells.per.mL" = "Cells.per.mL") %>%
+ rename("Total.BM.per.mL" = "BM.ug.per.mL") %>%
+ rename("Total.LCEFA.per.mL" = "LCEFA.per.mL") %>%
+ ungroup
+## `summarise()` has grouped output by 'Year', 'Month', 'DateTime', 'StationCode',
+## 'Region', 'WYType', 'FlowPulseType', 'FlowPulseCategory'. You can override
+## using the `.groups` argument.
+# Separate out data frames by density type and add zeros
+# Biovolume Density + Group
+df_phyto_FASTR_grp_BV <- df_phyto_FASTR_grp %>%
+ select(Year:Group,BV.um3.per.mL)
+
+# Add zeros for taxa that weren't detected
+temp <- pivot_wider(df_phyto_FASTR_grp_BV,
+ names_from = Group,
+ values_from = BV.um3.per.mL,
+ values_fill = 0)
+
+df_phyto_FASTR_grp_BV <- pivot_longer(temp,
+ cols = Cyanobacteria:last_col(),
+ names_to = "Group",
+ values_to = "BV.um3.per.mL")
+
+# Create Biovolume-only data frame at genus level
+df_phyto_FASTR_gen_BV <- df_phyto_FASTR_gen %>%
+ select(Year:Genus,BV.um3.per.mL)
+
+temp <- pivot_wider(df_phyto_FASTR_gen_BV,
+ names_from = Genus,
+ values_from = BV.um3.per.mL,
+ values_fill = 0)
+
+df_phyto_FASTR_gen_BV <- pivot_longer(temp,
+ cols = Chlorella:last_col(),
+ names_to = "Genus",
+ values_to = "BV.um3.per.mL")
+
+rm(temp)
+Use the vegan package in R to calculate the Bray-Curtis similarity +coefficients for phytoplankton community comparisons in each year (2014 +- 2019). A good rule of thumb: stress < 0.05 provides an excellent +representation in reduced dimensions, < 0.1 is great, < 0.2 is +good/ok, and stress < 0.3 provides a poor representation.
+## Square root transformation
+## Wisconsin double standardization
+## Run 0 stress 0.171739
+## Run 1 stress 0.1799737
+## Run 2 stress 0.1742501
+## Run 3 stress 0.1716342
+## ... New best solution
+## ... Procrustes: rmse 0.02073212 max resid 0.115767
+## Run 4 stress 0.1727315
+## Run 5 stress 0.1742501
+## Run 6 stress 0.1723086
+## Run 7 stress 0.1824937
+## Run 8 stress 0.1774677
+## Run 9 stress 0.1716779
+## ... Procrustes: rmse 0.01221018 max resid 0.05862438
+## Run 10 stress 0.1725571
+## Run 11 stress 0.171739
+## ... Procrustes: rmse 0.0207781 max resid 0.1158136
+## Run 12 stress 0.1716342
+## ... New best solution
+## ... Procrustes: rmse 0.0002862936 max resid 0.0005694585
+## ... Similar to previous best
+## Run 13 stress 0.1723087
+## Run 14 stress 0.1774723
+## Run 15 stress 0.1742608
+## Run 16 stress 0.1742503
+## Run 17 stress 0.1716342
+## ... New best solution
+## ... Procrustes: rmse 0.0002895222 max resid 0.0008556222
+## ... Similar to previous best
+## Run 18 stress 0.1726868
+## Run 19 stress 0.1740751
+## Run 20 stress 0.1716779
+## ... Procrustes: rmse 0.01197123 max resid 0.05782632
+## *** Best solution repeated 1 times
+
+## Square root transformation
+## Wisconsin double standardization
+## Run 0 stress 0.2131056
+## Run 1 stress 0.2131107
+## ... Procrustes: rmse 0.000519833 max resid 0.003952992
+## ... Similar to previous best
+## Run 2 stress 0.2137802
+## Run 3 stress 0.2130855
+## ... New best solution
+## ... Procrustes: rmse 0.003409788 max resid 0.03065725
+## Run 4 stress 0.2133165
+## ... Procrustes: rmse 0.008312165 max resid 0.07816092
+## Run 5 stress 0.2135564
+## ... Procrustes: rmse 0.01354384 max resid 0.1364567
+## Run 6 stress 0.213105
+## ... Procrustes: rmse 0.003334044 max resid 0.03023528
+## Run 7 stress 0.2135563
+## ... Procrustes: rmse 0.01369097 max resid 0.1373975
+## Run 8 stress 0.215335
+## Run 9 stress 0.2131049
+## ... Procrustes: rmse 0.003315144 max resid 0.03026994
+## Run 10 stress 0.2133171
+## ... Procrustes: rmse 0.00848766 max resid 0.0798828
+## Run 11 stress 0.2147763
+## Run 12 stress 0.2130858
+## ... Procrustes: rmse 0.0001543822 max resid 0.0006121159
+## ... Similar to previous best
+## Run 13 stress 0.2133171
+## ... Procrustes: rmse 0.008489425 max resid 0.07996174
+## Run 14 stress 0.213105
+## ... Procrustes: rmse 0.003335774 max resid 0.03027651
+## Run 15 stress 0.2135562
+## ... Procrustes: rmse 0.01363277 max resid 0.1372427
+## Run 16 stress 0.2131539
+## ... Procrustes: rmse 0.004738893 max resid 0.03763404
+## Run 17 stress 0.2131759
+## ... Procrustes: rmse 0.005476212 max resid 0.04365836
+## Run 18 stress 0.2130857
+## ... Procrustes: rmse 0.000138796 max resid 0.0005762975
+## ... Similar to previous best
+## Run 19 stress 0.213556
+## ... Procrustes: rmse 0.01367729 max resid 0.1373718
+## Run 20 stress 0.2140375
+## *** Best solution repeated 2 times
+
+## Square root transformation
+## Wisconsin double standardization
+## Run 0 stress 0.1944156
+## Run 1 stress 0.1944144
+## ... New best solution
+## ... Procrustes: rmse 0.000382865 max resid 0.001909889
+## ... Similar to previous best
+## Run 2 stress 0.1944468
+## ... Procrustes: rmse 0.002756107 max resid 0.01413661
+## Run 3 stress 0.1944155
+## ... Procrustes: rmse 0.0008671 max resid 0.00682672
+## ... Similar to previous best
+## Run 4 stress 0.1944328
+## ... Procrustes: rmse 0.001707608 max resid 0.0118279
+## Run 5 stress 0.1944164
+## ... Procrustes: rmse 0.001614454 max resid 0.009112283
+## ... Similar to previous best
+## Run 6 stress 0.1944153
+## ... Procrustes: rmse 0.001445943 max resid 0.009471164
+## ... Similar to previous best
+## Run 7 stress 0.1944172
+## ... Procrustes: rmse 0.001745561 max resid 0.01017501
+## Run 8 stress 0.1944148
+## ... Procrustes: rmse 0.001000083 max resid 0.006661262
+## ... Similar to previous best
+## Run 9 stress 0.1944155
+## ... Procrustes: rmse 0.000381632 max resid 0.001980546
+## ... Similar to previous best
+## Run 10 stress 0.1944144
+## ... Procrustes: rmse 0.0001141491 max resid 0.0006230301
+## ... Similar to previous best
+## Run 11 stress 0.1944161
+## ... Procrustes: rmse 0.001485528 max resid 0.008536007
+## ... Similar to previous best
+## Run 12 stress 0.1944147
+## ... Procrustes: rmse 0.0008702444 max resid 0.004735455
+## ... Similar to previous best
+## Run 13 stress 0.1944146
+## ... Procrustes: rmse 0.001234144 max resid 0.00872337
+## ... Similar to previous best
+## Run 14 stress 0.1944326
+## ... Procrustes: rmse 0.001719908 max resid 0.01188171
+## Run 15 stress 0.1954796
+## Run 16 stress 0.2068974
+## Run 17 stress 0.1944161
+## ... Procrustes: rmse 0.001550354 max resid 0.008565922
+## ... Similar to previous best
+## Run 18 stress 0.1944146
+## ... Procrustes: rmse 0.001106216 max resid 0.008095053
+## ... Similar to previous best
+## Run 19 stress 0.1944177
+## ... Procrustes: rmse 0.001782505 max resid 0.01102172
+## Run 20 stress 0.194415
+## ... Procrustes: rmse 0.001367209 max resid 0.009455656
+## ... Similar to previous best
+## *** Best solution repeated 13 times
+
+## Square root transformation
+## Wisconsin double standardization
+## Run 0 stress 0.1895903
+## Run 1 stress 0.1897947
+## ... Procrustes: rmse 0.009975627 max resid 0.05265113
+## Run 2 stress 0.1906022
+## Run 3 stress 0.1897305
+## ... Procrustes: rmse 0.01772661 max resid 0.133512
+## Run 4 stress 0.1933988
+## Run 5 stress 0.1895894
+## ... New best solution
+## ... Procrustes: rmse 0.01388348 max resid 0.1203614
+## Run 6 stress 0.1897621
+## ... Procrustes: rmse 0.009261005 max resid 0.06637024
+## Run 7 stress 0.1899092
+## ... Procrustes: rmse 0.01950227 max resid 0.1236913
+## Run 8 stress 0.195562
+## Run 9 stress 0.1914162
+## Run 10 stress 0.1931959
+## Run 11 stress 0.1907902
+## Run 12 stress 0.1921355
+## Run 13 stress 0.194268
+## Run 14 stress 0.1939197
+## Run 15 stress 0.1933385
+## Run 16 stress 0.1897273
+## ... Procrustes: rmse 0.009105963 max resid 0.05559605
+## Run 17 stress 0.1923635
+## Run 18 stress 0.1961133
+## Run 19 stress 0.191121
+## Run 20 stress 0.1930784
+## Run 21 stress 0.1896396
+## ... Procrustes: rmse 0.006644986 max resid 0.0390687
+## Run 22 stress 0.1906213
+## Run 23 stress 0.1928022
+## Run 24 stress 0.191636
+## Run 25 stress 0.1906089
+## Run 26 stress 0.1916571
+## Run 27 stress 0.1898253
+## ... Procrustes: rmse 0.01696262 max resid 0.1223649
+## Run 28 stress 0.1955546
+## Run 29 stress 0.1932549
+## Run 30 stress 0.1909511
+## Run 31 stress 0.1935664
+## Run 32 stress 0.1976113
+## Run 33 stress 0.1917427
+## Run 34 stress 0.1917284
+## Run 35 stress 0.1942728
+## Run 36 stress 0.1936999
+## Run 37 stress 0.1945568
+## Run 38 stress 0.1898387
+## ... Procrustes: rmse 0.01659778 max resid 0.1208494
+## Run 39 stress 0.1931959
+## Run 40 stress 0.1979693
+## Run 41 stress 0.1897235
+## ... Procrustes: rmse 0.01124413 max resid 0.06557715
+## Run 42 stress 0.1898238
+## ... Procrustes: rmse 0.01683725 max resid 0.1223217
+## Run 43 stress 0.1923242
+## Run 44 stress 0.1926558
+## Run 45 stress 0.1898334
+## ... Procrustes: rmse 0.01513478 max resid 0.1219036
+## Run 46 stress 0.1913842
+## Run 47 stress 0.1897298
+## ... Procrustes: rmse 0.008975949 max resid 0.05563052
+## Run 48 stress 0.1921273
+## Run 49 stress 0.1918316
+## Run 50 stress 0.1980422
+## Run 51 stress 0.1896155
+## ... Procrustes: rmse 0.003701401 max resid 0.02233224
+## Run 52 stress 0.1906273
+## Run 53 stress 0.1922303
+## Run 54 stress 0.1919807
+## Run 55 stress 0.1932551
+## Run 56 stress 0.1930746
+## Run 57 stress 0.1901826
+## Run 58 stress 0.1907895
+## Run 59 stress 0.1897901
+## ... Procrustes: rmse 0.008930043 max resid 0.04971957
+## Run 60 stress 0.1917658
+## Run 61 stress 0.1932894
+## Run 62 stress 0.1917414
+## Run 63 stress 0.1915189
+## Run 64 stress 0.1942447
+## Run 65 stress 0.1925553
+## Run 66 stress 0.1949433
+## Run 67 stress 0.1910727
+## Run 68 stress 0.1897211
+## ... Procrustes: rmse 0.007269772 max resid 0.05297184
+## Run 69 stress 0.1905817
+## Run 70 stress 0.1932645
+## Run 71 stress 0.1912897
+## Run 72 stress 0.1910341
+## Run 73 stress 0.1915931
+## Run 74 stress 0.1944675
+## Run 75 stress 0.1895674
+## ... New best solution
+## ... Procrustes: rmse 0.01274907 max resid 0.1206398
+## Run 76 stress 0.1923927
+## Run 77 stress 0.1897941
+## ... Procrustes: rmse 0.007347698 max resid 0.05013562
+## Run 78 stress 0.1896217
+## ... Procrustes: rmse 0.01346878 max resid 0.1172194
+## Run 79 stress 0.1961435
+## Run 80 stress 0.1906074
+## Run 81 stress 0.1911673
+## Run 82 stress 0.1924783
+## Run 83 stress 0.1908723
+## Run 84 stress 0.1908301
+## Run 85 stress 0.1932223
+## Run 86 stress 0.1933609
+## Run 87 stress 0.1916363
+## Run 88 stress 0.1897222
+## ... Procrustes: rmse 0.01764094 max resid 0.1343025
+## Run 89 stress 0.191584
+## Run 90 stress 0.193264
+## Run 91 stress 0.1898609
+## ... Procrustes: rmse 0.009457616 max resid 0.07291282
+## Run 92 stress 0.1911669
+## Run 93 stress 0.1970428
+## Run 94 stress 0.1896411
+## ... Procrustes: rmse 0.01448483 max resid 0.1189741
+## Run 95 stress 0.1921884
+## Run 96 stress 0.1915692
+## Run 97 stress 0.1918064
+## Run 98 stress 0.189788
+## ... Procrustes: rmse 0.005484163 max resid 0.03886502
+## Run 99 stress 0.1935432
+## Run 100 stress 0.1939257
+## Run 101 stress 0.1921067
+## Run 102 stress 0.1972544
+## Run 103 stress 0.1926558
+## Run 104 stress 0.1897223
+## ... Procrustes: rmse 0.01775078 max resid 0.1344298
+## Run 105 stress 0.1909636
+## Run 106 stress 0.1916958
+## Run 107 stress 0.1915198
+## Run 108 stress 0.1915205
+## Run 109 stress 0.1921865
+## Run 110 stress 0.189622
+## ... Procrustes: rmse 0.0139468 max resid 0.1273695
+## Run 111 stress 0.1977597
+## Run 112 stress 0.1917267
+## Run 113 stress 0.1908725
+## Run 114 stress 0.1943388
+## Run 115 stress 0.1945873
+## Run 116 stress 0.1905559
+## Run 117 stress 0.1934113
+## Run 118 stress 0.1896145
+## ... Procrustes: rmse 0.01331169 max resid 0.1181038
+## Run 119 stress 0.1915871
+## Run 120 stress 0.1944296
+## Run 121 stress 0.1949918
+## Run 122 stress 0.1927801
+## Run 123 stress 0.1924308
+## Run 124 stress 0.1906174
+## Run 125 stress 0.1968921
+## Run 126 stress 0.1923624
+## Run 127 stress 0.1910576
+## Run 128 stress 0.1927729
+## Run 129 stress 0.195346
+## Run 130 stress 0.1896208
+## ... Procrustes: rmse 0.01383497 max resid 0.1266143
+## Run 131 stress 0.196324
+## Run 132 stress 0.1929855
+## Run 133 stress 0.1906194
+## Run 134 stress 0.1931814
+## Run 135 stress 0.1955323
+## Run 136 stress 0.1896211
+## ... Procrustes: rmse 0.01344674 max resid 0.1176706
+## Run 137 stress 0.1898426
+## ... Procrustes: rmse 0.008459036 max resid 0.05082267
+## Run 138 stress 0.1897562
+## ... Procrustes: rmse 0.01672158 max resid 0.1353291
+## Run 139 stress 0.1954845
+## Run 140 stress 0.1932209
+## Run 141 stress 0.1935605
+## Run 142 stress 0.1940391
+## Run 143 stress 0.1898237
+## ... Procrustes: rmse 0.01025345 max resid 0.05480253
+## Run 144 stress 0.1935867
+## Run 145 stress 0.1910097
+## Run 146 stress 0.1997389
+## Run 147 stress 0.1974073
+## Run 148 stress 0.189834
+## ... Procrustes: rmse 0.007557327 max resid 0.05140839
+## Run 149 stress 0.1897327
+## ... Procrustes: rmse 0.01646526 max resid 0.1340014
+## Run 150 stress 0.1915852
+## Run 151 stress 0.1905832
+## Run 152 stress 0.1948827
+## Run 153 stress 0.1912945
+## Run 154 stress 0.189589
+## ... Procrustes: rmse 0.005443428 max resid 0.04631368
+## Run 155 stress 0.1933951
+## Run 156 stress 0.1931611
+## Run 157 stress 0.1897277
+## ... Procrustes: rmse 0.01643828 max resid 0.133551
+## Run 158 stress 0.1896212
+## ... Procrustes: rmse 0.01344496 max resid 0.1177917
+## Run 159 stress 0.1896401
+## ... Procrustes: rmse 0.01429524 max resid 0.1176958
+## Run 160 stress 0.193601
+## Run 161 stress 0.191121
+## Run 162 stress 0.1914283
+## Run 163 stress 0.1896397
+## ... Procrustes: rmse 0.01422834 max resid 0.1178875
+## Run 164 stress 0.193289
+## Run 165 stress 0.1911224
+## Run 166 stress 0.1915844
+## Run 167 stress 0.1905808
+## Run 168 stress 0.1898968
+## ... Procrustes: rmse 0.01563165 max resid 0.1139724
+## Run 169 stress 0.1896182
+## ... Procrustes: rmse 0.01364735 max resid 0.1239508
+## Run 170 stress 0.1931552
+## Run 171 stress 0.1929347
+## Run 172 stress 0.192174
+## Run 173 stress 0.1898339
+## ... Procrustes: rmse 0.01133763 max resid 0.07001905
+## Run 174 stress 0.1966668
+## Run 175 stress 0.1898248
+## ... Procrustes: rmse 0.01068642 max resid 0.06096265
+## Run 176 stress 0.1906053
+## Run 177 stress 0.1905322
+## Run 178 stress 0.1991652
+## Run 179 stress 0.1914442
+## Run 180 stress 0.191957
+## Run 181 stress 0.1896407
+## ... Procrustes: rmse 0.01429016 max resid 0.116931
+## Run 182 stress 0.193258
+## Run 183 stress 0.1948732
+## Run 184 stress 0.1908872
+## Run 185 stress 0.1915658
+## Run 186 stress 0.1912438
+## Run 187 stress 0.1931942
+## Run 188 stress 0.2000997
+## Run 189 stress 0.1897219
+## ... Procrustes: rmse 0.01753967 max resid 0.134269
+## Run 190 stress 0.1915386
+## Run 191 stress 0.193229
+## Run 192 stress 0.1906084
+## Run 193 stress 0.1992432
+## Run 194 stress 0.1917564
+## Run 195 stress 0.1915527
+## Run 196 stress 0.1978064
+## Run 197 stress 0.192074
+## Run 198 stress 0.1897476
+## ... Procrustes: rmse 0.01639169 max resid 0.1336338
+## Run 199 stress 0.1936178
+## Run 200 stress 0.19058
+## Run 201 stress 0.1926558
+## Run 202 stress 0.1935054
+## Run 203 stress 0.1924323
+## Run 204 stress 0.193636
+## Run 205 stress 0.1915838
+## Run 206 stress 0.1896659
+## ... Procrustes: rmse 0.01517555 max resid 0.1380996
+## Run 207 stress 0.1906179
+## Run 208 stress 0.1898433
+## ... Procrustes: rmse 0.01181096 max resid 0.07488563
+## Run 209 stress 0.1906112
+## Run 210 stress 0.1896407
+## ... Procrustes: rmse 0.01392701 max resid 0.1168902
+## Run 211 stress 0.1914431
+## Run 212 stress 0.1913669
+## Run 213 stress 0.1937428
+## Run 214 stress 0.1909157
+## Run 215 stress 0.1953611
+## Run 216 stress 0.1906066
+## Run 217 stress 0.1935581
+## Run 218 stress 0.1932052
+## Run 219 stress 0.193196
+## Run 220 stress 0.1956965
+## Run 221 stress 0.1943427
+## Run 222 stress 0.1923993
+## Run 223 stress 0.1923191
+## Run 224 stress 0.1895892
+## ... Procrustes: rmse 0.01259026 max resid 0.1190304
+## Run 225 stress 0.1896226
+## ... Procrustes: rmse 0.01354072 max resid 0.1189862
+## Run 226 stress 0.1905839
+## Run 227 stress 0.1931957
+## Run 228 stress 0.1931356
+## Run 229 stress 0.1955484
+## Run 230 stress 0.1938366
+## Run 231 stress 0.1897209
+## ... Procrustes: rmse 0.0175776 max resid 0.1343006
+## Run 232 stress 0.1916146
+## Run 233 stress 0.1906301
+## Run 234 stress 0.1942664
+## Run 235 stress 0.1895907
+## ... Procrustes: rmse 0.005650632 max resid 0.04769659
+## Run 236 stress 0.1896277
+## ... Procrustes: rmse 0.01348837 max resid 0.1159627
+## Run 237 stress 0.1895883
+## ... Procrustes: rmse 0.01255931 max resid 0.1190452
+## Run 238 stress 0.1961186
+## Run 239 stress 0.1927524
+## Run 240 stress 0.1909499
+## Run 241 stress 0.1915205
+## Run 242 stress 0.1908374
+## Run 243 stress 0.1907584
+## Run 244 stress 0.1906184
+## Run 245 stress 0.1940027
+## Run 246 stress 0.1931748
+## Run 247 stress 0.1898224
+## ... Procrustes: rmse 0.01030259 max resid 0.0554775
+## Run 248 stress 0.1905802
+## Run 249 stress 0.1898228
+## ... Procrustes: rmse 0.01041511 max resid 0.05705223
+## Run 250 stress 0.1947886
+## Run 251 stress 0.1928064
+## Run 252 stress 0.1907908
+## Run 253 stress 0.1898029
+## ... Procrustes: rmse 0.009338604 max resid 0.05185346
+## Run 254 stress 0.1897251
+## ... Procrustes: rmse 0.01658959 max resid 0.1336738
+## Run 255 stress 0.194801
+## Run 256 stress 0.1938735
+## Run 257 stress 0.1951146
+## Run 258 stress 0.1916244
+## Run 259 stress 0.1919925
+## Run 260 stress 0.1918572
+## Run 261 stress 0.191754
+## Run 262 stress 0.1934006
+## Run 263 stress 0.1936579
+## Run 264 stress 0.1936423
+## Run 265 stress 0.1899088
+## ... Procrustes: rmse 0.01574558 max resid 0.1137615
+## Run 266 stress 0.1934007
+## Run 267 stress 0.193876
+## Run 268 stress 0.1934597
+## Run 269 stress 0.1905544
+## Run 270 stress 0.1908205
+## Run 271 stress 0.195844
+## Run 272 stress 0.1914707
+## Run 273 stress 0.1944296
+## Run 274 stress 0.1985556
+## Run 275 stress 0.1917966
+## Run 276 stress 0.1908282
+## Run 277 stress 0.1965629
+## Run 278 stress 0.1896142
+## ... Procrustes: rmse 0.01339756 max resid 0.1192769
+## Run 279 stress 0.1898247
+## ... Procrustes: rmse 0.01051423 max resid 0.05952224
+## Run 280 stress 0.1910356
+## Run 281 stress 0.1930772
+## Run 282 stress 0.1934633
+## Run 283 stress 0.1931165
+## Run 284 stress 0.1962593
+## Run 285 stress 0.1965061
+## Run 286 stress 0.1974134
+## Run 287 stress 0.1959307
+## Run 288 stress 0.1917935
+## Run 289 stress 0.1903922
+## Run 290 stress 0.1923355
+## Run 291 stress 0.1960622
+## Run 292 stress 0.1938365
+## Run 293 stress 0.193181
+## Run 294 stress 0.1916963
+## Run 295 stress 0.1908124
+## Run 296 stress 0.1896154
+## ... Procrustes: rmse 0.01354413 max resid 0.1200708
+## Run 297 stress 0.1895905
+## ... Procrustes: rmse 0.01284672 max resid 0.1215462
+## Run 298 stress 0.1915522
+## Run 299 stress 0.1955251
+## Run 300 stress 0.1914669
+## Run 301 stress 0.1975519
+## Run 302 stress 0.1896406
+## ... Procrustes: rmse 0.01442246 max resid 0.1196683
+## Run 303 stress 0.189838
+## ... Procrustes: rmse 0.01156624 max resid 0.07248831
+## Run 304 stress 0.1908202
+## Run 305 stress 0.1917537
+## Run 306 stress 0.1897993
+## ... Procrustes: rmse 0.009007352 max resid 0.05184579
+## Run 307 stress 0.1922244
+## Run 308 stress 0.1909323
+## Run 309 stress 0.1909309
+## Run 310 stress 0.1920236
+## Run 311 stress 0.191624
+## Run 312 stress 0.1915207
+## Run 313 stress 0.1895661
+## ... New best solution
+## ... Procrustes: rmse 0.001269865 max resid 0.009959143
+## ... Similar to previous best
+## *** Best solution repeated 1 times
+
+## Square root transformation
+## Wisconsin double standardization
+## Run 0 stress 0.1948295
+## Run 1 stress 0.1972171
+## Run 2 stress 0.1956095
+## Run 3 stress 0.1943622
+## ... New best solution
+## ... Procrustes: rmse 0.02713387 max resid 0.1871345
+## Run 4 stress 0.1968831
+## Run 5 stress 0.195176
+## Run 6 stress 0.1964519
+## Run 7 stress 0.1983735
+## Run 8 stress 0.1965895
+## Run 9 stress 0.1949891
+## Run 10 stress 0.197387
+## Run 11 stress 0.1955686
+## Run 12 stress 0.1965345
+## Run 13 stress 0.1973584
+## Run 14 stress 0.1959554
+## Run 15 stress 0.1955709
+## Run 16 stress 0.1968558
+## Run 17 stress 0.1943747
+## ... Procrustes: rmse 0.01858822 max resid 0.1573974
+## Run 18 stress 0.1963404
+## Run 19 stress 0.1957003
+## Run 20 stress 0.1957245
+## Run 21 stress 0.1940382
+## ... New best solution
+## ... Procrustes: rmse 0.01970152 max resid 0.1768252
+## Run 22 stress 0.1971585
+## Run 23 stress 0.1942057
+## ... Procrustes: rmse 0.01719366 max resid 0.1564763
+## Run 24 stress 0.1969554
+## Run 25 stress 0.196278
+## Run 26 stress 0.1945957
+## Run 27 stress 0.1962854
+## Run 28 stress 0.1963868
+## Run 29 stress 0.1948691
+## Run 30 stress 0.1959674
+## Run 31 stress 0.1957724
+## Run 32 stress 0.1972058
+## Run 33 stress 0.1945045
+## ... Procrustes: rmse 0.01755659 max resid 0.1205744
+## Run 34 stress 0.1950518
+## Run 35 stress 0.1949894
+## Run 36 stress 0.1961915
+## Run 37 stress 0.195606
+## Run 38 stress 0.197221
+## Run 39 stress 0.1984455
+## Run 40 stress 0.1967174
+## Run 41 stress 0.1961556
+## Run 42 stress 0.1959828
+## Run 43 stress 0.1964993
+## Run 44 stress 0.1944356
+## ... Procrustes: rmse 0.02444921 max resid 0.1833991
+## Run 45 stress 0.1981101
+## Run 46 stress 0.1944676
+## ... Procrustes: rmse 0.02185838 max resid 0.1669605
+## Run 47 stress 0.1954688
+## Run 48 stress 0.197367
+## Run 49 stress 0.1977092
+## Run 50 stress 0.1969778
+## Run 51 stress 0.1974244
+## Run 52 stress 0.1954648
+## Run 53 stress 0.1943538
+## ... Procrustes: rmse 0.01797617 max resid 0.1598972
+## Run 54 stress 0.1958871
+## Run 55 stress 0.1952714
+## Run 56 stress 0.1967838
+## Run 57 stress 0.196318
+## Run 58 stress 0.1957168
+## Run 59 stress 0.1947616
+## Run 60 stress 0.1971381
+## Run 61 stress 0.1972116
+## Run 62 stress 0.1944656
+## ... Procrustes: rmse 0.04875559 max resid 0.2129848
+## Run 63 stress 0.1964095
+## Run 64 stress 0.1967164
+## Run 65 stress 0.1972541
+## Run 66 stress 0.198516
+## Run 67 stress 0.1954722
+## Run 68 stress 0.1994759
+## Run 69 stress 0.1963823
+## Run 70 stress 0.1983984
+## Run 71 stress 0.1957931
+## Run 72 stress 0.1953022
+## Run 73 stress 0.1980429
+## Run 74 stress 0.1969466
+## Run 75 stress 0.1948577
+## Run 76 stress 0.1966478
+## Run 77 stress 0.195981
+## Run 78 stress 0.1953723
+## Run 79 stress 0.1976135
+## Run 80 stress 0.1951707
+## Run 81 stress 0.1963438
+## Run 82 stress 0.1967419
+## Run 83 stress 0.1968166
+## Run 84 stress 0.1961797
+## Run 85 stress 0.1965432
+## Run 86 stress 0.1956322
+## Run 87 stress 0.1966701
+## Run 88 stress 0.1952317
+## Run 89 stress 0.1966245
+## Run 90 stress 0.200022
+## Run 91 stress 0.1958387
+## Run 92 stress 0.1978159
+## Run 93 stress 0.1964852
+## Run 94 stress 0.1959969
+## Run 95 stress 0.1973299
+## Run 96 stress 0.1972319
+## Run 97 stress 0.1975279
+## Run 98 stress 0.198421
+## Run 99 stress 0.1960575
+## Run 100 stress 0.1963403
+## Run 101 stress 0.1998174
+## Run 102 stress 0.1968793
+## Run 103 stress 0.1986071
+## Run 104 stress 0.1950377
+## Run 105 stress 0.1970601
+## Run 106 stress 0.196615
+## Run 107 stress 0.1955119
+## Run 108 stress 0.1969225
+## Run 109 stress 0.1968677
+## Run 110 stress 0.19411
+## ... Procrustes: rmse 0.01532347 max resid 0.1582207
+## Run 111 stress 0.1964922
+## Run 112 stress 0.1939686
+## ... New best solution
+## ... Procrustes: rmse 0.02016994 max resid 0.1771657
+## Run 113 stress 0.1949802
+## Run 114 stress 0.1959918
+## Run 115 stress 0.1987794
+## Run 116 stress 0.1980934
+## Run 117 stress 0.2022944
+## Run 118 stress 0.1946934
+## Run 119 stress 0.1971622
+## Run 120 stress 0.1976449
+## Run 121 stress 0.1949804
+## Run 122 stress 0.1968289
+## Run 123 stress 0.1968501
+## Run 124 stress 0.198245
+## Run 125 stress 0.1965472
+## Run 126 stress 0.1946124
+## Run 127 stress 0.1986995
+## Run 128 stress 0.1951192
+## Run 129 stress 0.1980334
+## Run 130 stress 0.1962358
+## Run 131 stress 0.1948764
+## Run 132 stress 0.2000849
+## Run 133 stress 0.1946658
+## Run 134 stress 0.1946852
+## Run 135 stress 0.1996634
+## Run 136 stress 0.1968009
+## Run 137 stress 0.1959079
+## Run 138 stress 0.1967445
+## Run 139 stress 0.1967626
+## Run 140 stress 0.1976519
+## Run 141 stress 0.1956726
+## Run 142 stress 0.1977417
+## Run 143 stress 0.1956638
+## Run 144 stress 0.1962754
+## Run 145 stress 0.1954058
+## Run 146 stress 0.1968376
+## Run 147 stress 0.1956515
+## Run 148 stress 0.1982759
+## Run 149 stress 0.1952148
+## Run 150 stress 0.1954805
+## Run 151 stress 0.1966938
+## Run 152 stress 0.1985255
+## Run 153 stress 0.1984637
+## Run 154 stress 0.1955645
+## Run 155 stress 0.1960745
+## Run 156 stress 0.1965791
+## Run 157 stress 0.1959729
+## Run 158 stress 0.1962307
+## Run 159 stress 0.1950916
+## Run 160 stress 0.1955537
+## Run 161 stress 0.195608
+## Run 162 stress 0.1976324
+## Run 163 stress 0.1959039
+## Run 164 stress 0.1968625
+## Run 165 stress 0.1979047
+## Run 166 stress 0.1948864
+## Run 167 stress 0.1969696
+## Run 168 stress 0.1969508
+## Run 169 stress 0.1963646
+## Run 170 stress 0.1986712
+## Run 171 stress 0.1955385
+## Run 172 stress 0.1957761
+## Run 173 stress 0.1963647
+## Run 174 stress 0.1965521
+## Run 175 stress 0.1971403
+## Run 176 stress 0.1945159
+## Run 177 stress 0.1976395
+## Run 178 stress 0.194625
+## Run 179 stress 0.1960141
+## Run 180 stress 0.1957703
+## Run 181 stress 0.1991184
+## Run 182 stress 0.197228
+## Run 183 stress 0.1959958
+## Run 184 stress 0.1970487
+## Run 185 stress 0.1975668
+## Run 186 stress 0.1958509
+## Run 187 stress 0.1959249
+## Run 188 stress 0.1990631
+## Run 189 stress 0.1967124
+## Run 190 stress 0.1949879
+## Run 191 stress 0.1963849
+## Run 192 stress 0.1967959
+## Run 193 stress 0.1962111
+## Run 194 stress 0.1940779
+## ... Procrustes: rmse 0.005855018 max resid 0.0680215
+## Run 195 stress 0.1970503
+## Run 196 stress 0.1946733
+## Run 197 stress 0.1946977
+## Run 198 stress 0.1957057
+## Run 199 stress 0.1950924
+## Run 200 stress 0.1961848
+## Run 201 stress 0.1980782
+## Run 202 stress 0.1952524
+## Run 203 stress 0.1963078
+## Run 204 stress 0.1974424
+## Run 205 stress 0.1959737
+## Run 206 stress 0.1961854
+## Run 207 stress 0.1981339
+## Run 208 stress 0.1950794
+## Run 209 stress 0.194148
+## ... Procrustes: rmse 0.01753445 max resid 0.1623211
+## Run 210 stress 0.1950373
+## Run 211 stress 0.1963604
+## Run 212 stress 0.1944121
+## ... Procrustes: rmse 0.02593139 max resid 0.1689554
+## Run 213 stress 0.1967486
+## Run 214 stress 0.1974959
+## Run 215 stress 0.196877
+## Run 216 stress 0.1982982
+## Run 217 stress 0.1950337
+## Run 218 stress 0.1985535
+## Run 219 stress 0.2023546
+## Run 220 stress 0.1955064
+## Run 221 stress 0.1977017
+## Run 222 stress 0.1977683
+## Run 223 stress 0.1951582
+## Run 224 stress 0.195171
+## Run 225 stress 0.1952659
+## Run 226 stress 0.197176
+## Run 227 stress 0.1963156
+## Run 228 stress 0.1956006
+## Run 229 stress 0.1965089
+## Run 230 stress 0.196408
+## Run 231 stress 0.1971956
+## Run 232 stress 0.1979329
+## Run 233 stress 0.1977415
+## Run 234 stress 0.1949302
+## Run 235 stress 0.1968043
+## Run 236 stress 0.1959782
+## Run 237 stress 0.1961475
+## Run 238 stress 0.1959541
+## Run 239 stress 0.1946896
+## Run 240 stress 0.1953088
+## Run 241 stress 0.1971324
+## Run 242 stress 0.1969595
+## Run 243 stress 0.1965559
+## Run 244 stress 0.1964831
+## Run 245 stress 0.1981458
+## Run 246 stress 0.1964857
+## Run 247 stress 0.1961299
+## Run 248 stress 0.1975404
+## Run 249 stress 0.1955814
+## Run 250 stress 0.1958579
+## Run 251 stress 0.1954394
+## Run 252 stress 0.1956357
+## Run 253 stress 0.1963502
+## Run 254 stress 0.1948209
+## Run 255 stress 0.1949208
+## Run 256 stress 0.1948465
+## Run 257 stress 0.1979883
+## Run 258 stress 0.1941473
+## ... Procrustes: rmse 0.01811321 max resid 0.1592738
+## Run 259 stress 0.1970578
+## Run 260 stress 0.195633
+## Run 261 stress 0.1963655
+## Run 262 stress 0.1954747
+## Run 263 stress 0.1979997
+## Run 264 stress 0.1950198
+## Run 265 stress 0.1963881
+## Run 266 stress 0.1951168
+## Run 267 stress 0.1955054
+## Run 268 stress 0.1949236
+## Run 269 stress 0.1965371
+## Run 270 stress 0.1955364
+## Run 271 stress 0.1973774
+## Run 272 stress 0.1951195
+## Run 273 stress 0.1971001
+## Run 274 stress 0.195303
+## Run 275 stress 0.195889
+## Run 276 stress 0.196243
+## Run 277 stress 0.1978837
+## Run 278 stress 0.1946756
+## Run 279 stress 0.1953685
+## Run 280 stress 0.1975812
+## Run 281 stress 0.1983464
+## Run 282 stress 0.196114
+## Run 283 stress 0.1969137
+## Run 284 stress 0.1968935
+## Run 285 stress 0.1965913
+## Run 286 stress 0.1952428
+## Run 287 stress 0.1949229
+## Run 288 stress 0.1975762
+## Run 289 stress 0.2027064
+## Run 290 stress 0.1988609
+## Run 291 stress 0.1967732
+## Run 292 stress 0.1961371
+## Run 293 stress 0.1941312
+## ... Procrustes: rmse 0.02578706 max resid 0.1836086
+## Run 294 stress 0.1967686
+## Run 295 stress 0.1945943
+## Run 296 stress 0.1960046
+## Run 297 stress 0.1992348
+## Run 298 stress 0.1963737
+## Run 299 stress 0.1970285
+## Run 300 stress 0.1963328
+## Run 301 stress 0.1959941
+## Run 302 stress 0.1955817
+## Run 303 stress 0.1964946
+## Run 304 stress 0.1965427
+## Run 305 stress 0.196068
+## Run 306 stress 0.1964755
+## Run 307 stress 0.1957655
+## Run 308 stress 0.195069
+## Run 309 stress 0.1963647
+## Run 310 stress 0.1962629
+## Run 311 stress 0.1956904
+## Run 312 stress 0.1949249
+## Run 313 stress 0.1977817
+## Run 314 stress 0.1954493
+## Run 315 stress 0.1996754
+## Run 316 stress 0.1981487
+## Run 317 stress 0.1969629
+## Run 318 stress 0.1948469
+## Run 319 stress 0.1955264
+## Run 320 stress 0.1970171
+## Run 321 stress 0.1980643
+## Run 322 stress 0.1948531
+## Run 323 stress 0.1959385
+## Run 324 stress 0.1969819
+## Run 325 stress 0.1955502
+## Run 326 stress 0.1989472
+## Run 327 stress 0.1969258
+## Run 328 stress 0.1975152
+## Run 329 stress 0.1961503
+## Run 330 stress 0.1962048
+## Run 331 stress 0.1980131
+## Run 332 stress 0.1954045
+## Run 333 stress 0.1970949
+## Run 334 stress 0.1962414
+## Run 335 stress 0.197135
+## Run 336 stress 0.1960453
+## Run 337 stress 0.198219
+## Run 338 stress 0.1954481
+## Run 339 stress 0.1951492
+## Run 340 stress 0.1964684
+## Run 341 stress 0.1980847
+## Run 342 stress 0.1959226
+## Run 343 stress 0.1978133
+## Run 344 stress 0.1964054
+## Run 345 stress 0.1966999
+## Run 346 stress 0.1964876
+## Run 347 stress 0.1963074
+## Run 348 stress 0.1973617
+## Run 349 stress 0.1961349
+## Run 350 stress 0.1993755
+## Run 351 stress 0.1949832
+## Run 352 stress 0.1973227
+## Run 353 stress 0.1967578
+## Run 354 stress 0.1950124
+## Run 355 stress 0.1982633
+## Run 356 stress 0.1961275
+## Run 357 stress 0.1961174
+## Run 358 stress 0.1962278
+## Run 359 stress 0.1952574
+## Run 360 stress 0.1982571
+## Run 361 stress 0.1964328
+## Run 362 stress 0.19679
+## Run 363 stress 0.196467
+## Run 364 stress 0.1974428
+## Run 365 stress 0.1953128
+## Run 366 stress 0.1967444
+## Run 367 stress 0.1947836
+## Run 368 stress 0.1954273
+## Run 369 stress 0.1953529
+## Run 370 stress 0.1962055
+## Run 371 stress 0.1964416
+## Run 372 stress 0.1974131
+## Run 373 stress 0.1959594
+## Run 374 stress 0.2022894
+## Run 375 stress 0.1998593
+## Run 376 stress 0.1950022
+## Run 377 stress 0.1982089
+## Run 378 stress 0.197859
+## Run 379 stress 0.1953776
+## Run 380 stress 0.1962358
+## Run 381 stress 0.1956814
+## Run 382 stress 0.193887
+## ... New best solution
+## ... Procrustes: rmse 0.006044063 max resid 0.05732527
+## Run 383 stress 0.1955629
+## Run 384 stress 0.1971861
+## Run 385 stress 0.1972543
+## Run 386 stress 0.1951017
+## Run 387 stress 0.1948975
+## Run 388 stress 0.1970578
+## Run 389 stress 0.1946798
+## Run 390 stress 0.1944855
+## Run 391 stress 0.1980281
+## Run 392 stress 0.1950626
+## Run 393 stress 0.1966963
+## Run 394 stress 0.1959456
+## Run 395 stress 0.1971511
+## Run 396 stress 0.1963851
+## Run 397 stress 0.1964475
+## Run 398 stress 0.1966047
+## Run 399 stress 0.1955375
+## Run 400 stress 0.1966444
+## Run 401 stress 0.1962763
+## Run 402 stress 0.1960914
+## Run 403 stress 0.1974672
+## Run 404 stress 0.1957151
+## Run 405 stress 0.1974397
+## Run 406 stress 0.197232
+## Run 407 stress 0.1980386
+## Run 408 stress 0.1960084
+## Run 409 stress 0.1957376
+## Run 410 stress 0.1961153
+## Run 411 stress 0.2028988
+## Run 412 stress 0.196619
+## Run 413 stress 0.1956764
+## Run 414 stress 0.1970012
+## Run 415 stress 0.1954659
+## Run 416 stress 0.1957025
+## Run 417 stress 0.1960459
+## Run 418 stress 0.1966594
+## Run 419 stress 0.1979222
+## Run 420 stress 0.197876
+## Run 421 stress 0.1963762
+## Run 422 stress 0.1971267
+## Run 423 stress 0.195163
+## Run 424 stress 0.19528
+## Run 425 stress 0.1952496
+## Run 426 stress 0.202851
+## Run 427 stress 0.1965226
+## Run 428 stress 0.1955928
+## Run 429 stress 0.1969169
+## Run 430 stress 0.1965022
+## Run 431 stress 0.1980185
+## Run 432 stress 0.1979916
+## Run 433 stress 0.1939413
+## ... Procrustes: rmse 0.00553856 max resid 0.0468134
+## Run 434 stress 0.1983388
+## Run 435 stress 0.1950978
+## Run 436 stress 0.1981827
+## Run 437 stress 0.1959512
+## Run 438 stress 0.1965433
+## Run 439 stress 0.1959386
+## Run 440 stress 0.1985888
+## Run 441 stress 0.1968039
+## Run 442 stress 0.1994619
+## Run 443 stress 0.195941
+## Run 444 stress 0.1961696
+## Run 445 stress 0.1966057
+## Run 446 stress 0.1966775
+## Run 447 stress 0.1960106
+## Run 448 stress 0.1966233
+## Run 449 stress 0.1968561
+## Run 450 stress 0.1954134
+## Run 451 stress 0.1963467
+## Run 452 stress 0.1959199
+## Run 453 stress 0.1973048
+## Run 454 stress 0.1982941
+## Run 455 stress 0.1987085
+## Run 456 stress 0.195694
+## Run 457 stress 0.1961945
+## Run 458 stress 0.1984781
+## Run 459 stress 0.1955714
+## Run 460 stress 0.1972831
+## Run 461 stress 0.1968898
+## Run 462 stress 0.196071
+## Run 463 stress 0.1965107
+## Run 464 stress 0.1948227
+## Run 465 stress 0.1966497
+## Run 466 stress 0.1957486
+## Run 467 stress 0.195559
+## Run 468 stress 0.1981373
+## Run 469 stress 0.1959647
+## Run 470 stress 0.194187
+## ... Procrustes: rmse 0.01854373 max resid 0.1560759
+## Run 471 stress 0.196129
+## Run 472 stress 0.1959289
+## Run 473 stress 0.1997089
+## Run 474 stress 0.1949986
+## Run 475 stress 0.1951098
+## Run 476 stress 0.1964261
+## Run 477 stress 0.1940928
+## ... Procrustes: rmse 0.02219909 max resid 0.1712116
+## Run 478 stress 0.1947124
+## Run 479 stress 0.1968976
+## Run 480 stress 0.1967509
+## Run 481 stress 0.1961956
+## Run 482 stress 0.1966074
+## Run 483 stress 0.1950424
+## Run 484 stress 0.1962915
+## Run 485 stress 0.1971847
+## Run 486 stress 0.1971154
+## Run 487 stress 0.1956938
+## Run 488 stress 0.1980714
+## Run 489 stress 0.196942
+## Run 490 stress 0.1982367
+## Run 491 stress 0.1958801
+## Run 492 stress 0.1959479
+## Run 493 stress 0.1956743
+## Run 494 stress 0.1955249
+## Run 495 stress 0.195383
+## Run 496 stress 0.1969338
+## Run 497 stress 0.1957301
+## Run 498 stress 0.1973942
+## Run 499 stress 0.1947377
+## Run 500 stress 0.1971661
+## *** Best solution was not repeated -- monoMDS stopping criteria:
+## 331: no. of iterations >= maxit
+## 169: stress ratio > sratmax
+
+## Square root transformation
+## Wisconsin double standardization
+## Run 0 stress 0.1812837
+## Run 1 stress 0.1792325
+## ... New best solution
+## ... Procrustes: rmse 0.09462102 max resid 0.4475736
+## Run 2 stress 0.1779978
+## ... New best solution
+## ... Procrustes: rmse 0.07691572 max resid 0.2447974
+## Run 3 stress 0.1786597
+## Run 4 stress 0.1859918
+## Run 5 stress 0.1799088
+## Run 6 stress 0.1747167
+## ... New best solution
+## ... Procrustes: rmse 0.04975444 max resid 0.2121499
+## Run 7 stress 0.1749613
+## ... Procrustes: rmse 0.01163383 max resid 0.07486832
+## Run 8 stress 0.184274
+## Run 9 stress 0.1752255
+## Run 10 stress 0.1758472
+## Run 11 stress 0.174718
+## ... Procrustes: rmse 0.0007670905 max resid 0.003794574
+## ... Similar to previous best
+## Run 12 stress 0.1758696
+## Run 13 stress 0.1747188
+## ... Procrustes: rmse 0.0008787697 max resid 0.003974277
+## ... Similar to previous best
+## Run 14 stress 0.1766647
+## Run 15 stress 0.1807411
+## Run 16 stress 0.1758473
+## Run 17 stress 0.1775005
+## Run 18 stress 0.1849466
+## Run 19 stress 0.1759661
+## Run 20 stress 0.1766347
+## *** Best solution repeated 2 times
+
+