In brouwern/DRmencia: Data and code to replicate analyses for Latta, Brouwer, et al. 2018. PeerJ.

knitr::opts_chunk$set(echo = TRUE)

NMDS

Fit NMDS

diet_NMDS <- metaMDS(trait_dat[,c(10:15)])

Check stress

stressplot(diet_NMDS)

Base plot of NMDS

ordiplot(diet_NMDS)
orditorp(diet_NMDS,display="species")
ordiellipse(diet_NMDS,
         groups=trait_dat$site,
         label=T,
         kind = "se")

GGPLOT NMDS

https://chrischizinski.github.io/rstats/vegan-ggplot2/

data.scores <- as.data.frame(scores(diet_NMDS))  #Using the scores function from vegan to extract the site scores and convert to a data.frame
data.scores$site <- rownames(data.scores)  # create a column of site names, from the rownames of data.scores
data.scores$grp <- trait_dat$site[-i.0]  #  add the grp variable created earlier
head(data.scores)  #look at the data

species.scores <- as.data.frame(scores(diet_NMDS, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

ggplot() + 
  geom_text(data=species.scores,
            aes(x=NMDS1,y=NMDS2,
                label=species),alpha=0.5) +  # add the species labels
  geom_point(data=data.scores,
             aes(x=NMDS1,y=NMDS2
                 ,shape=site,
                 colour=site),
             size=1) + # add the point markers
  #geom_text(data=data.scores,aes(x=NMDS1,y=NMDS2,label=site),size=6,vjust=0) +  # add the site labels
  #scale_colour_manual(values=c("A" = "red", "B" = "blue")) +
  coord_equal() +
  theme_bw()

Hulls

grp.brdlf <- data.scores[data.scores$grp == "Broadleaf", ][chull(data.scores[data.scores$grp == 
    "Broadleaf", c("NMDS1", "NMDS2")]), ]  # hull values for grp A
grp.corral <- data.scores[data.scores$grp == "El Corral", ][chull(data.scores[data.scores$grp == 
    "El Corral", c("NMDS1", "NMDS2")]), ]  # hull values for grp B
grp.caoba <- data.scores[data.scores$grp == "La Caoba", ][chull(data.scores[data.scores$grp == 
    "La Caoba", c("NMDS1", "NMDS2")]), ]  # hull values for grp B
grp.cueva <- data.scores[data.scores$grp == "La Cueva", ][chull(data.scores[data.scores$grp == 
    "La Cueva", c("NMDS1", "NMDS2")]), ]  # hull values for grp B
grp.morelia<- data.scores[data.scores$grp == "morelia", ][chull(data.scores[data.scores$grp == 
    "morelia", c("NMDS1", "NMDS2")]), ]  # hull values for grp B

hull.data <- rbind(grp.brdlf, 
                   grp.corral,
                   grp.caoba,
                   grp.cueva,
                   grp.morelia)  #combine grp.a and grp.b

ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=grp,group=grp),alpha=0.30) + # add the convex hulls
  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species),alpha=0.5) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=grp,colour=grp),size=4) + # add the point markers
  #scale_colour_manual(values=c("A" = "red", "B" = "blue")) +
  coord_equal() +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())

brouwern/DRmencia documentation built on May 6, 2019, 12:24 p.m.