# Annotated R Code for Reid et al. 2020. Annual understory plant recovery
# dynamics in a temperate woodland mosaic during a decade of ecological
# restoration. Natural Areas Journal.

#Last updated by JLR on 2019-12-09

###############################################################################
# SET UP 
###############################################################################

#Libraries
library(BiodiversityR) #Diversity tools
library(plyr) #Data management tools
library(reshape2) #Data management tools
library(shape) #Draws nice arrows

#Data
veg<-read.csv('Data_Veg.csv',header=TRUE)
pheno<-read.csv('Data_Phenology.csv',header=TRUE)
traits<-read.csv('Data_SpecList.csv',header=TRUE)

#Format data
veg<-subset(veg,L=='H')
veg<-merge(veg,traits,by='Species')
veg<-merge(veg,pheno,by='Species')

#Unique years
Years<-unique(veg$Y)[order(unique(veg$Y))]

#define spring & summer phenology
pheno$PHENO<-ifelse(pheno$Start<=6 & pheno$End<=7,"Spring",'Summer')

###############################################################################
# FUNCTIONS
###############################################################################

## Restoration intervention timeline bars
RxLines<-function(){
  par(xpd=FALSE)
  abline(v=1+(328/365),lty=1,col='black',lwd=1) #prescribed burns
  abline(v=4+(57/365),lty=1,col='black',lwd=1)
  abline(v=7+(74/365),lty=1,col='black',lwd=1)
  abline(v=9+(319/365),lty=1,col='black',lwd=1)
  abline(v=12+(10/365),lty=1,col='black',lwd=1)
  abline(v=6,lty=2,col='black',lwd=1) #Juniper removal
  box()
  axis(side=1,at=seq(2001,2013,1),labels=NA)
}

## Rx toptext
RxToptext<-function(){
  mtext(side=3,line=0.5,text='Fire',at=1+(328/365),las=2)
  mtext(side=3,line=0.5,text='Fire',at=4+(57/365),las=2)
  mtext(side=3,line=0.5,text='Fire',at=7+(74/365),las=2)
  mtext(side=3,line=0.5,text='Fire',at=9+(319/365),las=2)
  mtext(side=3,line=0.5,text='Fire',at=12+(10/365),las=2)
  mtext(side=3,line=0.5,text='Thinning',at=6,las=2)
}


###############################################################################
###############################################################################
# SUMMARY STATISTICS
###############################################################################

#Total number of stems sampled
print(stemCount<-sum(veg$Stems))

#Percent of stems identified to species level
print(stemCompleteIDCount<-sum(veg[veg$idSpecLevel==1,]$Stems))
print(percentIDcomplete<-(stemCompleteIDCount/stemCount)*100)

#Total no. species encountered
print(noSpecies<-length(unique(veg[veg$idSpecLevel==1,]$Species)))

#Total no. genera encountered
print(noGenus<-length(unique(veg[veg$idSpecLevel==1,]$Genus)))

#Total no. families encountered  
print(noFamily<-length(unique(veg[veg$idSpecLevel==1,]$Family)))

#Top 3 families for stems
FamStems<-setNames(aggregate(veg$Stems,by=list(veg$Family),sum),c('Family','Stems'))
FamStems<-FamStems[order(-FamStems$Stems),]
FamStems$Percent<-(FamStems$Stems/sum(FamStems$Stems))*100

#Top 3 families for diversity
FamDiv<-with(veg[veg$idSpecLevel==1,],setNames(aggregate(Species,by=list(Family),function(x)length(unique(x))),c('Family','Species')))
FamDiv<-FamDiv[order(-FamDiv$Species),]

#Family summary table
FamTable<-merge(FamStems,FamDiv,by='Family',all.x=TRUE,all.y=TRUE)
FamTable$Percent<-sprintf("%.3f",FamTable$Percent)
write.csv(FamTable,"Table_FamTable.csv",row.names=FALSE)

#Stem count by species, 2001-2012
aggAll<-with(veg,setNames(aggregate(Stems,by=list(Species,Y),sum),c('Species','Y','Stems')))
aggAll<-dcast(aggAll,Species~Y,value.var='Stems',sum)
firstthree<-(aggAll$`2001`+aggAll$`2002`+aggAll$`2003`)/3
lastthree<-(aggAll$`2010`+aggAll$`2011`+aggAll$`2012`)/3
aggAll$threeYrDelta<-sprintf("%.0f",((lastthree-firstthree)/firstthree)*100)
aggAll$threeYrDelta[aggAll$threeYrDelta=='NaN']<-0
write.csv(aggAll,'Table_GroundStemsSummary.csv',row.names=FALSE)

#Most abundance species in 2001
DomSp<-function(Hab,Yr){
  v<-subset(veg,idSpecLevel==1 & H==Hab & Y==Yr)
  v1<-setNames(aggregate(v$Q,by=list(v$Species),length),c('Species','Q'))
  v1<-v1[with(v1,order(-Q)),]
  print(head(v1))
  plot(v1$Q~seq(1:length(rownames(v1))))
}

DomSp('J',2001)
DomSp('W',2001)
DomSp('R',2001)

###############################################################################
# CALCULATE FLORISTIC QUALITY INDICES
###############################################################################

#Remove incomplete species identities
vegC<-veg[veg$idSpecLevel==1,]

#Subset to a single habitat & season (or not)
#Vs<-vegC[vegC$H=='J',]
Vs<-vegC

#Create a dataframe to store Species Richness results
FQA<-data.frame(
  Y=as.integer(), #Year
  S.obs=as.integer(), #No. species observed
  S.chao1=as.numeric(), #Chao1 estimate
  se.chao1=as.numeric(), #Chao1 standard error
  S.ACE=as.numeric(), #ACE estimate
  se.ACE=as.numeric(), #ACE standard error
  SW=as.numeric(), #Shannon-Weiner diversity
  FA=as.numeric(), #Fisher's alpha
  PJ=as.numeric(), #Pielou's J
  meanC=as.numeric(), #Mean C score (omitting non-natives)
  nn=as.numeric(), #Native species richness
  percentNative=as.numeric(), #Percent native species
  I=as.numeric(), #FQI (omitting non-natives)
  AI=as.numeric(), #Adjusted FQI
  Iy<-as.numeric(), #Cover-weighted FQI
  Checker=as.numeric() #Should equal zero, if not, species richness calculated wrong somewhere
)

#Pretty FQA
pFQA<-data.frame(
  Year=as.integer(), #Year
  ObservedSpecies=as.integer(), #No. species observed
  NativeSpecies=as.numeric(), #Native species richness
  percentNative=as.numeric(), #Percent native species
  Chao=as.numeric(), #Chao1 estimate
  ACE=as.numeric(), #ACE estimate
  Shannon=as.numeric(), #Shannon-Weiner diversity
  Fisher=as.numeric(), #Fisher's alpha
  Pielou=as.numeric(), #Pielou's J
  I=as.numeric(), #FQI (omitting non-natives)
  meanC=as.numeric(), #Mean C score (omitting non-natives)
  AI=as.numeric(), #Adjusted FQI
  Iy<-as.numeric() #Cover-weighted FQI
)


