#Prep BRF data and workspace for analyses
#Clear workspace
rm(list=ls())

#Setup
#### Packages
library(permute)
library(lattice)
library(vegan)
library(vegetarian)
library(car)
library(reshape)
library(reshape2)
library(lubridate)
library(dplyr)
#Note that a package not published on the CRAN server "biostats.r" is sourced later. For more information please contact Kevin McGarigal of the University of Massachusetts, Amherst

#### Read in Data 
source('biostats.R')
BRF1<-read.csv("hf097-01.csv") #Read in Black Rock Forest data (hf-097-01 from the Harvard Forest data archive)

#Select Hand & litter collections
antshand <- subset(BRF1, trap.type=="Hand")
antslitter <- subset(BRF1, trap.type=="Litter")
ants1 <- rbind(antshand, antslitter)

##Clean up data
ants1<-ants1[ants1$plot!="A5",]
ants1<-ants1[ants1$plot!="A6",]
ants1<-ants1[ants1$plot!="B6",]
ants1<-ants1[ants1$plot!="B7",]
ants1<-ants1[ants1$plot!="B5",]
ants1<-ants1[ants1$plot!="Z5",]
ants1<-ants1[ants1$plot!="GX",]
ants1<-ants1[ants1$plot!="GC",]
year <- year(ants1$date)
ants1<-cbind(year,ants1)

# Create a vector of ones to append to each record, this will be eventually summed for each species.
incidence <- rep(1,dim(ants1)[1])
# append incidence vector to data frame
ants1<-cbind(ants1, incidence)

# Create wide format abundance matrix
ants2 <- matrix(NA, nrow=96, ncol=(5+length(unique(ants1$spec.code))))
colnames(ants2) <- c("year","block","plot","treatment","deer",as.vector(unique(ants1$spec.code)))
blocks <- c(rep("A",4), rep("B",4), rep("C",4))
treatments <- c("NO","O50","C","OG","O50","OG","NO","C","OG","C","O50","NO")
pretreatment <- c(2006,2007)
posttreatment <- c(2009,2011,2015)
plots <- unique(ants1$plot)
species <-as.vector(unique(ants1$spec.code))
  
# Calculate number of ants collected per plot per species by year for pre-treatment years without exclosures
for(i in 1:length(pretreatment)){# year 
  for(j in 1:length(as.vector(unique(ants1$plot)))){
      for(k in 1:length(species)){
        temp <- ants1[ants1$year==pretreatment[i],]
        ants2[(12*(i-1))+j, 1] <- unique(pretreatment[i])
        ants2[(12*(i-1))+j, 2] <- blocks[j]
        ants2[(12*(i-1))+j, 3] <- as.vector(unique(ants1$plot))[j]
        ants2[(12*(i-1))+j, 4] <- treatments[j]
        ants2[(12*(i-1))+j, 5] <- "NA"  
        temp <- subset(temp, plot==plots[j])
        temp_species <- as.vector(temp$spec.code)
        if(any(species[k]==temp_species)){
          ants2[(12*(i-1))+j, 5+k] <- sum(temp_species==species[k])
        }else{
          ants2[(12*(i-1))+j, 5+k] <- 0
        }
      }
  }
  rm("temp")
}

# Calculate number of ants collected per plot per species by year for post-treatment years without exclosures (main)
for(i in 1:length(posttreatment)){# year 
  for(j in 1:length(as.vector(unique(ants1$plot)))){
      for(k in 1:length(species)){
        temp <- ants1[ants1$year==posttreatment[i],]
        ants2[24+(12*(i-1))+j, 1] <- unique(posttreatment[i])
        ants2[24+(12*(i-1))+j, 2] <- blocks[j]
        ants2[24+(12*(i-1))+j, 3] <- as.vector(unique(ants1$plot))[j]
        ants2[24+(12*(i-1))+j, 4] <- treatments[j]
        ants2[24+(12*(i-1))+j, 5] <- "main"  
        temp <- subset(temp, plot==plots[j])
        temp <- subset(temp, deer=="main")
        temp_species <- as.vector(temp$spec.code)
        if(any(species[k]==temp_species)){
          ants2[24+(12*(i-1))+j, 5+k] <- sum(temp_species==species[k])
        }else{
          ants2[24+(12*(i-1))+j, 5+k] <- 0
        }
      }
  }
  rm("temp")
}

# Calculate number of ants collected per plot per species by year for post-treatment years with exclosures (exclosure)
for(i in 1:length(posttreatment)){# year
  for(j in 1:length(as.vector(unique(ants1$plot)))){
      for(k in 1:length(species)){
        temp <- ants1[ants1$year==posttreatment[i],]
        ants2[60+(12*(i-1))+j, 1] <- unique(posttreatment[i])
        ants2[60+(12*(i-1))+j, 2] <- blocks[j]
        ants2[60+(12*(i-1))+j, 3] <- as.vector(unique(ants1$plot))[j]
        ants2[60+(12*(i-1))+j, 4] <- treatments[j]
        ants2[60+(12*(i-1))+j, 5] <- "exclosure"  
        temp <- subset(temp, plot==plots[j])
        temp <- subset(temp, deer=="exclosure")
        temp_species <- as.vector(temp$spec.code)
        if(any(species[k]==temp_species)){
          ants2[60+(12*(i-1))+j, 5+k] <- sum(temp_species==species[k])
        }else{
          ants2[60+(12*(i-1))+j, 5+k] <- 0
        }
      }
  }
  rm("temp")
}

write.csv(ants2, file="BRF_data_for_analysis.csv")
rm(list=c("ants","ants2"))


##Prep HFR data for analysis  
HFR<-read.csv("hf118-01.csv") #Read in HFR data from Harvard Forest data archive

HFR$abundance[is.na(HFR$abundance)==TRUE]<-1 #Change 2007's NA's into occurences

occur <- rep(1,dim(HFR)[1]) # Create a vector from which to calculate the frequency of occurence for each species
# Create wide format abundance matrix

# Create wide format abundance matrix
ants2 <- matrix(NA, nrow=136, ncol=(5+length(unique(ants$code)))) # nrow = # of plots times # of years
colnames(ants2) <- c("year","block","plot","treatment","moose",as.vector(unique(ants$code)))
blocks <- c(rep("Ridge",3), rep("Valley",4),"Ridge")
treatments <- c("G","L","He","L","G","He","Hw","Hw")
preexclosure <- c(2003:2011)
postexclosure <- c(2012:2015)
plots <- unique(ants$plot)
species <-as.vector(unique(ants$code))

# Calculate number of ants collected per plot per species by year for years without exclosures
for(i in 1:length(preexclosure)){# year 
  for(j in 1:length(as.vector(unique(ants$plot)))){
      for(k in 1:length(species)){
        temp <- ants[ants$year==preexclosure[i],]
        ants2[(8*(i-1))+j, 1] <- preexclosure[i]
        ants2[(8*(i-1))+j, 2] <- blocks[j]
        ants2[(8*(i-1))+j, 3] <- as.vector(unique(ants$plot))[j]
        ants2[(8*(i-1))+j, 4] <- treatments[j]
        ants2[(8*(i-1))+j, 5] <- "NA"  
        temp <- subset(temp, plot==plots[j])
        temp_species <- as.vector(temp$code)
        if(any(species[k]==temp_species)){
          ants2[(8*(i-1))+j, 5+k] <- sum(temp_species==species[k])
        }else{
          ants2[(8*(i-1))+j, 5+k] <- 0
        }
      }
  }
  rm("temp")
}

# Calculate number of ants collected per plot per species by year for years without exclosures 
for(i in 1:length(postexclosure)){# year 
  for(j in 1:length(as.vector(unique(ants$plot)))){
      for(k in 1:length(species)){
        temp <- ants[ants$year==postexclosure[i],]
        ants2[72+(8*(i-1))+j, 1] <- unique(postexclosure[i])
        ants2[72+(8*(i-1))+j, 2] <- blocks[j]
        ants2[72+(8*(i-1))+j, 3] <- as.vector(unique(ants$plot))[j]
        ants2[72+(8*(i-1))+j, 4] <- treatments[j]
        ants2[72+(8*(i-1))+j, 5] <- "exterior"  
        temp <- subset(temp, plot==plots[j])
        temp <- subset(temp, moose.cage=="Exterior")
        temp_species <- as.vector(temp$code)
        if(any(species[k]==temp_species)){
          ants2[72+(8*(i-1))+j, 5+k] <- sum(temp_species==species[k])
        }else{
          ants2[72+(8*(i-1))+j, 5+k] <- 0
        }
      }
  }
  rm("temp")
}

# Calculate number of ants collected per plot per species by year for post-treatment years with exclosures (exclosure)
for(i in 1:length(postexclosure)){# year 
  for(j in 1:length(as.vector(unique(ants$plot)))){
      for(k in 1:length(species)){
        temp <- ants[ants$year==postexclosure[i],]
        ants2[104+(8*(i-1))+j, 1] <- unique(postexclosure[i])
        ants2[104+(8*(i-1))+j, 2] <- blocks[j]
        ants2[104+(8*(i-1))+j, 3] <- as.vector(unique(ants$plot))[j]
        ants2[104+(8*(i-1))+j, 4] <- treatments[j]
        ants2[104+(8*(i-1))+j, 5] <- "interior"  
        temp <- subset(temp, plot==plots[j])
        temp <- subset(temp, moose.cage=="Interior")
        temp_species <- as.vector(temp$code)
        if(any(species[k]==temp_species)){
          ants2[104+(8*(i-1))+j, 5+k] <- sum(temp_species==species[k])
        }else{
          ants2[104+(8*(i-1))+j, 5+k] <- 0
        }
      }
  }
  rm("temp")
}

write.csv(ants2, file="Simes_data_for_analysis.csv")

rm(list=c("ants1","ants2"))

##Rarefaction by treatment based on different orders of diversity (see Jost 2007)

BRF2 <- data.frame(read.csv(file="BRF_data_for_analysis.csv",header=TRUE, stringsAsFactors=FALSE))
HFR2 <- data.frame(read.csv(file="Simes_data_for_analysis.csv",header=TRUE,stringsAsFactors=FALSE))

# Source iNEXT package (no longer available on CRAN)
source("C:/Users/srecord/Dropbox/Sime-BRF_ms/analyses/iNEXT/R/iNEXT.r")
source("C:/Users/srecord/Dropbox/Sime-BRF_ms/analyses/iNEXT/R/EstIndex.r")
source("C:/Users/srecord/Dropbox/Sime-BRF_ms/analyses/iNEXT/R/JADE.r")
source("C:/Users/srecord/Dropbox/Sime-BRF_ms/analyses/iNEXT/R/invChat.r")
library(ggplot2)

# BRF Rarefaction
treatmentsums_BRF <- data.frame(matrix(NA, nrow=4,ncol=dim(BRF2)[2]-4))
colnames(treatmentsums_BRF) <- c("treatment", colnames(BRF2[6:length(colnames(BRF2))]))

for(i in 1:4){
  treatmentsums_BRF[i,1] <- unique(BRF2$treatment)[i]
  temp <- subset(BRF2, treatment==unique(BRF2$treatment)[i])
  treatmentsums_BRF[i,2:dim(treatmentsums_BRF)[2]] <- apply(temp[,6:dim(temp)[2]], MARGIN=2, FUN=sum)
}

#Turn treatmentsums into a list for iNEXT functions
NO <- treatmentsums_BRF[1,2:dim(treatmentsums_BRF)[2]]
O50 <- treatmentsums_BRF[2,2:dim(treatmentsums_BRF)[2]]
Control <- treatmentsums_BRF[3,2:dim(treatmentsums_BRF)[2]]
OG <- treatmentsums_BRF[4,2:dim(treatmentsums_BRF)[2]]

treatmentsums_BRF_list <- list(NO=NO, O50=O50, Control=Control, OG=OG)

out_BRF <- iNEXT(treatmentsums_BRF_list, q=c(0,1,2), datatype="abundance", endpoint=2000)

# HFR Rarefaction
treatmentsums_HFR <- data.frame(matrix(NA, nrow=4,ncol=dim(HFR2)[2]-4))
colnames(treatmentsums_HFR) <- c("treatment", colnames(HFR2[6:length(colnames(HFR2))]))

for(i in 1:4){
  treatmentsums_HFR[i,1] <- unique(HFR2$treatment)[i]
  temp <- subset(HFR2, treatment==unique(HFR2$treatment)[i])
  treatmentsums_HFR[i,2:dim(treatmentsums_HFR)[2]] <- apply(temp[,6:dim(temp)[2]], MARGIN=2, FUN=sum)
}

#Turn treatmentsums into a list for iNEXT functions
G <- treatmentsums_HFR[1,2:dim(treatmentsums_HFR)[2]]
L <- treatmentsums_HFR[2,2:dim(treatmentsums_HFR)[2]]
He <- treatmentsums_HFR[3,2:dim(treatmentsums_HFR)[2]]
Hw <- treatmentsums_HFR[4,2:dim(treatmentsums_HFR)[2]]

treatmentsums_HFR_list <- list(G=G, L=L, He=He, Hw=Hw)

out_HFR <- iNEXT(treatmentsums_HFR_list, q=c(0,1,2), datatype="abundance", endpoint=2000)

# Color plots
p1 <- ggiNEXT(out_BRF, type=1, facet.var="order", se=TRUE, color.var='site') 
p1 + labs(title="Black Rock Forest") 
p2 <- ggiNEXT(out_HFR, type=1, facet.var="order", se=TRUE)
p2 + labs(title="Harvard Forest")+ scale_colour_manual(values=c('goldenrod4','darkgreen','purple','red')) + scale_fill_manual(values=c('goldenrod4','darkgreen','purple','red'))

# Black and white plots
p1 <- ggiNEXT(out_BRF, type=1, facet.var="order", se=TRUE, grey=TRUE)
p1 + labs(title="Black Rock Forest")
p2 <- ggiNEXT(out, type=1, facet.var="order", se=TRUE, grey=TRUE) 
p2 + labs(title="Harvard Forest")

##Species accumulation curves in Appendix S2
# Split graphics space for 2 plots
par(mfrow=c(1,2))
#BRF
spec<-as.vector(BRF1$spec.code)
yeaR<-year(BRF1$date)
antsR<-cbind(spec,yeaR) #Need as.vector in there otherwise it returns factor numbers. 

