Analysis for Climbing impact manuscript

NMDS Components

Code by Georgia Harrison

March 11, 2021

library packages

library(vegan)
library(ggplot2)
library(readr)
library(dplyr)
library(tidyr)
library(ggpubr) #for combining multiple ggplots into one figure

Import the data

data <- read.csv("~/Research/Manuscript/diff_HARRISON_LGWA_MATRIX.csv", header=TRUE)
data <- data %>%  group_by(T_LOCATION) %>% 
  filter(T_LOCATION %in% c("HB", "TR")| n() == 1)
data <- data %>% group_by(P_LOCATION) %>% 
  filter(P_LOCATION %in% c("FACE")| n() == 1)

note I am using "diff" dataset has the weird two plots removed. something weird happening with two plots for community and lichen nmds. They have NMDS values of 18 and 20 when all of the rest were within much closer range. I looked at the original data and they seemed normal-ish, def low species richness, mostly vascular plants and lacking lichens. I removed these weird plots:

1. Cave Route _ 2_ climbed _ 2
2. Jim Dandy 8 Climbed 3

tutorials I used For example NMDS, NMDS figures and envfit, another NMDS example, ANDONIS, and dealing with missing values in matrix

First, need to smush surface heterogeneity down into one or two values per plot.

create a matrix of only surface features

which(colnames(data)=="LEDGE1_L") #23
## [1] 23
which(colnames(data)=="POCKET1_D") #52
## [1] 52
scom = data[,c(1, 23:52)] #only surface heterogeneity 

Issue 1: this data set contains NA vlues (4 plots without surface heterogenity values)

Get rid of those

na.rows <- scom[rowSums(is.na(scom)) > 0,]
head(na.rows)
## # A tibble: 4 x 31
##   NAME   LEDGE1_L LEDGE1_W LEDGE2_L LEDGE2_W LEDGE3_L LEDGE3_W LEDGE4_L LEDGE4_W
##   <chr>     <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 BONGO~       NA       NA       NA       NA       NA       NA       NA       NA
## 2 BONGO~       NA       NA       NA       NA       NA       NA       NA       NA
## 3 UNCLI~       NA       NA       NA       NA       NA       NA       NA       NA
## 4 UNCLI~       NA       NA       NA       NA       NA       NA       NA       NA
## # ... with 22 more variables: LEDGE5_L <dbl>, LEDGE5_W <dbl>, LEDGE6_L <dbl>,
## #   LEDGE6_W <dbl>, LEDGE7_L <dbl>, LEDGE7_W <dbl>, LEDGE8_L <int>,
## #   LEDGE8_W <dbl>, LEDGE9_L <int>, LEDGE9_W <int>, CRACK1_L <dbl>,
## #   CRACK1_W <dbl>, CRACK1_D <dbl>, CRACK2_L <dbl>, CRACK2_W <dbl>,
## #   CRACK2_D <dbl>, CRACK3_L <int>, CRACK3_W <int>, CRACK3_D <dbl>,
## #   POCKET1_L <int>, POCKET1_W <dbl>, POCKET1_D <dbl>

#there are 4 rows that contain NA values 
which((scom$NAME)=="BONGO FURY_7_UNCLIMBED_1") #24
## [1] 24
which((scom$NAME)=="BONGO FURY_8_UNCLIMBED_1") #26
## [1] 26
which((scom$NAME)=="UNCLIMBED_4_7_UNCLIMBED_4") #265
## [1] 265
which((scom$NAME)=="UNCLIMBED_4_8_UNCLIMBED_4") #266
## [1] 266

newscom <- scom[ !(scom$NAME %in% c("BONGO FURY_7_UNCLIMBED_1", "BONGO FURY_8_UNCLIMBED_1", "UNCLIMBED_4_7_UNCLIMBED_4", "UNCLIMBED_4_8_UNCLIMBED_4")), ]

which(is.na(newscom)) #double check that you actually removed the NAs
## integer(0)

#remove the same rows from the data
data1 <- data[ !(data$NAME %in% c("BONGO FURY_7_UNCLIMBED_1", "BONGO FURY_8_UNCLIMBED_1", "UNCLIMBED_4_7_UNCLIMBED_4", "UNCLIMBED_4_8_UNCLIMBED_4")), ]

Issue 2: Some plots (n=27) have 0 surface features.

I will add 1 to everything to get around this.