for(i in 1:11){ #Loop to estimate species richness for each year
  Y=Years[i] #Select year_i
  
  #Create a vegan matrix for year i
  matrixi<-with(Vs[Vs$Y==Y,],setNames(aggregate(Stems,by=list(Species),sum),c('Species','Abun')))
  
  #Calculate estimated & observed species richness
  estR<-estimateR(matrixi[,2]) #Estimate species richness
  S.obs<-estR[1] #Extract observed species richness
  S.chao1<-estR[2] #Extract Chao1 estimate
  se.chao1<-estR[3] #Extract Chao1 standard error
  S.ACE<-estR[4] #Extract ACE estimate
  se.ACE<-estR[5] #Extract ACE standard error
  
  #Some diversity functions
  SW<-diversity(matrixi[,2]) #Calculate Shannon-Weiner diversity
  FA<-fisher.alpha(matrixi[,2]) #Calculate Fisher's Alpha
  PJ<-SW/log(S.obs) #Calculate Pielou's J
  
  #Some FQA functions
  syi<-merge(matrixi,traits,by='Species',all.x=TRUE,all.y=FALSE)
  syi$C[syi$C==99]<-0 #Change non-native C to zero
  syi<-syi[!is.na(syi$C),] #remove NAs
  meanC<-mean(syi$C) #Calculate meanC
  nn<-length(syi$Species[syi$NT=='Nt']) #Calculate native species richness (nn)
  nt<-length(syi$Species) #Calculate total species richness (nt)
  percentNative<-(nn/nt)*100 #Calculate percent native species
  I<-meanC*sqrt(nn) #Calculate FQI
  AI<-100*(meanC/10)*(sqrt(nn)/sqrt(nt)) #Calculate adjusted FQI
  
  #Cover-weighted FQI
  library(plyr)
  Vs$Q10<-round_any(Vs$Q,10,ceiling)   #Create a variable combining 5 & 10 m quadrat observations
  castSprI<-acast(Vs[Vs$Y==Y,],Species~T+Q10,value.var='Cov',max) #Take the max value for each species to prevent double-counting between SP & SU
  castSprI[is.infinite(castSprI)]<-0 #Fill in zeros
  meanCov<-setNames(as.data.frame(rowMeans(castSprI)),'meanCov')
  meanCov$Species<-rownames(meanCov)
  meanCov<-merge(meanCov,traits,by='Species',all.x=TRUE,all.y=FALSE)
  meanCov$C[meanCov$C==99]<-0 #Make non-native C values 0
  meanCov$MCXC<-meanCov$meanCov*meanCov$C
  Cy<-sum(meanCov$MCXC)
  Iy<-Cy/length(meanCov$Species)  
  
  #Add a checker - species richness calculated two different ways should be the same
  Checker<-nt-S.obs
  
  #FQA Pretty
  chaoPretty<-paste(sprintf('%.1f',S.chao1),sprintf('%.1f',se.chao1),sep= " ± ")
  acePretty<-paste(sprintf('%.1f',S.ACE),sprintf('%.1f',se.ACE),sep= " ± ")
  fqaPretty<-as.data.frame(cbind(Year=Y,
                                 S.obs,nn,percentNative=sprintf('%.1f',percentNative),
                                 chaoPretty,acePretty,
                                 Shannon=sprintf('%.2f',SW),Fisher=sprintf('%.2f',FA),
                                 Pielou=sprintf('%.2f',PJ),
                                 I=sprintf('%.2f',I),meanC=sprintf('%.2f',meanC),
                                 AI=sprintf('%.2f',AI),Iy=sprintf('%.2f',Iy)))  
  
  #FQA Stat
  fqaStat<-as.data.frame(cbind(Y,
                               S.obs,nn,percentNative,
                               S.chao1,S.ACE,
                               SW,FA,PJ,
                               I,meanC,AI,Iy))
  
  #Add results to Species Richness output
  FQA<-rbind(FQA,fqaStat)
  pFQA<-rbind(pFQA,fqaPretty)
}


#Format output
FQA<-as.data.frame(t(FQA[,2:(length(colnames(FQA)))]))
colnames(FQA)<-Years

#Simple linear regression for each variable
Time<-c(0,1,2,3,4,6,7,8,9,10,11)

#Create a dataframe to store Species Richness results
LMResults<-data.frame(
  adjR<-as.character(),
  Beta=as.character(),
  t=as.character()
)

Pnone<-data.frame(pnone=as.numeric())

for(i in 1:length(rownames(FQA))){
  M<-summary(lm(as.numeric(FQA[i,1:11])~Time))
  adjR<-sprintf("%.2f",M$adj.r.squared)
  Beta<-paste(sprintf("%.3f",M$coefficients[2,1]),sprintf("%.3f",M$coefficients[2,2]),sep=" ± ")
  t<-sprintf("%.2f",M$coefficients[2,3])
  p<-as.numeric(M$coefficients[2,4])
  results<-cbind(adjR,Beta,t)
  LMResults<-rbind(LMResults,results)
  Pnone=rbind(Pnone,p)
}

PBon<-p.adjust(Pnone[,1],method='bonferroni')

LMResults$Padj<-paste(sprintf('%.4f',PBon),
                      ifelse(PBon<=0.0001,'****',
                             ifelse(PBon<=0.001,'***',
                                    ifelse(PBon<=0.01,'**',
                                           ifelse(PBon<=0.05,'*',
                                                  ifelse(PBon<=0.1,'.',''))))),sep="")

PrettyFQA<-t(pFQA[,2:length(colnames(pFQA))])
colnames(PrettyFQA)<-Years
PrettyFQA<-cbind(PrettyFQA,LMResults)

#Save results
#write.csv(PrettyFQA,'Table_FQA_Forest.csv',row.names=TRUE)
#write.csv(t(FQA),'Table_FQA_DBW_Raw.csv',row.names=TRUE)

###############################################################################
# ANALYZE FLORISTIC QUALITY DYNAMICS
###############################################################################

#Read FQA data from exported tables
FS<-read.csv("Table_FQA_DBW_Raw.csv")
  FS$Y<-as.numeric(FS$X)+0.5
FG<-read.csv("Table_FQA_Glade_Raw.csv")
  FG$Y<-as.numeric(FG$X)+0.5
FW<-read.csv("Table_FQA_Woodland_Raw.csv")
  FW$Y<-as.numeric(FW$X)+0.5
FF<-read.csv("Table_FQA_Forest_Raw.csv")
  FF$Y<-as.numeric(FF$X)+0.5

#Simple years for analysis
Ys<-c(1,2,3,4,5,7,8,9,10,11,12)
  
#Did species density increase linearly or logistically?
plot(FS[,3]~Ys)
LinMod<-lm(FS[,3]~Ys);summary(LinMod)
plot(resid(LinMod))

#Linear model for glade species density
plot(FG[,3]~Ys)
LinMod<-lm(FG[,3]~Ys);summary(LinMod)
abline(lm(FG[,3]~Ys))
plot(resid(LinMod))

#Linear model for woodland species density
plot(FW[,3]~Ys)
LinMod<-lm(FW[,3]~Ys);summary(LinMod)
abline(lm(FW[,3]~Ys))
plot(resid(LinMod))

#Log model for woodland species density
plot(FW[,3]~Ys)
LogMod<-lm(FW[,3]~log(Ys));summary(LogMod)
lines(Ys,predict(LogMod),col='red')