numspec<-c(length(unique(subset(antsR, yeaR=="2006", select=spec))),length(unique(subset(antsR, yeaR<="2007", select=spec))),length(unique(subset(antsR, yeaR<="2009", select=spec))),length(unique(subset(antsR, yeaR<="2011", select=spec))),length(unique(subset(antsR, yeaR<="2015", select=spec))))

plot(numspec~unique(yeaR),main= "Black Rock Forest",xlab="Year",ylab="Total Number of Species", cex.axis=1.2, cex.lab=1.2, pch=16, ylim=c(10,50), cex=1.3)
rm(list=c("spec","yeaR", "antsR"))

#Determine total number of occurence samples from BRF from 2006-2015
sum(apply(BRF2[,6:dim(BRF2)[2]], MARGIN=1, FUN=sum))

#HFR
yearnames<-c("2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015")

#Determine total number of occurence samples from HFR from 2003-2015
sum(apply(HFR2[,6:dim(HFR2)[2]], MARGIN=1, FUN=sum))

plot(numspec~unique(yeaR), main= "Harvard Forest",xlab="Year",ylab="Total Number of Species", cex.axis=1.2, cex.lab=1.2, pch=16, ylim=c(10,50), cex=1.3)

######### Univariate Shannon diversity #####################
library(vegetarian)
library(nlme)

shannonBRF <- vector('numeric',dim(BRF2)[1])
for(i in 1:dim(BRF2)[1]){
  shannonBRF[i] <- d(BRF2[i,6:47], q=1)
}

shannonBRF <- cbind(BRF2[,1:5], shannonBRF)

Shannonpvals <- NULL
Shannondf <- NULL

# BRF pre-treatment SHannon diversity analysis
shannonBRFpre <- shannonBRF[shannonBRF$year==2006,]
mod <- lme(shannonBRF ~ treatment , random = ~ 1|block, data=shannonBRFpre)
Anova(mod, type='II')
plot(resid(mod)) # model residuals do not indicate issues with non-normality or variance

pout <- matrix(NA, ncol=4, nrow=1)
names(pout) <- c("test", "predictor", "F", "rawP")
pout[,1] <- rep("BRFShannonpre", 1)
pout[,2] <- c("treatment")
pout[,3] <- Anova(mod, type='II')[[1]]
pout[,4] <- Anova(mod, type='II')[[3]]

Shannonpvals <- rbind(Shannonpvals,pout)

dfout <- matrix(NA, ncol=3, nrow=2)
names(dfout) <- c("test", "term", "df")
dfout[,1] <- rep("BRFShannonpre", 2)
dfout[,2] <- c("treatment","Residuals")
dfout[1,3] <- Anova(mod, type='II')[[2]]
dfout[2,3] <- summary(mod)$tTable[1,3]
Shannondf <- rbind(Shannondf, dfout)

# clean up
rm(list=c("pout", "dfout", "mod"))

# Post-treatment analysis
shannonBRFpost <- shannonBRF[shannonBRF$year>2006,]
mod <- lme(shannonBRF ~ treatment/deer + year + treatment*year, random = ~ 1|block, data=shannonBRFpost)
Anova(mod, type='II')
plot(resid(mod)) # model residuals do not indicate issues with non-normality or variance

pout <- matrix(NA, ncol=4, nrow=4)
names(pout) <- c("test", "predictor", "F", "rawP")
pout[,1] <- rep("BRFShannonpost", 4)
pout[,2] <- c("treatment","year",'treatment:deer','treatment:year')
pout[,3] <- Anova(mod, type='II')[[1]]
pout[,4] <- Anova(mod, type='II')[[3]]

Shannonpvals <- rbind(Shannonpvals,pout)

dfout <- matrix(NA, ncol=3, nrow=5)
names(dfout) <- c("test", "term", "df")
dfout[,1] <- rep("BRFShannonpost", 5)
dfout[,2] <- c("treatment", "year",'treatment:deer','treatment:year',"Residuals")
dfout[1:4,3] <- Anova(mod, type='II')[[2]]
dfout[5,3] <- summary(mod)$tTable[1,3]
Shannondf <- rbind(Shannondf, dfout)

# clean up
rm(list=c("pout", "dfout", "mod"))

# Contrasts for deer exclosure 
OG <- shannonBRFpost[shannonBRFpost$treatment=="OG",]
NO <- shannonBRFpost[shannonBRFpost$treatment=="NO",]
O50 <- shannonBRFpost[shannonBRFpost$treatment=="O50",]
C<- shannonBRFpost[shannonBRFpost$treatment=="C",]

mod <- lme(shannonBRF ~ deer, random=~1|block, data=OG)
Anova(mod, type='II')
mod <- lme(shannonBRF ~ deer, random=~1|block, data=O50)
Anova(mod, type='II')
mod <- lme(shannonBRF ~ deer, random=~1|block, data=NO)
Anova(mod, type='II')
mod <- lme(shannonBRF ~ deer, random=~1|block, data=C)
Anova(mod, type='II')

shannonHF <- vector('numeric',dim(HFR2)[1])
for(i in 1:dim(HFR2)[1]){
  shannonHF[i] <- d(HFR2[i,7:52], q=1)
}

shannonHF <- cbind(HFR2[,1:6], shannonHF)

# HF pre-treatment SHannon diversity analysis
shannonHFpre <- shannonHF[shannonHF$year<2005,]
mod <- lme(shannonHF ~ treatment , random = ~ 1|block, data=shannonHFpre)
Anova(mod, type='II')
plot(resid(mod)) # model residuals do not indicate issues with non-normality or variance

pout <- matrix(NA, ncol=4, nrow=1)
names(pout) <- c("test", "predictor", "F", "rawP")
pout[,1] <- rep("HFShannonpre", 1)
pout[,2] <- c("treatment")
pout[,3] <- Anova(mod, type='II')[[1]]
pout[,4] <- Anova(mod, type='II')[[3]]

Shannonpvals <- rbind(Shannonpvals,pout)

dfout <- matrix(NA, ncol=3, nrow=2)
names(dfout) <- c("test", "term", "df")
dfout[,1] <- rep("HFShannonpre", 2)
dfout[,2] <- c("treatment","Residuals")
dfout[1,3] <- Anova(mod, type='II')[[2]]
dfout[2,3] <- summary(mod)$tTable[1,3]
Shannondf <- rbind(Shannondf, dfout)

# clean up
rm(list=c("pout", "dfout", "mod"))

# Pre-treatment contrasts
HeGL <- shannonHFpre[shannonHFpre$treatment!='Hw',]

mod <- lme(shannonHF ~ treatment, random=~1|block, HeGL)
Anova(mod, type='II')

# Post-treatment analysis all years, no exclosure
shannonHFpost <- shannonHF[shannonHF$year>2004,]
shannonHFinterior <- subset(shannonHFpost, moose!='interior')
shannonHFpost <- rbind(shannonHFpost[shannonHFpost$year<2012,],shannonHFinterior)
infest <- c(rep('preinfest',40),rep('postinfest',48))
shannonHFpost <- cbind(shannonHFpost,infest)

mod <- lme(shannonHF ~ treatment + year + treatment*year + infest, random = ~ 1|block, data=shannonHFpost)
Anova(mod, type='II')
plot(resid(mod)) # model residuals do not indicate issues with non-normality or variance

pout <- matrix(NA, ncol=4, nrow=4)
names(pout) <- c("test", "predictor", "F", "rawP")
pout[,1] <- rep("HFShannonpost", 4)
pout[,2] <- c("treatment","year",'infest','treatment:year')
pout[,3] <- Anova(mod, type='II')[[1]]
pout[,4] <- Anova(mod, type='II')[[3]]

Shannonpvals <- rbind(Shannonpvals,pout)

dfout <- matrix(NA, ncol=3, nrow=5)
names(dfout) <- c("test", "term", "df")
dfout[,1] <- rep("HFShannonpost", 5)
dfout[,2] <- c("treatment", "year",'infest','treatment:year',"Residuals")
dfout[1:4,3] <- Anova(mod, type='II')[[2]]
dfout[5,3] <- summary(mod)$tTable[1,3]
Shannondf <- rbind(Shannondf, dfout)

# clean up
rm(list=c("pout", "dfout", "mod"))

#Contrasts
Hepost <- subset(shannonHFpost, treatment=='He')
Hwpost <- subset(shannonHFpost, treatment=='Hw')
Gpost <- subset(shannonHFpost, treatment=='G')
Lpost <- subset(shannonHFpost, treatment=='L')

HeG <- rbind(Hepost, Gpost)
mod <- lme(shannonHF ~ treatment, random = ~ 1|block, data=HeG)
Anova(mod, type='II')

HeHw <- rbind(Hepost, Hwpost)
mod <- lme(shannonHF ~ treatment, random = ~ 1|block, data=HeHw)
Anova(mod, type='II')

HeL <- rbind(Hepost, Lpost)
mod <- lme(shannonHF ~ treatment, random = ~ 1|block, data=HeL)
Anova(mod, type='II')

GL <- rbind(Gpost, Lpost)
mod <- lme(shannonHF ~ treatment, random = ~ 1|block, data=GL)
Anova(mod, type='II')

HwL <- rbind(Hwpost, Lpost)
mod <- lme(shannonHF ~ treatment, random = ~ 1|block, data=HwL)
Anova(mod, type='II')

HwG <- rbind(Hwpost, Gpost)
mod <- lme(shannonHF ~ treatment, random = ~ 1|block, data=HwG)
Anova(mod, type='II')

# Post-treatment analysis exclosure
shannonHFdeer <- shannonHF[shannonHF$year>2011,]

mod <- lme(shannonHF ~ treatment + year + treatment/moose + treatment*year, random = ~ 1|block, data=shannonHFdeer)
Anova(mod, type='II')
plot(resid(mod)) # model residuals do not indicate issues with non-normality or variance

pout <- matrix(NA, ncol=4, nrow=4)
names(pout) <- c("test", "predictor", "F", "rawP")
pout[,1] <- rep("HFShannondeer", 4)
pout[,2] <- c("treatment","year",'treatment:moose','treatment:year')
pout[,3] <- Anova(mod, type='II')[[1]]
pout[,4] <- Anova(mod, type='II')[[3]]

Shannonpvals <- rbind(Shannonpvals,pout)

dfout <- matrix(NA, ncol=3, nrow=5)
names(dfout) <- c("test", "term", "df")
dfout[,1] <- rep("HFShannondeer", 5)
dfout[,2] <- c("treatment", "year",'treatment:deer','treatment:year',"Residuals")
dfout[1:4,3] <- Anova(mod, type='II')[[2]]
dfout[5,3] <- summary(mod)$tTable[1,3]
Shannondf <- rbind(Shannondf, dfout)

# clean up
rm(list=c("pout", "dfout", "mod"))

write.csv(Shannonpvals, "Shannonpvals.csv", row.names=FALSE)
write.csv(Shannondf, "Shannondf.csv", row.names=FALSE)

# Plot Shannon diversity
# sum dbh for individuals with multiple stems
meanShannonBRF <- summarise(group_by(shannonBRFpost, year, treatment, deer), 
                            meanShannon = mean(shannonBRF)
) 
shannonHFpost <- shannonHFpost[shannonHFpost$year<2012,]
meanShannonHF <- summarise(group_by(shannonHFpost, year, treatment, moose), 
                           meanShannon = mean(shannonHF)
) 

meanShannonHFdeer <- summarise(group_by(shannonHFdeer, year, treatment, moose), 
                               meanShannon = mean(shannonHF)
) 

OG <- meanShannonBRF[meanShannonBRF$treatment=="OG",]
OGex <- OG[OG$deer=='exclosure',]
OGmain <- OG[OG$deer=='main',]

O50 <- meanShannonBRF[meanShannonBRF$treatment=="O50",]
O50ex <- O50[O50$deer=='exclosure',]
O50main <- O50[O50$deer=='main',]

NO <- meanShannonBRF[meanShannonBRF$treatment=="NO",]
NOex <- NO[NO$deer=='exclosure',]
NOmain <- NO[NO$deer=='main',]

C <- meanShannonBRF[meanShannonBRF$treatment=="C",]
Cex <- C[C$deer=='exclosure',]
Cmain <- C[C$deer=='main',]

He <- meanShannonHF[meanShannonHF$treatment=='He',]
Hw <- meanShannonHF[meanShannonHF$treatment=='Hw',]
G <- meanShannonHF[meanShannonHF$treatment=='G',]
L <- meanShannonHF[meanShannonHF$treatment=='L',]

Hedeer <- meanShannonHFdeer[meanShannonHFdeer$treatment=='He',]
Heex <- Hedeer[Hedeer$moose=='interior',]
Hemain <- Hedeer[Hedeer$moose=='exterior',]
Hwdeer <- meanShannonHFdeer[meanShannonHFdeer$treatment=='Hw',]
Hwex <- Hwdeer[Hwdeer$moose=='interior',]
Hwmain <- Hwdeer[Hwdeer$moose=='exterior',]
Ldeer <- meanShannonHFdeer[meanShannonHFdeer$treatment=='L',]
Lex <- Ldeer[Ldeer$moose=='interior',]
Lmain <- Ldeer[Ldeer$moose=='exterior',]
Gdeer <- meanShannonHFdeer[meanShannonHFdeer$treatment=='G',]
Gex <- Gdeer[Gdeer$moose=='interior',]
Gmain <- Gdeer[Gdeer$moose=='exterior',]

par(mfrow=c(2,2))

plot(He$year, He$meanShannon, type='line', ylim=c(0,max(meanShannonHFdeer$meanShannon)), xlim=c(2005,2015), ylab="Shannon diversity", xlab="Year", col='darkgreen', cex.lab=1.2, cex.axis=1.2, cex.main=1.2, main='HF-HeRE', lwd=3)
lines(Hemain$year, Hemain$meanShannon, col='darkgreen', lwd=3)
lines(Heex$year, Heex$meanShannon, col='darkgreen', lty='dashed', lwd=3)
lines(G$year, G$meanShannon, col='goldenrod', lwd=3)
lines(Gmain$year, Gmain$meanShannon, col='goldenrod', lwd=3)
lines(Gex$year, Gex$meanShannon, col='goldenrod', lty='dashed', lwd=3)
lines(L$year, L$meanShannon, col='red', lwd=3)
lines(Lmain$year, Lmain$meanShannon, col='red', lwd=3)
lines(Lex$year, Lex$meanShannon, col='red', lty='dashed', lwd=3)
lines(Hw$year, Hw$meanShannon, col='purple', lwd=3)
lines(Hwmain$year, Hwmain$meanShannon, col='purple', lwd=3)
lines(Hwex$year, Hwex$meanShannon, col='purple', lty='dashed', lwd=3)
mtext('a)',3,adj=0, cex=1.5, line=2)
plot.new()
legend(locator(1), legend=c("Hemlock","Girdled", "Logged","Hardwood"), fill=c("darkgreen","goldenrod","red", "purple"), bty='n' , cex=1.3)