newscom <- newscom[,2:31] #remove the name column
newscom <- newscom + 1 #add 1 to all values
m_scom = as.matrix(newscom) #convert com to a matrix

NMDS of just surface features

First, let's try with one dimension

set.seed(123)
s1nmds = metaMDS(m_scom, k=1)

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.3301436 
## Run 1 stress 0.4936894 
## Run 2 stress 0.4648671 
## Run 3 stress 0.4930178 
## Run 4 stress 0.5738235 
## Run 5 stress 0.4923481 
## Run 6 stress 0.429893 
## Run 7 stress 0.4482455 
## Run 8 stress 0.49192 
## Run 9 stress 0.5737158 
## Run 10 stress 0.5736427 
## Run 11 stress 0.4929222 
## Run 12 stress 0.4281314 
## Run 13 stress 0.489482 
## Run 14 stress 0.4399156 
## Run 15 stress 0.4896801 
## Run 16 stress 0.4997382 
## Run 17 stress 0.5746972 
## Run 18 stress 0.5740741 
## Run 19 stress 0.4905896 
## Run 20 stress 0.4484604 
## *** No convergence -- monoMDS stopping criteria:
##      2: stress ratio > sratmax
##     18: scale factor of the gradient < sfgrmin

stressplot(s1nmds)

s1nmds

## 
## Call:
## metaMDS(comm = m_scom, k = 1) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(m_scom)) 
## Distance: bray 
## 
## Dimensions: 1 
## Stress:     0.3301436 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(m_scom))'

#pretty high stress (over 0.3)

k1surface = as.data.frame(scores(s1nmds))
k1surface$CL_UNCL = data1$CL_UNCL
k1surface$T_LOCATION = data1$T_LOCATION
k1surface$name = data1$NAME
k1surface$s_1 <- k1surface$NMDS1

And with two dimensions

s2nmds = metaMDS(m_scom, k=2)

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.157139 
## Run 1 stress 0.1868757 
## Run 2 stress 0.1861934 
## Run 3 stress 0.1708245 
## Run 4 stress 0.1827327 
## Run 5 stress 0.1738167 
## Run 6 stress 0.2015736 
## Run 7 stress 0.216196 
## Run 8 stress 0.1879059 
## Run 9 stress 0.1934866 
## Run 10 stress 0.2201841 
## Run 11 stress 0.2075779 
## Run 12 stress 0.1718933 
## Run 13 stress 0.2059365 
## Run 14 stress 0.1983714 
## Run 15 stress 0.1616603 
## Run 16 stress 0.1810809 
## Run 17 stress 0.1982177 
## Run 18 stress 0.1901671 
## Run 19 stress 0.180665 
## Run 20 stress 0.1774621 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     13: stress ratio > sratmax
##      6: scale factor of the gradient < sfgrmin

stressplot(s2nmds)

s2nmds

## 
## Call:
## metaMDS(comm = m_scom, k = 2) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(m_scom)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.157139 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(m_scom))'

#low stress (0.15)
#better r2

k2surface = as.data.frame(scores(s2nmds))
k2surface$CL_UNCL = data1$CL_UNCL
k2surface$T_LOCATION = data1$T_LOCATION
k2surface$name = data1$NAME
k2surface$s_1 <- k2surface$NMDS1
k2surface$s_2 <- k2surface$NMDS2

NMDS to community data with envfit

with one dimensional surface heterogeneity

com = data1[,54:212] #community data - all species
m_com = as.matrix(com) #convert com to a matrix

#create an environment data set.
env <- data1[,c("SLOPE", "ASPECT_NORTH", "ASPECT_EAST")]
env$s_1 <- k1surface$s_1
head(env)

## # A tibble: 6 x 4
##   SLOPE ASPECT_NORTH ASPECT_EAST     s_1
##   <int>        <dbl>       <dbl>   <dbl>
## 1   130        -0.91        0.42 -0.252 
## 2   140        -0.91        0.42 -0.0735
## 3   130        -0.91        0.42 -0.176 
## 4    90        -0.91        0.42 -0.242 
## 5   110        -0.91        0.42  0.0114
## 6   140        -0.91        0.42 -0.0739

#rename columns so it will look nice in the graph
colnames(env)[1]  <- "Slope"
colnames(env)[2]  <- "Northness"
colnames(env)[3]  <- "Eastness"
colnames(env)[4]  <- "SH"


#nmds code - same from above 
set.seed(123)
en2nmds = metaMDS(m_com, k=3, distance = "bray")

