Nonmetric Multidimensional Scaling- Top 50 biomass species

Data

#load packages
library(tidyverse)
library(vegan)
library(ggrepel)
library(ggforce)
library(ggnewscale)
library(ggthemes)
library(here)
library(rmarkdown)

#set up ggplot theme
theme_set(theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"),
        axis.text.x = element_text(colour = "black", face = "bold", size = 12),
        legend.text = element_text(size = 12, face ="bold", colour ="black"),
        legend.position = "right", 
        axis.title.y = element_text(face = "bold", size = 14, angle=90),
        axis.title.x = element_text(face = "bold", size = 12, colour = "black"),
        legend.title = element_text(size = 14, colour = "black", face = "bold"),
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
        legend.key=element_blank()))

#load in data
trawl_data_arrange<-read.csv(here("Data/species_biomass_matrix.csv"))[-1]
paged_table(head(trawl_data_arrange))
set.seed(123)

Set up data for NMDS

  • split community matrix into two dataframes- one for grouping variables and one for species biomass
  • calculate dissimilarity matrix with Bray-Curtis distances
#set up final grouping data into dataframe
ME_group_data<-trawl_data_arrange[, c(1,2,3,55,56,57,58)]
ME_NMDS_data<-as.matrix(trawl_data_arrange[,4:53])

#calculate distance matrix
ME_NMDS_distance<- vegdist(ME_NMDS_data, method="bray")

Run the NMDS and extract scores

  • change in community composition
  • uses rank order
  • stress < 0.2 is good, < 0.1 is great, <0.05 is excellent representation in reduced dimensions
ME_NMDS=metaMDS(ME_NMDS_distance, # Our community-by-species matrix
                k=2, # The number of reduced dimensions
                method="bray",
                trymax=200) #increase iterations
