Nonmetric Multidimensional Scaling- Functional group abundance

Data

functional groups based on NOAA IEA feeding guilds
average across depth strata using the NOAA IEA procedure

#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/group_abundance_matrix.csv"))[-1]
trawl_data_arrange[is.na(trawl_data_arrange)]<-0

paged_table(head(trawl_data_arrange))

ABCDEFGHIJ0123456789

	Season <chr>	Region <int>	Year <int>	benthivore <dbl>	benthos <dbl>	piscivore <dbl>	planktivore <dbl>	undefined <dbl>	YEAR_GROUPS <chr>
1	Fall	1	2000	222.5042	51.845679	537.3266	169.7228	253.31852	2000-2004
2	Fall	1	2001	337.9533	41.111111	975.7028	348.6578	124.90123	2000-2004
3	Fall	1	2002	136.8917	3.019231	664.2088	281.9691	46.25153	2000-2004
4	Fall	1	2003	430.0938	4.500000	1057.9583	528.3438	865.65625	2000-2004
5	Fall	1	2004	534.5424	2.000000	1272.9108	1443.0556	614.95066	2000-2004
6	Fall	1	2005	372.2325	3.707602	800.3662	605.3618	2289.39474	2005-2009

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,9,10,11,12,13)]
ME_NMDS_data<-as.matrix(trawl_data_arrange[,4:8])

#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.1375323 
## Run 1 stress 0.1517229 
## Run 2 stress 0.1404538 
## Run 3 stress 0.1375321 
## ... New best solution
## ... Procrustes: rmse 0.0001374314  max resid 0.001853758 
## ... Similar to previous best
## Run 4 stress 0.1554313 
## Run 5 stress 0.1517228 
## Run 6 stress 0.1457243 
## Run 7 stress 0.137532 
## ... New best solution
## ... Procrustes: rmse 2.066821e-05  max resid 0.0001872769 
## ... Similar to previous best
## Run 8 stress 0.1676676 
## Run 9 stress 0.1404536 
## Run 10 stress 0.1375321 
## ... Procrustes: rmse 1.782204e-05  max resid 0.0001659615 
## ... Similar to previous best
## Run 11 stress 0.151723 
## Run 12 stress 0.1544327 
## Run 13 stress 0.1554312 
## Run 14 stress 0.1375323 
## ... Procrustes: rmse 0.0002670889  max resid 0.003317331 
## ... Similar to previous best
## Run 15 stress 0.1375321 
## ... Procrustes: rmse 1.485174e-05  max resid 0.0001391151 
## ... Similar to previous best
## Run 16 stress 0.1404533 
## Run 17 stress 0.1404533 
## Run 18 stress 0.1375323 
## ... Procrustes: rmse 0.0002888453  max resid 0.003483537 
## ... Similar to previous best
## Run 19 stress 0.1457244 
## Run 20 stress 0.1404533 
## *** Solution reached

#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

Plots

Region

ellipses can be all data points or 95% confidence level

#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()+
  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()+
  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_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")+
    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