## Wisconsin double standardization
## Run 0 stress 0.1788156 
## Run 1 stress 0.1764567 
## ... New best solution
## ... Procrustes: rmse 0.02650077  max resid 0.1733628 
## Run 2 stress 0.1749737 
## ... New best solution
## ... Procrustes: rmse 0.009823129  max resid 0.132702 
## Run 3 stress 0.1784228 
## Run 4 stress 0.1756506 
## Run 5 stress 0.1761569 
## Run 6 stress 0.1778584 
## Run 7 stress 0.1762312 
## Run 8 stress 0.1768725 
## Run 9 stress 0.1755505 
## Run 10 stress 0.1756811 
## Run 11 stress 0.1806595 
## Run 12 stress 0.1782912 
## Run 13 stress 0.1749562 
## ... New best solution
## ... Procrustes: rmse 0.004491965  max resid 0.04580302 
## Run 14 stress 0.1751377 
## ... Procrustes: rmse 0.007414258  max resid 0.07901636 
## Run 15 stress 0.1781558 
## Run 16 stress 0.1770615 
## Run 17 stress 0.1763083 
## Run 18 stress 0.1752335 
## ... Procrustes: rmse 0.01141235  max resid 0.1663677 
## Run 19 stress 0.1774953 
## Run 20 stress 0.1783089 
## *** No convergence -- monoMDS stopping criteria:
##     18: no. of iterations >= maxit
##      2: stress ratio > sratmax

stressplot(en2nmds)

en2nmds

## 
## Call:
## metaMDS(comm = m_com, distance = "bray", k = 3) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(m_com) 
## Distance: bray 
## 
## Dimensions: 3 
## Stress:     0.1749562 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(m_com)'

en2 = envfit(en2nmds, env, permutations = 999, na.rm=TRUE)
en2

## 
## ***VECTORS
## 
##              NMDS1    NMDS2     r2 Pr(>r)    
## Slope      0.24298  0.97003 0.0562  0.001 ***
## Northness  0.77145 -0.63629 0.5066  0.001 ***
## Eastness  -0.98486  0.17338 0.3479  0.001 ***
## SH         0.99802  0.06290 0.0638  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## 10 observations deleted due to missingness

Figures for community NMDS in general and envfit

#extract values to plot with ggplot
e2.data.scores = as.data.frame(scores(en2nmds))
e2.data.scores$Site = data1$T_LOCATION
e2.data.scores$Climbing = data1$CL_UNCL
en_coord_cont = as.data.frame(scores(en2, "vectors")) * ordiArrowMul(en2)
en_coord_cat = as.data.frame(scores(en2, "factors")) * ordiArrowMul(en2)

theme_set(
  theme(axis.title = element_text(size = 10, face = "bold", colour = "black"),
        panel.background = element_blank(), 
        panel.border = element_rect(fill = NA, colour = "black"), 
        axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.key = element_blank(), 
        legend.title = element_text(size = 10, face = "bold", colour = "black"),
        legend.text = element_text(size = 9, colour = "black")))

#just dots - colored by climbing
gg <- ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2)) 
gg <- gg + geom_point(data = e2.data.scores, aes(colour = Climbing), size = 2, alpha = 0.6) + scale_colour_manual(values = c("orange", "steelblue"))
gg <- gg + labs(colour = "Climbing")
gg

#by site
gg1 <- ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2, col=Site)) 
gg1 <- gg1 + geom_point(data = e2.data.scores, aes(colour = Site), size = 2, alpha = 0.6) + 
  scale_colour_manual(values = c("orange", "steelblue")) 
gg1 <- gg1 + stat_ellipse(size=1.5)
gg1 <- gg1+ labs(colour = "Site")
gg1

#by climbing
gg2 <- ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2, col=Climbing)) 
gg2 <- gg2 + geom_point(data = e2.data.scores, aes(colour = Climbing), size = 2, alpha = 0.6) + 
  scale_colour_manual(values = c("orange", "steelblue")) 
gg2 <- gg2 + stat_ellipse(size=1.5)
gg2 <- gg2 + labs(colour = "Climbing")
gg2

community_fig <- ggarrange(gg1, gg2, ncol=2, legend = "top")
community_fig 

#now lets add in environmental vectors
#only vectors
gg = ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2))
gg <- gg + geom_point(data = e2.data.scores, aes(colour = Climbing), alpha = 0)
gg <- gg + theme(legend.position = "none")
gg = gg + geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
         data = en_coord_cont, size =1, alpha = 0.7, colour = "black") 
