####figures 

##input pre-code
rm(list=ls())
setwd("~/sfuvault/")
library("readxl")
pc9 <- read_excel("chapter3_completedata.xlsx", sheet= "PC2019") 
pc9<-as.data.frame(pc9)
pc9.species<-subset.data.frame(pc9[,30:72])
as.data.frame(pc9.env<-subset.data.frame(pc9[,1:12]))
pc9$treatment<-factor(pc9$treatment, levels=c("CF","C", "N", "P", "K","NK","PK", "NP","NPK","NPKF"))
library(BiodiversityR)
richness <- diversityresult(pc9.species, index='richness', method='each site', sortit=FALSE, digits=3)
#and now we can add these two columns to our dataset 'pc9[year]'
pc9$richness <- richness[,1]
pc9.heath<-subset(pc9,pc9$habitat=="heath")
pc9.melur<-subset(pc9,pc9$habitat=="melur")
#melur richness
number<-c(0,1,0,0,0,0,0,1,0,0)
mydata <- data.frame(pc9.melur$richness, pc9.melur$treatment,number)
names(mydata) <- c("value", "group","number")
# function for computing mean, DS, max and min values
min.mean.sd.max <- function(x) {
  r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

p1 <- ggplot(aes(y = value, x = factor(group), fill=factor(number)), data = mydata)
p1 <- p1 + stat_summary(fun.data = min.mean.sd.max, geom="boxplot") 
p1<- p1 + geom_jitter(position=position_jitter(width=0), size=3) + xlab("treatment") + ylab("species richness index")+theme_classic(base_size=15)
p1<- p1+ scale_fill_discrete(name="fence")
p1<-p1+coord_cartesian(ylim=c(5,23))
p1


##heath richness
number<-c(0,1,0,0,0,0,0,1,0,0)
mydata <- data.frame(pc9.heath$richness, pc9.heath$treatment,number)
names(mydata) <- c("value", "group","number")
# function for computing mean, DS, max and min values
min.mean.sd.max <- function(x) {
  r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

p1 <- ggplot(aes(y = value, x = factor(group), fill=factor(number)), data = mydata)
p1 <- p1 + stat_summary(fun.data = min.mean.sd.max, geom="boxplot") 
p1<- p1 + geom_jitter(position=position_jitter(width=0.3), size=3)+ xlab("treatment") + ylab("species richness index")+theme_classic(base_size=15)
p1<- p1+ scale_fill_discrete(name="fence")
p1<-p1+coord_cartesian(ylim=c(5,23))
p1








###code for biomass
setwd("~/sfuvault/")
library("readxl")
nutnet <- read_excel("chapter3_completedata.xlsx", sheet= "2019 biomass") 
nutnet<- as.data.frame(nutnet)

data.spp= nutnet[,c(9:43)]
data.spp<-as.data.frame(data.spp)
rownames(data.spp)=nutnet$plot #just to make sure that each row has the name of the corresponding plot
data.env= nutnet[,c(1:8,44)]
data.env=as.data.frame(data.env)
rownames(data.env)=nutnet$plot 
data.spp.UF <- subset(data.spp[2:35])  #this dataset should have 60 observations and 34var (removed moss)
data.UF <- subset(nutnet) 
data.UF$treatment=factor(data.UF$treatment, levels=c("CF","C", "N", "P", "K","NK","PK", "NP","NPK","NPKF")) #use control as the reference level
data.UF.heath<-subset(data.UF, data.UF$habitat=="heath")
data.UF.melur<-subset(data.UF, data.UF$habitat=="melur")

##biomass HEATH
number<-c(0,1,0,0,0,0,0,1,0,0)
mydata <- data.frame(data.UF.heath$biomass, data.UF.heath$treatment, number)
names(mydata) <- c("value", "group", "number")
# function for computing mean, DS, max and min values
names(mydata) <- c("value", "group","number")
# function for computing mean, DS, max and min values
min.mean.sd.max <- function(x) {
  r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

p1 <- ggplot(aes(y = value, x = factor(group), fill=factor(number)), data = mydata)
p1 <- p1 + stat_summary(fun.data = min.mean.sd.max, geom="boxplot") 
p1<- p1 + geom_jitter(position=position_jitter(width=0), size=3) + xlab("treatment") + ylab("biomass g/m^-2")+theme_classic(base_size=15)
p1<- p1+ scale_fill_discrete(name="fenced")
p1<-p1+coord_cartesian(ylim=c(0,1700))
p1

###ADD LETTERING... 

###MELUR
number<-c(0,1,0,0,0,0,0,1,0,0)
mydata <- data.frame(data.UF.melur$biomass, data.UF.melur$treatment, number)
names(mydata) <- c("value", "group", "number")
# function for computing mean, DS, max and min values
names(mydata) <- c("value", "group","number")
# function for computing mean, DS, max and min values
min.mean.sd.max <- function(x) {
  r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

p1 <- ggplot(aes(y = value, x = factor(group), fill=factor(number)), data = mydata)
p1 <- p1 + stat_summary(fun.data = min.mean.sd.max, geom="boxplot") 
p1<- p1 + geom_jitter(position=position_jitter(width=0), size=3) + xlab("treatment") + ylab("biomass g/m^-2")+theme_classic(base_size=15)
p1<- p1+ scale_fill_discrete(name="fenced")
p1<-p1+coord_cartesian(ylim=c(0,1700))
p1
##figure for biomass combined
number<-c(0,0,1,1,2,2,3,3,1,2)
mydata <- data.frame(data.UF$biomass, data.UF$treatment, number)
names(mydata) <- c("value", "group", "number")
# function for computing mean, DS, max and min values 
names(mydata) <- c("value", "group","number")
# function for computing mean, DS, max and min values
min.mean.sd.max <- function(x) {
  r <- c(min(x), mean(x) - sd(x), mean(x), mean(x) + sd(x), max(x))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

p1 <- ggplot(aes(y = value, x = factor(group), fill=factor(number)), data = mydata)
p1 <- p1 + stat_summary(fun.data = min.mean.sd.max, geom="boxplot") 
p1<- p1 + geom_jitter(position=position_jitter(width=0), size=3) + xlab("treatment") + ylab("biomass g/m^-2")+theme_classic()
p1<- p1+ scale_fill_discrete(name="number of nutrients")
p1