plot(OGmain$year, OGmain$meanShannon, type='line', ylim=c(0,max(meanShannonBRF$meanShannon)), xlim=c(2009,2015), ylab="Shannon diversity", xlab="Year", col='darkorchid2', cex.lab=1.2, cex.axis=1.2, cex.main=1.2, main='BRF-FOFE', lwd=3)
lines(OGex$year, OGex$meanShannon, col='darkorchid2', lwd=3, lty='dashed')
lines(O50main$year, O50main$meanShannon, col="cyan", lwd=3)
lines(O50ex$year, O50ex$meanShannon, col="cyan", lwd=3, lty='dashed')
lines(NOmain$year, NOmain$meanShannon, col="darkolivegreen3", lwd=3)
lines(NOex$year, NOex$meanShannon, col="darkolivegreen3", lwd=3, lty='dashed')
lines(Cmain$year, Cmain$meanShannon, col="salmon", lwd=3)
lines(Cex$year, Cex$meanShannon, col="salmon", lwd=3, lty='dashed')
mtext('b)',3,adj=0, cex=1.5, line=2)
plot.new()
legend(locator(1), legend=c("Control","50% Oak girdled", "100% Oak girdled","Non-oaks girdled"), fill=c("salmon","cyan","darkorchid2", "darkolivegreen3"), bty='n' , cex=1.3)

############ Plot ant genera over time for each experiment ###############
ants<-HFR

ants$abundance[is.na(ants$abundance)==TRUE]<-1 #Change 2007's NA's into occurences
incidence <- rep(1,dim(ants)[1]) # Create a vector from which to calculate the frequency of occurence for each species
ants <- cbind(ants, incidence)

HFgenus <- summarise(group_by(ants, year, treatment, moose.cage, genus), 
                     incidence = sum(incidence)
) 

# determine top ten genera
HFtopten <- summarise(group_by(HFgenus, genus), 
                      incidence = sum(incidence)
) 
HFgenus1 <- HFgenus
HFgenus <- HFgenus[HFgenus$year<2012,]
HFgenus <- HFgenus[HFgenus$year>2004,]
HFgenus2 <-HFgenus1[HFgenus1$year>2011,]

#Get incidence totals for different treatments to calculate relative abundances
G <- HFgenus[HFgenus$treatment=='Girdled',]
G <- summarise(group_by(G, year), incidence = sum(incidence)) 
G2<- HFgenus2[HFgenus2$treatment=='Girdled',]
Gex <- G2[G2$moose.cage=='Interior',]
Gex <- summarise(group_by(Gex, year), incidence = sum(incidence)) 
Gmain <- G2[G2$moose.cage=='Exterior',]
Gmain <- summarise(group_by(Gmain, year), incidence = sum(incidence)) 

He <- HFgenus[HFgenus$treatment=='HemlockControl',]
He <- summarise(group_by(He, year), incidence = sum(incidence)) 
He2<- HFgenus2[HFgenus2$treatment=='HemlockControl',]
Heex <- He2[He2$moose.cage=='Interior',]
Heex <- summarise(group_by(Heex, year), incidence = sum(incidence)) 
Hemain <- He2[He2$moose.cage=='Exterior',]
Hemain <- summarise(group_by(Hemain, year), incidence = sum(incidence)) 

Hw <- HFgenus[HFgenus$treatment=='HardwoodControl',]
Hw <- summarise(group_by(Hw, year), incidence = sum(incidence)) 
Hw2<- HFgenus2[HFgenus2$treatment=='HardwoodControl',]
Hwex <- Hw2[Hw2$moose.cage=='Interior',]
Hwex <- summarise(group_by(Hwex, year), incidence = sum(incidence)) 
Hwmain <- Hw2[Hw2$moose.cage=='Exterior',]
Hwmain <- summarise(group_by(Hwmain, year), incidence = sum(incidence)) 

L <- HFgenus[HFgenus$treatment=='Logged',]
L <- summarise(group_by(L, year), incidence = sum(incidence)) 
L2<- HFgenus2[HFgenus2$treatment=='Logged',]
Lex <- L2[L2$moose.cage=='Interior',]
Lex <- summarise(group_by(Lex, year), incidence = sum(incidence)) 
Lmain <- L2[L2$moose.cage=='Exterior',]
Lmain <- summarise(group_by(Lmain, year), incidence = sum(incidence)) 

# Select data for top ten genera
aph <- HFgenus[HFgenus$genus=='Aphaenogaster',]
aphHe <- aph[aph$treatment=='HemlockControl',]
aphG <- aph[aph$treatment=='Girdled',]
aphL <- aph[aph$treatment=='Logged',]
aphHw <- aph[aph$treatment=='HardwoodControl',]
aphaenogaster <- HFgenus2[HFgenus2$genus=='Aphaenogaster',]
aphex <- aphaenogaster[aphaenogaster$moose.cage=='Interior',]
aphexHe <- aphex[aphex$treatment=='HemlockControl',]
aphexG <- aphex[aphex$treatment=='Girdled',]
aphexL <- aphex[aphex$treatment=='Logged',]
aphexHw <- aphex[aphex$treatment=='HardwoodControl',]
aphmain <- aphaenogaster[aphaenogaster$moose.cage=='Exterior',]
aphmainHe <- aphmain[aphmain$treatment=='HemlockControl',]
aphmainG <- aphmain[aphmain$treatment=='Girdled',]
aphmainL <- aphmain[aphmain$treatment=='Logged',]
aphmainHw <- aphmain[aphmain$treatment=='HardwoodControl',]

cam <- HFgenus[HFgenus$genus=='Camponotus',]
camHe <- cam[cam$treatment=='HemlockControl',]
camG <- cam[cam$treatment=='Girdled',]
camL <- cam[cam$treatment=='Logged',]
camHw <- cam[cam$treatment=='HardwoodControl',]
camponotus <- HFgenus2[HFgenus2$genus=='Camponotus',]
camex <- camponotus[camponotus$moose.cage=='Interior',]
camexHe <- camex[camex$treatment=='HemlockControl',]
camexG <- camex[camex$treatment=='Girdled',]
camexL <- camex[camex$treatment=='Logged',]
camexHw <- camex[camex$treatment=='HardwoodControl',]
cammain <- camponotus[camponotus$moose.cage=='Exterior',]
cammainHe <- cammain[cammain$treatment=='HemlockControl',]
cammainG <- cammain[cammain$treatment=='Girdled',]
cammainL <- cammain[cammain$treatment=='Logged',]
cammainHw <- cammain[cammain$treatment=='HardwoodControl',]

cre <- HFgenus[HFgenus$genus=='Crematogaster',]
creHe <- cre[cre$treatment=='HemlockControl',]
creG <- cre[cre$treatment=='Girdled',]
creL <- cre[cre$treatment=='Logged',]
creHw <- cre[cre$treatment=='HardwoodControl',]
Crematogaster <- HFgenus2[HFgenus2$genus=='Crematogaster',]
creex <- Crematogaster[Crematogaster$moose.cage=='Interior',]
creexHe <- creex[creex$treatment=='HemlockControl',]
creexG <- creex[creex$treatment=='Girdled',]
creexL <- creex[creex$treatment=='Logged',]
creexHw <- creex[creex$treatment=='HardwoodControl',]
cremain <- Crematogaster[Crematogaster$moose.cage=='Exterior',]
cremainHe <- cremain[cremain$treatment=='HemlockControl',]
cremainG <- cremain[cremain$treatment=='Girdled',]
cremainL <- cremain[cremain$treatment=='Logged',]
cremainHw <- cremain[cremain$treatment=='HardwoodControl',]

form <- HFgenus[HFgenus$genus=='Formica',]
formHe <- form[form$treatment=='HemlockControl',]
formG <- form[form$treatment=='Girdled',]
formL <- form[form$treatment=='Logged',]
formHw <- form[form$treatment=='HardwoodControl',]
formica <- HFgenus2[HFgenus2$genus=='Formica',]
formex <- formica[formica$moose.cage=='Interior',]
formexHe <- formex[formex$treatment=='HemlockControl',]
formexG <- formex[formex$treatment=='Girdled',]
formexL <- formex[formex$treatment=='Logged',]
formexHw <- formex[formex$treatment=='HardwoodControl',]
formmain <- formica[formica$moose.cage=='Exterior',]
formmainHe <- formmain[formmain$treatment=='HemlockControl',]
formmainG <- formmain[formmain$treatment=='Girdled',]
formmainL <- formmain[formmain$treatment=='Logged',]
formmainHw <- formmain[formmain$treatment=='HardwoodControl',]

las <- HFgenus[HFgenus$genus=='Lasius',]
lasHe <- las[las$treatment=='HemlockControl',]
lasG <- las[las$treatment=='Girdled',]
lasL <- las[las$treatment=='Logged',]
lasHw <- las[las$treatment=='HardwoodControl',]
lasius <- HFgenus2[HFgenus2$genus=='Lasius',]
lasex <- lasius[lasius$moose.cage=='Interior',]
lasexHe <- lasex[lasex$treatment=='HemlockControl',]
lasexG <- lasex[lasex$treatment=='Girdled',]
lasexL <- lasex[lasex$treatment=='Logged',]
lasexHw <- lasex[lasex$treatment=='HardwoodControl',]
lasmain <- lasius[lasius$moose.cage=='Exterior',]
lasmainHe <- lasmain[lasmain$treatment=='HemlockControl',]
lasmainG <- lasmain[lasmain$treatment=='Girdled',]
lasmainL <- lasmain[lasmain$treatment=='Logged',]
lasmainHw <- lasmain[lasmain$treatment=='HardwoodControl',]

myr <- HFgenus[HFgenus$genus=='Myrmecina',]
myrHe <- myr[myr$treatment=='HemlockControl',]
myrG <- myr[myr$treatment=='Girdled',]
myrL <- myr[myr$treatment=='Logged',]
myrHw <- myr[myr$treatment=='HardwoodControl',]
Myrmecina <- HFgenus2[HFgenus2$genus=='Myrmecina',]
myrex <- Myrmecina[Myrmecina$moose.cage=='Interior',]
myrexHe <- myrex[myrex$treatment=='HemlockControl',]
myrexG <- myrex[myrex$treatment=='Girdled',]
myrexL <- myrex[myrex$treatment=='Logged',]
myrexHw <- myrex[myrex$treatment=='HardwoodControl',]
myrmain <- Myrmecina[Myrmecina$moose.cage=='Exterior',]
myrmainHe <- myrmain[myrmain$treatment=='HemlockControl',]
myrmainG <- myrmain[myrmain$treatment=='Girdled',]
myrmainL <- myrmain[myrmain$treatment=='Logged',]
myrmainHw <- myrmain[myrmain$treatment=='HardwoodControl',]

myrmic <- HFgenus[HFgenus$genus=='Myrmica',]
myrmicHe <- myrmic[myrmic$treatment=='HemlockControl',]
myrmicG <- myrmic[myrmic$treatment=='Girdled',]
myrmicL <- myrmic[myrmic$treatment=='Logged',]
myrmicHw <- myrmic[myrmic$treatment=='HardwoodControl',]
myrmica <- HFgenus2[HFgenus2$genus=='Myrmica',]
myrmicex <- myrmica[myrmica$moose.cage=='Interior',]
myrmicexHe <- myrmicex[myrmicex$treatment=='HemlockControl',]
myrmicexG <- myrmicex[myrmicex$treatment=='Girdled',]
myrmicexL <- myrmicex[myrmicex$treatment=='Logged',]
myrmicexHw <- myrmicex[myrmicex$treatment=='HardwoodControl',]
myrmicmain <- myrmica[myrmica$moose.cage=='Exterior',]
myrmicmainHe <- myrmicmain[myrmicmain$treatment=='HemlockControl',]
myrmicmainG <- myrmicmain[myrmicmain$treatment=='Girdled',]
myrmicmainL <- myrmicmain[myrmicmain$treatment=='Logged',]
myrmicmainHw <- myrmicmain[myrmicmain$treatment=='HardwoodControl',]

ste <- HFgenus[HFgenus$genus=='Stenamma',]
steHe <- ste[ste$treatment=='HemlockControl',]
steG <- ste[ste$treatment=='Girdled',]
steL <- ste[ste$treatment=='Logged',]
steHw <- ste[ste$treatment=='HardwoodControl',]
Stenamma <- HFgenus2[HFgenus2$genus=='Stenamma',]
stemicex <- Stenamma[Stenamma$moose.cage=='Interior',]
stemicexHe <- stemicex[stemicex$treatment=='HemlockControl',]
stemicexG <- stemicex[stemicex$treatment=='Girdled',]
stemicexL <- stemicex[stemicex$treatment=='Logged',]
stemicexHw <- stemicex[stemicex$treatment=='HardwoodControl',]
stemicmain <- Stenamma[Stenamma$moose.cage=='Exterior',]
stemicmainHe <- stemicmain[stemicmain$treatment=='HemlockControl',]
stemicmainG <- stemicmain[stemicmain$treatment=='Girdled',]
stemicmainL <- stemicmain[stemicmain$treatment=='Logged',]
stemicmainHw <- stemicmain[stemicmain$treatment=='HardwoodControl',]