## Run 0 stress 0.2121976 
## Run 1 stress 0.2057908 
## ... New best solution
## ... Procrustes: rmse 0.03260853  max resid 0.3907521 
## Run 2 stress 0.207885 
## Run 3 stress 0.2066022 
## Run 4 stress 0.2126567 
## Run 5 stress 0.2117096 
## Run 6 stress 0.2090653 
## Run 7 stress 0.2127132 
## Run 8 stress 0.2061185 
## ... Procrustes: rmse 0.008997804  max resid 0.07750829 
## Run 9 stress 0.2120618 
## Run 10 stress 0.2121856 
## Run 11 stress 0.2061184 
## ... Procrustes: rmse 0.008992806  max resid 0.07749652 
## Run 12 stress 0.2064777 
## Run 13 stress 0.2058398 
## ... Procrustes: rmse 0.003094323  max resid 0.03803567 
## Run 14 stress 0.2078695 
## Run 15 stress 0.2135671 
## Run 16 stress 0.2061183 
## ... Procrustes: rmse 0.008977382  max resid 0.07733821 
## Run 17 stress 0.2118637 
## Run 18 stress 0.2065612 
## Run 19 stress 0.2124997 
## Run 20 stress 0.2065748 
## Run 21 stress 0.2133277 
## Run 22 stress 0.2057111 
## ... New best solution
## ... Procrustes: rmse 0.00395844  max resid 0.03825867 
## Run 23 stress 0.2120847 
## Run 24 stress 0.2063852 
## Run 25 stress 0.2061811 
## ... Procrustes: rmse 0.01033487  max resid 0.08282011 
## Run 26 stress 0.2129884 
## Run 27 stress 0.20681 
## Run 28 stress 0.2139782 
## Run 29 stress 0.2057848 
## ... Procrustes: rmse 0.002541256  max resid 0.0229569 
## Run 30 stress 0.2123711 
## Run 31 stress 0.2120959 
## Run 32 stress 0.2090882 
## Run 33 stress 0.2065496 
## Run 34 stress 0.2314908 
## Run 35 stress 0.2078696 
## Run 36 stress 0.2078851 
## Run 37 stress 0.2065612 
## Run 38 stress 0.2064649 
## Run 39 stress 0.2068109 
## Run 40 stress 0.2068101 
## Run 41 stress 0.2066023 
## Run 42 stress 0.2228289 
## Run 43 stress 0.2117409 
## Run 44 stress 0.2346884 
## Run 45 stress 0.2065999 
## Run 46 stress 0.2065853 
## Run 47 stress 0.2065871 
## Run 48 stress 0.20787 
## Run 49 stress 0.4165635 
## Run 50 stress 0.206812 
## Run 51 stress 0.2063172 
## Run 52 stress 0.212956 
## Run 53 stress 0.2273767 
## Run 54 stress 0.2121424 
## Run 55 stress 0.2131374 
## Run 56 stress 0.2090879 
## Run 57 stress 0.2286814 
## Run 58 stress 0.2063667 
## Run 59 stress 0.2220723 
## Run 60 stress 0.2063955 
## Run 61 stress 0.2063136 
## Run 62 stress 0.229283 
## Run 63 stress 0.2078751 
## Run 64 stress 0.2064627 
## Run 65 stress 0.2126292 
## Run 66 stress 0.2061264 
## ... Procrustes: rmse 0.007105145  max resid 0.07700621 
## Run 67 stress 0.2120917 
## Run 68 stress 0.2129648 
## Run 69 stress 0.206487 
## Run 70 stress 0.2061524 
## ... Procrustes: rmse 0.00734878  max resid 0.07782809 
## Run 71 stress 0.2060866 
## ... Procrustes: rmse 0.00739737  max resid 0.0739193 
## Run 72 stress 0.416565 
## Run 73 stress 0.2120858 
## Run 74 stress 0.2302658 
## Run 75 stress 0.2065616 
## Run 76 stress 0.2145424 
## Run 77 stress 0.2065613 
## Run 78 stress 0.2117192 
## Run 79 stress 0.2065547 
## Run 80 stress 0.2078697 
## Run 81 stress 0.2138939 
## Run 82 stress 0.2064631 
## Run 83 stress 0.2068108 
## Run 84 stress 0.2120437 
## Run 85 stress 0.2283524 
## Run 86 stress 0.2057105 
## ... New best solution
## ... Procrustes: rmse 0.0002644395  max resid 0.00239894 
## ... Similar to previous best
## *** Solution reached
ME_NMDS
## 
## Call:
## metaMDS(comm = ME_NMDS_distance, k = 2, trymax = 200, method = "bray") 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     ME_NMDS_distance 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2057105 
## Stress type 1, weak ties
## Two convergent solutions found after 86 tries
## Scaling: centring, PC rotation 
## Species: scores missing
#extract NMDS scores for ggplot
data.scores = as.data.frame(scores(ME_NMDS))
#add columns to data frame 
data.scores$Stratum = trawl_data_arrange$Stratum
data.scores$Region = trawl_data_arrange$Region
data.scores$Year = trawl_data_arrange$Year
data.scores$Season= trawl_data_arrange$Season
data.scores$Year_groups= trawl_data_arrange$YEAR_GROUPS
data.scores$Year_decades= trawl_data_arrange$YEAR_DECADES
data.scores$Region_new=trawl_data_arrange$REGION_NEW
data.scores$Region_year=trawl_data_arrange$REGION_YEAR
data.scores$Season_year=trawl_data_arrange$SEASON_YEAR


#calculate center of ellipse for each region/year grouping
data.scores_2<-group_by(data.scores, Region_new, Year_groups, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Region_new, Year_groups,Season, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:24)
rownames(d)<-NULL
centers<-full_join(d,id)

Plots

Region

#both seasons
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")

#calculate center of ellipse for each region/year grouping
data.scores_2<-group_by(data.scores_region, Region_new, Year_groups)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Region_new, Year_groups, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

p <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2, shape=factor(Year_groups),color=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Region_new)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Region_new),level = 0.95,size=1)+
  scale_color_colorblind()+
 # geom_text_repel(data= data.scores_region, aes(x= NMDS1, y=NMDS2, label= Year))+
  labs(x = "NMDS1", colour = "Region", y = "NMDS2", shape = "Year")+
  geom_point(data=centers, aes(NMDS1, NMDS2,shape=Year_groups, color=Region_new), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, color=Region_new), size=2)


#spring
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")%>%
  filter(Season=="Spring")

#calculate center of ellipse for each region/year grouping
data.scores_2<-group_by(data.scores_region, Region_new, Year_groups, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Region_new, Year_groups,Season, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

p1 <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2, shape=factor(Year_groups),color=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Region_new)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Region_new),level = 0.95,size=1)+
  scale_color_colorblind()+
 # geom_text_repel(data= data.scores_region, aes(x= NMDS1, y=NMDS2, label= Year))+
  labs(x = "NMDS1", colour = "Region", y = "NMDS2", shape = "Year", title="Spring")+
  geom_point(data=centers, aes(NMDS1, NMDS2,shape=Year_groups, color=Region_new), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, color=Region_new), size=2)