#Linear model for forest species density
plot(FF[,3]~Ys)
LinMod<-lm(FF[,3]~Ys);summary(LinMod)
abline(lm(FW[,3]~Ys))
plot(resid(LinMod))

#Linear model for Site-level I
plot(FS$I~Ys)
LinMod<-lm(FS$I~Ys);summary(LinMod)
LogMod<-lm(FS$I~log(Ys));summary(LogMod)
lines(Ys,predict(LinMod),col='red')

#Linear model for glade I
plot(FW$I~Ys)
LinMod<-lm(FW$I~Ys);summary(LinMod)
LogMod<-lm(FW$I~log(Ys));summary(LogMod)

#model for DBW Iy
plot(FS$Iy~Ys)
LinMod<-lm(FS$Iy~Ys);summary(LinMod)
LogMod<-lm(FS$Iy~log(Ys));summary(LogMod)
lines(Ys,predict(LogMod),col='red')

#Model for glade Iy
plot(FG$Iy~Ys)
LinMod<-lm(FG$Iy~Ys);summary(LinMod)
LogMod<-lm(FG$Iy~log(Ys));summary(LogMod)
lines(Ys,predict(LinMod),col='red')

#Model for woodland Iy
plot(FW$Iy~Ys)
LinMod<-lm(FW$Iy~Ys);summary(LinMod)
LogMod<-lm(FW$Iy~log(Ys));summary(LogMod)
lines(Ys,predict(LinMod),col='red')

#Model for forest Iy
plot(FF$Iy~Ys)
LogMod<-lm(FF$Iy~log(Ys));summary(LogMod)
lines(Ys,predict(LogMod),col='red')

#Broken stick model for forest Iy
library(segmented)
FIY<-FF$Iy
LinMod2<-lm(FIY~Ys)
SegMod<-segmented(LinMod2,seg.Z=~Ys);summary(SegMod)
lines(Ys,predict(SegMod),col='red')

#Graphics function for lines and shaded area
fqaGraphics<-function(x){
  abline(v=2001+(328/365),lty=2,col='gray85') #prescribed burns
  abline(v=2004+(57/365),lty=2,col='gray85')
  abline(v=2007+(74/365),lty=2,col='gray85')
  abline(v=2009+(319/365),lty=2,col='gray85')
  abline(v=2012+(10/365),lty=2,col='gray85')
  polygon(c(2005.8,2005.8,2007,2007),c(-1000,3000,3000,-1000),col='gray85',border=NA) #Shade gray - shrubs
  box()
  par(xpd=NA)
  text(2002,x,pos=3,offset=1.3,'fire',srt=45,col='gray75',adj=0)
  text(2004,x,pos=3,offset=1.3,'fire',srt=45,col='gray75',adj=0)
  text(2007.7,x,pos=3,offset=1.3,'fire',srt=45,col='gray75',adj=0)
  text(2010,x,pos=3,offset=1.3,'fire',srt=45,col='gray75',adj=0)
  text(2012,x,pos=3,offset=1.3,'fire',srt=45,col='gray75',adj=0)
  text(2006.5,x,pos=3,offset=1.7,'juniper\nremoval',srt=45,col='gray75',adj=0)
  par(new=TRUE)
  par(xpd=FALSE)
}

#Axis function to prevent overplotting
axisFQA<-function(x){
  axis(side=1,at=seq(2002,2010,4))
  axis(side=1,at=seq(2004,2012,4))
}

DBWpalCOL<-c('darkgreen','black','maroon')
DBWpalBG<-c('yellow','white','pink')

#FQA Scatterplot
tiff('Fig_FQA_Scatter_Site.tif',res=300,height=80,width=180,units='mm')
par(mfrow=c(1,3),mai=c(1,0.5,0.5,0.1),omi=c(0,0,0,0),mgp=c(2.5,0.7,0))

plot(FS$S.obs~FS$Y,type='n',ylim=c(150,210),las=1,ylab='Obs. native species density',xlab=NA,xaxt='n')
fqaGraphics(210)
points(FS$nn~FS$Y,pch=16)
abline(lm(FS$nn~FS$Y))
axisFQA()

par(mgp=c(2,0.7,0))
plot(FS$I~FS$Y,type='n',las=1,ylim=c(49,56),ylab=expression(paste('Floristic quality (',italic(I),')')),xlab=NA,xaxt='n')
fqaGraphics(56)
points(FS$I~FS$Y,pch=16)
abline(lm(FS$I~FS$Y))
axisFQA()

plot(FS$Iy~FS$Y,type='n',las=1,
     ylab=expression(paste('Cover-weighted floristic quality ( ',italic(I[gamma]),')')),
     xlab=NA,ylim=c(1,2.6),xaxt='n')
fqaGraphics(2.6)
points(FS$Iy~FS$Y,pch=16)
lines(Ys+2000.5,predict(lm(FS$Iy~log(Ys))))
axisFQA()

dev.off()

#FQA Scatterplot - Community level
tiff('Fig_FQA_Scatter_Comm.tif',res=300,height=80,width=180,units='mm')
par(mfrow=c(1,3),mai=c(1,0.5,0.5,0.1),omi=c(0,0,0,0),mgp=c(2.5,0.7,0))

plot(FS$S.obs~FS$Y,type='n',ylim=c(90,150),las=1,ylab='Obs. native species density',xlab=NA,xaxt='n')
fqaGraphics(150)
lines(FF$nn~FF$Y,pch=22,col=DBWpalCOL[3],bg=DBWpalBG[3],type='b',lty=2)
points(FG$nn~FG$Y,pch=21,bg=DBWpalBG[1],col=DBWpalCOL[1])
abline(lm(FG$nn~FG$Y),col=DBWpalCOL[1])
points(FW$nn~FW$Y,pch=24,bg=DBWpalBG[2],col=DBWpalCOL[2])
lines(Ys+2000.5,predict(lm(FW$nn~log(Ys))),lty=3,col=DBWpalCOL[2])
axisFQA()
legend('topleft',legend=c('Juniper','Woodland','Forest'),
       pch=c(21,24,22),pt.bg=DBWpalBG,col=DBWpalCOL,lty=c(1,3,2),bty='n')

plot(FS$I~FS$Y,type='n',las=1,ylim=c(35,45),ylab=expression(paste('Floristic quality (',italic(I),')')),xlab=NA,xaxt='n')
fqaGraphics(45)
lines(FW$I~FW$Y,pch=24,bg=DBWpalBG[2],col=DBWpalCOL[2],type='b',lty=3)
lines(FF$I~FF$Y,pch=22,bg=DBWpalBG[3],col=DBWpalCOL[3],type='b',lty=2)
lines(FG$I~FG$Y,pch=21,bg=DBWpalBG[1],col=DBWpalCOL[1],type='b')
axisFQA()

plot(FS$Iy~FS$Y,type='n',las=1,
     ylab=expression(paste('Cover-weighted floristic quality ( ',italic(I[gamma]),')')),xlab=NA,ylim=c(1,5),xaxt='n')