gg = gg + geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), 
                 check_overlap = TRUE, colour = "black", 
                 fontface = "bold", label = row.names(en_coord_cont), vjust="outward")
gg

#add vectors to gg1 - by site
gg3 <- gg1 + geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
         data = en_coord_cont, size =1, alpha = 0.7, colour = "black") 
gg3 <- gg3 + geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), 
                 check_overlap = TRUE, colour = "black", 
                 fontface = "bold", label = row.names(en_coord_cont), vjust="inward")
gg3

#add vectors to gg2 - by climbing
gg4 <- gg2 + geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
         data = en_coord_cont, size =1, alpha = 0.7, colour = "black") 
gg4 <- gg4 + geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), 
                 check_overlap = TRUE, colour = "black", 
                 fontface = "bold", label = row.names(en_coord_cont), vjust="inward")
gg4

env_fig <- ggarrange(gg3, gg4, ncol=2, legend = "top")
env_fig 

Envfit to community data

with two dimensional surface heterogeneity. I think only using 1-D surface feature is better

com = data1[,54:212] #community data - all species
m_com = as.matrix(com) #convert com to a matrix

#create an environment data set.
env <- data1[,c("SLOPE", "ASPECT_NORTH", "ASPECT_EAST")]
env$s_1 <- k2surface$s_1
env$s_2 <- k2surface$s_2
#rename columns so it will look nice in the graph
colnames(env)[1]  <- "Slope"
colnames(env)[2]  <- "Northness"
colnames(env)[3]  <- "Eastness"
colnames(env)[4]  <- "SH1"
colnames(env)[5]  <- "SH2"

#nmds code - same from above 
set.seed(123)
en2nmds = metaMDS(m_com, k=3, distance = "bray")

## Wisconsin double standardization
## Run 0 stress 0.1788156 
## Run 1 stress 0.1764567 
## ... New best solution
## ... Procrustes: rmse 0.02650077  max resid 0.1733628 
## Run 2 stress 0.1749737 
## ... New best solution
## ... Procrustes: rmse 0.009823129  max resid 0.132702 
## Run 3 stress 0.1784228 
## Run 4 stress 0.1756506 
## Run 5 stress 0.1761569 
## Run 6 stress 0.1778584 
## Run 7 stress 0.1762312 
## Run 8 stress 0.1768725 
## Run 9 stress 0.1755505 
## Run 10 stress 0.1756811 
## Run 11 stress 0.1806595 
## Run 12 stress 0.1782912 
## Run 13 stress 0.1749562 
## ... New best solution
## ... Procrustes: rmse 0.004491965  max resid 0.04580302 
## Run 14 stress 0.1751377 
## ... Procrustes: rmse 0.007414258  max resid 0.07901636 
## Run 15 stress 0.1781558 
## Run 16 stress 0.1770615 
## Run 17 stress 0.1763083 
## Run 18 stress 0.1752335 
## ... Procrustes: rmse 0.01141235  max resid 0.1663677 
## Run 19 stress 0.1774953 
## Run 20 stress 0.1783089 
## *** No convergence -- monoMDS stopping criteria:
##     18: no. of iterations >= maxit
##      2: stress ratio > sratmax

stressplot(en2nmds)

en2nmds

## 
## Call:
## metaMDS(comm = m_com, distance = "bray", k = 3) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(m_com) 
## Distance: bray 
## 
## Dimensions: 3 
## Stress:     0.1749562 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(m_com)'

en2 = envfit(en2nmds, env, permutations = 999, na.rm=TRUE)
en2

## 
## ***VECTORS
## 
##              NMDS1    NMDS2     r2 Pr(>r)    
## Slope      0.24298  0.97003 0.0562  0.001 ***
## Northness  0.77145 -0.63629 0.5066  0.001 ***
## Eastness  -0.98486  0.17338 0.3479  0.001 ***
## SH1        0.99980 -0.01993 0.0720  0.001 ***
## SH2       -0.93824 -0.34598 0.0062  0.378    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## 10 observations deleted due to missingness

Figures

#extract values to plot with ggplot
e2.data.scores = as.data.frame(scores(en2nmds))
e2.data.scores$Site = data1$T_LOCATION
e2.data.scores$Climbing = data1$CL_UNCL
en_coord_cont = as.data.frame(scores(en2, "vectors")) * ordiArrowMul(en2)
en_coord_cat = as.data.frame(scores(en2, "factors")) * ordiArrowMul(en2)