tapmic <- HFgenus[HFgenus$genus=='Tapinoma',]
tapHe <- tap[tap$treatment=='HemlockControl',]
tapG <- tap[tap$treatment=='Girdled',]
tapL <- tap[tap$treatment=='Logged',]
tapHw <- tap[tap$treatment=='HardwoodControl',]
Tapinoma <- HFgenus2[HFgenus2$genus=='Tapinoma',]
tapmicex <- Tapinoma[Tapinoma$moose.cage=='Interior',]
tapmicexHe <- tapmicex[tapmicex$treatment=='HemlockControl',]
tapmicexG <- tapmicex[tapmicex$treatment=='Girdled',]
tapmicexL <- tapmicex[tapmicex$treatment=='Logged',]
tapmicexHw <- tapmicex[tapmicex$treatment=='HardwoodControl',]
tapmicmain <- Tapinoma[Tapinoma$moose.cage=='Exterior',]
tapmicmainHe <- tapmicmain[tapmicmain$treatment=='HemlockControl',]
tapmicmainG <- tapmicmain[tapmicmain$treatment=='Girdled',]
tapmicmainL <- tapmicmain[tapmicmain$treatment=='Logged',]
tapmicmainHw <- tapmicmain[tapmicmain$treatment=='HardwoodControl',]

temmic <- HFgenus[HFgenus$genus=='Temnothorax',]
temHe <- tem[tem$treatment=='HemlockControl',]
temG <- tem[tem$treatment=='Girdled',]
temL <- tem[tem$treatment=='Logged',]
temHw <- tem[tem$treatment=='HardwoodControl',]
Temnothorax <- HFgenus2[HFgenus2$genus=='Temnothorax',]
temmicex <- Temnothorax[Temnothorax$moose.cage=='Interior',]
temmicexHe <- temmicex[temmicex$treatment=='HemlockControl',]
temmicexG <- temmicex[temmicex$treatment=='Girdled',]
temmicexL <- temmicex[temmicex$treatment=='Logged',]
temmicexHw <- temmicex[temmicex$treatment=='HardwoodControl',]
temmicmain <- Temnothorax[Temnothorax$moose.cage=='Exterior',]
temmicmainHe <- temmicmain[temmicmain$treatment=='HemlockControl',]
temmicmainG <- temmicmain[temmicmain$treatment=='Girdled',]
temmicmainL <- temmicmain[temmicmain$treatment=='Logged',]
temmicmainHw <- temmicmain[temmicmain$treatment=='HardwoodControl',]

years <- c(2005:2011)
years2 <- c(2012:2015)

# Begin plotting genera over time
par(mfrow=c(3,4))
plot(years, (aphHe$incidence/He$incidence)*100, type='line', ylim=c(0,100), xlim=c(2005,2015), xlab='Year', ylab='Relative Abundance', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='Hemlock')
lines(years2, (aphexHe$incidence/Heex$incidence)*100,lty='dashed', lwd=3, col='red')
lines(years2, (aphmainHe$incidence/Hemain$incidence)*100, lwd=3, col='red')
lines(years, (camHe$incidence/He$incidence)*100,lwd=3, col='orange')
lines(years2, (c(0,0,cammainHe$incidence)/Hemain$incidence)*100, lwd=3, col='orange')
lines(years2, (camexHe$incidence/Heex$incidence)*100,lty='dashed', lwd=3, col='orange')
lines(years, (rep(0,7)/He$incidence)*100,lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='skyblue')
lines(years, (c(formHe$incidence,rep(0,6))/He$incidence)*100,lwd=3, col='green')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='green')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='green')
lines(years, (c(1,3,2,1,0,0,2)/He$incidence)*100,lwd=3, col='purple')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='purple')
lines(years2, (c(0,1,2,0)/Heex$incidence)*100,lty='dashed', lwd=3, col='purple')
lines(years, (rep(0,7)/He$incidence)*100,lwd=3, col='magenta')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='magenta')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='magenta')
lines(years, (c(1,rep(0,5),4)/He$incidence)*100,lwd=3, col='brown')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='brown')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='brown')
lines(years, (c(2,0,2,rep(0,4))/He$incidence)*100,lwd=3, col='darkgreen')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='darkgreen')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='darkgreen')
lines(years, (c(2,0,2,rep(0,4))/He$incidence)*100,lwd=3, col='black')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='black')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='black')
lines(years, (c(2,0,2,rep(0,4))/He$incidence)*100,lwd=3, col='pink')
lines(years2, (rep(0,4)/Hemain$incidence)*100, lwd=3, col='pink')
lines(years2, (rep(0,4)/Heex$incidence)*100,lty='dashed', lwd=3, col='pink')
mtext('a)',3,adj=0, cex=1.5, line=2)

plot(years, (aphG$incidence/G$incidence)*100, type='line', ylim=c(0,100), xlim=c(2005,2015), xlab='Year', ylab='', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='Girdled')
lines(years2, (aphexG$incidence/Gex$incidence)*100,lty='dashed', lwd=3, col='red')
lines(years2, (aphmainG$incidence/Gmain$incidence)*100, lwd=3, col='red')
lines(years, (camG$incidence/G$incidence)*100,lwd=3, col='orange')
lines(years2, (cammainG$incidence/Gmain$incidence)*100, lwd=3, col='orange')
lines(years2, (camexG$incidence/Gex$incidence)*100,lty='dashed', lwd=3, col='orange')
lines(years, (rep(0,7)/G$incidence)*100,lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Gmain$incidence)*100, lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Gex$incidence)*100,lty='dashed', lwd=3, col='skyblue')
lines(years, (c(formG$incidence,rep(0,6))/G$incidence)*100,lwd=3, col='green')
lines(years2, (formmainG$incidence/Gmain$incidence)*100, lwd=3, col='green')
lines(years2, (formexG$incidence/Gex$incidence)*100,lty='dashed', lwd=3, col='green')
lines(years, (c(1,1,3,0,2,0,0)/G$incidence)*100,lwd=3, col='purple')
lines(years2, (c(0,2,4,0)/Gmain$incidence)*100, lwd=3, col='purple')
lines(years2, (c(0,1,3,2)/Gex$incidence)*100,lty='dashed', lwd=3, col='purple')
lines(years, (rep(0,7)/G$incidence)*100,lwd=3, col='magenta')
lines(years2, (rep(0,4)/Gmain$incidence)*100, lwd=3, col='magenta')
lines(years2, (c(rep(0,3),1)/Gex$incidence)*100,lty='dashed', lwd=3, col='magenta')
lines(years, (c(0,1,0,1,2,0,4)/G$incidence)*100,lwd=3, col='brown')
lines(years2, (myrmicmainG$incidence/Gmain$incidence)*100, lwd=3, col='brown')
lines(years2, (myrmicexG$incidence/Gex$incidence)*100,lty='dashed', lwd=3, col='brown')
lines(years, (c(steG$incidence,rep(0,3))/G$incidence)*100,lwd=3, col='darkgreen')
lines(years2, (c(1,rep(0,3))/Gmain$incidence)*100, lwd=3, col='darkgreen')
lines(years2, (c(1,0,0,1)/Gex$incidence)*100,lty='dashed', lwd=3, col='darkgreen')
lines(years, (rep(0,7)/G$incidence)*100,lwd=3, col='black')
lines(years2, (rep(0,4)/Gmain$incidence)*100, lwd=3, col='black')
lines(years2, (rep(0,4)/Gex$incidence)*100,lty='dashed', lwd=3, col='black')
lines(years, (c(0,3,0,1,0,0,0)/G$incidence)*100,lwd=3, col='pink')
lines(years2, (c(3,4,0,4)/Gmain$incidence)*100, lwd=3, col='pink')
lines(years2, (c(2,rep(0,3))/Gex$incidence)*100,lty='dashed', lwd=3, col='pink')
mtext('b)',3,adj=0, cex=1.5, line=2)

plot(years, (aphL$incidence/L$incidence)*100, type='line', ylab='Relative abundnace', ylim=c(0,100), xlim=c(2005,2015), xlab='Year', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='Logged')
lines(years2, (aphexL$incidence/Lex$incidence)*100,lty='dashed', lwd=3, col='red')
lines(years2, (aphmainL$incidence/Lmain$incidence)*100, lwd=3, col='red')
lines(years, (camL$incidence/L$incidence)*100,lwd=3, col='orange')
lines(years2, (cammainL$incidence/Lmain$incidence)*100, lwd=3, col='orange')
lines(years2, (camexL$incidence/Lex$incidence)*100,lty='dashed', lwd=3, col='orange')
lines(years, (rep(0,7)/L$incidence)*100,lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Lmain$incidence)*100, lwd=3, col='skyblue')
lines(years2, (c(0,1,0,1)/Lex$incidence)*100,lty='dashed', lwd=3, col='skyblue')
lines(years, (formL$incidence/L$incidence)*100,lwd=3, col='green')
lines(years2, (formmainL$incidence/Lmain$incidence)*100, lwd=3, col='green')
lines(years2, (formexL$incidence/Lex$incidence)*100,lty='dashed', lwd=3, col='green')
lines(years, (lasL$incidence/L$incidence)*100,lwd=3, col='purple')
lines(years2, (lasmainL$incidence/Lmain$incidence)*100, lwd=3, col='purple')
lines(years2, (lasexL$incidence/Lex$incidence)*100,lty='dashed', lwd=3, col='purple')
lines(years, (rep(0,7)/L$incidence)*100,lwd=3, col='magenta')
lines(years2, (rep(0,4)/Lmain$incidence)*100, lwd=3, col='magenta')
lines(years2, (c(0,2,0,1)/Lex$incidence)*100,lty='dashed', lwd=3, col='magenta')
lines(years, (c(0,1,1,0,1,2,6)/L$incidence)*100,lwd=3, col='brown')
lines(years2, (myrmicmainL$incidence/Lmain$incidence)*100, lwd=3, col='brown')
lines(years2, (myrmicexL$incidence/Lex$incidence)*100,lty='dashed', lwd=3, col='brown')
lines(years, (c(1,rep(0,5),2)/L$incidence)*100,lwd=3, col='darkgreen')
lines(years2, (rep(0,4)/Lmain$incidence)*100, lwd=3, col='darkgreen')
lines(years2, (c(0,3,2,0)/Lex$incidence)*100,lty='dashed', lwd=3, col='darkgreen')
lines(years, (c(rep(0,5),1,1)/L$incidence)*100,lwd=3, col='black')
lines(years2, (c(2,2,0,0)/Lmain$incidence)*100, lwd=3, col='black')
lines(years2, (rep(0,4)/Lex$incidence)*100,lty='dashed', lwd=3, col='black')
lines(years, (c(1,2,0,0,1,1,0)/L$incidence)*100,lwd=3, col='pink')
lines(years2, (c(0,3,0,2)/Lmain$incidence)*100, lwd=3, col='pink')
lines(years2, (c(0,1,0,0)/Lex$incidence)*100,lty='dashed', lwd=3, col='pink')
mtext('c)',3,adj=0, cex=1.5, line=2)

plot(years, (c(aphHw$incidence)/Hw$incidence)*100, type='line', ylim=c(0,100), xlim=c(2005,2015), xlab='Year', ylab='', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='Hardwood')
lines(years2, (aphexHw$incidence/Hwex$incidence)*100,lty='dashed', lwd=3, col='red')
lines(years2, (aphmainHw$incidence/Hwmain$incidence)*100, lwd=3, col='red')
lines(years, (c(camHw$incidence,0)/Hw$incidence)*100,lwd=3, col='orange')
lines(years2, (cammainHw$incidence/Hwmain$incidence)*100, lwd=3, col='orange')
lines(years2, (camexHw$incidence/Hwex$incidence)*100,lty='dashed', lwd=3, col='orange')
lines(years, (rep(0,7)/Hw$incidence)*100,lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Hwmain$incidence)*100, lwd=3, col='skyblue')
lines(years2, (rep(0,4)/Hwex$incidence)*100,lty='dashed', lwd=3, col='skyblue')
lines(years, (c(formHw$incidence,0)/Hw$incidence)*100,lwd=3, col='green')
lines(years2, (formmainHw$incidence/Hwmain$incidence)*100, lwd=3, col='green')
lines(years2, (c(formexHw$incidence,1)/Hwex$incidence)*100,lty='dashed', lwd=3, col='green')
lines(years, (c(lasHw$incidence,0)/Hw$incidence)*100,lwd=3, col='purple')
lines(years2, (lasmainHw$incidence/Hwmain$incidence)*100, lwd=3, col='purple')
lines(years2, (c(0,0,lasexHw$incidence)/Hwex$incidence)*100,lty='dashed', lwd=3, col='purple')
lines(years, (rep(0,7)/Hw$incidence)*100,lwd=3, col='magenta')
lines(years2, (rep(0,4)/Hwmain$incidence)*100, lwd=3, col='magenta')
lines(years2, (rep(0,4)/Hwex$incidence)*100,lty='dashed', lwd=3, col='magenta')
lines(years, (c(myrmicHw$incidence,0)/Hw$incidence)*100,lwd=3, col='brown')
lines(years2, (c(3,0,2,1)/Hwmain$incidence)*100, lwd=3, col='brown')
lines(years2, (myrmicexHw$incidence/Hwex$incidence)*100,lty='dashed', lwd=3, col='brown')
lines(years, (c(1,2,0,2,1,0,0)/Hw$incidence)*100,lwd=3, col='darkgreen')
lines(years2, (rep(0,4)/Hwmain$incidence)*100, lwd=3, col='darkgreen')
lines(years2, (c(0,1,0,0)/Hwex$incidence)*100,lty='dashed', lwd=3, col='darkgreen')
lines(years, (c(0,1,rep(0,5))/Hw$incidence)*100,lwd=3, col='black')
lines(years2, (c(0,0,0,0)/Hwmain$incidence)*100, lwd=3, col='black')
lines(years2, (rep(0,4)/Hwex$incidence)*100,lty='dashed', lwd=3, col='black')
lines(years, (c(1,1,0,1,0,0,0)/Hw$incidence)*100,lwd=3, col='pink')
lines(years2, (c(0,0,0,0)/Hwmain$incidence)*100, lwd=3, col='pink')
lines(years2, (c(1,0,0,0)/Hwex$incidence)*100,lty='dashed', lwd=3, col='pink')
mtext('d)',3,adj=0, cex=1.5, line=2)

# plot genera of BRF ants
#Select Hand & litter collections
antshand <- subset(BRF1, trap.type=="Hand")
antslitter <- subset(BRF1, trap.type=="Litter")
ants1 <- rbind(antshand, antslitter)

