########################################################################
#        Principal response curves and species scores for thesis chapter 3 

#                              22-Mar-2019
########################################################################
###included is PRC plots for heath and melur as well as the species scores divided by 
##species groups 

rm(list=ls())
setwd("~/sfuvault/R/PRC")
library("readxl")
#add columns year to each data set, keep the species and treatment, and subset for habitat

pc6 <- read_excel("Chapter3_completedata copy.xlsx", sheet= "PC2016") 
pc6.heath<-subset(pc6, pc6$habitat=="heath")
pc6.melur<-subset(pc6, pc6$habitat=="melur")
pc6.heath$year<-paste("2016")
pc6.melur$year<-paste("2016")
pc6$year<-paste("2016")
pc6.heath$year<-paste("2016")
pc6.heath<-pc6.heath[c(1:30),c(1,5,67,28:66)]
pc6.melur<-pc6.melur[c(1:30),c(1,5,67,28:66)]
PC6<-pc6[c(1:60),c(1,5,67,28:66)]

pc7 <- read_excel("chapter3_completedata copy.xlsx", sheet= "PC2017") 
pc7.heath<-subset(pc7, pc7$habitat=="heath")
pc7.melur<-subset(pc7, pc7$habitat=="melur")
pc7$year<-paste("2017")
pc7.heath$year<-paste("2017")
pc7.melur$year<-paste("2017")
pc7.heath$year<-paste("2017")
PC7<-pc7[,c(1,4,56,19:55)]
pc7.heath<-pc7.heath[,c(1,4,56,19:55)]
pc7.melur<-pc7.melur[,c(1,4,56,19:55)]

pc8 <- read_excel("chapter3_completedata copy.xlsx", sheet= "PC2018")
pc8.heath<-subset(pc8, pc8$habitat=="heath")
pc8.melur<-subset(pc8, pc8$habitat=="melur")
pc8$year<-paste("2018")
pc8.heath$year<-paste("2018")
pc8.melur$year<-paste("2018")
PC8<-pc8[,c(1,4,71,29:70)]
pc8.heath<-pc8.heath[,c(1,4,71,29:70)]
pc8.melur<-pc8.melur[,c(1,4,71,29:70)]

pc9 <- read_excel("chapter3_completedata copy.xlsx", sheet= "PC2019") 
pc9.heath<-subset(pc9, pc9$habitat=="heath")
pc9.melur<-subset(pc9, pc9$habitat=="melur")
pc9$year<-paste("2019")
pc9.heath$year<-paste("2019")
pc9.melur$year<-paste("2019")
PC9<-pc9[,c(1,4,73,30:72)]
pc9.heath<-pc9.heath[,c(1,4,73,30:72)]
pc9.melur<-pc9.melur[,c(1,4,73,30:72)]

##run for melur
PC.all.melur<-merge(pc6.melur, pc7.melur, all=TRUE)
PC.all.melur2<-merge(PC.all.melur, pc8.melur, all=TRUE)
PC.all.melur3<-merge(PC.all.melur2, pc9.melur, all=TRUE)


PC.all.melur3[is.na(PC.all.melur3)] <- 0

row.names(PC.all.melur3)<-paste(PC.all.melur3$plot, PC.all.melur3$year, sep=".")

data.spp<-PC.all.melur3[,c(4:52)]

##code for species scores grouped by species groups (i.e. poa glauce and alpina grouped into poa species)
data.spp.melur<-read.csv("data_spp_melur.csv")  
data.spp<-data.spp.melur[,c(2:32)]

data.env<-PC.all.melur3[,c(1:3)]
data.env$fTTM <- relevel(factor(data.env$treatment, labels=c("C","CF","K","N","NK", "NP","NPK","NPKF","P","PK")), ref="C")
data.env$fYEAR <- factor(data.env$year) 

library(vegan)
data.env$fYEAR= as.factor(data.env$year)
log.data.spp <- log(data.spp+1) #need to log transform prior to analyses
dev.off()

prc.all.melur <- prc(log.data.spp, treatment=data.env$fTTM, time=data.env$fYEAR)

plot(prc.all.melur, type="b", species=T, pch=19, cex=2, col=c("black","deepskyblue2","darksalmon","Yellow","darkorange3","chartreuse1","chartreuse2","red","darkorchid3"), 
     legpos=NA, ylim=c(-1,1), las=1, xlab="year")
legend(x="topleft", y=1, pch=c(7,10,10,NA), col=c("black","deepskyblue2","darksalmon","Yellow","darkorange3","chartreuse1","chartreuse2","red","darkorchid3"), lty=c(1,3,2,1),
       legend=c("CF","K","N","NK", "NP","NPK","NPKF","P","PK"))


##summary of species coef
sum_prcmelur<-summary(prc.all.melur, scaling=1)
sum.species<-sum_prcmelur$sp
scores<-scores(sum.species)
range(scores)