fqaGraphics(5)
points(FG$Iy~FG$Y,pch=21,bg=DBWpalBG[1],col=DBWpalCOL[1])
abline(lm(FG$Iy~FG$Y),col=DBWpalCOL[1])
points(FF$Iy~FF$Y,pch=22,bg=DBWpalBG[3],col=DBWpalCOL[3])
lines(FF$Y,predict(segmented(lm(FIY~Ys))),col=DBWpalCOL[3],lty=2)
points(FW$Iy~FW$Y,pch=24,bg=DBWpalBG[2],col=DBWpalCOL[2])
abline(lm(FW$Iy~FW$Y),lty=3)
axisFQA()

dev.off()
###############################################################################
#POLYGON PLOT - GUILD S + Cover ~ TIME
###############################################################################

#Unique habitat list
Habitats<-c('J','W','R')

#Subset to drop incomplete IDs
v<-droplevels(subset(veg,idSpecLevel==1))

#Calculate species richness for each guild per habitat per year
GuildRich<-function(x){
vx<-droplevels(subset(v,H==Habitats[x]))
vx2<-setNames(aggregate(vx$Species,by=list(vx$Y,vx$PHYS),function(x)length(unique(x))),c('Y','Guild','Species'))
vx3<-dcast(vx2,Y~Guild)
vx3$S<-rowSums(vx3[,2:length(colnames(vx3))])
assign('vx3',vx3,envir=.GlobalEnv)
}

#Draw polygons for species richness per guild per year
GuildPoly<-function(vx3){
polygon(c(vx3$Y,rev(vx3$Y)),c(rep(0,11),rev(vx3$S)),col='blue') #Forbs
polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$FORB,rev(vx3$S)),col='pink') #TREE
polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$FORB+vx3$TREE,rev(vx3$S)),col='green') #Grass
polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$FORB+vx3$TREE+vx3$GRASS,rev(vx3$S)),col='seagreen') #Shrub
polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$FORB+vx3$TREE+vx3$GRASS+vx3$SHRUB,rev(vx3$S)),col='yellow') #Sedge
polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$FORB+vx3$TREE+vx3$GRASS+vx3$SHRUB+vx3$SEDGE,rev(vx3$S)),col='gray80') #Woody vine
polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$FORB+vx3$TREE+vx3$GRASS+vx3$SHRUB+vx3$SEDGE+vx3$WVINE,rev(vx3$S)),col='red') #Fern
}

#Calculate mean percent cover per 0.5 m2 for each guild per habitat per year (taking the max of 5 & 10 quads)
GuildCov<-function(x){
  vx<-droplevels(subset(v,H==Habitats[x]))
  vx$Q10<-paste(vx$T,round_any(vx$Q,10,ceiling))
  vx2<-setNames(aggregate(vx$Cov,by=list(vx$Y,vx$PHYS,vx$Q10),max),c('Y','PHYS','Q10','Cov'))
  vx3<-dcast(vx2,Y+PHYS~Q10)
  vx3[is.na(vx3)]<-0
  vx3$Cov<-rowMeans(vx3[,3:length(colnames(vx3))])
  vx4<-dcast(vx3,Y~PHYS,value.var='Cov')
  vx4$Cov<-rowSums(vx4[,2:length(colnames(vx4))])
  assign('vx4',vx4,envir=.GlobalEnv)
}

#Draw polygons for mean percent cover per guild per year
GuildCovPoly<-function(vx4){
  polygon(c(vx4$Y,rev(vx4$Y)),c(rep(0,11),rev(vx4$Cov)),col='blue') #Forbs
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$FORB,rev(vx4$Cov)),col='pink') #TREE
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$FORB+vx4$TREE,rev(vx4$Cov)),col='green') #Grass
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$FORB+vx4$TREE+vx4$GRASS,rev(vx4$Cov)),col='seagreen') #Shrub
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$FORB+vx4$TREE+vx4$GRASS+vx4$SHRUB,rev(vx4$Cov)),col='yellow') #Sedge
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$FORB+vx4$TREE+vx4$GRASS+vx4$SHRUB+vx4$SEDGE,rev(vx4$Cov)),col='gray80') #Woody vine
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$FORB+vx4$TREE+vx4$GRASS+vx4$SHRUB+vx4$SEDGE+vx4$WVINE,rev(vx4$Cov)),col='red') #Fern
}

#Analyze guild richness by treatment
GuildRich(1)
vx3$ForbPC<-with(vx3,FORB/S)
with(vx3,plot(ForbPC~Ys))
with(vx3,summary(lm(ForbPC~Ys)))
vx3$BA<-c(rep('B',5),rep('A',6))
with(vx3,summary(lm(ForbPC~BA)))

#Assign C categories
v$catC<-ifelse(v$C<=3,'Ruderal',
               ifelse(v$C<=6,'Matrix',
                      ifelse(v$C<=10,'Conservative','Non-native')))

#Calculate species richness for each C per habitat per year
cRich<-function(x){
  vx<-droplevels(subset(v,H==Habitats[x]))
  vx2<-setNames(aggregate(vx$Species,by=list(vx$Y,vx$catC),function(x)length(unique(x))),c('Y','catC','Species'))
  vx3<-dcast(vx2,Y~catC)
  vx3$S<-rowSums(vx3[,2:length(colnames(vx3))])
  assign('vx3',vx3,envir=.GlobalEnv)
}


cols<-c('black','gray40','gray80','white')

#Draw polygons for species richness per C per year
cPoly<-function(vx3){
  polygon(c(vx3$Y,rev(vx3$Y)),c(rep(0,11),rev(vx3$S)),col=cols[1]) #Non-native
  polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$`Non-native`,rev(vx3$S)),col=cols[2]) #Ruderal
  polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$`Non-native`+vx3$Ruderal,rev(vx3$S)),col=cols[3]) #Matrix
  polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$`Non-native`+vx3$Ruderal+vx3$Matrix,rev(vx3$S)),col=cols[4]) #Conservative
}

#Calculate mean percent cover per 0.5 m2 for each C per habitat per year (taking the max of 5 & 10 quads)
cCov<-function(x){
  vx<-droplevels(subset(v,H==Habitats[x]))
  vx$Q10<-paste(vx$T,round_any(vx$Q,10,ceiling))
  vx2<-setNames(aggregate(vx$Cov,by=list(vx$Y,vx$catC,vx$Q10),max),c('Y','catC','Q10','Cov'))
  vx3<-dcast(vx2,Y+catC~Q10)
  vx3[is.na(vx3)]<-0
  vx3$Cov<-rowMeans(vx3[,3:length(colnames(vx3))])
  vx4<-dcast(vx3,Y~catC,value.var='Cov')
  vx4$Cov<-rowSums(vx4[,2:length(colnames(vx4))])
  assign('vx4',vx4,envir=.GlobalEnv)
}

#Draw polygons for mean percent cover per C per year
cCovPoly<-function(vx4){
  polygon(c(vx4$Y,rev(vx4$Y)),c(rep(0,11),rev(vx4$Cov)),col=cols[1]) #Non-native
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$`Non-native`,rev(vx4$Cov)),col=cols[2]) #Ruderal
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$`Non-native`+vx4$Ruderal,rev(vx4$Cov)),col=cols[3]) #Matrix
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$`Non-native`+vx4$Ruderal+vx4$Matrix,rev(vx4$Cov)),col=cols[4]) #Conservative
}

#Plot change in mean cover by (A) physiognomy and (B) conservatism
tiff('Fig_PolygonsMeanCov.tif',res=300,height=140,width=180,units='mm')
par(mfrow=c(2,3),mai=c(0.1,0.1,0,0),omi=c(0.8,0.5,0.5,1.0),mgp=c(2.5,0.7,0))