theme_set(
  theme(axis.title = element_text(size = 10, face = "bold", colour = "black"),
        panel.background = element_blank(), 
        panel.border = element_rect(fill = NA, colour = "black"), 
        axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.key = element_blank(), 
        legend.title = element_text(size = 10, face = "bold", colour = "black"),
        legend.text = element_text(size = 9, colour = "black")))


#just dots - colored by climbing
gg <- ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2)) 
gg <- gg + geom_point(data = e2.data.scores, aes(colour = Climbing), size = 2, alpha = 0.6) + scale_colour_manual(values = c("orange", "steelblue"))
gg <- gg + labs(colour = "Climbing")
gg

#by site
gg1 <- ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2, col=Site)) 
gg1 <- gg1 + geom_point(data = e2.data.scores, aes(colour = Site), size = 2, alpha = 0.6) + 
  scale_colour_manual(values = c("orange", "steelblue")) 
gg1 <- gg1 + stat_ellipse(size=1.5)
gg1 <- gg1+ labs(colour = "Site")
gg1

#by climbing
gg2 <- ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2, col=Climbing)) 
gg2 <- gg2 + geom_point(data = e2.data.scores, aes(colour = Climbing), size = 2, alpha = 0.6) + 
  scale_colour_manual(values = c("orange", "steelblue")) 
gg2 <- gg2 + stat_ellipse(size=1.5)
gg2 <- gg2 + labs(colour = "Climbing")
gg2

community_fig <- ggarrange(gg1, gg2, ncol=2, legend = "top")
community_fig 

#now lets add in environmental vectors
#only vectors
gg = ggplot(data = e2.data.scores, aes(x = NMDS1, y = NMDS2))
gg = gg + geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
         data = en_coord_cont, size =1, alpha = 0.7, colour = "black") 
gg = gg + geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), 
                 check_overlap = TRUE, colour = "black", 
                 fontface = "bold", label = row.names(en_coord_cont), vjust="inward")
gg

#add vectors to gg1 - by site
gg3 <- gg1 + geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
         data = en_coord_cont, size =1, alpha = 0.7, colour = "black") 
gg3 <- gg3 + geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), 
                 check_overlap = TRUE, colour = "black", 
                 fontface = "bold", label = row.names(en_coord_cont), vjust="inward")
gg3

#add vectors to gg2 - by climbing
gg4 <- gg2 + geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
         data = en_coord_cont, size =1, alpha = 0.7, colour = "black") 
gg4 <- gg4 + geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), 
                 check_overlap = TRUE, colour = "black", 
                 fontface = "bold", label = row.names(en_coord_cont), vjust="inward")
gg4

env_fig <- ggarrange(gg3, gg4, ncol=2, legend = "top")
env_fig 

#ggsave(env_fig, device = "jpeg", plot = env_fig, dpi=700, limitsize = TRUE)



ggsave("nmds.jpeg", plot = env_fig,  device = "jpeg", width = 8, height = 4, dpi = 300, scale = 1)

ADONIS

A PERMANOVA

env$Climbing <- data1$CL_UNCL
env$Site <- data1$T_LOCATION
adonis(m_com ~ Climbing*Site, data=env) 

## 
## Call:
## adonis(formula = m_com ~ Climbing * Site, data = env) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Climbing        1     3.452  3.4520  15.973 0.04153  0.001 ***
## Site            1     9.171  9.1713  42.436 0.11033  0.001 ***
## Climbing:Site   1     0.479  0.4786   2.214 0.00576  0.026 *  
## Residuals     324    70.024  0.2161         0.84238           
## Total         327    83.125                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple linear regressions

Between: slope, northnes, eastness, and SH to richness, diverstiy and cover

First, recall the indic data set (or create again)

Doing this all with data1 file - which excluded the four plot with blank SH values

total.abundance.matrix <- data1[,54:212]
lichen.abundance.matrix <- data1[,54:133]
moss.abundance.matrix <- data1 [,134:174]
plant.abundance.matrix <- data1 [,175:212]

indic <- data1[,c("CL_UNCL", "T_LOCATION", "SLOPE", "ASPECT_NORTH", "ASPECT_EAST")]
indic$SH1 <- k2surface$s_1
indic$SH2 <- k2surface$s_2
indic$T_LOCATION <- as.factor((indic$T_LOCATION))
indic$CL_UNCL <- as.factor((indic$CL_UNCL))

