---
title: "Pascoe P. P., Houghton M., Jones H. P, Weldrick C., Trebilco R & Shaw, J. D. The effect of seabird presence and seasonality on ground-active spider communities across temperate islands"
author: "Penelope P Pascoe"
date: "3/11/2022"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

Load libraries
```{r}
library(tidyverse)
library(DescTools)
library(lme4)
library(vegan)
library(nlme)
library(spaa)
library(multcomp)
library(pairwiseAdonis)
```


Read in data
```{r}
sp<-read.csv("C:/Users/pppascoe/OneDrive - University of Tasmania/Phd/Chapters/Tasmanian Spiders/For_submission/Revison 1/Pascoe_etal_2022_spider_data.csv")
sp <- sp %>% mutate_if(is.character,as.factor)
sp$round.F<-as.factor(sp$round)
```

# Basic Results

How many spiders captured? = 1592
```{r}
sum(sp$no.spiders)
```

How many days were traps deployed?
* across the five islands, five sampling rounds and 15 traps deployed per island and sampling round, trapping occurred for 3075 days
```{r}
days<-sp%>%group_by(island, round)%>%summarise(deployment.days=max(days))%>%mutate(trap.days=deployment.days*15)%>%summarise(total.trap.days=sum(trap.days))
sum(days$total.trap.days)
```

How many spiders caught per trapping day?
* 0.52 spiders caught per trapping day
```{r}
sum(sp$no.spiders) /sum(days$total.trap.days)
```

age and sex class breakup of total spiders caught
* Of the spiders caught, 1075 were adults and 490 were juveniles and 27 could not be classed by age
* Of the adults 58.6% were male and 41.4% female
```{r}
sp%>% group_by(age)%>%summarise(n=sum(no.spiders))

sp%>% group_by(age,sex)%>%filter(age=="adult")%>%summarise(n=sum(no.spiders),X.n=(sum(no.spiders)/1075)*100	)

```

How many were identified to family?
* 99.63% of adults were IDd to family
```{r}
'%ni%' <- function(x,y)!('%in%'(x,y))
sp%>%filter(family%ni%c("to.check", "unk", "unk.damaged", "unk.juv","unk_family") & age=="adult")%>%summarise(n=sum(no.spiders)) # 1071 adults were identified to family

(1071/1075)*100 # % of adults IDd to family
```


How many families were found?
* 26 families were found
```{r}
no.fam<-sp%>%filter(family%ni%c("to_check", "unk", "unk.damaged", "unk.juv","unk_family" ) )%>%group_by(family)%>%summarise(n=n())
dim(no.fam)

```

Which were the most common families?
* Lycosidae and Linyphiidae comprised 24.7% and 19.7% of the spiders caught, respectively.
```{r}
no.fam %>% dplyr::arrange(desc(n))
sp%>%group_by(family)%>%summarise(n=sum(no.spiders))%>%filter(family%in%c("Lycosidae", "Linyphiidae"))%>%mutate(per = n/1592*100 ) 
```
What spiders were the most widespread?
* These spiders were found on all 5 islands
```{r}
widespread<-sp%>%group_by(island, family)%>%summarise(n=sum(no.spiders))
widespread<-as.data.frame(table(widespread$family))
widespread%>%filter(Freq==5)

```

How widespread were the two most abundant families across the islands and sites?
* Linyphiidae were found across all five islands and all locations within an island, and Lycosidae were found at all locations except for at the colony and outside-colony locations on Wedge Island
```{r}
site<-sp%>%group_by(island, seabirds_present, family)%>%summarise(n=sum(no.spiders))
site1<-site%>%filter(family%in% c("Lycosidae","Linyphiidae"))
site2<-as.data.frame(table(site1$island,site1$family, site1$seabirds_present))
site2%>%filter(Freq==1) 
```


How correlated are habitat class and seabird colony status?
* very correlated (0.90)
```{r}
ContCoef(sp$seabirds_present, sp$v.code.s, correct = TRUE) 
```