GuildCov(1)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n')
mtext(side=2,outer=TRUE,'Mean cover (%)',line=2)
fqaGraphics(100)
GuildCovPoly(vx4)
legend('bottomleft','Juniper',text.col='white',bty='n')

GuildCov(2)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n',ylab='',yaxt='n')
fqaGraphics(100)
GuildCovPoly(vx4)
legend('bottomleft','Woodland',text.col='white',bty='n')

GuildCov(3)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n',ylab='',yaxt='n')
fqaGraphics(100)
GuildCovPoly(vx4)
legend('bottomleft','Forest',text.col='white',bty='n')

par(xpd=NA)
legend(2012.3,100,legend=c('Forb','Tree','Grass','Shrub','Sedge','Woody vine','Fern'),
       fill=c('blue','pink','green','seagreen','yellow','gray80','red'),bty='n',title='Physiognomy')
par(xpd=FALSE)

cCov(1)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',ylab='Mean cover (%)',xaxt='n')
fqaGraphics(1000)
cCovPoly(vx4)
axisFQA()
legend('bottomleft','Juniper',text.col='white',bty='n')

cCov(2)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n',ylab='',yaxt='n')
fqaGraphics(1000)
cCovPoly(vx4)
axisFQA()
legend('bottomleft','Woodland',text.col='white',bty='n')

cCov(3)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n',ylab='',yaxt='n')
fqaGraphics(1000)
cCovPoly(vx4)
axisFQA()
legend('bottomleft','Forest',text.col='white',bty='n')

par(xpd=NA)
legend(2012.3,100,legend=c('Non-native','Ruderal','Matrix','Conservative'),fill=c(cols[1:4]),
       title="Conservatism",bty='n')
par(xpd=FALSE)

dev.off()

#Plot change in species density by (A) physiognomy and (B) conservatism
tiff('Fig_PolygonsSpecDens.tif',res=300,height=140,width=180,units='mm')
par(mfrow=c(2,3),mai=c(0.1,0.1,0,0),omi=c(0.8,0.5,0.5,1.0),mgp=c(2.5,0.7,0))

GuildRich(1)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',xaxt='n',ylab='')
fqaGraphics(160)
GuildPoly(vx3)
legend('bottomleft','Juniper',text.col='white',bty='n')

GuildRich(2)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n',yaxt='n')
fqaGraphics(160)
GuildPoly(vx3)
legend('bottomleft','Woodland',text.col='white',bty='n')

GuildRich(3)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n',yaxt='n')
fqaGraphics(160)
GuildPoly(vx3)
legend('bottomleft','Forest',text.col='white',bty='n')

par(xpd=NA)
legend(2012.3,160,legend=c('Forb','Tree','Grass','Shrub','Sedge','Woody vine','Fern'),
       fill=c('blue','pink','green','seagreen','yellow','gray80','red'),bty='n',title='Physiognomy')
par(xpd=FALSE)

cRich(1)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n')
fqaGraphics(1600)
cPoly(vx3)
axisFQA()
legend('bottomleft','Juniper',text.col='white',bty='n')

cRich(2)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n',yaxt='n')
fqaGraphics(1600)
cPoly(vx3)
axisFQA()
legend('bottomleft','Woodland',text.col='white',bty='n')

cRich(3)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n',yaxt='n')
fqaGraphics(1600)
cPoly(vx3)
axisFQA()
legend('bottomleft','Forest',text.col='white',bty='n')

par(xpd=NA)
legend(2012.3,160,legend=c('Non-native','Ruderal','Matrix','Conservative'),fill=c(cols[1:4]),
       title="Conservatism",bty='n')
par(xpd=FALSE)

mtext(side=2,outer=TRUE,'Species density',line=2)

dev.off()

###############################################################################
#POLYGON PLOT - C ~ TIME
###############################################################################

#Subset to drop incomplete IDs
v<-droplevels(subset(veg,idSpecLevel==1))

#Assign C categories
v$catC<-ifelse(v$C<=3,'Ruderal',
               ifelse(v$C<=6,'Matrix',
                      ifelse(v$C<=10,'Conservative','Non-native')))

#Calculate species richness for each C per habitat per year
cRich<-function(x){
  vx<-droplevels(subset(v,H==Habitats[x]))
  vx2<-setNames(aggregate(vx$Species,by=list(vx$Y,vx$catC),function(x)length(unique(x))),c('Y','catC','Species'))
  vx3<-dcast(vx2,Y~catC)
  vx3$S<-rowSums(vx3[,2:length(colnames(vx3))])
  assign('vx3',vx3,envir=.GlobalEnv)
}


cols<-c('black','gray40','gray80','white')

#Draw polygons for species richness per C per year
cPoly<-function(vx3){
  polygon(c(vx3$Y,rev(vx3$Y)),c(rep(0,11),rev(vx3$S)),col=cols[1]) #Non-native
  polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$`Non-native`,rev(vx3$S)),col=cols[2]) #Ruderal
  polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$`Non-native`+vx3$Ruderal,rev(vx3$S)),col=cols[3]) #Matrix
  polygon(c(vx3$Y,rev(vx3$Y)),c(vx3$`Non-native`+vx3$Ruderal+vx3$Matrix,rev(vx3$S)),col=cols[4]) #Conservative
}

#Calculate mean percent cover per 0.5 m2 for each C per habitat per year (taking the max of 5 & 10 quads)
cCov<-function(x){
  vx<-droplevels(subset(v,H==Habitats[x]))
  vx$Q10<-paste(vx$T,round_any(vx$Q,10,ceiling))
  vx2<-setNames(aggregate(vx$Cov,by=list(vx$Y,vx$catC,vx$Q10),max),c('Y','catC','Q10','Cov'))
  vx3<-dcast(vx2,Y+catC~Q10)
  vx3[is.na(vx3)]<-0
  vx3$Cov<-rowMeans(vx3[,3:length(colnames(vx3))])
  vx4<-dcast(vx3,Y~catC,value.var='Cov')
  vx4$Cov<-rowSums(vx4[,2:length(colnames(vx4))])
  assign('vx4',vx4,envir=.GlobalEnv)
}

#Draw polygons for mean percent cover per C per year
cCovPoly<-function(vx4){
  polygon(c(vx4$Y,rev(vx4$Y)),c(rep(0,11),rev(vx4$Cov)),col=cols[1]) #Non-native
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$`Non-native`,rev(vx4$Cov)),col=cols[2]) #Ruderal
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$`Non-native`+vx4$Ruderal,rev(vx4$Cov)),col=cols[3]) #Matrix
  polygon(c(vx4$Y,rev(vx4$Y)),c(vx4$`Non-native`+vx4$Ruderal+vx4$Matrix,rev(vx4$Cov)),col=cols[4]) #Conservative
}

#Plot C frequency of S & Cover by habitat & year
tiff('Fig_CPoly.tif',res=300,height=160,width=180,units='mm')
par(mfrow=c(2,3),mai=c(0.3,0.5,0.5,0),omi=c(0.8,0,0,0.2),mgp=c(2.5,0.7,0))