#Fall
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")%>%
  filter(Season=="Fall")

#calculate center of ellipse for each region/year grouping
data.scores_2<-group_by(data.scores_region, Region_new, Year_groups, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Region_new, Year_groups,Season, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)


p2 <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2, shape=factor(Year_groups),color=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Region_new)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Region_new),level = 0.95,size=1)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Region", y = "NMDS2", shape = "Year",title="Fall")+
  geom_point(data=centers, aes(NMDS1, NMDS2, shape=Year_groups, color=Region_new), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, color=Region_new), size=2)


#Spring
data.scores_east<-filter(data.scores, Region_new=="East of Penobscot Bay")%>%
  filter(Season=="Spring")
data.scores_west<-filter(data.scores, Region_new=="West of Penobscot Bay")%>%
  filter(Season=="Spring")

p3<- ggplot() + 
  geom_point(data=data.scores_east, aes(x = NMDS1, y = NMDS2,color=factor(Region_year), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_east,aes(x=NMDS1,y=NMDS2,color=factor(Region_year)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_east, aes(x=NMDS1, y=NMDS2, color=Region_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("grey", "black", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_east, aes(x=NMDS1,y=NMDS2,group=factor(Region_new)))+
  #geom_text_repel(data=data.scores_east,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold", color="black")+
  ggnewscale::new_scale_color()+
  geom_point(data=data.scores_west, aes(x = NMDS1, y = NMDS2,color=factor(Region_year), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_west,aes(x=NMDS1,y=NMDS2,color=factor(Region_year)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_west, aes(x=NMDS1, y=NMDS2, color=Region_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("light blue", "navy", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_west, aes(x=NMDS1,y=NMDS2,group=factor(Region_new)))+
  #geom_text_repel(data=data.scores_west,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold", color="black")+
  labs(x = "NMDS1", colour = "Region/Year", y = "NMDS2", shape = "Region", fill="Region/Year", title="Spring")


#Fall
data.scores_east<-filter(data.scores, Region_new=="East of Penobscot Bay")%>%
  filter(Season=="Fall")
data.scores_west<-filter(data.scores, Region_new=="West of Penobscot Bay")%>%
  filter(Season=="Fall")

p4<- ggplot() + 
  geom_point(data=data.scores_east, aes(x = NMDS1, y = NMDS2,color=factor(Region_year), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_east,aes(x=NMDS1,y=NMDS2,color=factor(Region_year)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_east, aes(x=NMDS1, y=NMDS2, color=Region_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("grey", "black", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_east, aes(x=NMDS1,y=NMDS2,group=factor(Region_new)))+
  #geom_text_repel(data=data.scores_east,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold",   color="black")+
  ggnewscale::new_scale_color()+
  geom_point(data=data.scores_west, aes(x = NMDS1, y = NMDS2,color=factor(Region_year),      shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_west,aes(x=NMDS1,y=NMDS2,color=factor(Region_year)), size=1,   alpha=0.15)+
  stat_ellipse(data=data.scores_west, aes(x=NMDS1, y=NMDS2, color=Region_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("light blue", "navy", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_west, aes(x=NMDS1,y=NMDS2,group=factor(Region_new)))+
  #geom_text_repel(data=data.scores_west,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold",   color="black")+
  labs(x = "NMDS1", colour = "Region/Year", y = "NMDS2", shape = "Year", fill="Region/Year", title="Fall")

p

p1

p2

p3

p4

Time

#both seasons
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")

#centers for year plots
data.scores_2<-group_by(data.scores_region, Year_groups)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:4)
rownames(d)<-NULL
centers<-full_join(d,id)

p <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2,color=factor(Year_groups), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Year_groups)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Year_groups),level = 0.95,size=1)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Year", y = "NMDS2", shape = "Region")+
  geom_point(data=centers, aes(NMDS1, NMDS2, color=Year_groups), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2), size=2)


#Spring
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")%>%
  filter(Season=="Spring")

#centers for year plots
data.scores_2<-group_by(data.scores_region, Year_groups, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:4)
rownames(d)<-NULL
centers<-full_join(d,id)

p1 <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2,color=factor(Year_groups), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Year_groups)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Year_groups),level = 0.95,size=1)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Year", y = "NMDS2", shape = "Region", title="Spring")+
  geom_point(data=centers, aes(NMDS1, NMDS2, color=Year_groups), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2), size=2)

#Fall
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")%>%
  filter(Season=="Fall")

#centers for year plots
data.scores_2<-group_by(data.scores_region, Year_groups, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:4)
rownames(d)<-NULL
centers<-full_join(d,id)

p2 <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2,color=factor(Year_groups), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Year_groups)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Year_groups),level = 0.95,size=1)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Year", y = "NMDS2", shape = "Region", title="Fall")+
  geom_point(data=centers, aes(NMDS1, NMDS2, color=Year_groups), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2), size=2)


#Spring
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")%>%
  filter(Season=="Spring")

#centers for year plots
data.scores_2<-group_by(data.scores_region, Year_groups,Region_new, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season,Region_new, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

p3 <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2,color=factor(Year_groups), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Year_groups)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Year_groups),level = 0.95,size=1)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Year", y = "NMDS2", shape = "Region", title="Spring")+
  geom_point(data=centers, aes(NMDS1, NMDS2, color=Year_groups, shape=Region_new), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, group=Region_new), size=2)


