1 Packages and data

Load packages and read in datasets for class exercises:

lfish <- read.csv(file =  "data/fishcoms_lizardisland.csv")

2 Exercise 1

We will continue using the “fishcoms_lizardisland.csv” dataset, but this time, we will explore differences in fish communities between the two years, 2011 and 2015.

2.1 Part I

  1. Calculate fish abundance per transect, site, and year for each species.
  2. Put the dataset in wide format.
  3. Split the metadata from the numeric data.
lizard.sum <- lfish %>%
  group_by(Site, transect, year, species) %>%
  summarize(abundance = sum(abundance))
## `summarise()` has grouped output by 'Site', 'transect', 'year'. You can
## override using the `.groups` argument.
lizard.wide <- lizard.sum %>%
  spread(key = species, value = abundance, fill = 0)

lizard.meta <- lizard.wide[1:3]

lizard.raw <- lizard.wide %>%
  ungroup() %>%
  select(-c(1:3))

2.2 Part II

  1. Run a PCoA with Bray-Curtis dissimilarity.
  2. Check eigenvalues. Add a correction if needed.
  3. Plot the PCoA.
lizard.d.bray <- vegdist(lizard.raw, method = "bray")

lizard.pcoa.bray <-  pcoa(lizard.d.bray)

lizard.pcoa.bray <-  pcoa(lizard.d.bray, correction = "cailliez")

lizard.scores.b <- as.data.frame(lizard.pcoa.bray$vectors) %>%
  bind_cols(lizard.meta)

lizard.hulls.b <- lizard.scores.b %>%
  group_by(year) %>%
  slice(chull(Axis.1, Axis.2))

lizard.plot.b <- ggplot() +
  geom_point(data = lizard.scores.b, aes(x = Axis.1, y = Axis.2, color = as.factor(year)), size = 2) +
  geom_polygon(data = lizard.hulls.b, aes(x = Axis.1, Axis.2, fill = as.factor(year)), alpha = 0.5) +
  scale_fill_fish_d(option = "Stegastes_nigricans", end = 0.8) +
  scale_color_fish_d(option = "Stegastes_nigricans", end = 0.8) +
  theme_bw() +
  xlab("PCo1 (12.5%)") +
  ylab("PCo2 (8.4%)") +
  ggtitle("PCoA") +
  guides(color = guide_legend(title = "Year"), fill = "none")
lizard.plot.b

2.3 Part III

  1. Run an nMDS with Bray-Curtis dissimilarity.
  2. Check stressplot.
  3. Plot the nMDS.
lizard.mds <- metaMDS(comm = lizard.raw, k = 2, trymax = 1000, distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2318191 
## Run 1 stress 0.2612653 
## Run 2 stress 0.2559672 
## Run 3 stress 0.2318425 
## ... Procrustes: rmse 0.003716228  max resid 0.02701714 
## Run 4 stress 0.2318197 
## ... Procrustes: rmse 0.003266044  max resid 0.0271231 
## Run 5 stress 0.2328584 
## Run 6 stress 0.2318197 
## ... Procrustes: rmse 0.003266679  max resid 0.02712859 
## Run 7 stress 0.2768707 
## Run 8 stress 0.2318197 
## ... Procrustes: rmse 0.003265132  max resid 0.02712373 
## Run 9 stress 0.2328682 
## Run 10 stress 0.2328509 
## Run 11 stress 0.2328584 
## Run 12 stress 0.2318198 
## ... Procrustes: rmse 0.003268323  max resid 0.02712856 
## Run 13 stress 0.2798722 
## Run 14 stress 0.2328583 
## Run 15 stress 0.2318197 
## ... Procrustes: rmse 0.003265576  max resid 0.02712908 
## Run 16 stress 0.2328635 
## Run 17 stress 0.2318198 
## ... Procrustes: rmse 0.003269891  max resid 0.027155 
## Run 18 stress 0.2328461 
## Run 19 stress 0.2586182 
## Run 20 stress 0.2318702 
## ... Procrustes: rmse 0.005120638  max resid 0.03843087 
## Run 21 stress 0.2326943 
## Run 22 stress 0.2328509 
## Run 23 stress 0.23187 
## ... Procrustes: rmse 0.00606659  max resid 0.03860568 
## Run 24 stress 0.232851 
## Run 25 stress 0.2318724 
## ... Procrustes: rmse 0.005265016  max resid 0.03847786 
## Run 26 stress 0.2318192 
## ... Procrustes: rmse 6.72684e-05  max resid 0.0003216106 
## ... Similar to previous best
## *** Best solution repeated 1 times
lizard.mds
## 
## Call:
## metaMDS(comm = lizard.raw, distance = "bray", k = 2, trymax = 1000) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(lizard.raw)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2318191 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 26 tries
## The best solution was from try 0 (metric scaling or null solution)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(lizard.raw))'
stressplot(lizard.mds)

lizard.scores <- as.data.frame(lizard.mds$points) %>%
  bind_cols(lizard.meta)

lizard.hulls <- lizard.scores %>%
  group_by(year) %>%
  slice(chull(MDS1, MDS2))

lizard.plot.mds <- ggplot() +
  geom_point(data = lizard.scores, aes(x = MDS1, y = MDS2, color = as.factor(year)), size = 2) +
  geom_polygon(data = lizard.hulls, aes(x = MDS1, MDS2, fill = as.factor(year)), alpha = 0.5) +
  scale_fill_fish_d(option = "Stegastes_nigricans", end = 0.8) +
  scale_color_fish_d(option = "Stegastes_nigricans", end = 0.8) +
  theme_bw() +
  xlab("MDS1") +
  ylab("MDS2") +
  ggtitle("nMDS") +
  guides(color = guide_legend(title = "Year"), fill = "none")
lizard.plot.mds