cRich(1)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='Obs. species density',xaxt='n')
fqaGraphics(160)
cPoly(vx3)
axisFQA()
legend('bottomleft','Glade',text.col='white',bty='n')

cRich(3)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n')
fqaGraphics(160)
cPoly(vx3)
axisFQA()
legend('bottomleft','Woodland',text.col='white',bty='n')

cRich(2)
plot(vx3$S~vx3$Y,type='l',ylim=c(0,160),las=1,xlab='',ylab='',xaxt='n')
fqaGraphics(160)
cPoly(vx3)
axisFQA()
legend('bottomleft','Forest',text.col='white',bty='n')

cCov(1)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',ylab='Mean cover (%)',xaxt='n')
fqaGraphics(100)
cCovPoly(vx4)
axisFQA()
legend('bottomleft','Glade',text.col='white',bty='n')

cCov(3)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n',ylab='')
fqaGraphics(100)
cCovPoly(vx4)
axisFQA()
legend('bottomleft','Woodland',text.col='white',bty='n')

par(xpd=NA)
legend(1990,-30,ncol=4,legend=c('Non-native','Ruderal species','Matrix species','Conservative species'),
       fill=c(cols[1:4]))
par(xpd=FALSE)

cCov(2)
plot(vx4$Cov~vx4$Y,type='l',ylim=c(0,100),las=1,xlab='',xaxt='n',ylab='')
fqaGraphics(100)
cCovPoly(vx4)
axisFQA()
legend('bottomleft','Forest',text.col='white',bty='n')

dev.off()

###############################################################################
# NMDS
###############################################################################

#Pairwise adonis Function - Paula Martins/Robert Danczak rendition of original function by Pedro Martinez Arbizu
pairwise.adonis <- function(x,factors, sim.method,binary=FALSE, p.adjust.m)
{
  library(vegan)
  co = as.matrix(combn(unique(factors),2))
  pairs = c()
  F.Model =c()
  R2 = c()
  p.value = c()
  
  for(elem in 1:ncol(co)){
    ad = adonis(x[factors %in% c(as.character(co[1,elem]),as.character(co[2,elem])),] ~ 
                  factors[factors %in% c(as.character(co[1,elem]),as.character(co[2,elem]))] , method =sim.method,binary=FALSE);
    pairs = c(pairs,paste(co[1,elem],'vs',co[2,elem]));
    F.Model =c(F.Model,ad$aov.tab[1,4]);
    R2 = c(R2,ad$aov.tab[1,5]);
    p.value = c(p.value,ad$aov.tab[1,6])
  }
  p.adjusted = p.adjust(p.value,method=p.adjust.m)
  pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted)
  return(pairw.res)
}

#Subset veg to included only complete species identifications
V<-droplevels(subset(veg,idSpecLevel==1))

#Add a Year_Transect variable
V$YT<-with(V,paste(Y,T,sep='_'))

#Cast into a vegan-style matrix, calculating average percent cover per quadrat
VegMat<-acast(V,YT~ACRON,value.var='Cov',function(x)round_any(sum(x)/10,accuracy=1,f=ceiling))

#Remove species that account for <5% of records (i.e., those present in <16.5 transect years)
SpecCount<-apply(VegMat,2,function(x)length(x[x>0]))
VegMat<-VegMat[,SpecCount>17]

#Environmental dataframe
env<-setNames(as.data.frame(rownames(VegMat)),'YT')
env$H<-substr(env$YT,6,6)
env$HY<-paste(env$H,substr(env$YT,1,4),sep="_")

#Abbreviate species names for plotting
#colnames(DM)<-make.cepnames(colnames(DM)) #No need to do this since we already have abbreviations

#Relativize by species maximum percent cover
#VegMat<-apply(VegMat,2,function(x)x/max(x)) #Better fit w/o relativization

#Set species plotting prioritization
priSpp<-colSums(VegMat)

#Run NMDS
MD1<-metaMDS(VegMat,distance='bray',k=3,trymax=1000)

#Check stress
stressplot(MD1)

#PERMANOVA
adonis(VegMat~HY,data=env,method='bray',binary=FALSE,perm=999) #Regular permanova

#Pairwise PERMANOVA
ad1<-pairwise.adonis(x=VegMat,factors=env$HY,sim.method='bray',binary=FALSE,p.adjust.m='bonferroni');ad1 #Pairwise permanova

#PERMDISP (test for dispersion)
bd1<-betadisper(vegdist(VegMat,method='bray'),group=env$HY,type='centroid',bias.adjust=TRUE)
anova(bd1)
tk1<-as.data.frame(TukeyHSD(bd1)$group);tk1
tk1$pairsForward<-gsub("-"," vs ",rownames(tk1)) #Add a merging column
#Add a reversed merging column (for some reason, the pairs are backwards in PERMDISP vs. PERMANOVA)
tk1$pairsReversed<-paste(substr(tk1$pairsForward,11,16),substr(tk1$pairsForward,1,6),sep=' vs ')

#Merge PERMANOVA & PERMDISP results
permResults<-merge(ad1,tk1,by.x='pairs',by.y='pairsReversed',all.x=TRUE)

#Pairs I want
PairsIWant<-c('J_2001 vs R_2001',
              'J_2001 vs W_2001',
              'R_2001 vs W_2001',
              'J_2012 vs R_2012',
              'J_2012 vs W_2012',
              'R_2012 vs W_2012',
              'J_2001 vs J_2012',
              'R_2001 vs R_2012',
              'W_2001 vs W_2012')

#Subset to include only selected pairs
permResults<-permResults[permResults$pairs %in% PairsIWant,]

#Format PERMANOVA/PERMDISP results table
permResults<-subset(permResults,select=c('pairs','F.Model','R2','p.value','diff','lwr','upr','p adj'))
permResults$diff<-permResults$diff*-1
permResults$lwr<-permResults$lwr*-1
permResults$upr<-permResults$upr*-1
permResults$L<-pmin(permResults$lwr,permResults$upr)
permResults$U<-pmax(permResults$lwr,permResults$upr)
permResults$diff<-with(permResults,paste(sprintf("%.3f",diff)," (",sprintf("%.3f",L),"-",sprintf("%.3f",U),")",sep=""))
permResults<-subset(permResults,select=c('pairs','F.Model','R2','p.value','diff','p adj'))
colnames(permResults)<-c('Contrast','F','R2','P1','Dispersion','P2')
permResults$F<-round(permResults$F,2)
permResults$R2<-round(permResults$R2,2)
permResults$P1<-round(permResults$P1,4)
permResults$P2<-round(permResults$P2,4)
permResults

#Save pairwise PERMANOVA table
write.csv(permResults,'Table_GroundPERMANOVA.csv',row.names=FALSE)

#Export species scores
SpecScores<-as.data.frame(MD1$species)
SpecScores$ACRON<-rownames(SpecScores)
SpecScores<-merge(SpecScores,traits,all.x=TRUE,all.y=FALSE)
SpecScores$PCH<-with(SpecScores,
                     ifelse(PHYS=='FORB',21,
                            ifelse(PHYS=='FERN',21,
                                   ifelse(PHYS=='SEDGE',21,
                                          ifelse(PHYS=='GRASS',21,25)))))
SpecScores$BG<-with(SpecScores,
                     ifelse(NT=='Ad','black',
                            ifelse(C<=3,'gray40',
                                   ifelse(C>=7,'white','gray80'))))