# Results Fig 2: % Aundance by family
```{r}
sp$family<-as.factor(sp$family)


(Figure2_percent<-sp
  %>% mutate(family = recode_factor(family, `unk.juv` = "unknown immatures", 'unk.damaged' = "unknown damaged", 'unk_family' = "unknown family"))%>%
    group_by(family)%>%summarise(n.family = sum(no.spiders))%>%
    mutate(X.std.spi=(n.family/(1592))*100)%>%
    ggplot(., aes(x=reorder(family, -X.std.spi), y=X.std.spi, fill="black"))+geom_col(fill="darkgreen")+theme_classic(base_size = 25) + theme(axis.text.x = element_text( size = 15, angle = 45, hjust = 1))+guides(fill="none")+xlab("")+ylab("% of total spiders captured"))
```


# What islands had average highest spider activity density each round?
* Across the sampling rounds, spider activity density was greatest on Whale Bone Point (11.8 ± 4.64 SD spiders captured island wide per trapping day) and        lowest on Wedge Island (4.06 ± 1.65 SD spiders captured island wide per trapping day), with the differences between island statistically significant (ANOVA,   F4,20 = 4.51, P = 0.009) 

```{r}
t<-sp%>%group_by(island, round.F)%>%summarise(n.sp=sum(no.spiders), t.days = max(days))%>%mutate(trap.days = t.days*5)%>%mutate(std.sp=n.sp/t.days) 

t%>%group_by(island)%>%summarise(
  sd.mean.n=mean(std.sp),
  sd.sd=sd(std.sp))

fit_island<-lm(std.sp~island, data=t)
anova(fit_island)
summary(fit_island)
```


# Round and collony effects on activity density
```{r}
#Standardise spider captures by trap days for each colony location each round
df_1<-sp%>%group_by(island, round.F,seabirds_present)%>%summarise(n.sp=sum(no.spiders), t.days = max(days))%>%mutate(trap.days = t.days*5)%>%mutate(std.sp=n.sp/t.days)
df_1$round.F<-as.factor(df_1$round.F)


# There was no support for an effect of island and sampling round interacting  
top_F1<-lmer(std.sp~(seabirds_present*round.F)+(1|island), data=df_1, na.action=na.fail)
drop1(top_F1,test="Chisq")

# There was no support for an effect of island and sampling round at all
top_F2<-lmer(std.sp~(seabirds_present+round.F)+(1|island), data=df_1, na.action=na.fail)
drop1(top_F2,test="Chisq")

# Visual inspection of the data suggested an interaction between colony and island. On Maatsuyker Island spider activity-density was highest in the outside seabird colony location compared to either of the close to or colony locations with 5.60 ± 2.01 and 2.66 ± 1.65 spiders captured at the outside and colony locations each sampling round, respectively (TukeyHSD t=2.807, P = 0.0391). 

df_1%>%filter(island=="Maatsuyker")%>%group_by(seabirds_present)%>%
  summarise(mean.n=mean(n.sp),
            sd=sd(n.sp),
            sd.mean.n=mean(std.sp),
            sd.sd=sd(std.sp))
tt_MI<-df_1%>%filter(island=="Maatsuyker")%>%group_by(seabirds_present)
fit_col_MI<-lm(std.sp~seabirds_present,data=tt_MI)
anova(fit_col_MI)
summary(glht(fit_col_MI, linfct = mcp(seabirds_present = "Tukey")))

# Across the other islands spider activity density was greater at colony sites that at outside-colony sites, with an average of 3.16 ± 2.43 SD and 1.60 ± 1.05 SD spiders captured per trap day each sampling round in each location, respectively (TukeyHSD t=-2.537, P = 0.037).

df_1%>%filter(island!="Maatsuyker")%>%group_by(seabirds_present)%>%
  summarise(mean.n=mean(n.sp),
            sd=sd(n.sp),
            sd.mean.n=mean(std.sp),
            sd.sd=sd(std.sp))
tt<-df_1%>%filter(island!="Maatsuyker")%>%group_by(seabirds_present)
fit_col<-lm(std.sp~seabirds_present,data=tt)
anova(fit_col)
summary(glht(fit_col, linfct = mcp(seabirds_present = "Tukey")))


# Adding island into the model as a fix effect interacting with colony status confirmed the significance of this interaction (LRT, G8=27.12, P <0.001) 

top_mod<-lmer(std.sp~(seabirds_present*island)+(1|island), data=df_1, na.action=na.fail)
drop1(top_mod,test="Chisq")
summary(top_mod)


```