##Clean up data
ants1<-ants1[ants1$plot!="A5",]
ants1<-ants1[ants1$plot!="A6",]
ants1<-ants1[ants1$plot!="B6",]
ants1<-ants1[ants1$plot!="B7",]
ants1<-ants1[ants1$plot!="B5",]
ants1<-ants1[ants1$plot!="Z5",]
ants1<-ants1[ants1$plot!="GX",]
ants1<-ants1[ants1$plot!="GC",]
year <- year(ants1$date)
ants1<-cbind(year,ants1)

# Create a vector of ones to append to each record, this will be eventually summed for each species.
incidence <- rep(1,dim(ants1)[1])
# append incidence vector to data frame
ants1<-cbind(ants1, incidence)
antspost <- ants1[ants1$year>2008,]

BRFgenus <- summarise(group_by(antspost, year, treatment, deer, genus), 
                      incidence = sum(incidence)
) 

# determine top ten genera
BRFtopten <- summarise(group_by(BRFgenus, genus), 
                       incidence = sum(incidence)
) 
BRFgenus <- BRFgenus[BRFgenus$year>2007,]

#Get incidence totals for different treatments to calculate relative abundances
OG <- BRFgenus[BRFgenus$treatment=='O',]
OGex <- OG[OG$deer=='exclosure',]
OGex <- summarise(group_by(OGex, year), incidence = sum(incidence)) 
OGmain <- OG[OG$deer=='main',]
OGmain <- summarise(group_by(OGmain, year), incidence = sum(incidence)) 
NO <- BRFgenus[BRFgenus$treatment=='N',]
NOex <- NO[NO$deer=='exclosure',]
NOex <- summarise(group_by(NOex, year), incidence = sum(incidence)) 
NOmain <- NO[NO$deer=='main',]
NOmain <- summarise(group_by(NOmain, year), incidence = sum(incidence)) 
C <- BRFgenus[BRFgenus$treatment=='C',]
Cex <- C[C$deer=='exclosure',]
Cex <- summarise(group_by(Cex, year), incidence = sum(incidence)) 
Cmain <- C[C$deer=='main',]
Cmain <- summarise(group_by(Cmain, year), incidence = sum(incidence)) 
O50 <- BRFgenus[BRFgenus$treatment=='O50',]
O50ex <- O50[O50$deer=='exclosure',]
O50ex <- summarise(group_by(O50ex, year), incidence = sum(incidence)) 
O50main <- O50[O50$deer=='main',]
O50main <- summarise(group_by(O50main, year), incidence = sum(incidence)) 

# Select data for top ten genera
aphaenogaster <- BRFgenus[BRFgenus$genus=='Aphaenogaster',]
aphex <- aphaenogaster[aphaenogaster$deer=='exclosure',]
aphexOG <- aphex[aphex$treatment=='O',]
aphexC <- aphex[aphex$treatment=='C',]
aphexO50 <- aphex[aphex$treatment=='O50',]
aphexNO <- aphex[aphex$treatment=='N',]
aphmain <- aphaenogaster[aphaenogaster$deer=='main',]
aphmainOG <- aphmain[aphmain$treatment=='O',]
aphmainC <- aphmain[aphmain$treatment=='C',]
aphmainO50 <- aphmain[aphmain$treatment=='O50',]
aphmainNO <- aphmain[aphmain$treatment=='N',]

brachymyrmex <- BRFgenus[BRFgenus$genus=='Brachymyrmex',]
braex <- brachymyrmex[brachymyrmex$deer=='exclosure',]
braexOG <- braex[braex$treatment=='O',]
braexC <- braex[braex$treatment=='C',]
braexO50 <- braex[braex$treatment=='O50',]
braexNO <- braex[braex$treatment=='N',]
bramain <- brachymyrmex[brachymyrmex$deer=='main',]
bramainOG <- bramain[bramain$treatment=='O',]
bramainC <- bramain[bramain$treatment=='C',]
bramainO50 <- bramain[bramain$treatment=='O50',]
bramainNO <- bramain[bramain$treatment=='N',]

camponotus <- BRFgenus[BRFgenus$genus=='Camponotus',]
camex <- camponotus[camponotus$deer=='exclosure',]
camexOG <- camex[camex$treatment=='O',]
camexC <- camex[camex$treatment=='C',]
camexO50 <- camex[camex$treatment=='O50',]
camexNO <- camex[camex$treatment=='N',]
cammain <- camponotus[camponotus$deer=='main',]
cammain <- camponotus[camponotus$deer=='main',]
cammainOG <- cammain[cammain$treatment=='O',]
cammainC <- cammain[cammain$treatment=='C',]
cammainO50 <- cammain[cammain$treatment=='O50',]
cammainNO <- cammain[cammain$treatment=='N',]

formica <- BRFgenus[BRFgenus$genus=='Formica',]
forex <- formica[formica$deer=='exclosure',]
forexOG <- forex[forex$treatment=='O',]
forexC <- forex[forex$treatment=='C',]
forexO50 <- forex[forex$treatment=='O50',]
forexNO <- forex[forex$treatment=='N',]
formain <- formica[formica$deer=='main',]
formainOG <- formain[formain$treatment=='O',]
formainC <- formain[formain$treatment=='C',]
formainO50 <- formain[formain$treatment=='O50',]
formainNO <- formain[formain$treatment=='N',]

lasius <- BRFgenus[BRFgenus$genus=='Lasius',]
lasex <- lasius[lasius$deer=='exclosure',]
lasexOG <- lasex[lasex$treatment=='O',]
lasexC <- lasex[lasex$treatment=='C',]
lasexO50 <- lasex[lasex$treatment=='O50',]
lasexNO <- lasex[lasex$treatment=='N',]
lasmain <- lasius[lasius$deer=='main',]
lasmainOG <- lasmain[lasmain$treatment=='O',]
lasmainC <- lasmain[lasmain$treatment=='C',]
lasmainO50 <- lasmain[lasmain$treatment=='O50',]
lasmainNO <- lasmain[lasmain$treatment=='N',]

myrmica <- BRFgenus[BRFgenus$genus=='Myrmica',]
myrex <- myrmica[myrmica$deer=='exclosure',]
myrexOG <- myrex[myrex$treatment=='O',]
myrexC <- myrex[myrex$treatment=='C',]
myrexO50 <- myrex[myrex$treatment=='O50',]
myrexNO <- myrex[myrex$treatment=='N',]
myrmain <- myrmica[myrmica$deer=='main',]
myrmainOG <- myrmain[myrmain$treatment=='O',]
myrmainC <- myrmain[myrmain$treatment=='C',]
myrmainO50 <- myrmain[myrmain$treatment=='O50',]
myrmainNO <- myrmain[myrmain$treatment=='N',]

ponera <- BRFgenus[BRFgenus$genus=='Ponera',]
ponex <- ponera[ponera$deer=='exclosure',]
ponexOG <- ponex[ponex$treatment=='O',]
ponexC <- ponex[ponex$treatment=='C',]
ponexO50 <- ponex[ponex$treatment=='O50',]
ponexNO <- ponex[ponex$treatment=='N',]
ponmain <- ponera[ponera$deer=='main',]
ponmainOG <- ponmain[ponmain$treatment=='O',]
ponmainC <- ponmain[ponmain$treatment=='C',]
ponmainO50 <- ponmain[ponmain$treatment=='O50',]
ponmainNO <- ponmain[ponmain$treatment=='N',]

pre <- BRFgenus[BRFgenus$genus=='Prenolepis',]
preex <- pre[pre$deer=='exclosure',]
preexOG <- preex[preex$treatment=='O',]
preexC <- preex[preex$treatment=='C',]
preexO50 <- preex[preex$treatment=='O50',]
preexNO <- preex[preex$treatment=='N',]
premain <- pre[pre$deer=='main',]
premainOG <- premain[premain$treatment=='O',]
premainC <- premain[premain$treatment=='C',]
premainO50 <- premain[premain$treatment=='O50',]
premainNO <- premain[premain$treatment=='N',]

tap <- BRFgenus[BRFgenus$genus=='Tapinoma',]
tapex <- tap[tap$deer=='exclosure',]
tapexOG <- tapex[tapex$treatment=='O',]
tapexC <- tapex[tapex$treatment=='C',]
tapexO50 <- tapex[tapex$treatment=='O50',]
tapexNO <- tapex[tapex$treatment=='N',]
tapmain <- tap[tap$deer=='main',]
tapmainOG <- tapmain[tapmain$treatment=='O',]
tapmainC <- tapmain[tapmain$treatment=='C',]
tapmainO50 <- tapmain[tapmain$treatment=='O50',]
tapmainNO <- tapmain[tapmain$treatment=='N',]

tem <- BRFgenus[BRFgenus$genus=='Temnothorax',]
temex <- tem[tem$deer=='exclosure',]
temexOG <- temex[temex$treatment=='O',]
temexC <- temex[temex$treatment=='C',]
temexO50 <- temex[temex$treatment=='O50',]
temexNO <- temex[temex$treatment=='N',]
temmain <- tem[tem$deer=='main',]
temmainOG <- temmain[temmain$treatment=='O',]
temmainC <- temmain[temmain$treatment=='C',]
temmainO50 <- temmain[temmain$treatment=='O50',]
temmainNO <- temmain[temmain$treatment=='N',]

range(BRFgenus$incidence)

BRFtoptengen <- c("Aphaenogaster","Brachymyrmex","Camponotus","Formica","Lasius","Myrmica","Ponera", "Prenolepis","Tapinoma","Temnothorax")
years <- c(2009,2011,2015)