SpecScores$COL<-with(SpecScores,
                     ifelse(PHYS=='FORB','blue',
                            ifelse(PHYS=='FERN','red',
                                   ifelse(PHYS=='SEDGE','yellow',
                                          ifelse(PHYS=='GRASS','green',
                                                 ifelse(PHYS=='TREE','pink',
                                                        ifelse(PHYS=='SHRUB','seagreen','gray80')))))))

#write.csv(SpecScores,"Table_NMDSGroundSpecScores.csv",row.names=TRUE)

#Extract site scores for plotting
MDsite<-as.data.frame(MD1$points)
MDsite$HY<-paste(substr(rownames(MDsite),6,6),substr(rownames(MDsite),3,4),sep='')

#Calculate site centroids for each habitat_x_year
MDcen<-setNames(aggregate(list(MDsite$MDS1,MDsite$MDS2,MDsite$MDS3),by=list(MDsite$HY),mean),c('HY','MDS1','MDS2','MDS3'))
MDcen$H<-substr(MDcen$HY,1,1)
MDcen$Y<-substr(MDcen$HY,2,3)

#Set plot controls
plotCols<-c('darkgreen','red','black')

#Plot NMDS result
tiff('Fig_GroundNMDSBray_Pub_v23.tif',res=300,height=80,width=160,units='mm')
par(mfrow=c(1,2),mai=c(0,0.4,0.3,0),omi=c(1,0.3,0,1.2),mgp=c(1.5,0.3,0),
    tcl=-0.2,las=1,cex.axis=0.8)
