Veg Structure Ordination

setwd("C:/Users/micayla.lakey/Documents/Project_Info/Veg Structure")
veg.all <- read.csv("./veg.all.csv")
  
 veg.str <-
  veg.all %>% 
            mutate(Litter_2=as.numeric(Litter_2),
              Litter_4=as.numeric(as.character(Litter_4))) %>% 
            mutate(VOR=rowMeans(select(., VOR_N:VOR_W)),
              Litter=rowSums(select(., starts_with("Litter")), na.rm = TRUE)) %>%
           select(year, management, Location, Pasture, Patch, SubPatch, Transect, VOR, Litter) %>% 
            group_by(year, management, Location, Pasture, Patch, SubPatch) %>% 
            summarise(VOR=mean(VOR),
              Litter=mean(Litter)) %>%  
    unite("yearmanloc", c("year", "management", "Location"))
            
veg.co <-  
  veg.all %>%
  #filter(year==2018) %>%
  
  select(year, management, Location, Pasture, Patch, SubPatch, Transect:IntWdy, -Quadrat) %>% 
  
  group_by(year, management, Location, Pasture, Patch, SubPatch) %>% 
  summarise_all(list(mean)) %>%
  ungroup() %>% 
  mutate(sum=select(., POPR:IntWdy) %>% rowSums()) %>% 
  select(-Transect) %>% 
  unite("yearmanloc", c("year", "management", "Location")) 
  
  
  veg.co <- veg.co %>% mutate(NatLeg = replace_na(NatLeg, 0))

#metric mds 

veg.spp <- veg.co %>% select(., POPR:IntWdy)

veg.cap <- capscale(veg.spp~1, distance = "euc", metaMDSdist=TRUE,
                    engine=monoMDS, autotransform=TRUE)

## Square root transformation
## Wisconsin double standardization

var.view(veg.cap, 10)

##                       MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8 MDS9 MDS10
## Eigenvalue            2.23 1.16 0.89 0.74 0.59 0.44 0.35 0.30 0.22  0.13
## Proportion Explained  0.31 0.16 0.12 0.10 0.08 0.06 0.05 0.04 0.03  0.02
## Cumulative Proportion 0.31 0.47 0.60 0.70 0.78 0.84 0.89 0.94 0.97  0.99

sppscores(veg.cap) <- veg.co %>% select(., POPR:IntWdy)

#use veg.cap below

plot(veg.cap, display = "species")

spp.sc <- scores(veg.cap, display = "species") %>%
          as.data.frame() %>%
          rownames_to_column(var="functgrp") 


veg.str<- veg.str %>%
  separate(yearmanloc, c("year", "management", "location"), sep="_") %>%
  mutate(management=as.factor(management), 
         year=as.factor(year)) %>% 
         unite("LocationYear", c("location", "year"), remove=FALSE)



#test env vars

fit.mds <- envfit(veg.cap~management+VOR+Litter, veg.str, choices=c(1:4), strata = veg.str$LocationYear)
fit.mds

## 
## ***VECTORS
## 
##            MDS1     MDS2     MDS3     MDS4     r2 Pr(>r)    
## VOR     0.79560 -0.52657  0.23555  0.18510 0.4536  0.001 ***
## Litter  0.93432  0.01970 -0.33985 -0.10561 0.3419  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                         MDS1    MDS2    MDS3    MDS4
## managementburned     -0.1724 -0.0889  0.0732  0.1644
## managementcontinuous  0.4388  0.1617  0.0134 -0.1036
## managementtwice-over -0.0822 -0.4915 -0.1236 -0.0067
## managementunburned   -0.1863  0.0637  0.0095  0.0227
## 
## Goodness of fit:
##                r2 Pr(>r)  
## management 0.1783  0.091 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999

pairwise.factorfit(veg.cap, veg.str$management, 
                   nperm = 999, p.method = "fdr")

## 
##  Pairwise comparisons using factor fitting to an ordination 
## 
## data:  veg.cap by veg.str$management
## 999 permutations 
## 
##            burned continuous twice-over
## continuous 0.0012 -          -         
## twice-over 0.0012 0.0012     -         
## unburned   0.1160 0.0012     0.0012    
## 
## P value adjustment method: fdr

litfit <- as.data.frame(round(scores(fit.mds, "vectors"), 3)) %>%
          rownames_to_column() #%>%
          #filter(rowname != "VOR")

#gg ordiplot

gg_ordiplot(veg.cap, groups = veg.str$management, 
            spiders=TRUE, ellipse=FALSE, plot=TRUE)

man.gg <- gg_ordiplot(veg.cap, groups = veg.str$management, 
                       spiders=TRUE, ellipse=FALSE, plot=FALSE) 
vegord.gg <-
man.gg$df_spiders %>% 
  plyr::rename(c("x"="X", 
                 "y"="Y") )%>%
  ggplot() + theme_ord(12) +
  geom_vline(xintercept = 0, lty=3, color="darkgrey") +
  geom_hline(yintercept = 0, lty=3, color="darkgrey") +
  labs(x="MDS 1", y="MDS 2") + #took out caption = "Species scores <|0.3| excluded
  geom_segment(aes(x=cntr.x, y=cntr.y,
                   xend=X, yend=Y, color=Group), 
               size=1.0, show.legend = FALSE, alpha = 0.5) +
  geom_segment(data = litfit, aes(x=0, y=0, xend = MDS1*.75, yend=MDS2*.75), color = "dark red",
               lwd = 1.25, arrow = arrow(length = unit(0.02, "npc")))+ #litter and VOR vectors
  geom_text(data = litfit, aes(x=MDS1*.9, y=MDS2*.5, label = rowname), color = "dark red", 
             fontface = "italic", size = 5)+ #litter and VOR labels
  geom_text(data = filter(spp.sc, functgrp!="IntWdy", functgrp!="IntC4"),
            aes(x=MDS1*.2, y=MDS2*.2, label=functgrp), fontface="bold")+ #species labels
  geom_label(aes(x=cntr.x, y=cntr.y, 
                 label=Group, color=Group), 
             fontface="bold", size=4,
             label.size = 0, 
             label.r = unit(0.5, "lines"), 
             show.legend = FALSE)+
  
  #geom_text(data=filter(spp.sc, abs(MDS1)>0.1, abs(MDS2)>0.1),
         #   aes(x=MDS1*0.5, y=MDS2*0.5, label=functgrp), fontface="bold")+
          scale_color_manual(values = cbPal5)

vegord.gg