#rename columns so it will look nice in the graph
colnames(indic)[1]  <- "Climbing"
colnames(indic)[2]  <- "Site"
colnames(indic)[3]  <- "Slope"
colnames(indic)[4]  <- "Northness"
colnames(indic)[5]  <- "Eastness"




## Determine speciess richnes, shannon and total abundance for each 
#all
indic$AllRichness <- rowSums(total.abundance.matrix>0)
indic$AllShannon <- diversity(total.abundance.matrix) # shannon is default
indic$AllCover <- rowSums(total.abundance.matrix) 

#lichens
indic$LRichness <- rowSums(lichen.abundance.matrix>0)
indic$LShannon <- diversity(lichen.abundance.matrix) # shannon is default
indic$LCover <- rowSums(lichen.abundance.matrix)

#mosses
indic$MRichness <- rowSums(moss.abundance.matrix>0)
indic$MShannon <- diversity(moss.abundance.matrix) # shannon is default
indic$MCover <- rowSums(moss.abundance.matrix)

#plants
indic$PRichness <- rowSums(plant.abundance.matrix>0)
indic$PShannon <- diversity(plant.abundance.matrix) # shannon is default
indic$PCover <- rowSums(plant.abundance.matrix)

linear models

richness

lmr <- lm(AllRichness ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmr)

## 
## Call:
## lm(formula = AllRichness ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2107 -1.6881 -0.1715  1.4255  6.6388 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.545946   0.561474  11.658  < 2e-16 ***
## Slope        0.012028   0.005315   2.263 0.024325 *  
## Northness    0.378621   0.258084   1.467 0.143372    
## Eastness     1.095402   0.280050   3.911 0.000113 ***
## SH1         -4.971966   1.098437  -4.526 8.54e-06 ***
## SH2          0.637720   1.238958   0.515 0.607111    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.305 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.1665, Adjusted R-squared:  0.1532 
## F-statistic: 12.47 on 5 and 312 DF,  p-value: 4.807e-11

lmr <- lm(LRichness ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmr)

## 
## Call:
## lm(formula = LRichness ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7215 -1.3392 -0.1028  1.3471  6.9851 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.090269   0.516087  11.801  < 2e-16 ***
## Slope        0.007450   0.004886   1.525  0.12829    
## Northness    0.521987   0.237222   2.200  0.02851 *  
## Eastness     0.830402   0.257412   3.226  0.00139 ** 
## SH1         -2.853067   1.009644  -2.826  0.00502 ** 
## SH2         -1.127661   1.138806  -0.990  0.32284    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.119 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.08262,    Adjusted R-squared:  0.06792 
## F-statistic:  5.62 on 5 and 312 DF,  p-value: 5.618e-05

lmr <- lm(MRichness ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmr)

## 
## Call:
## lm(formula = MRichness ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5036 -0.7071 -0.0647  0.2396  4.0244 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.5965502  0.2195005   2.718  0.00694 **
## Slope        0.0002952  0.0020780   0.142  0.88711   
## Northness   -0.1478197  0.1008944  -1.465  0.14390   
## Eastness     0.2602608  0.1094817   2.377  0.01805 * 
## SH1         -1.2906738  0.4294184  -3.006  0.00287 **
## SH2          0.8767494  0.4843534   1.810  0.07124 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9012 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.1288, Adjusted R-squared:  0.1149 
## F-statistic: 9.227 on 5 and 312 DF,  p-value: 3.316e-08

lmr <- lm(PRichness ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmr)

## 
## Call:
## lm(formula = PRichness ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9527 -0.3703 -0.2198 -0.0406  6.1897 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.140873   0.191305  -0.736   0.4621  
## Slope        0.004283   0.001811   2.365   0.0186 *
## Northness    0.004454   0.087934   0.051   0.9596  
## Eastness     0.004740   0.095419   0.050   0.9604  
## SH1         -0.828225   0.374259  -2.213   0.0276 *
## SH2          0.888631   0.422138   2.105   0.0361 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7854 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.05933,    Adjusted R-squared:  0.04425 
## F-statistic: 3.935 on 5 and 312 DF,  p-value: 0.001785

Diversity

lmd <- lm(AllShannon ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmd)