plot(MD1,choices=c(1,2),type='n',ylab="",xlab="",xlim=c(-1.1,1.4))
#abline(h=0,lty=2,col='gray70')
#abline(v=0,lty=2,col='gray70')
with(MDcen[MDcen$Y=='01',],arrows(MDS1[1],MDS2[1],MDS1[3],MDS2[3],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='01',],arrows(MDS1[3],MDS2[3],MDS1[2],MDS2[2],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='01',],arrows(MDS1[2],MDS2[2],1.0,MDS2[2],length=0,lty=3,col='gray50'))
text(1.0,MDcen[MDcen$Y=='01',]$MDS2[2],'2001',adj=0,col='gray50',cex=0.8)
with(MDcen[MDcen$Y=='05',],arrows(MDS1[1],MDS2[1],MDS1[3],MDS2[3],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='05',],arrows(MDS1[3],MDS2[3],MDS1[2],MDS2[2],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='05',],arrows(MDS1[2],MDS2[2],1.0,MDS2[2],length=0,lty=3,col='gray50'))
text(1.0,MDcen[MDcen$Y=='05',]$MDS2[2],'2005',adj=0,col='gray50',cex=0.8)
with(MDcen[MDcen$Y=='12',],arrows(MDS1[1],MDS2[1],MDS1[3],MDS2[3],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='12',],arrows(MDS1[3],MDS2[3],MDS1[2],MDS2[2],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='12',],arrows(MDS1[2],MDS2[2],1.0,MDS2[2],length=0,lty=3,col='gray50'))
text(1.0,MDcen[MDcen$Y=='12',]$MDS2[2],'2012',adj=0,col='gray50',cex=0.8)
for(i in 1:11){with(MDcen[MDcen$H=='J',],Arrows(MDS1[i],MDS2[i],MDS1[i+1],MDS2[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[1]))}
for(i in 1:11){with(MDcen[MDcen$H=='J',],points(MDS2[i]~MDS1[i],pch=21,bg=DBWpalBG[1],col=DBWpalCOL[1]))}
for(i in 1:11){with(MDcen[MDcen$H=='W',],Arrows(MDS1[i],MDS2[i],MDS1[i+1],MDS2[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[2]))}
for(i in 1:11){with(MDcen[MDcen$H=='W',],points(MDS2[i]~MDS1[i],pch=24,bg=DBWpalBG[2],col=DBWpalCOL[2]))}
for(i in 1:11){with(MDcen[MDcen$H=='R',],Arrows(MDS1[i],MDS2[i],MDS1[i+1],MDS2[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[3]))}
for(i in 1:11){with(MDcen[MDcen$H=='R',],points(MDS2[i]~MDS1[i],pch=22,bg=DBWpalBG[3],col=DBWpalCOL[3]))}
text(-0.4,0.4,'Juniper',adj=1,cex=0.8,col=DBWpalCOL[1])
text(0.5,-0.45,'Forest',adj=0,cex=0.8,col=DBWpalCOL[3])
text(-0.2,-0.6,'Woodland',adj=0,cex=0.8,col=DBWpalCOL[2])
mtext(side=3,adj=0.03,line=-1.2,'A')
mtext(side=2,line=2,'NMDS Axis 2',las=3,cex=0.8)
mtext(side=1,line=2,'NMDS Axis 1',cex=0.8)

plot(MD1,choices=c(1,2),type='n',ylab="",xlab="",xlim=c(-1.1,1.4))
#abline(h=0,lty=2,col='gray70')
#abline(v=0,lty=2,col='gray70')
with(SpecScores,points(MDS2~MDS1,pch=PCH,bg=BG,col='black'))
mtext(side=3,adj=0.03,line=-1.2,'B')
mtext(side=1,line=2,'NMDS Axis 1',cex=0.8)

par(xpd=NA)
legend(1.5,1.2,legend=c('Woody','Herbaceous'),pch=c(25,21),bty='n',title='Physiognamy',
       cex=0.8,pt.cex=1)
legend(1.5,0.1,legend=c('Non-native','Ruderal','Matrix','Conservative'),
       fill=c('black','gray40','gray80','white'),bty='n',title='Conservatism',
       cex=0.8,pt.cex=1)
par(xpd=FALSE)

dev.off()


#Plot NMDS result - full figure for supplement
tiff('Fig_GroundNMDSBray_SM_v23.tif',res=300,height=150,width=160,units='mm')
par(mfrow=c(2,2),mai=c(0,0.4,0.3,0),omi=c(1,0.3,0,1.2),mgp=c(1.5,0.3,0),
    tcl=-0.2,las=1,cex.axis=0.8)

plot(MD1,choices=c(1,2),type='n',ylab="",xlab="",xlim=c(-1.1,1.4))
with(MDcen[MDcen$Y=='01',],arrows(MDS1[1],MDS2[1],MDS1[3],MDS2[3],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='01',],arrows(MDS1[3],MDS2[3],MDS1[2],MDS2[2],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='01',],arrows(MDS1[2],MDS2[2],1.0,MDS2[2],length=0,lty=3,col='gray50'))
text(1.0,MDcen[MDcen$Y=='01',]$MDS2[2],'2001',adj=0,col='gray50',cex=0.8)
with(MDcen[MDcen$Y=='05',],arrows(MDS1[1],MDS2[1],MDS1[3],MDS2[3],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='05',],arrows(MDS1[3],MDS2[3],MDS1[2],MDS2[2],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='05',],arrows(MDS1[2],MDS2[2],1.0,MDS2[2],length=0,lty=3,col='gray50'))
text(1.0,MDcen[MDcen$Y=='05',]$MDS2[2],'2005',adj=0,col='gray50',cex=0.8)
with(MDcen[MDcen$Y=='12',],arrows(MDS1[1],MDS2[1],MDS1[3],MDS2[3],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='12',],arrows(MDS1[3],MDS2[3],MDS1[2],MDS2[2],length=0,lty=3,col='gray50'))
with(MDcen[MDcen$Y=='12',],arrows(MDS1[2],MDS2[2],1.0,MDS2[2],length=0,lty=3,col='gray50'))
text(1.0,MDcen[MDcen$Y=='12',]$MDS2[2],'2012',adj=0,col='gray50',cex=0.8)
for(i in 1:11){with(MDcen[MDcen$H=='J',],Arrows(MDS1[i],MDS2[i],MDS1[i+1],MDS2[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[1]))}
for(i in 1:11){with(MDcen[MDcen$H=='J',],points(MDS2[i]~MDS1[i],pch=21,bg=DBWpalBG[1],col=DBWpalCOL[1]))}
for(i in 1:11){with(MDcen[MDcen$H=='W',],Arrows(MDS1[i],MDS2[i],MDS1[i+1],MDS2[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[2]))}
for(i in 1:11){with(MDcen[MDcen$H=='W',],points(MDS2[i]~MDS1[i],pch=24,bg=DBWpalBG[2],col=DBWpalCOL[2]))}
for(i in 1:11){with(MDcen[MDcen$H=='R',],Arrows(MDS1[i],MDS2[i],MDS1[i+1],MDS2[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[3]))}
for(i in 1:11){with(MDcen[MDcen$H=='R',],points(MDS2[i]~MDS1[i],pch=22,bg=DBWpalBG[3],col=DBWpalCOL[3]))}
text(-0.4,0.4,'Juniper',adj=1,cex=0.8,col=DBWpalCOL[1])
text(0.5,-0.45,'Forest',adj=0,cex=0.8,col=DBWpalCOL[3])
text(-0.2,-0.6,'Woodland',adj=0,cex=0.8,col=DBWpalCOL[2])
mtext(side=3,adj=0.03,line=-1.2,'A')
mtext(side=2,line=2,'NMDS Axis 2',las=3,cex=0.8)

plot(MD1,choices=c(1,2),type='n',ylab="",xlab="",xlim=c(-1.1,1.4))
orditorp(MD1, display = "species", choices=c(1,2), priority = priSpp,pcol='gray70',cex=0.8)
mtext(side=3,adj=0.01,line=-1.2,'B')

plot(MD1,choices=c(1,3),type='n',xlab="",ylab="",xlim=c(-1.1,1.4))
for(i in 1:11){with(MDcen[MDcen$H=='J',],Arrows(MDS1[i],MDS3[i],MDS1[i+1],MDS3[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[1]))}
for(i in 1:11){with(MDcen[MDcen$H=='J',],points(MDS3[i]~MDS1[i],pch=21,bg=DBWpalBG[1],col=DBWpalCOL[1]))}
for(i in 1:11){with(MDcen[MDcen$H=='R',],Arrows(MDS1[i],MDS3[i],MDS1[i+1],MDS3[i+1],arr.type='triangle',arr.length=0,col=DBWpalCOL[3]))}
for(i in 1:11){with(MDcen[MDcen$H=='R',],points(MDS3[i]~MDS1[i],pch=22,bg=DBWpalBG[3],col=DBWpalCOL[3]))}
for(i in 1:11){with(MDcen[MDcen$H=='W',],Arrows(MDS1[i],MDS3[i],MDS1[i+1],MDS3[i+1],arr.type='triangle',arr.length=0.2,col=DBWpalCOL[2]))}
for(i in 1:11){with(MDcen[MDcen$H=='W',],points(MDS3[i]~MDS1[i],pch=24,bg=DBWpalBG[2],col=DBWpalCOL[2]))}
mtext(side=3,adj=0.03,line=-1.2,'C')
mtext(side=1,line=2,'NMDS Axis 1',cex=0.8)
mtext(side=2,line=2,'NMDS Axis 3',las=3,cex=0.8)

plot(MD1,choices=c(1,3),type='n',ylab="",xlim=c(-1.1,1.4),xlab='')
#abline(h=0,lty=2,col='gray70')
#abline(v=0,lty=2,col='gray70')
orditorp(MD1, display = "species", choices=c(1,3), priority = priSpp,pcol='gray70',cex=0.8)
mtext(side=3,adj=0.03,line=-1.2,'D')
mtext(side=1,line=2,'NMDS Axis 1',cex=0.8)

dev.off()

###############################################################################  
#TABLE SHOWING INCREASERS, DECREASERS, AND STABLE SPECIES
###############################################################################

vIDS<-veg[veg$idSpecLevel=='1',]#Select fully identified species
vIDS_aggStems<-setNames(aggregate(vIDS$Stems,by=list(vIDS$Species,vIDS$Y),sum),c('Species','Y','Stems')) #Aggregate stem count
vIDS_aggStems<-subset(vIDS_aggStems,Y=='2001' | Y>2009) #Select only 2001, 2010, 2011, and 2012
vIDS_aggStems<-setNames(dcast(vIDS_aggStems,Species~Y,sum),c('Species','Y01','Y10','Y11','Y12')) #Sum stems
vIDS_aggStems$lastThree<-round(with(vIDS_aggStems,(Y10+Y11+Y12)/3),1) #Take mean of last three years
vIDS_aggStems$Delta<-vIDS_aggStems$lastThree-vIDS_aggStems$Y01 #Get change over time 01-last three
vx<-vIDS_aggStems #Give a short name for ease below

#New species
newSpecies<-vx[vx$Y01==0,][with(vx[vx$Y01==0,],order(-lastThree)),]
newSpecies<-merge(newSpecies,traits,by='Species')
newSpecies<-merge(newSpecies,pheno,by='Species')
write.csv(newSpecies,'Table_newSpecies.csv',row.names=FALSE)

#Species that disappeared
disappeared<-vx[vx$Y01>0 & vx$Y10<1 & vx$Y11<1 & vx$Y12<1,]
disappeared<-disappeared[with(disappeared,order(Delta)),]
disappeared<-merge(disappeared,traits,by='Species')
disappeared<-merge(disappeared,pheno,by='Species')
write.csv(disappeared,'Table_disappeared.csv',row.names=FALSE)  

#Biggest increasers
increasers<-head(vx[with(vx,order(-Delta)),],20)
increasers<-merge(increasers,traits,by='Species')
increasers<-merge(increasers,pheno,by='Species')
write.csv(increasers,'Table_increasers.csv',row.names=FALSE)

#Biggest decreasers
decreasers<-head(vx[with(vx,order(Delta)),],20)
decreasers<-merge(decreasers,traits,by='Species')
decreasers<-merge(decreasers,pheno,by='Species')
write.csv(decreasers,'Table_decreasers.csv',row.names=FALSE)

#Most stable species
stable<-vx[with(vx,order(abs(Delta))),]
stable<-subset(stable,Y01>=10)
stable<-head(stable,20)
stable<-merge(stable,traits,by='Species')
stable<-merge(stable,pheno,by='Species')
write.csv(stable,'Table_stable.csv',row.names=FALSE)

vx2<-merge(vx,traits,by='Species')
vx2<-merge(vx2,pheno,by='Species')
write.csv(vx2,'Table_IDS.csv',row.names=FALSE)