#Plot the modeled output - Figure 3: predicted spider activity density for each island and colony status
```{r}
predicted_values<-df_1%>%dplyr::select("std.sp", "seabirds_present", "island")%>%
  modelr::add_predictions(top_mod)

predicted_values%>%mutate(island_reorder = fct_relevel(island, "Cape Queen Elizabeth", "Whale Bone Point", "Courts", "Wedge", "Maatsuyker")) %>% filter(round.F==1)%>%group_by(island)%>%summarise(n=sum(pred))
top_island<-lm(std.sp~island, data=df_1)
summary(top_island)


(Fig3<-predicted_values%>%mutate(island_reorder = fct_relevel(island, "Cape Queen Elizabeth", "Whale Bone Point", "Courts", "Wedge", "Maatsuyker"))%>%mutate(colony_reorder = fct_relevel(seabirds_present, "colony", "close","no")) %>% filter(round.F==1)%>%dplyr::select(colony_reorder, island_reorder, pred)%>%
    ggplot(., aes(x=island_reorder, y=pred, fill=colony_reorder,width=.8 ))+geom_col(col="black")+theme_classic(base_size = 27)+xlab("")+ylab("Predicted number of spiders caught per round")+ theme(axis.text.x = element_text( size = 27, angle = 45, hjust = 1))+
    scale_fill_brewer(palette = "Greens") + guides(fill=guide_legend(title="Colony Status")))
```


# Rarefied fmily-level richness

Make data frames
```{r}
## Standardise spider captures
d.col<-sp%>%
  filter(family%ni% c("unk.damaged", "unk.juv", "unk_family"))%>%
  dplyr::group_by( island, family, seabirds_present, round.F)%>%
  summarise(n = sum(no.spiders))%>%
  pivot_wider(names_from = family, values_from=n, values_fn = sum) 

d.col[is.na(d.col)] = 0
d.col<-as.data.frame(d.col)

(d.col.table<-d.col %>%
    mutate(rich = rowSums(d.col[, -1:-3]!=0))%>%rowwise()%>%mutate(abund=sum(c_across("Amaurobiidae" :"Stiphidiidae")))%>%dplyr::select(island, seabirds_present, round.F, rich, abund))
# d.col.table<-d.col.table%>%unite(ID, c("island", "seabirds_present", "round.F"))
d.col.table$ID <- paste(d.col$island, d.col$seabirds_present, d.col$round.F, sep="_")
d.col.table<-as.data.frame(d.col.table)
row.names(d.col.table) <- d.col.table$ID

rs <- d.col.table$abund
quantile(rs) # The smallest number of spiders

Srar <- rarefy(d.col[, -1:-3], min(rs)) # the expected species richness in each sample
head(Srar)
df_srar<-as.data.frame(Srar)
rar.col<-cbind(d.col.table,df_srar )%>% dplyr::select(island, seabirds_present, round.F, ID,abund, rich,Srar)

```


Model rarefied richness
```{r}
rar.col$round.F<-as.factor(rar.col$round.F)

# Rarefied family richness did not differ between sampling events (LRT, G4=4.80, P=0.309), but did differ between colony status locations however (LRT, G2=6.88, P=0.032).
F.4_rar<-lmer(Srar~(seabirds_present+round.F)+(1|island), data=rar.col, na.action = na.fail)
drop1(F.4_rar, test="Chisq")


#Richness greater at outside-colony locations (1.79 ± 0.144) than colony locations (1.68 ± 0.185) (Tukey’s, P=0.044). Near-colony locations did not differ significantly from the other two locations.
Rar.top<-lmer(Srar~(seabirds_present)+(1|island), data=rar.col, na.action = na.fail)
drop1(Rar.top, test="Chisq")
summary(Rar.top)
summary(glht(Rar.top, linfct = mcp(seabirds_present = "Tukey")))

rar.col%>%group_by(seabirds_present)%>%summarise(mean=mean(Srar), sd=sd(Srar))
```


# Age and sex structure

