Bat_data_ordination

Set up stuff

rm(list=ls())
pacman::p_load(plyr,dplyr,ggplot2,tidyr, vegan)
setwd("C:/Users/Katy/Desktop/thesisr")

Organize data

# Read in the .csv file and save as a new object.
data.allparks <- read.csv("./MixedModel_allparks.csv")

# Change the Year column to a factor class.
data.allparks$Year <- as.factor(data.allparks$Year)

# Add a new column for Calls per Night.
data.allparks$CallsPN <- round(data.allparks$Calls/data.allparks$Nights, 3)

# Temporarily subset the data by year. Remove the Park,Year, Nights, Calls columns.
# Rearrange data so each species has a column for each year's Calls Per Night value.
# Rename Calls Per Night columns to include year information. Then merge all data
# back to one data frame, with rows for each Site.  
data.allparks.2016 <- data.allparks %>% 
        subset(Year == "2016") %>%
        subset(select = -c(Park,Year,Nights,Calls)) %>%
        spread(key = Species, value = CallsPN) %>%
        rename(EPFU.cpn.2016 = EPFU, LABO.cpn.2016 = LABO, LACI.cpn.2016 = LACI,
               LANO.cpn.2016 = LANO, MYLU.cpn.2016 = MYLU, MYSE.cpn.2016 = MYSE,
               MYSO.cpn.2016 = MYSO, NYHU.cpn.2016 = NYHU, PESU.cpn.2016 = PESU)

data.allparks.2017 <- data.allparks %>% 
        subset(Year == "2017") %>%
        subset(select = -c(Park,Year,Nights,Calls)) %>%
        spread(key = Species, value = CallsPN) %>%
        rename(EPFU.cpn.2017 = EPFU, LABO.cpn.2017 = LABO, LACI.cpn.2017 = LACI,
               LANO.cpn.2017 = LANO, MYLU.cpn.2017 = MYLU, MYSE.cpn.2017 = MYSE,
               MYSO.cpn.2017 = MYSO, NYHU.cpn.2017 = NYHU, PESU.cpn.2017 = PESU)

data.ord <- merge(data.allparks.2016, data.allparks.2017, by = "Site") 

# Because the data have to be all numeric for the ordination function, remove the 
# column containing the site names. 
data.ord.num <- data.ord %>%   
                subset(select = -c(Site))

# Because the ordination function can't handle "n/a" values (or at least, I don't
# know how to make it take these), remove the columns for the species that weren't
# an option at all parks - MYSO, NYHU, and PESU. 
data.ord.num2 <- data.ord.num %>% 
                 subset(select= -c(MYSO.cpn.2016, NYHU.cpn.2016, PESU.cpn.2016, 
                                   MYSO.cpn.2017, NYHU.cpn.2017, PESU.cpn.2017))

# Check that it looks right before doing the MDS
str(data.ord.num2)

## 'data.frame':    185 obs. of  12 variables:
##  $ EPFU.cpn.2016: num  2.429 0.4 0.429 1.545 0.091 ...
##  $ LABO.cpn.2016: num  13 0.1 7.143 2.636 0.455 ...
##  $ LACI.cpn.2016: num  5 0.4 12 3.727 0.182 ...
##  $ LANO.cpn.2016: num  20.571 3.3 3.143 9.545 0.727 ...
##  $ MYLU.cpn.2016: num  340 27.6 13.9 18.5 66.9 ...
##  $ MYSE.cpn.2016: num  2.29 1.9 1.14 0 8.91 ...
##  $ EPFU.cpn.2017: num  10.833 0.889 0.111 7.143 0 ...
##  $ LABO.cpn.2017: num  7.667 0.222 0 15.571 0 ...
##  $ LACI.cpn.2017: num  7.917 0.778 0.222 80.429 0 ...
##  $ LANO.cpn.2017: num  124.667 5.222 0.889 17.714 0 ...
##  $ MYLU.cpn.2017: num  192.583 32.222 0.556 41.429 0 ...
##  $ MYSE.cpn.2017: num  0.083 2 0.889 0.571 0 ...

head(data.ord.num2)

##   EPFU.cpn.2016 LABO.cpn.2016 LACI.cpn.2016 LANO.cpn.2016 MYLU.cpn.2016
## 1         2.429        13.000         5.000        20.571       340.000
## 2         0.400         0.100         0.400         3.300        27.600
## 3         0.429         7.143        12.000         3.143        13.857
## 4         1.545         2.636         3.727         9.545        18.545
## 5         0.091         0.455         0.182         0.727        66.909
## 6         0.167         0.417         0.167         0.500        28.083
##   MYSE.cpn.2016 EPFU.cpn.2017 LABO.cpn.2017 LACI.cpn.2017 LANO.cpn.2017
## 1         2.286        10.833         7.667         7.917       124.667
## 2         1.900         0.889         0.222         0.778         5.222
## 3         1.143         0.111         0.000         0.222         0.889
## 4         0.000         7.143        15.571        80.429        17.714
## 5         8.909         0.000         0.000         0.000         0.000
## 6         7.500         0.000         0.000         0.071         0.071
##   MYLU.cpn.2017 MYSE.cpn.2017
## 1       192.583         0.083
## 2        32.222         2.000
## 3         0.556         0.889
## 4        41.429         0.571
## 5         0.000         0.000
## 6         1.500         1.786

Run the metaMDS

# Run the metaMDS function using mostly default values. 
mds.1 <- metaMDS(data.ord.num2, distance = "bray", k = 2, 
                   try = 20, trymax = 20, engine = "monoMDS", 
                   autotransform =TRUE, noshare = FALSE, 
                   wascores = TRUE, expand = TRUE, trace = 1, plot = FALSE)

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1585962 
## Run 1 stress 0.175029 
## Run 2 stress 0.1662595 
## Run 3 stress 0.1758025 
## Run 4 stress 0.1733935 
## Run 5 stress 0.1725098 
## Run 6 stress 0.1660093 
## Run 7 stress 0.1782417 
## Run 8 stress 0.1656276 
## Run 9 stress 0.1765813 
## Run 10 stress 0.1904926 
## Run 11 stress 0.1744928 
## Run 12 stress 0.1615381 
## Run 13 stress 0.1657807 
## Run 14 stress 0.173871 
## Run 15 stress 0.1811771 
## Run 16 stress 0.1870154 
## Run 17 stress 0.1833334 
## Run 18 stress 0.1701852 
## Run 19 stress 0.1814895 
## Run 20 stress 0.1853286 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     19: stress ratio > sratmax

# Display the default plot
plot.1 <- plot(mds.1, type = "t")

# Declutter the plot
plot.2 <- ordiplot(mds.1, type = "n")
points(plot.2, "sites", pch=21, col="red", cex=0.75)
points(plot.2, "species", pch="+", col="blue", cex=1.2)
text(plot.2, "species", col="blue", cex=0.5)