plot(years, (aphmainC$incidence/Cmain$incidence)*100, type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year', ylab='Relative Abundance', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main="Oak Control")
lines(years, (bramainC$incidence/Cmain$incidence)*100, lwd=3, col='blue')
lines(years, (cammainC$incidence/Cmain$incidence)*100, lwd=3, col='orange')
lines(years, (formainC$incidence/Cmain$incidence)*100, lwd=3, col='green')
lines(years, (lasmainC$incidence/Cmain$incidence)*100, lwd=3, col='purple')
lines(years, (c(1,0,1)/Cmain$incidence)*100, lwd=3, col='brown')
lines(years, (c(ponmainC$incidence,0,0)/Cmain$incidence)*100, lwd=3, col='yellow')
lines(years, (c(0,premainC$incidence)/Cmain$incidence)*100, lwd=3, col='grey')
lines(years, (c(0,0,tapmainC$incidence)/Cmain$incidence)*100, lwd=3, col='black')
lines(years, (c(0,temmainC$incidence)/Cmain$incidence)*100, lwd=3, col='pink')
mtext('e)',3,adj=0, cex=1.5, line=2)

plot(years, (aphmainO50$incidence/O50main$incidence)*100, type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year',  cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main="Oak 50% Girdled", ylab='')
lines(years, (bramainO50$incidence/O50main$incidence)*100, lwd=3, col='blue')
lines(years, (c(0,cammainO50$incidence)/O50main$incidence)*100, lwd=3, col='orange')
lines(years, (formainO50$incidence/O50main$incidence)*100, lwd=3, col='green')
lines(years, (lasmainO50$incidence/O50main$incidence)*100, lwd=3, col='purple')
lines(years, (c(0,myrmainO50$incidence)/O50main$incidence)*100, lwd=3, col='brown')
lines(years, (ponmainO50$incidence/O50main$incidence)*100, lwd=3, col='yellow')
lines(years, (c(premainO50$incidence,0)/O50main$incidence)*100, lwd=3, col='grey')
lines(years, (rep(0,3)/O50main$incidence)*100, lwd=3, col='black')
lines(years, (c(0,temmainOG$incidence)/O50main$incidence)*100, lwd=3, col='pink')
mtext('f)',3,adj=0, cex=1.5, line=2)

plot(years, (aphmainOG$incidence/OGmain$incidence)*100, type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year',  main='100% Oak Girdled', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', ylab='')
lines(years, (bramainOG$incidence/OGmain$incidence)*100, lwd=3, col='blue')
lines(years, (cammainOG$incidence/OGmain$incidence)*100, lwd=3, col='orange')
lines(years, (formainOG$incidence/OGmain$incidence)*100, lwd=3, col='green')
lines(years, (lasmainOG$incidence/OGmain$incidence)*100, lwd=3, col='purple')
lines(years, (c(0,myrmainOG$incidence)/OGmain$incidence)*100, lwd=3, col='brown')
lines(years, (ponmainOG$incidence/OGmain$incidence)*100, lwd=3, col='yellow')
lines(years, (c(0,premainOG$incidence)/OGmain$incidence)*100, lwd=3, col='grey')
lines(years, (c(0,tapmainOG$incidence,0)/OGmain$incidence)*100, lwd=3, col='black')
lines(years, (c(0,temmainOG$incidence)/OGmain$incidence)*100, lwd=3, col='pink')
mtext('g)',3,adj=0, cex=1.5, line=2)

plot(years, (aphmainNO$incidence/NOmain$incidence)*100, type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year',  main='Non-Oaks Girdled', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', ylab='')
lines(years, (bramainNO$incidence/NOmain$incidence)*100, lwd=3, col='blue')
lines(years, (cammainNO$incidence/NOmain$incidence)*100, lwd=3, col='orange')
lines(years, (formainNO$incidence/NOmain$incidence)*100, lwd=3, col='green')
lines(years, (lasmainNO$incidence/NOmain$incidence)*100, lwd=3, col='purple')
lines(years, (c(0,0,myrmainNO$incidence)/NOmain$incidence)*100, lwd=3, col='brown')
lines(years, (ponmainNO$incidence/NOmain$incidence)*100, lwd=3, col='yellow')
lines(years, (c(premainNO$incidence,0)/NOmain$incidence)*100, lwd=3, col='grey')
lines(years, (tapmainNO$incidence/NOmain$incidence)*100, lwd=3, col='black')
lines(years, (c(0,temmainNO$incidence)/NOmain$incidence)*100, lwd=3, col='pink')
mtext('h)',3,adj=0, cex=1.5, line=2)
######

plot(years, (aphexC$incidence/Cex$incidence)*100, type='line', lty='dashed', ylim=c(0,70), xlim=c(2009,2015), xlab='Year', ylab='Relative Abundance', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='Oak Control')
lines(years, (braexC$incidence/Cex$incidence)*100, lty='dashed', lwd=3, col='blue')
lines(years, (c(camexC$incidence,0)/Cex$incidence)*100, lty='dashed', lwd=3, col='orange')
lines(years, (forexC$incidence/Cex$incidence)*100,lty='dashed', lwd=3, col='green')
lines(years, (lasexC$incidence/Cex$incidence)*100, lty='dashed',lwd=3, col='purple')
lines(years, (c(0,0,myrexC$incidence)/Cex$incidence)*100,lty='dashed', lwd=3, col='brown')
lines(years, (c(0,0,0)/Cex$incidence)*100, lwd=3, col='yellow',lty='dashed')
lines(years, (c(0,preexC$incidence,0)/Cex$incidence)*100,lty='dashed', lwd=3, col='grey')
lines(years, (c(0,0,tapexC$incidence)/Cex$incidence)*100, lty='dashed',lwd=3, col='black')
lines(years, (c(0,temexC$incidence)/Cex$incidence)*100, lty='dashed',lwd=3, col='pink')
mtext('i)',3,adj=0, cex=1.5, line=2)

plot(years, (aphexO50$incidence/O50ex$incidence)*100, lty='dashed', type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year', ylab='Incidences', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='50% Oak Girdled')
lines(years, (braexO50$incidence/O50ex$incidence)*100, lty='dashed', lwd=3, col='blue')
lines(years, (c(0,0,0)/O50ex$incidence)*100, lty='dashed', lwd=3, col='orange')
lines(years, (forexO50$incidence/O50ex$incidence)*100, lty='dashed', lwd=3, col='green')
lines(years, (lasexO50$incidence/O50ex$incidence)*100, lty='dashed', lwd=3, col='purple')
lines(years, (c(0,myrexO50$incidence)/O50ex$incidence)*100, lty='dashed', lwd=3, col='brown')
lines(years, (ponexO50$incidence/O50ex$incidence)*100, lty='dashed', lwd=3, col='yellow')
lines(years, (c(preexO50$incidence,0,0)/O50ex$incidence)*100, lty='dashed', lwd=3, col='grey')
lines(years, (rep(0,3)/O50ex$incidence)*100, lty='dashed', lwd=3, col='black')
lines(years, (c(0,0,0)/O50ex$incidence)*100, lty='dashed', lwd=3, col='pink')
mtext('j)',3,adj=0, cex=1.5, line=2)

plot(years, (aphexOG$incidence/OGex$incidence)*100, lty='dashed', type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year', ylab='' , cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='100% Oak Girdled')
lines(years, (braexOG$incidence/OGex$incidence)*100, lty='dashed', lwd=3, col='blue')
lines(years, (c(1,0,1)/OGex$incidence)*100, lty='dashed', lwd=3, col='orange')
lines(years, (forexOG$incidence/OGex$incidence)*100, lty='dashed', lwd=3, col='green')
lines(years, (lasexOG$incidence/OGex$incidence)*100, lty='dashed', lwd=3, col='purple')
lines(years, (c(0,myrexOG$incidence,0)/OGex$incidence)*100, lty='dashed', lwd=3, col='brown')
lines(years, (ponexOG$incidence/OGex$incidence)*100, lty='dashed', lwd=3, col='yellow')
lines(years, (c(0,0,0)/OGex$incidence)*100, lty='dashed', lwd=3, col='grey')
lines(years, (c(0,0,0)/OGex$incidence)*100, lty='dashed', lwd=3, col='black')
lines(years, (c(0,0,temexOG$incidence)/OGex$incidence)*100, lty='dashed', lwd=3, col='pink')
mtext('k)',3,adj=0, cex=1.5, line=2)

plot(years, (aphexNO$incidence/NOex$incidence)*100, lty='dashed', type='line', ylim=c(0,70), xlim=c(2009,2015), xlab='Year', ylab='', cex.axis=1.5, cex.main=1.5, cex.lab=1.5, lwd=3, col='red', main='Non-Oaks Girdled')
lines(years, (braexNO$incidence/NOex$incidence)*100, lty='dashed', lwd=3, col='blue')
lines(years, (c(camexNO$incidence,0)/NOex$incidence)*100, lty='dashed', lwd=3, col='orange')
lines(years, (forexNO$incidence/NOex$incidence)*100, lty='dashed', lwd=3, col='green')
lines(years, (lasexNO$incidence/NOex$incidence)*100, lty='dashed', lwd=3, col='purple')
lines(years, (c(0,0,0)/NOex$incidence)*100, lty='dashed', lwd=3, col='brown')
lines(years, (rep(0,3)/NOex$incidence)*100, lty='dashed', lwd=3, col='yellow')
lines(years, (c(2,0,2)/NOex$incidence)*100, lty='dashed', lwd=3, col='grey')
lines(years, (tapexNO$incidence/NOex$incidence)*100, lty='dashed', lwd=3, col='black')
lines(years, (c(1,0,1)/NOex$incidence)*100, lty='dashed', lwd=3, col='pink')
mtext('l)',3,adj=0, cex=1.5, line=2)

par(mfrow=c(1,1))
plot.new()
legend(locator(1), legend=c('Aphaenogaster','Brachmyrmex','Camponotus', 'Crematogaster','Formica','Lasius', 'Myrmecina','Myrmica','Ponera','Prenolepis', 'Stenamma','Tapinoma','Temnothorax'), fill=c("red","blue","orange","skyblue","green", "purple",'magenta',"brown","yellow", "grey",'darkgreen', "black","pink"), bty='n' , cex=1.3)

#### Prep BRF data for multivariate analyses that will be performed in PERMANOVA add-on package of PRIMER
# Check data 
testBRF <- BRF2[rowSums(BRF2[,-c(1:6)])!=0,]# Check for and remove any rows without any species found. No all zero rows found.
dim(testBRF)
dim(BRF2)
rm("testBRF")
BRF2pre1 <- subset(BRF2, year<2008) # Pre-treatment data

# Relativize frequencies of species across rows
BRF2pre <- data.stand(BRF2pre1[,6:dim(BRF2pre1)[2]], method='total', margin='row') 

# Check for multivariate outliers
mv.outliers(BRF2pre[,6:dim(BRF2pre)[2]], method='bray')# No outliers

# Check assumptions of PERMANCOVA (check for non-significant multivariate dispersion) 
anova(betadisper(vegdist(BRF2pre1[,6:dim(BRF2pre1)[2]], method='bray'),BRF2pre1$treatment)) 

# BRF post-treatment data analysis
BRF2post1 <- subset(BRF2, year>2007)# Removes 2006 & 2007, which are the control years. 
testBRFpost <- BRF2post1[rowSums(BRF2post1[,-c(1:6)])!=0,]# Check for and remove any rows without any species found. No all zero rows found.
dim(BRF2post1)
dim(testBRFpost)

BRF2postrel <- data.stand(BRF2post1[,6:dim(BRF2post1)[2]], method='total', margin='row')
mv.outliers(BRF2postrel, method='bray')# No outliers
BRF2postrel <- cbind(BRF2post1[,1:5], BRF2postrel)

# Check assumptions of PERMANCOVA (check for non-significant multivariate dispersion)
anova(betadisper(vegdist(BRF2postrel[,6:dim(BRF2postrel)[2]], method='bray'),BRF2postrel$year))
anova(betadisper(vegdist(BRF2postrel[,6:dim(BRF2postrel)[2]], method='bray'),BRF2postrel$treatment)) 
anova(betadisper(vegdist(BRF2postrel[,6:dim(BRF2postrel)[2]], method='bray'),BRF2postrel$block))
anova(betadisper(vegdist(BRF2postrel[,6:dim(BRF2postrel)[2]], method='bray'),BRF2postrel$deer)) 

#### Prep HFR data for multivariate analyses of taxonomic composition

# Check pre-treatment years to make sure there were no pre-existing difference between treatments
HFR2pre <- subset(HFR2, year<2005) # Pre-treatment data
HFR2pre <- subset(HFR2pre, treatment!='Hw') # Exclude hardwood controls as we expect them to be different pre-treatment

# Check data 
test <- HFR2pre[rowSums(HFR2pre[,6:dim(HFR2pre)[2]])==0,]# Check for and remove any rows without any species found. No all zero rows found
dim(test)

# Relativize frequencies of species across rows
HFR2prerel <- data.stand(HFR2pre[,6:dim(HFR2pre)[2]], method='total', margin='row') 
test <- HFR2prerel[rowSums(HFR2prerel)==0,]# Check for and remove any rows without any species found. No all zero rows found
dim(test)


# Make sure pre-treatment years did not differ across treatments
mv.outliers(HFR2prerel, method='bray')# No outliers

# Check assumptions of PERMANCOVA (check for non-significant multivariate dispersion) 
anova(betadisper(vegdist(HFR2prerel, method='bray'),HFR2pre$treatment)) 

# HFR post-treatment data analysis
HFR2post <- subset(HFR2, year>2004)# Removes 2003 & 2004, which are the control years. 
HFR2post <- susbset(HFR2post, moose.cage=='exterior')
HFR2preinfest <- subset(HFR2post, year<2010)
HFR2postinfest <- subset(HFR2post, year>=2010)
infest <- c(rep(0, dim(HFR2preinfest)[1]), rep(1, dim(HFR2postinfest)[1]))

#Row relativize species frequencies
HFR2postrel <- data.stand(HFR2post[,6:dim(HFR2post)[2]], method='total', margin='row')

# Check for multivariate outliers
mv.outliers(HFR2postrel, method='bray')# outlier on row 81 
HFR2postno <- cbind(HFR2post[,1:5],infest,HFR2postrel)
HFR2postno <- HFR2postno[-81,] # Remove outliers

# Check assumptions of PERMANCOVA (check for non-significant multivariate dispersion)
anova(betadisper(vegdist(HFR2postno[,6:dim(HFR2postno)[2]], method='bray'),HFR2postno$infest)) 
anova(betadisper(vegdist(HFR2postno[,6:dim(HFR2postno)[2]], method='bray'),HFR2postno$treatment)) 
anova(betadisper(vegdist(HFR2postno[,6:dim(HFR2postno)[2]], method='bray'),HFR2postno$block)) 

# Multivariate taonomic large herbivore exclosure analysis for HFR
HFRex <- subset(HFR2postno, year>2011)

# Check assumptions of PERMANCOVA (check for non-significant multivariate dispersion)
anova(betadisper(vegdist(HFRex[,6:dim(HFRex)[2]], method='bray'),HFRex$moose)) 

##NMDS for canopy treatments (Appendix S2)

#BRF ordination

# Create scree plot to find number of axes needed to reduce stress to 0.05
nmds.scree(BRF2postrel, distance='bray', k=10, trymax=200, autotransform=FALSE)
BRF_nmds <- metaMDS(BRF2postrel, distance='bray', k=8, trymax=200, autotransform=FALSE)
stressplot(BRF_nmds)

#HFR ordination 
# Create scree plot to find number of axes needed to reduce stress to 0.05
nmds.scree(HFR2postrel, distance='bray', k=10, trymax=200, autotransform=FALSE)
HFR_nmds <- metaMDS(HFR2postrel, distance='bray', k=5, trymax=200, autotransform=FALSE)
stressplot(HFR_nmds) 

par(mfrow=c(1,2))

p<- ordiplot(HFR_nmds, choices=c(1,2), display='sites',type='n', main="Harvard Forest Canopy Treatments", xlab="NMDS1", ylab="NMDS2", xaxt='n', yaxt='n', cex.lab=1.2)

points(p$sites[1:29,],col="goldenrod4", pch=16, cex=1.2)
points(p$sites[30:59,], col="darkgreen", pch=16, cex=1.2)
points(p$sites[60:89,], col="purple", pch=16, cex=1.2)
points(p$sites[90:119,], col="red", pch=16, cex=1.2)

ordiellipse(p, groups=HFR2postno$treatment, lty=1, col="darkgreen", show.groups="He", kind='se', lwd=3) 
ordiellipse(p, groups=HFR2postno$treatment, lty=1, col="goldenrod4", show.groups="G", kind='se', lwd=3) 
ordiellipse(p, groups=HFR2postno$treatment, lty=1, col="red", show.groups="L", kind='se', lwd=3) 
ordiellipse(p, groups=HFR2postno$treatment, lty=1, col="purple",show.groups="Hw", kind='se', lwd=3) 

legend(locator(1), legend=c("Hardwood","Hemlock", "Girdled","Logged"), fill=c("purple","darkgreen", "goldenrod4","red"), bty='n' , cex=1.3)
mtext("a)", side=3, at=-1, line=1, cex=2)

l<- ordiplot(BRF_nmds, cex.lab=1.2,choices=c(1,2), display='sites',type='n', main="Black Rock Forest Canopy Treatments", xlab="NMDS1", ylab="NMDS2", xaxt='n', yaxt='n')
points(l$sites[1:18,],col="salmon", pch=16, cex=1.2)
points(l$sites[19:36,], col="darkolivegreen3", pch=16, cex=1.2)
points(l$sites[37:54,], col="cyan", pch=16, cex=1.2)
points(l$sites[55:72,], col="darkorchid2", pch=16, cex=1.2)

ordiellipse(l, groups=BRF2post1$treatment, lty=1, col="darkolivegreen3", show.groups="NO", kind='se', lwd=3) 
ordiellipse(l, groups=BRF2post1$treatment, lty=1, col="salmon", show.groups="C", kind='se', lwd=3) 
ordiellipse(l, groups=BRF2post1$treatment, lty=1, col="darkorchid2", show.groups="OG", kind='se', lwd=3) 
ordiellipse(l, groups=BRF2post1$treatment, lty=1, col="cyan",show.groups="O50", kind='se', lwd=3) 

legend(locator(1), legend=c("Control","50% Oak girdled", "100% Oak girdled","Non-oaks girdled"), fill=c("salmon","cyan","darkorchid2", "darkolivegreen3"), bty='n' , cex=1.3)
mtext("b)", side=3, at=-1.2, line=1, cex=2)

##Deer Exclosure NMDS (Appendix S2)

#Fit HFRex NMDS
# Create scree plot to find number of axes needed to reduce stress to 0.05
nmds.scree(HFRex[,6:dim(HFRex)[2]], distance='bray', k=10, trymax=200, autotransform=FALSE)#  k=8 for 0.05 stress
HFRex_nmds <- metaMDS(HFRex[,6:dim(HFRex)[2]], distance='bray', k=7, trymax=200, autotransform=FALSE)
stressplot(HFRex_nmds) #The linear fit is only 0.984, non-metric fit is .998 

#Fit BRF NMDS
nmds.scree(BRF2post1exrel, distance='bray', k=10, trymax=200, autotransform=FALSE)#  k=8 for 0.05 stress
BRFex_nmds <- metaMDS(BRF2post1exrel, distance='bray', k=8, trymax=200, autotransform=FALSE)
stressplot(BRFex_nmds) #The linear fit is only 0.974, non-metric fit is .998 

par(mfrow=c(1,2))

p<- ordiplot(HFRex_nmds, choices=c(1,2), display='sites',type='n', main="Harvard Forest Exclosure Treatments", cex.lab=1.2, ylab="NMDS2", xlab="NMDS1", xaxt='n', yaxt='n')
points(p$sites[1:31,],col="gray", pch=16, cex=1.2)
points(p$sites[32:63,], col="black", pch=16, cex=1.2)

ordiellipse(p, groups=HFRex$moose, lty=1, col="gray", show.groups="exterior", kind='se', lwd=3) 
ordiellipse(p, groups=HFRex$moose, lty=1, col="black", show.groups="interior", kind='se', lwd=3) 

legend(locator(1), legend=c("Exterior", "Interior"), fill=c("gray","black"), bty='n' , cex=1.3)
mtext("a)", side=3, at=-.7, line=1, cex=2)

l<- ordiplot(BRFex_nmds, choices=c(1,2), display='sites',type='n', main="Black Rock Forest Exclosure Treatments")
points(l$sites[1:36,],col="black", pch=16, cex=1.2)
points(l$sites[37:72,], col="gray", pch=16, cex=1.2)

ordiellipse(l, groups=BRF2post1ex$deer, lty=1, col="black", show.groups="exclosure", kind='se', lwd=3) 
ordiellipse(l, groups=BRF2post1ex$deer, lty=1, col="gray", show.groups="main", kind='se', lwd=3) 

legend(locator(1), legend=c("Exterior", "Interior"), fill=c("gray","black"), bty='n' , cex=1.3)
mtext("b)", side=3, at=-1.25, line=1, cex=2)

########## Functional Traits Analyses

library(GUniFrac)##
library(ape)##
library(picante)##
library(vegan)##
library(geiger)
library(SDMTools)
library(cluster)
library(abind)

##simes species data
simes<-read.csv("Simes_data_for_analysis.csv", header=T)
str(simes)
dim(simes)
colnames(simes)

#remove "myrsmi" and "lepcan" because we don't have traits for them
simes_red<-simes[,-which(names(simes) %in% c("myrsmi","lepcan" ))]

dim(simes_red)

#remove non-ant species columns
simes_red<-simes_red[, 6:length(simes_red)]

#transform species counts to relatve abindance

simes_r.a<-simes_red/rowSums(simes_red)

##trait data
all_traits<-read.csv("FuncTrait.csv", header=T, na.strings=c("","NA"))
str(all_traits)
summary(all_traits)

#extract traits for Simes
simes_traits<-all_traits[all_traits$Species_Code%in% colnames(simes_r.a) ,] 
dim(simes_traits)



#remove non-traits from trait data matrix and give row names

simes_traits2<-simes_traits[,5:16]


row.names(simes_traits2)<-simes_traits[,1]


#turn traits in distance matrix using the gower index.
trait.dist.mat<-as.matrix(daisy(simes_traits2,metric="gower"))


#calculate mpwd for all comparisons
obs<-comdist(simes_r.a,trait.dist.mat, abundance.weighted=T)



##### null model

#dendrogram method not using
#function to calculate mpwd on a randomized dendrogram
#comdist.shuff<-function(x){
#  as.matrix(comdist(simes_r.a,
#	cophenetic(tipShuffle(trait_dendrogram)), abundance.weighted=T))
#}


#shuffle names
comdist.shuff<-function(x){
rownames(x)<-sample(rownames(x))
as.matrix(comdist(simes_r.a,
as.matrix(daisy(x[colnames(simes_r.a),],metric="gower")), abundance.weighted=T))
}


#runs randomization to get NULL distribution of MPWD 
nulls<-replicate(999,comdist.shuff(simes_traits2))


##means for all pairS
nulls.means<-apply(nulls,c(1:2), mean, na.rm=T)

###sd for all pairwise
nulls.sds<-apply(nulls,c(1:2), sd, na.rm=T)



#STANDARDIZE observed values
ses<-(as.matrix(obs)-nulls.means)/nulls.sds

write.csv(ses, file="ses_pw_All_simes.csv")



#extract the comparisons we want (i.e., 2003 vs each year for each treatment)
#not sure what you want to do with the double count for 2012-2015 (interior/exterior)

standardized_betafd=matrix(NA, 8,16)

for(i in 1:8){
for(j in 1:16){
standardized_betafd[i,j]<-ses[i,i+j*8]
}
}

#rows are the treatments and columns are the years
standardized_betafd

write.csv(standardized_betafd, file="ses_pw_simes.csv")

raw_betafd=matrix(NA, 8,16)

for(i in 1:8){
for(j in 1:16){
raw_betafd[i,j]<-as.matrix(obs)[i,i+j*8]
}
}


#rows are the treatments and columns are the years
raw_betafd

write.csv(raw_betafd, file="raw_pw_simes.csv")





########################################################################################
########################################################################################
#BRF ANALYSIS


##BRF species data
BRF<-read.csv("BRF_data_for_analysis.csv", header=T)
str(BRF_red)
dim(BRF_red)


#remove "myrsmi" and "lepcan" because we don't have traits for them
BRF_red<-BRF[,-which(names(BRF) %in% c("myrsmi","lepcan" ))]

#remove non ant species columns
BRF_red<-BRF_red[, 6:length(BRF_red)]
dim(BRF_red)
#transform species counts to relatve abindance

BRF_r.a<-BRF_red/rowSums(BRF_red)

##trait data
all_traits<-read.csv("FuncTraitv2.csv", header=T, na.strings=c("","NA"))
str(all_traits)


#extract traits for BRF
BRF_traits<-all_traits[all_traits$Species_Code %in% colnames(BRF_r.a) ,] 

dim(BRF_traits)
#extract traits for BRF




#remove non-traits from trait data matrix and give row names

BRF_traits2<-BRF_traits[,5:16]

row.names(BRF_traits2)<-BRF_traits[,1]


#turn traits in distance matrix using the gower index.

trait.dist.mat=as.matrix(daisy(BRF_traits2,metric="gower"))


#calculate mpwd for all comparisons
obs<-comdist(BRF_r.a,trait.dist.mat, abundance.weighted=T)



##### null model
#dendrogram method, not using

#trait_dendrogram<-as.phylo(hclust(daisy(BRF_traits2,metric="gower")))


#function to calculate mpwd on a randomized dendrogram
#comdist.shuff<-function(x){
#	as.matrix(comdist(BRF_r.a,
#	cophenetic(tipShuffle(trait_dendrogram)), abundance.weighted=T))
#}


#shuffle names
comdist.shuff<-function(x){
rownames(x)<-sample(rownames(x))
as.matrix(comdist(BRF_r.a,
as.matrix(daisy(x[colnames(BRF_r.a),],metric="gower")), abundance.weighted=T))
}



#runs randomization to get NULL distribution of MPWD 
nulls<-replicate(999,comdist.shuff(BRF_traits2))


##means for all pairS
nulls.means<-apply(nulls,c(1:2), mean, na.rm=T)

###sd for all pairwise
nulls.sds<-apply(nulls,c(1:2), sd, na.rm=T)



#STANDARDIZE observed values
ses<-(as.matrix(obs)-nulls.means)/nulls.sds

write.csv(ses, file="ses_pw_All_BRF.csv")



#extract the comparisons we want (i.e., 2003 vs each year for each treatment)
#not sure what you want to do with the double count for 2012-2015 (interior/exterior)

standardized_betafd=matrix(NA, 12,7)

for(i in 1:12){
for(j in 1:7){
standardized_betafd[i,j]<-ses[i,i+j*12]
}
}


#rows are the treatments and columns are the years

standardized_betafd
write.csv(standardized_betafd, file="ses_pw_BRF.csv")

raw_betafd=matrix(NA, 12,7)

for(i in 1:12){
for(j in 1:7){
raw_betafd[i,j]<-as.matrix(obs)[i,i+j*12]
}
}



#rows are the treatments and columns are the years
raw_betafd


write.csv(raw_betafd, file="raw_pw_BRF.csv")


####################################################################################################
####################################################################################################
###################################################################################################
#NEAREST NEIGHBOR SDISTANCE


##simes species data
simes<-read.csv("Simes_data_for_analysis.csv", header=T)
str(simes)
dim(simes)
colnames(simes)

#remove "myrsmi" and "lepcan" because we don't have traits for them
simes_red<-simes[,-which(names(simes) %in% c("myrsmi","lepcan" ))]

dim(simes_red)
#remove non-ant species columns
simes_red<-simes_red[, 6:length(simes_red)]

#transform species counts to relatve abindance

simes_r.a<-simes_red/rowSums(simes_red)

##trait data
all_traits<-read.csv("FuncTraitv2.csv", header=T, na.strings=c("","NA"))
str(all_traits)
summary(all_traits)

#extract traits for Simes
simes_traits<-all_traits[all_traits$Species_Code%in% colnames(simes_r.a) ,] 
dim(simes_traits)



#remove non-traits from trait data matrix and give row names

simes_traits2<-simes_traits[,5:16]


row.names(simes_traits2)<-simes_traits[,1]


#turn traits in distance matrix using the gower index.
trait.dist.mat<-as.matrix(daisy(simes_traits2,metric="gower"))


#calculate mpwd for all comparisons
obs<-comdistnt(simes_r.a,trait.dist.mat, abundance.weighted=T)



##### null model

#dendrogram method not using
#function to calculate mpwd on a randomized dendrogram
#comdist.shuff<-function(x){
#	as.matrix(comdist(simes_r.a,
#	cophenetic(tipShuffle(trait_dendrogram)), abundance.weighted=T))
#}


#shuffle names
comdist.shuff<-function(x){
rownames(x)<-sample(rownames(x))
as.matrix(comdistnt(simes_r.a,
as.matrix(daisy(x[colnames(simes_r.a),],metric="gower")), abundance.weighted=T))
}


#runs randomization to get NULL distribution of MPWD 
nulls<-replicate(999,comdist.shuff(simes_traits2))


##means for all pairS
nulls.means<-apply(nulls,c(1:2), mean, na.rm=T)

###sd for all pairwise
nulls.sds<-apply(nulls,c(1:2), sd, na.rm=T)



#STANDARDIZE observed values
ses<-(as.matrix(obs)-nulls.means)/nulls.sds

write.csv(ses, file="ses_nn_All_simes.csv")



#extract the comparisons we want (i.e., 2003 vs each year for each treatment)
#not sure what you want to do with the double count for 2012-2015 (interior/exterior)

standardized_betafd=matrix(NA, 8,16)

for(i in 1:8){
for(j in 1:16){
standardized_betafd[i,j]<-ses[i,i+j*8]
}
}

#rows are the treatments and columns are the years
standardized_betafd

write.csv(standardized_betafd, file="ses_nn_simes.csv")

raw_betafd=matrix(NA, 8,16)

for(i in 1:8){
for(j in 1:16){
raw_betafd[i,j]<-as.matrix(obs)[i,i+j*8]
}
}


#rows are the treatments and columns are the years
raw_betafd

write.csv(raw_betafd, file="raw_nn_simes.csv")





########################################################################################
########################################################################################
#BRF ANALYSIS


##BRF species data
BRF<-read.csv("BRF_data_for_analysis.csv", header=T)
str(BRF_red)
dim(BRF_red)


#remove "myrsmi" and "lepcan" because we don't have traits for them
BRF_red<-BRF[,-which(names(BRF) %in% c("myrsmi","lepcan" ))]

#remove non ant species columns
BRF_red<-BRF_red[, 6:length(BRF_red)]
dim(BRF_red)
#transform species counts to relatve abindance

BRF_r.a<-BRF_red/rowSums(BRF_red)

##trait data
all_traits<-read.csv("FuncTraitv2.csv", header=T, na.strings=c("","NA"))
str(all_traits)


#extract traits for BRF
BRF_traits<-all_traits[all_traits$Species_Code %in% colnames(BRF_r.a) ,] 

dim(BRF_traits)
#extract traits for BRF




#remove non-traits from trait data matrix and give row names

BRF_traits2<-BRF_traits[,5:16]

row.names(BRF_traits2)<-BRF_traits[,1]


#turn traits in distance matrix using the gower index.

trait.dist.mat=as.matrix(daisy(BRF_traits2,metric="gower"))


#calculate mpwd for all comparisons
obs<-comdistnt(BRF_r.a,trait.dist.mat, abundance.weighted=T)



##### null model
#dendrogram method, not using

#trait_dendrogram<-as.phylo(hclust(daisy(BRF_traits2,metric="gower")))


#function to calculate mpwd on a randomized dendrogram
#comdist.shuff<-function(x){
#	as.matrix(comdist(BRF_r.a,
#	cophenetic(tipShuffle(trait_dendrogram)), abundance.weighted=T))
#}


#shuffle names
comdist.shuff<-function(x){
rownames(x)<-sample(rownames(x))
as.matrix(comdistnt(BRF_r.a,
as.matrix(daisy(x[colnames(BRF_r.a),],metric="gower")), abundance.weighted=T))
}



#runs randomization to get NULL distribution of MPWD 
nulls<-replicate(999,comdist.shuff(BRF_traits2))


##means for all pairS
nulls.means<-apply(nulls,c(1:2), mean, na.rm=T)

###sd for all pairwise
nulls.sds<-apply(nulls,c(1:2), sd, na.rm=T)



#STANDARDIZE observed values
ses<-(as.matrix(obs)-nulls.means)/nulls.sds

write.csv(ses, file="ses_nn_All_BRF.csv")



#extract the comparisons we want (i.e., 2003 vs each year for each treatment)
#not sure what you want to do with the double count for 2012-2015 (interior/exterior)

standardized_betafd=matrix(NA, 12,7)

for(i in 1:12){
for(j in 1:7){
standardized_betafd[i,j]<-ses[i,i+j*12]
}
}


#rows are the treatments and columns are the years

standardized_betafd
write.csv(standardized_betafd, file="ses_nn_BRF.csv")

raw_betafd=matrix(NA, 12,7)

for(i in 1:12){
for(j in 1:7){
raw_betafd[i,j]<-as.matrix(obs)[i,i+j*12]
}
}



#rows are the treatments and columns are the years
raw_betafd


write.csv(raw_betafd, file="raw_nn_BRF.csv")

# Combine files into separate Simes and BRF beta_FD.csvs

##Taxonomic and functional trait Change Analyses ANCOVAs
# Read in output files
HFR_trait <- read.csv("Simes_beta_FD.csv", header=TRUE, stringsAsFactors=FALSE)
BRF_trait <- read.csv("BRF_beta_FD.csv", header=TRUE, stringsAsFactors=FALSE)

# Load car library for extracting Type II test results from Anovas
library(car)
library(multcomp)
library(lme4)

# ANCOVA GLMM analyses for HFR
# Select post-treatment years
HFR_trait_post <- subset(HFR_trait, Year>2004)
HFR_trait_post <- subset(HFR_trait, moose.cage=='exterior')
HFRtrait_preinfest <- subset(HFR_trait_post, Year<2010)
HFRtrait_postinfest <- subset(HFR_trait_post, Year>=2010)
infest <- c(rep(0, dim(HFRtrait_preinfest)[1]), rep(1, dim(HFRtrait_postinfest)[1]))

# Graph response variables to inspect normality
mod <- lme(PW.SES ~  treatment + year + treatment*year + infest, random = ~ 1|block, data=HFR_trait_post)
plot(resid(mod)) 

#Exclosure analysis for HFR
HFRex_trait <- subset(HFR_trait_post, Year>2011)

#Compare deer exclosure treatments
mod <- lme(PW.SES ~ treatment + year + treatment/moose + treatment*year, random = ~ 1|block, data=HFRex_trait)
plot(resid(mod)) 

# HFR PWRAW
# Graph response variables to inspect normality
mod <- lme(PW.RAW ~  treatment + year + treatment*year + infest, random = ~ 1|block, data=HFR_trait_post)
plot(resid(mod)) 

#Compare deer exclosure treatments
mod <- lme(PW.RAW ~ treatment*moose, random = ~ 1|block, data=HFRex_trait)
plot(resid(mod)) 

# HFR NNRAW
# Graph response variables to inspect normality
mod <- lme(NN.RAW ~  treatment + year + treatment*year + infest, random = ~ 1|block, data=HFR_trait_post)
plot(resid(mod)) 

#Compare deer exclosure treatments
mod <- lme(NN.RAW ~ treatment*moose, random = ~ 1|block, data=HFRex_trait)
plot(resid(mod))

# HFR NNSES
# Graph response variables to inspect normality
mod <- lme(NN.SES ~  treatment + year + treatment*year + infest, random = ~ 1|block, data=HFR_trait_post)
plot(resid(mod)) 

#Compare deer exclosure treatments
mod <- lme(NN.SES ~ treatment*moose, random = ~ 1|block, data=HFRex_trait)
plot(resid(mod)) 

##### ANCOVA GLMM analyses for BRF traits
# select only post treatment data
BRF_trait <- subset(BRF_trait, year>2007)

#PW.SES
# Graph response variables to inspect normality
mod <- lme(PW.SES ~ treatment + year + treatment/moose + treatment*year, random = ~ 1|block, data=BRF_trait)
plot(resid(mod))

#PW.RAW
# Graph response variables to inspect normality
mod <- lme(PW.RAW ~ treatment + year + treatment/moose + treatment*year, random = ~ 1|block, data=BRF_trait)
plot(resid(mod))

#NN.SES
# Graph response variables to inspect normality
mod <- lme(NN.SES ~ treatment + year + treatment/moose + treatment*year, random = ~ 1|block, data=BRF_trait)
plot(resid(mod))

#PW.RAW
# Graph response variables to inspect normality
mod <- lme(NN.RAW ~ treatment + year + treatment/moose + treatment*year, random = ~ 1|block, data=BRF_trait)
plot(resid(mod))

###################### Correction for multiple tests
pvals <- read.csv("pvals.csv", header=TRUE, stringsAsFactors=FALSE) # This file incorporates p-values from PRIMER analyses, too.

# adjust pvals with holmes correction
holmpvals <- p.adjust(pvals[,4],method='holm')
pvals <- cbind(pvals, holmpvals)
write.csv(pvals, "pvals_corrected.csv")

#################### Boxplots for traits
# Read in output files
HFR_trait <- read.csv("Simes_beta_FD.csv", header=TRUE, stringsAsFactors=FALSE)
BRF_trait <- read.csv("BRF_beta_FD.csv", header=TRUE, stringsAsFactors=FALSE)

# Select post-treatment years
HFR_trait_post <- subset(HFR_trait, Year>2004)
HFR_trait_post <- subset(HFR_trait_post, moose.cage=='exterior')
HFRtrait_preinfest <- subset(HFR_trait_post, Year<2010)
HFRtrait_postinfest <- subset(HFR_trait_post, Year>=2010)
infest <- c(rep(0, dim(HFRtrait_preinfest)[1]), rep(1, dim(HFRtrait_postinfest)[1]))

par(mfrow=c(2,4))

boxplot(HFR_trait$PW.SES ~ HFR_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['pw'])), cex.lab=1.5)
mtext("a)", 3, at=0, line=1, cex=1.2)
mtext("*",3,at=2, cex=1.2)
mtext("**", 3, at=1, line=0, cex=1.2)
mtext("**", 3, at=3, line=0, cex=1.2)
mtext("***", 3, at=4, line=1.2, cex=1.2)
mtext("***", 3, at=3, line=1.2, cex=1.2)

boxplot(HFR_trait$PW.RAW ~ HFR_trait$treatment,  cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['pw'])), cex.lab=1.5)
mtext("b)", 3, at=0, line=1, cex=1.2)
mtext("*", 3, at=1, line=0,cex=1.2)
mtext("*", 3, at=2, line=0,cex=1.2)
mtext("**", 3, at=3, line=0,cex=1.2)
mtext("**", 3, at=4, line=0,cex=1.2)