##Age class structure
Proportion of adults
```{r}
df_age<-sp%>% filter(age%in%c("adult", "juv"))%>%
  group_by(island, round.F,seabirds_present, age)%>%
  summarise(n.sp=sum(no.spiders), t.days = max(days))%>%mutate(trap.days = t.days*5)%>%mutate(std.sp=n.sp/t.days)

df_age$age<-factor(df_age$age, levels = c("adult", "juv", "unk", "damaged"))
df_age$age<-factor(df_age$age, levels = c("juv", "adult")) 

##calculate proportions
ag2.2<- df_age%>%dplyr::select(-c(n.sp))%>%pivot_wider(names_from = "age", values_from="std.sp")
ag2.2[is.na(ag2.2)] <- 0
ag2.3<-ag2.2%>%mutate(total.n=adult+juv)
ag2.3<-ag2.3%>%mutate("prop.ad"=adult/total.n*100)
ag2.3$round.F<-as.factor(ag2.3$round.F)
```


modeling age class structure against sampling event and colony status
```{r}
#Full model - no interaction effect
age.model.1<-lmer(prop.ad ~ (round.F*seabirds_present)+(1|island), data=ag2.3, na.action=na.fail)
drop1(age.model.1, test="Chisq")

# There is an effect of round and colony status on age class structure
age.model.2<-lmer(prop.ad ~ (round.F+seabirds_present)+(1|island), data=ag2.3, na.action=na.fail)
drop1(age.model.2, test="Chisq")

#Age class structure varied with sampling event (LRT, G4 = 36.48, P < 0.001) and to a lesser extent colony status (LRT, G2 = 4.96, P = 0.08) 
drop1(age.model.2, test="Chisq")

#A greater proportion of spiders captured at outside-colony locations were adults than in colony locations (Figure 4), although the difference was not statistically significant (Tukey’s, z = 2.06, P = 0.098).
summary(glht(age.model.2, linfct = mcp(seabirds_present = "Tukey")))
summary(glht(age.model.2, linfct = mcp(round.F = "Tukey")))
```

##Sex structure (adults only)
Proportion of adult males
```{r}
(sex.std<-sp%>% filter(sex%in%c("M", "F"), age=="adult") %>% group_by(island, round.F, seabirds_present,days, sex)%>%summarise(n=sum(no.spiders))%>%mutate(n.std=n/(days*5)))
t.sp1<-sp%>% filter(sex%in%c("M", "F"), age=="adult") %>% group_by(island, round.F, seabirds_present,days, sex)%>%summarise(n=sum(no.spiders))

sex._wide<-t.sp1%>%pivot_wider(names_from = "sex", values_from="n")

sex._wide[which(is.na(sex._wide$F)),]$F<-0
sex._wide[which(is.na(sex._wide$M)),]$M<-0

sex._wide<-sex._wide%>%mutate(tot=F+M)
sex.std_wide<-sex._wide%>%mutate(std.n.F=F/(days*5), std.n.M=M/(days*5),std.n.tot=tot/(days*5))%>% mutate(std.prop.M= (std.n.M/std.n.tot)*100)

sex.std_wide$round.F<-as.factor(sex.std_wide$round.F)
```

Modelling sex structure against round and colony status
```{r}
#Full model - no interaction effect
Fit.sex.std1<-lmer(std.prop.M ~ (round.F*seabirds_present)+(1|island), data=sex.std_wide, na.action=na.fail)
drop1(Fit.sex.std1, test="Chisq")

#The proportion of adult males to females was most affected by sampling event (LRT, G4 = 52.084, P < 0.001)
Fit.sex.std2<-lmer(std.prop.M ~ (round.F+seabirds_present)+(1|island), data=sex.std_wide, na.action=na.fail)
drop1(Fit.sex.std2, test="Chisq")
 
```