#Fall
data.scores_region<-filter(data.scores,Region_new!="Penobscot Bay")%>%
  filter(Season=="Fall")

#centers for year plots
data.scores_2<-group_by(data.scores_region, Year_groups,Region_new, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season,Region_new, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

p4 <- ggplot() + 
  geom_point(data=data.scores_region, aes(x = NMDS1, y = NMDS2,color=factor(Year_groups), shape=factor(Region_new)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_region,aes(x=NMDS1,y=NMDS2,color=factor(Year_groups)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_region, aes(x=NMDS1, y=NMDS2, color=Year_groups),level = 0.95,size=1)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Year", y = "NMDS2", shape = "Region", title="Fall")+
  geom_point(data=centers, aes(NMDS1, NMDS2, color=Year_groups, shape=Region_new), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, group=Region_new), size=2)

p

p1

p2

p3

p4

Season

#all regions
#centers for year plots
data.scores_2<-group_by(data.scores, Year_groups, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

p <- ggplot(data.scores, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 3, aes(shape=factor(Year_groups), color=factor(Season)))+ 
  stat_ellipse(data=data.scores, aes(x=NMDS1, y=NMDS2, color=Season),level = 0.95,size=1)+
  #geom_mark_ellipse(aes(color=factor(Season)), size=1, alpha=0.15)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Season", y = "NMDS2", shape = "Year", fill="Season")+
    geom_point(data=centers, aes(NMDS1, NMDS2, color=Season, shape=Year_groups), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, color=Season), size=2)



#ellipse season and year plot

data.scores_west<-filter(data.scores, Region_new=="West of Penobscot Bay")
#centers for year plots
data.scores_2<-group_by(data.scores_west, Year_groups,Region_new, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season,Region_new, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

#West
p1 <- ggplot(data.scores_west, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 3, aes(shape=factor(Year_groups), color=factor(Season)))+ 
  stat_ellipse(data=data.scores_west, aes(x=NMDS1, y=NMDS2, color=Season),level = 0.95,size=1)+
  #geom_mark_ellipse(aes(color=factor(Season)), size=1, alpha=0.15)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Season", y = "NMDS2", shape = "Year", fill="Season", title="West")+
    geom_point(data=centers, aes(NMDS1, NMDS2, color=Season, shape=Year_groups), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, color=Season), size=2)


#East
data.scores_east<-filter(data.scores, Region_new=="East of Penobscot Bay")

#centers for year plots
data.scores_2<-group_by(data.scores_east, Year_groups,Region_new, Season)%>%
   mutate(indicator=cur_group_id())
id<-dplyr::select(data.scores_2, Year_groups,Season,Region_new, indicator)
id<-unique(id)

ctr<-NULL
Ncombo<-length(unique(data.scores_2$indicator))

for (i in 1:Ncombo) {
 combos<- filter(data.scores_2, indicator==i)%>%
   ungroup()%>%
   dplyr::select(1:2)
 
  ctr[[i]]<-MASS::cov.trob(combos)$center
  
}

d<-as.data.frame(ctr)
d<-as.data.frame(t(d))
d$indicator<-seq(1:8)
rownames(d)<-NULL
centers<-full_join(d,id)

p2 <- ggplot(data.scores_east, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 3, aes(shape=factor(Year_groups), color=factor(Season)))+ 
  stat_ellipse(data=data.scores_east, aes(x=NMDS1, y=NMDS2, color=Season),level = 0.95,size=1)+
  #geom_mark_ellipse(aes(color=factor(Season)), size=1, alpha=0.15)+
  scale_color_colorblind()+
  labs(x = "NMDS1", colour = "Season", y = "NMDS2", shape = "Year", fill="Season", title="East")+
  geom_point(data=centers, aes(NMDS1, NMDS2, color=Season, shape=Year_groups), size=10)+
  geom_path(data=centers, aes(NMDS1, NMDS2, color=Season), size=2)




#West
data.scores_spring<-filter(data.scores, Season=="Spring")%>%
  filter(Region_new=="West of Penobscot Bay")
data.scores_fall<-filter(data.scores, Season=="Fall")%>%
  filter(Region_new=="West of Penobscot Bay")

p3 <- ggplot() + 
  geom_point(data=data.scores_fall, aes(x = NMDS1, y = NMDS2,color=factor(Season_year), shape=factor(Season)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_fall,aes(x=NMDS1,y=NMDS2,color=factor(Season_year)), size=1, alpha=0.15)+
stat_ellipse(data=data.scores_fall, aes(x=NMDS1, y=NMDS2, color=Season_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("grey", "black", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_fall, aes(x=NMDS1,y=NMDS2,group=factor(Season)))+
  #geom_text_repel(data=data.scores_fall,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold", color="black")+
  ggnewscale::new_scale_color()+
  geom_point(data=data.scores_spring, aes(x = NMDS1, y = NMDS2,color=factor(Season_year), shape=factor(Season)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_spring,aes(x=NMDS1,y=NMDS2,color=factor(Season_year)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_spring, aes(x=NMDS1, y=NMDS2, color=Season_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("light blue", "navy", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_spring, aes(x=NMDS1,y=NMDS2,group=factor(Season)),color="navy")+
  #geom_text_repel(data=data.scores_spring,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold", color="black")+
  labs(x = "NMDS1", colour = "Season/Year", y = "NMDS2", shape = "Season", fill="Season/Year", title="West")


#East
data.scores_spring<-filter(data.scores, Season=="Spring")%>%
  filter(Region_new=="East of Penobscot Bay")
data.scores_fall<-filter(data.scores, Season=="Fall")%>%
  filter(Region_new=="East of Penobscot Bay")

p4 <- ggplot() + 
  geom_point(data=data.scores_fall, aes(x = NMDS1, y = NMDS2,color=factor(Season_year), shape=factor(Season)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_fall,aes(x=NMDS1,y=NMDS2,color=factor(Season_year)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_fall, aes(x=NMDS1, y=NMDS2, color=Season_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("grey", "black", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_fall, aes(x=NMDS1,y=NMDS2,group=factor(Season)))+
  #geom_text_repel(data=data.scores_fall,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold", color="black")+
  ggnewscale::new_scale_color()+
  geom_point(data=data.scores_spring, aes(x = NMDS1, y = NMDS2,color=factor(Season_year), shape=factor(Season)),size = 3)+ 
  #geom_mark_ellipse(data=data.scores_spring,aes(x=NMDS1,y=NMDS2,color=factor(Season_year)), size=1, alpha=0.15)+
  stat_ellipse(data=data.scores_spring, aes(x=NMDS1, y=NMDS2, color=Season_year),level = 0.95,size=1)+
  scale_color_manual(values=scales::seq_gradient_pal("light blue", "navy", "Lab")(seq(0,1,length.out=4)))+
  #geom_path(data=data.scores_spring, aes(x=NMDS1,y=NMDS2,group=factor(Season)),color="navy")+
  #geom_text_repel(data=data.scores_spring,aes(x=NMDS1,y=NMDS2,label=Year), size=4, fontface="bold", color="black")+
  labs(x = "NMDS1", colour = "Season/Year", y = "NMDS2", shape = "Year", fill="Season/Year", title="East")

p

p1

p2

p3

p4