##filter coef values that are less than 0 
summary<-Filter(function(x) x > 0.1, sum.species)

is.numeric(sum.species)
#plot scores for abundant species only (>1000 hits)
logabu <- colSums(log.data.spp)

plot(prc.all.melur, type="b",logabu>20, species=F, col=c("black","deepskyblue2","darksalmon","Yellow","darkorange3","chartreuse1","chartreuse2","red","darkorchid3"), 
     pch=c(12,3,3,3,3,3,12,3,3), scaling=3, lty=c(1,5,5,5,5,5,1,5,5), ylim=c(-.1,.7), cex=1.4, lwd=2, legpos=NA, las=1, xlab="year", main="gravelly desert >20% cover")

legend(x="left", cex=1, y=-1, pch=c(12,3,3,3,3,3,3,3,12),lty=c(1,5,5,5,5,5,5,5,1) , col=c("black","darksalmon","red","deepskyblue2","Yellow","darkorchid3","darkorange3","chartreuse2", "chartreuse1"),
       legend=c("CF","N","P","K","NK","PK","NP","NPK","NPKF"))

##linestack filtered 
linestack(summary,  cex = 1, side = "right", hoff = 2, air = 1.7, at = 0, axis = TRUE, ylim=c(-7,7))

#full linestack
linestack(sum.species,  cex = 1, side = "right", hoff = 2, air = 1.7, at = 0, axis = TRUE, ylim=c(-7,7))
scores(summary)

###run for heath
PC.all.heath<-merge(pc6.heath, pc7.heath, all=TRUE)
PC.all.heath2<-merge(PC.all.heath, pc8.heath, all=TRUE)
PC.all.heath3<-merge(PC.all.heath2, pc9.heath, all=TRUE)


PC.all.heath3[is.na(PC.all.heath3)] <- 0

row.names(PC.all.heath3)<-paste(PC.all.heath3$plot, PC.all.heath3$year, sep=".")


data.spp<-PC.all.heath3[,c(4:52)]

data.spp.heath<-read.csv("data_spp_heath.csv")  
data.spp<-data.spp.heath[,c(2:35)]
data.env<-PC.all.heath3[,c(1:3)]
data.env$fTTM <- relevel(factor(data.env$treatment, labels=c("C","CF","K","N","NK", "NP","NPK","NPKF","P","PK")), ref="C")
data.env$fYEAR <- factor(data.env$year) 

library(vegan)
data.env$fYEAR= as.factor(data.env$year)
log.data.spp <- log(data.spp+1) #need to log transform prior to analyses
dev.off()

prc.all.heath <- prc(log.data.spp, treatment=data.env$fTTM, time=data.env$fYEAR)

plot(prc.all.heath, type="b", species=T, pch=19, cex=1.4, col=c("orange","Brown","GreenYellow","DarkGreen","Black","Pink","Yellow","Blue","HotPink"), 
     legpos=NA, ylim=c(-.2,.7), las=1, xlab="year")
legend(x="bottomleft", cex=1, y=-1, pch=c(10,10,10,NA), col=c("orange","GreenYellow","Blue","Brown","DarkGreen","HotPink","Black","Pink","Yellow"), lty=c(1,3,2,1),
       legend=c("CF","N","P","K","NK","PK","NP","NPK","NPKF"))

##summary of species coefficients
sum_prcheath<-summary(prc.all.heath, scaling=1)
sum_prcheath$sp

#plot scores for abundant species only (>50 hits)

logabu <- colSums(data.spp)
plot(prc.all.heath, type="b",logabu>20, species=F, col=c("black","deepskyblue2","darksalmon","darkgoldenrod1","darkorange3","chartreuse1","chartreuse2","deeppink2","darkorchid3"), 
     pch=c(12,3,3,3,3,3,12,3,3), lwd=2, scaling=3, lty=c(1,5,5,5,5,5,1,5,5), ylim=c(-.1,.7), cex=1.4, legpos=NA, las=1, xlab="year", main="dwarf-shrub heath")

legend(x="topleft", cex=1, y=-1, pch=c(12,3,3,3,3,3,3,3,12),lty=c(1,5,5,5,5,5,5,5,1) , col=c("black","darksalmon","deeppink2","deepskyblue2","Yellow","darkorange3","darkorchid3","chartreuse2", "chartreuse1"),
       legend=c("CF","N","P","K","NK","PK","NP","NPK","NPKF"))

linestack(sum_prcheath$sp,  cex = 1, side = "right", hoff = 2, air = 1.7, at = 0, axis = TRUE, ylim=c(-20,20))

sum_prcheath<-summary(prc.all.heath, scaling=1)
sum.species<-sum_prcheath$sp
summary<-Filter(function(x)  x < -.1 | x > .1, sum.species)
scores<-scores(sum.species)
range(scores)
linestack(summary,  cex = 1, side = "right", hoff = 2, air = 1.7, at = 0, axis = TRUE, ylim=c(-7,7))
scores(summary)