boxplot(HFR_trait$NN.SES ~ HFR_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['nn'])), cex.lab=1.5)
mtext("c)", 3, at=0, line=1,cex=1.2)
mtext("*", 3, at=1, line=0,cex=1.2)
mtext("***", 3, at=1, line=1.2,cex=1.2)
mtext("**", 3, at=2, line=0,cex=1.2)
mtext("**", 3, at=3, line=0,cex=1.2)
mtext("***", 3, at=3, line=1.2,cex=1.2)
mtext("**", 3, at=4, line=0,cex=1.2)

boxplot(HFR_trait$NN.RAW ~ HFR_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['nn'])), cex.lab=1.5)
mtext("d)", 3, at=0, line=1, cex=1.2)
mtext("*", 3, at=1, line=0,cex=1.2)
mtext("**", 3, at=2, line=0,cex=1.2)
mtext("**", 3, at=3, line=0,cex=1.2)
mtext("**", 3, at=4, line=0,cex=1.2)
mtext("***", 3, at=1, line=1.2, cex=1.2)
mtext("***", 3, at=4, line=1.2, cex=1.2)
mtext("****", 3, at=1, line=2.4, cex=1.2)
mtext("****", 3, at=3, line=2.4, cex=1.2)

boxplot(BRF_trait$PW.SES~BRF_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['pw'])), cex.lab=1.5)
mtext("e)", 3, at=0, line=1, cex=1.2)