## 
## Call:
## lm(formula = AllShannon ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13179 -0.18699  0.03158  0.18017  0.67064 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.6646003  0.0711031  23.411  < 2e-16 ***
## Slope        0.0012507  0.0006731   1.858 0.064103 .  
## Northness    0.0165220  0.0326829   0.506 0.613548    
## Eastness     0.1207083  0.0354646   3.404 0.000752 ***
## SH1         -0.6372191  0.1391021  -4.581  6.7e-06 ***
## SH2          0.0761305  0.1568973   0.485 0.627857    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2919 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.1655, Adjusted R-squared:  0.1522 
## F-statistic: 12.38 on 5 and 312 DF,  p-value: 5.754e-11

lmd <- lm(LShannon ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmd)

## 
## Call:
## lm(formula = LShannon ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.66129 -0.18110  0.03664  0.21449  0.72491 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.5975967  0.0787207  20.294  < 2e-16 ***
## Slope        0.0007810  0.0007452   1.048  0.29547    
## Northness    0.0527646  0.0361844   1.458  0.14579    
## Eastness     0.1010849  0.0392641   2.574  0.01050 *  
## SH1         -0.4280712  0.1540048  -2.780  0.00577 ** 
## SH2         -0.1679662  0.1737065  -0.967  0.33432    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3232 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.06739,    Adjusted R-squared:  0.05245 
## F-statistic: 4.509 on 5 and 312 DF,  p-value: 0.0005536

lmd <- lm(MShannon ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmd)

## 
## Call:
## lm(formula = MShannon ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38344 -0.16358 -0.09692  0.00505  1.28685 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0849413  0.0735156   1.155 0.248803    
## Slope        0.0003022  0.0006960   0.434 0.664381    
## Northness    0.0260490  0.0337918   0.771 0.441369    
## Eastness     0.0777445  0.0366679   2.120 0.034775 *  
## SH1         -0.4785682  0.1438218  -3.328 0.000981 ***
## SH2          0.2399494  0.1622207   1.479 0.140108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3018 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.07884,    Adjusted R-squared:  0.06408 
## F-statistic: 5.341 on 5 and 312 DF,  p-value: 1e-04

lmd <- lm(PShannon ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmd)

## 
## Call:
## lm(formula = PShannon ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26315 -0.08281 -0.04469 -0.00845  1.56847 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.0680690  0.0549872  -1.238   0.2167  
## Slope        0.0011901  0.0005205   2.286   0.0229 *
## Northness   -0.0238424  0.0252752  -0.943   0.3463  
## Eastness    -0.0275238  0.0274264  -1.004   0.3164  
## SH1         -0.1832873  0.1075740  -1.704   0.0894 .
## SH2          0.2625735  0.1213358   2.164   0.0312 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2258 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.05929,    Adjusted R-squared:  0.04421 
## F-statistic: 3.933 on 5 and 312 DF,  p-value: 0.001795

Cover

lmc <- lm(AllCover ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmc)

## 
## Call:
## lm(formula = AllCover ~ Slope + Northness + Eastness + SH1 + 
##     SH2, data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.095  -7.311   0.149   6.504  40.233 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  29.08100    2.54157  11.442  < 2e-16 ***
## Slope         0.04198    0.02406   1.745  0.08198 .  
## Northness    -0.49683    1.16825  -0.425  0.67093    
## Eastness      3.38657    1.26768   2.671  0.00795 ** 
## SH1         -10.66936    4.97219  -2.146  0.03266 *  
## SH2          -6.11274    5.60828  -1.090  0.27658    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.43 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.0919, Adjusted R-squared:  0.07735 
## F-statistic: 6.315 on 5 and 312 DF,  p-value: 1.334e-05

lmc <- lm(LCover ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmc)

## 
## Call:
## lm(formula = LCover ~ Slope + Northness + Eastness + SH1 + SH2, 
##     data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.150  -7.183   0.207   6.882  44.136 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 27.89059    2.77446  10.053   <2e-16 ***
## Slope        0.02954    0.02627   1.125   0.2615    
## Northness    0.83225    1.27529   0.653   0.5145    
## Eastness     2.93613    1.38384   2.122   0.0346 *  
## SH1         -7.08733    5.42780  -1.306   0.1926    
## SH2         -8.18394    6.12217  -1.337   0.1823    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.39 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.03822,    Adjusted R-squared:  0.0228 
## F-statistic:  2.48 on 5 and 312 DF,  p-value: 0.03195