# Figure 4 - age and sex structure for the variables that the models found were significant for each
```{r}
# Make labels
roundlabs <- c("1. pre-fledging","2. start vacant colony","3. end vacant colony ","4. incubation", "5. pre-fledging")
roundlabs_season <- c("1. Autumn","2. Winter","3. Spring","4. Summer", "5. Autumn")
collabs <- c("outside", "near", "colony")

(Fig4_age_round<-ggplot(df_age, aes(x=round.F, y= std.sp, fill = age, width=.9))+geom_col(position = "fill")+scale_y_continuous(labels = scales::percent)+scale_fill_brewer(palette = "Greens", direction = 1, labels = c("Juvenile", "Adult"), name="Age Class")+theme_classic(base_size = 20)+xlab("Sampling Event")+ylab("% of each age class caught")+ theme(axis.text.x = element_text( size = 15, angle = 45, hjust = 1))+ scale_x_discrete(labels= roundlabs_season))


(Fig4_age_colony<-ggplot(df_age%>%mutate(colony_reorder = fct_relevel(seabirds_present, "no","close","colony"))%>%group_by(colony_reorder), aes(x=colony_reorder, y= std.sp, fill = age, width=.9))+geom_col(position = "fill")+scale_y_continuous(labels = scales::percent)+scale_fill_brewer(palette = "Greens", direction = 1, labels = c("Juvenile", "Adult"), name="Age Class")+theme_classic(base_size = 20)+xlab("Colony Status")+ylab("% of each age class caught")+ scale_x_discrete(labels= collabs)) 

(Fig4_sex_round<-ggplot(sex.std, aes(x=round.F, y= n.std, fill = sex, width=.9))+geom_col(position = "fill")+scale_y_continuous(labels = scales::percent)+scale_fill_brewer(palette = "Blues", direction = 1, name="Sex", labels=c("Female","Male"))+theme_classic(base_size = 20)+xlab("Sampling Event")+ylab("% of each sex")+theme(axis.text.x = element_text( size = 15, angle = 45, hjust = 1))+ scale_x_discrete(labels= roundlabs_season))

```


Which spiders made the winter sampling event different?
```{r}
# Standardise spider captures
sp1<-sp%>%mutate(n.std=no.spiders/days)

# in round2, which family dominates, and what % of the catch does it make up?
tt<-sp1%>%filter(age=="adult")%>%group_by(family, round)%>%summarise(n.family = sum(n.std))%>%pivot_wider(names_from = family, values_from=n.family)
tt[is.na(tt)] = 0
tt <- tt %>% rowwise(round)
tt %>% mutate(TOTAL = sum(c_across(Amaurobiidae:Zodariidae)))%>%pivot_longer(cols = Amaurobiidae:Zodariidae, names_to = "family",values_to = "proportion")%>%mutate(X.fam=(proportion/TOTAL)*100)%>%filter(round==2,X.fam==max(X.fam))

# 42.5% of adult spiders in round 2 are Linyphiiidaes. What proportion of them are Male
tt %>% mutate(TOTAL = sum(c_across(Amaurobiidae:Zodariidae)))%>%pivot_longer(cols = Amaurobiidae:Zodariidae, names_to = "family",values_to = "proportion")%>%mutate(X.fam=(proportion/TOTAL)*100)%>%filter(family=="Linyphiidae")

# 80.5% of round 2 adult linyphiidae are M
t<-sp1%>%filter(age=="adult")%>%group_by(family, round,sex)%>%summarise(n.family = sum(n.std))%>%pivot_wider(names_from = sex, values_from=n.family)%>%dplyr::select(female='F', male='M')
t[is.na(t)] = 0
t%>%mutate(total=sum(female + male))%>%mutate(X.M = (male/total)*100)%>%dplyr::filter(round==2, family=="Linyphiidae")
```

# Community Composition

Make a species matrix
```{r}
sp1<-sp%>%mutate(n.std=no.spiders/days)
df.spider.std<- sp1%>%
  filter(family%ni% c("unk.damaged", "unk.juv"))%>%
  mutate(loc.id = paste(island,seabirds_present, round.F,v.code.s,sep="_"))%>%
  dplyr::group_by(loc.id, island, seabirds_present, round.F,family, days, veg=v.code.s) %>%
  summarise(n.std = sum(n.std))%>%
  pivot_wider(names_from = family, values_from=n.std,values_fn = sum)


df.spider.std[is.na(df.spider.std)] = 0
mx.spider.std.df<-as.data.frame(df.spider.std)
row.names(mx.spider.std.df) <- mx.spider.std.df$loc.id
mx.spider.std<-mx.spider.std.df%>%dplyr::select(c(-loc.id, -island, -seabirds_present, -round.F, -days, -veg))
```

## PCA (unconstrained redundancy analysis)
```{r}
ord <- rda(decostand(mx.spider.std, method="hellinger"), retx=FALSE, method="euclidean", scale=FALSE) 
plot(ord)
scores(ord)
biplot(ord, scaling = "symmetric") 
```