boxplot(BRF_trait$PW.RAW~BRF_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['pw'])), cex.lab=1.5)
mtext("f)", 3, at=0, line=1, cex=1.2)

boxplot(BRF_trait$NN.SES~BRF_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['nn'])), cex.lab=1.5)
mtext("g)", 3, at=0, line=1, cex=1.2)

boxplot(BRF_trait$NN.RAW~BRF_trait$treatment, cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['nn'])), cex.lab=1.5)
mtext("h)", 3, at=0, line=1, cex=1.2)


# Exclosure boxplots
#Exclosure analysis for HFR (28-29 cm wide, full height, desktop)
HFRex_trait <- subset(HFR_trait, Year>2011)

par(mfrow=c(2,4))

boxplot(HFRex_trait$PW.SES~HFRex_trait$moose, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['pw'])), cex.lab=1.5)
mtext("a)", 3, at=0.1, line=1, cex=1.2)

boxplot(HFRex_trait$PW.RAW~HFRex_trait$moose, cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['pw'])), cex.lab=1.5)
mtext("b)", 3, at=0, line=1, cex=1.2)

boxplot(HFRex_trait$NN.SES~HFRex_trait$moose, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['nn'])), cex.lab=1.5)
mtext("c)", 3, at=0, line=1, cex=1.2)

boxplot(HFRex_trait$NN.RAW~HFRex_trait$moose, cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['nn'])), cex.lab=1.5)
mtext("d)", 3, at=0, line=1, cex=1.2)

boxplot(BRF_trait$PW.SES~BRF_trait$deer2, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['pw'])), cex.lab=1.5)
mtext("e)", 3, at=0.1, line=1, cex=1.2)

boxplot(BRF_trait$PW.RAW~BRF_trait$deer2, cex=1.5, cex.axis=1.5, xlab=expression(italic('D'['pw'])), cex.lab=1.5)
mtext("f)", 3, at=0, line=1, cex=1.2)

boxplot(BRF_trait$NN.SES~BRF_trait$deer2, cex=1.5, cex.axis=1.5, xlab=expression(italic('SES D'['nn'])), cex.lab=1.5)
mtext("g)", 3, at=0, line=1, cex=1.2)

boxplot(BRF_trait$NN.RAW~BRF_trait$deer2, cex=1.5, cex.axis=1.5, xlab=expression(italic('*D'['nn'])), cex.lab=1.5)
mtext("h)", 3, at=0, line=1, cex=1.2)

########## Appendix S3 Plots #####################
library(plyr)
library(dplyr)
library(tidyr)
library(ggplot2)

##simes species data
simes<-read.csv("Simes_data_for_analysis.csv", header=T)
str(simes)
dim(simes)
colnames(simes)
dat = unite(simes, bpt, block, plot, treatment, remove = T ) %>% 
  tbl_df() %>% 
  select(-myrsmi, -lepcan, -moose)#remove two species without trait data

dat = gather(dat, key = "sp", value = "freq", -year, -bpt)

dat = group_by(dat, year, bpt) %>% 
  summarise(total_freq = sum(freq, na.rm = T)) %>% 
  left_join(dat) %>% 
  mutate(rel_abund = freq/total_freq) %>% 
  select(-total_freq)

##trait data
all_traits<-read.csv("FuncTraitv2.csv", header=T, na.strings=c("","NA"), stringsAsFactors = F)
str(all_traits)
summary(all_traits)


#get continuous traits

traits = select(all_traits, Species_Code, HL:RLL) %>% 
  filter(Species_Code %in% unique(dat$sp)) %>% na.omit() %>% 
  rename(sp = Species_Code)

dat2 = left_join(dat, traits, by = "sp") %>% 
  gather(key = "cont_traits", value = "value", HL:RLL)

#calculate community weighted mean
dat3 = group_by(dat2, year, bpt, cont_traits) %>% 
  mutate(rel_v = rel_abund * value) %>% 
  summarise(cwm = sum(rel_v, na.rm = T)) %>% 
  arrange(bpt, cont_traits)

#subtract the relative community weighted mean from 2003 from each subsequent year
dat3_diff = group_by(dat3, bpt, cont_traits) %>% 
  mutate(cwm_diff = (cwm/max(cwm))-(cwm[year == 2003]/max(cwm)))

## get traits that are factors 
traits2 = select(all_traits, -HL, -REL, -RLL, -Subfamily, -Genus, -Species) %>% 
  filter(Species_Code %in% unique(dat$sp)) %>%
  rename(sp = Species_Code)

dat2_disc = left_join(dat, traits2, by = "sp") %>% 
  gather(key = "disc_traits", value = "value", Colony_size:Biogeographic_Affinity) %>% 
  mutate(value = stringr::str_trim(value))

#calculate community weighted mean
dat3_disc = group_by(dat2_disc, year, bpt, disc_traits, value) %>% 
  summarise(cwm = sum(rel_abund)) %>% 
  arrange(bpt, disc_traits)

#subtract the community weighted mean from 2003 from each subsequent year
dat3_diff_disc = group_by(dat3_disc, bpt, disc_traits, value) %>% 
  mutate(cwm_diff = cwm - cwm[year == 2003])


#Combine continuous and factor cwm diff
all_diff = bind_rows(rename(dat3_diff, traits = cont_traits) %>% 
                       mutate(status = "continous"), 
                     rename(dat3_diff_disc, traits = disc_traits) %>% 
                       mutate(status = "discrete")) %>% 
  rename(level = value)

# look at the average departure in functional similarity for continuous traits each plot from the 2003 pre-treatment baseline 
time_mean=filter(all_diff, status == "continous", year != 2003) %>% 
  group_by(bpt, traits) %>% 
  summarise(ave_diff = mean(cwm_diff, na.rm = T)) %>% 
  arrange(traits,bpt) %>% 
  separate(bpt, into = c("block","plot","treat"),sep = "_")

p<-ggplot(time_mean, aes(x=treat, y=ave_diff, colour=block))+
  geom_point()+
  xlab("Treatment")+
  ylab("Abundance weighted mean trait difference")+
  ggtitle("Continuous Traits", subtitle = NULL)+
  theme(plot.title = element_text(hjust = 0.5))+
  facet_grid(.~traits)
p


# # look at the average departure in functional similarity for factor traits for each plot from the 2003 pre-treatment baseline  
time_mean.disc=filter(all_diff, status == "discrete", year != 2003) %>% 
  group_by(bpt, traits,level) %>% 
  summarise(ave_diff = mean(cwm_diff, na.rm = T)) %>% 
  arrange(level,bpt) %>% 
  separate(bpt, into = c("block","plot","treat"),sep = "_")


# to plot each trait use the trait name from the FuncTraitv2.csv file in the "traits == "trait name"
cat<-filter(time_mean.disc, traits == "Behavior",level!="NA")
p<-ggplot(cat, aes(x=treat, y=ave_diff, colour=block))+
  geom_point()+
  xlab("Treatment")+
  ylab("Abundance weighted mean trait difference")+
  ggtitle("Behavior", subtitle = NULL)+
  theme(plot.title = element_text(hjust = 0.5))+
  facet_grid(.~level)
p