lmc <- lm(MCover ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmc)

## 
## Call:
## lm(formula = MCover ~ Slope + Northness + Eastness + SH1 + SH2, 
##     data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4453 -2.6242 -0.4988  1.1355 12.8761 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.5793377  0.7977194   1.980 0.048603 *  
## Slope        0.0007903  0.0075518   0.105 0.916720    
## Northness   -1.3681669  0.3666756  -3.731 0.000226 ***
## Eastness     0.4502733  0.3978838   1.132 0.258642    
## SH1         -1.5134237  1.5606134  -0.970 0.332916    
## SH2          0.7736467  1.7602611   0.440 0.660599    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.275 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.1114, Adjusted R-squared:  0.0972 
## F-statistic: 7.826 on 5 and 312 DF,  p-value: 5.885e-07

lmc <- lm(PCover ~ Slope + Northness + Eastness + SH1 + SH2, data = indic)
summary(lmc)

## 
## Call:
## lm(formula = PCover ~ Slope + Northness + Eastness + SH1 + SH2, 
##     data = indic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3072 -0.9424 -0.6395 -0.2375 13.0944 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.3889280  0.5491375  -0.708   0.4793  
## Slope        0.0116513  0.0051985   2.241   0.0257 *
## Northness    0.0390913  0.2524137   0.155   0.8770  
## Eastness     0.0001668  0.2738969   0.001   0.9995  
## SH1         -2.0685994  1.0743017  -1.926   0.0551 .
## SH2          1.2975529  1.2117360   1.071   0.2851  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.255 on 312 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.04131,    Adjusted R-squared:  0.02594 
## F-statistic: 2.689 on 5 and 312 DF,  p-value: 0.02136

Lichen only nmds

lichen output is not that different from all community. Just use all community

lcom = data[,54:133] #lichen only data 
m_lcom = as.matrix(lcom) #convert com to a matrix

#create an environment data set. But I only want climbing and site in
env = data[,c(20,6)] #same env data set as above, only with site and climbing

#nmds code
set.seed(123)
l_nmds = metaMDS(m_lcom, k=3, distance = "bray")
#k of 3 gives better stress (0.17) vs (0.25 from k=2)
stressplot(l_nmds)
l_nmds

en = envfit(l_nmds, env, permutations = 999, na.rm = TRUE)
en #only factors so far

plot(l_nmds)
plot(en)

Figures for lichen only NMDS

#extract values to plot with ggplot
l.data.scores = as.data.frame(scores(l_nmds))
l.data.scores$CL_UNCL = data$CL_UNCL
l.data.scores$T_LOCATION = data$T_LOCATION


theme_set(
  theme(axis.title = element_text(size = 10, face = "bold", colour = "black"),
        panel.background = element_blank(), 
        panel.border = element_rect(fill = NA, colour = "black"), 
        axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.key = element_blank(), 
        legend.title = element_text(size = 10, face = "bold", colour = "black"),
        legend.text = element_text(size = 9, colour = "black")))

gg <- ggplot(data = l.data.scores, aes(x = NMDS1, y = NMDS2)) 
gg <- gg + geom_point(data = l.data.scores, aes(colour = CL_UNCL), size = 2, alpha = 0.6) + scale_colour_manual(values = c("orange", "steelblue"))
gg <- gg + labs(colour = "Climbing")
gg

#by site
gg3 <- ggplot(data = l.data.scores, aes(x = NMDS1, y = NMDS2, col=T_LOCATION)) 
gg3 <- gg3 + geom_point(data = l.data.scores, aes(colour = T_LOCATION), size = 2, alpha = 0.6) + 
  scale_colour_manual(values = c("orange", "steelblue")) 
gg3 <- gg3 + stat_ellipse(size=1.5)
gg3 <- gg3+ labs(colour = "Site")
gg3

#by climbing
gg4 <- ggplot(data = l.data.scores, aes(x = NMDS1, y = NMDS2, col=CL_UNCL)) 
gg4 <- gg4 + geom_point(data = l.data.scores, aes(colour = CL_UNCL), size = 2, alpha = 0.6) + 
  scale_colour_manual(values = c("orange", "steelblue")) 
gg4 <- gg4 + stat_ellipse(size=1.5)
gg4 <- gg4 + labs(colour = "Climbing")
gg4


lichen_fig <- ggarrange(gg3, gg4, ncol=2, legend = "top")
lichen_fig