## Environmental data
```{r}
env2<-sp%>%filter(family%ni% c("unk.damaged", "unk.juv", "unk_family"))%>%dplyr::select(island, locality, round,seabirds_present)%>%
  mutate(loc.id=paste(island, seabirds_present,round, sep="_"))%>%
  dplyr::select(-locality, )%>%
  group_by(loc.id, island, round, colony=seabirds_present)%>%summarise()

env2$colony<-as.factor(env2$colony)
env2$island<-factor(env2$island, levels=c("Wedge", "Cape Queen Elizabeth" ,"Courts" ,"Maatsuyker" ,"Whale Bone Point"))
env2$round.F<-as.factor(env2$round)
env2<-env2[ with(env2, order(island)),]
```


## ellipsoid hulls & standard error & deviation ellipses

### By island
```{r}
#make the order of islands in env2 match the order in mx.spider.std
levels(env2$island) <- c("Cape Queen Elizabeth" ,"Courts" ,"Maatsuyker" ,"Wedge", "Whale Bone Point")

col_vec_is <- c("red", "green", "gold", "blue","purple")
col_vec_is_long<- rep(col_vec_is, 1, each=15) # in mx.spider.std, island  goes Cape Queen Elizabeth*15, Courts*15 etc

t<-biplot(ord, scaling = "species",
          display = "species",  col = "black", type=c("points"), ylim=c(-0.7, .7),xlim=c(-.9, .9),lwd = 5, cex=2)

## standard error of centroid  ellipse
ordiellipse(ord, groups = env2$island, draw = "polygon", col = col_vec_is,alpha=.4,
            scaling = "species", lwd = 2)

#add site (loc.id) points
points(ord, display = "site", scaling = "species",
       pch = 19, cex=1.5, col =alpha(col_vec_is_long, ), bg = alpha(col_vec_is_long, 1))


# #add species labels
orditorp(t,display = "site", cex=1, col="black", pcol=NA)

#Add legend
env2$island
lvl <- with(env2, levels(island))
legend("bottomright", legend = lvl,
       bty = "n", col = col_vec_is, pch = 19, cex=1)
```


### By colony
```{r}
#make the order of colony in env2 match the order in mx.spider.std
levels(env2$colony) <- c("near", "colony","outside")

col_vec_col <- c("blue", "green", "red")
col_vec_col_long<- rep(col_vec_col, 5, each=5) # in mx.spider.std, colony just goes closex5, colonyx5, nox5 (5 times)

t<-biplot(ord, scaling = "species",
          display = "species",  col = "black", type=c("points"), ylim=c(-0.7, .7),xlim=c(-.9, .9),lwd = 5, cex=2)

ordiellipse(ord, groups = env2$colony,
            draw = "polygon", col = col_vec_col,
            scaling = "species", lwd = 2)

points(ord, display = "site", scaling = "species",
       pch = 19, cex=1.5, col =alpha(col_vec_col_long, 1), bg = alpha(col_vec_col_long, 1))

#add species labels
orditorp(t,display = "species", cex=1, col="black", pcol=NA)

#Add legend
lvl <- with(env2, levels(colony))
legend("bottomright", legend = lvl,
       bty = "n", col = col_vec_col, pch = 19, cex=1)

```




### By round
```{r}
#make the order of colony in env2 match the order in mx.spider.std
env2$round.F<-as.factor(env2$round)
levels(env2$round.F) <- c("1. Autumn", "2. Winter", "3. Spring", "4. Summer", "5. Autumn")

col_vec <- c("red", "blue", "green", "yellow", "pink")
col_vec3<- rep(col_vec,15 ) # in mx.spider.std, round just goes 12345,12345,12345 (15 times)


t<-biplot(ord, scaling = "species",
          display = "species",  col = "black", type=c("points"), ylim=c(-0.7, .7),xlim=c(-0.9, .9))

## standard error of centroid  ellipse
ordiellipse(ord, groups = env2$round.F,
            draw = "polygon", col = col_vec,alpha=.4,
            scaling = "species", lwd = 2)

#add site (loc.id) points
points(ord, display = "site", scaling = "species",
       pch = 19, cex=1.5, col =alpha(col_vec3, 1), bg = alpha(col_vec3, 1))

#add species labels
orditorp(t,display = "species", cex=1, col="black", pcol="White")

#Add legend
lvl <- with(env2, levels(round.F))
legend("bottomright", legend = lvl,
       bty = "n", col = col_vec, pch = 19, cex=1)
```




# NMDS - for supplimentary materials
```{r}
(sol <- metaMDS(mx.spider.std[, -1], trymax = 1000, trace=TRUE))
plot(sol)
stressplot(sol)
```



## plotting NMDS

### Island
```{r}
levels(env2$island) <- c("Cape Queen Elizabeth" ,"Courts" ,"Maatsuyker" ,"Wedge", "Whale Bone Point")

col_vec_is <- c("red", "green", "gold", "blue","purple")
col_vec_is_long<- rep(col_vec_is, 1, each=15) # in mx.spider.std, island  goes wedge*15, CQE*15 etc

t<-ordiplot(sol, type = "points", display="sites")
points(sol, pch=19, cex=1.5, col=col_vec_is_long)
ordiellipse(sol,draw="polygon",groups=env2$island, col=col_vec_is,label=F)
orditorp(sol,display="species",col="black",air=0.01, cex=1)

#Add legend
env2$island
lvl <- with(env2, levels(island))
legend("bottomright", legend = lvl,
       bty = "n", col = col_vec_is, pch = 19, cex=1)
```



### colony
```{r}
levels(env2$colony) <- c("near", "colony","outside")

col_vec_col <- c("blue", "green", "red")
col_vec_col_long<- rep(col_vec_col, 5, each=5) # in mx.spider.std, round just goes 12345,12345,12345 (15 times)

ordiplot(sol, type = "points")
ordiellipse(sol,draw="polygon",groups=env2$colony, col=col_vec_col,label=F)
orditorp(sol,display="species",col="black",air=0.01, cex=1)
points(sol, pch=19, cex=1.5, col=col_vec_col_long)#,       pch = 19, cex=1.5, col =alpha(col_vec3, 1), bg = alpha(col_vec3, 1))

#Add legend
env2$colony
lvl <- with(env2, levels(colony))
legend("bottomright", legend = lvl,
       bty = "n", col = col_vec_col, pch = 19, cex=1)
```


### round
```{r}
env2$round.F<-as.factor(env2$round)
levels(env2$round.F) <- c("1. Autumn", "2. Winter", "3. Spring", "4. Summer", "5. Autumn")

col_vec <- c("red", "blue", "green", "yellow", "pink")
col_vec3<- rep(col_vec,15 ) # in mx.spider.std, round just goes 12345,12345,12345 (15 times)
disp <- "site"
scl <- "species"

ordiplot(sol, type = "points", display = "sites")
ordiellipse(sol,draw="polygon",groups=env2$round.F, col=col_vec,label=F)
orditorp(sol,display="species",col="black",air=0.01, cex=1)
points(sol, pch=19, cex=1.5, col=col_vec3)

#Add legend
env2$round.F
lvl <- with(env2, levels(round.F))
legend("bottomright", legend = lvl,
       bty = "n", col = col_vec, pch = 19, cex=1)
```


## PERMANOVA to test results of biplots
```{r}
dis_hell<- vegdist(decostand(mx.spider.std, method = "hellinger"), method = "euclidean") # this is the hellinger distance (even though it says euclidean too)

adonis2(dis_hell~round.F, data=env2, method = "hellinger", permutations = 10000, by="margin")
adonis2(dis_hell~colony, data=env2, method = "hellinger", permutations = 10000, by="margin")
adonis2(dis_hell~island, data=env2, method = "hellinger", permutations = 10000, by="margin")
```


## two-way comparisons
```{r}
(i2<-pairwise.adonis(dis_hell,factors=mx.spider.std.df$island,p.adjust.m='bonferroni', perm=10000)%>%arrange(., desc(p.adjusted)))
(c2<-pairwise.adonis(dis_hell,factors=mx.spider.std.df$seabirds_present,p.adjust.m='bonferroni',perm=10000)%>%arrange(., desc(p.adjusted)))
(r2<-pairwise.adonis(dis_hell,factors=mx.spider.std.df$round.F,p.adjust.m='bonferroni', perm=10000)%>%arrange(., desc(p.adjusted)))

pair_com<-rbind(i2, c2, r2)
```
