########################################################################
#                        BIOMASS AND BARE 2017
#              data analysis for chapter 2 of thesis
#                         run for figure below 
#                           22-Mar-2018
########################################################################
#note: code before we measured the effect of grazing. This has general habitat discriptions and figure 3 at the end of the code

rm(list=ls())
#first thing is to specify your working directory
setwd("~/sfuvault")
#this should specify the path to the folder where your data are stored
download.packages(rlang)
#your data should be saved as tab delimited text file
#to load your data, you have to create an 'object' that later you will be able to call
#here we're loading your data to an object named 'nutnet'
 library("readxl")
nutnet <- read_excel("chapter1_completedata.xls")  
#sites in rows, species in columns, no blank cells
#the file has 84 observations of 76 variables plus 1 for bare
View(nutnet)

#this is the whole dataset for 2017, but for now we want to focus on the
#experiment crossing fertilisers and grazing, and we are only interested in the first
#7 columns (indicating the plot, treatment and so on) and the last two (which are actually
#the variables we are interested in: total biomass and percent cover of bare ground)
grazing <- subset(nutnet, nutnet$chapter=="grazing")[,c(1:8,76,77)] 
#this dataset should have 48 observations of 9 variables plus one for bare


#a bit of organizing (just making sure the different levels make sense)
names(grazing) #just checking here what things are called :)
#and what those columns mean:
#(writing down what the variables are here is good practice, 
#so that you can recall exactly what your variables are,
#the possible values they can take, etc  when you look at your file again)
#"site": the two study locations, Audkuluheidi and Theistareykir
#"habitat": melur or heath
#"block": id of NutNET block (3 blocks per habitat: AH1,AH2,AH3,AM1,AM2...)
#"treatment": describes the crossed effect of fertilizer and grazing (4 treatments 
#             per block: C control unfenced, CF control fenced, 
#             NPK fertilized unfenced, NPKF fertilized fenced)
#"plot": specific ID code for each plot
#"NPK": fertilization treatment, either fertilized (NPK) or not (C)
#"fence": if plots are fenced (1) or not (0)
#"biomass": biomass (g/m2) harvested in the plots -- converted to g/m2

#other ways of quickly ckecking your data
summary(grazing)
str(grazing)

grazing$fTTM=relevel(factor(grazing$treatment), ref="C") #use control as the reference level
grazing$fTTMcf=relevel(factor(grazing$treatment), ref="CF") #use control-fenced as the reference level
grazing$fTTMnpk=relevel(factor(grazing$treatment), ref="NPK") #use NPK as the reference level

grazing$fFENCE=relevel(factor(grazing$fence), ref="0") #use unfenced as the reference level
grazing$fNPK=relevel(factor(grazing$NPK), ref="C") #use control as the reference level

grazing$fSITE=factor(grazing$site) #make sure factors are factors
grazing$fBLOCK=factor(grazing$block) #make sure factors are factors
grazing$fHAB=factor(grazing$habitat) #make sure factors are factors

grazing$fTTM2=paste(grazing$fHAB,".",grazing$fTTM, sep="") 
#this is a combination of habitat and treatment
grazing$fSITE2=paste(grazing$fSITE,".",grazing$fHAB, sep="") 
#this is a combination of site and habitat



########################################################################
#                 UNIVARIATE ANALYSIS
########################################################################
#with these analyses we want to see the effect of fertilisation and grazing
#in the two sites and habitats, on:
#a) the amount of plant biomass 
#b) the amount of bare ground

#some basic data exploration and analyses to see what things look like :)



#################################################################
#PLANT BIOMASS
mean(grazing$biomass); sd(grazing$biomass)

#just to see what the variable looks like (this affects the choice of model we make)
par(mfrow=c(1,3)) #just telling that we want three charts arranged in a row
boxplot(grazing$biomass)
hist(grazing$biomass, breaks=25, main= "plant biomass g/m2)")
dotchart(grazing$biomass)

#how is biomass related to other variables?
par(mfrow=c(1,1)) #we just want one chart now
boxplot<-boxplot(grazing$biomass~grazing$fTTM)

#we have there the median values of biomass in each treatment 
#(but this is not separating the different habitats)
par(mfrow=c(1,1)) 
n<-boxplot(grazing$biomass~grazing$fHAB)
n$stats
#there are definitively some difference between heath and melur :)
par(mfrow=c(1,1)) 
b<-boxplot(grazing$biomass~grazing$fSITE)
b$stats
#biomass seems to be more variable in Theistareykir


#to see the treatments in the different site/habitats we can use
library(lattice)
bwplot(grazing$biomass~grazing$fTTM|grazing$fSITE2)


##let's build our model :)
#we want to see the effects of treatment and habitat on biomass, 
#and we need to tell the model that our plots are grouped in blocks
#we will do that with a linear mixed effect model (LMM)
library(nlme)
mod.graz=lme(biomass~fFENCE*fNPK*fHAB*fSITE, random=~1|fBLOCK, method="REML", data=grazing)
summary(mod.graz)

#but this summary is messy, as there are many levels in the combinations of the interaction
#it is easier if we get just one value of significance for the whole interaction term
#to do that, we can compare the models with and without the interaction term
#to do this comparison we have to use method="ML" in the formula (just believe me ;P)
#if you want to read more about models and model selection, there is a very good book by Alain Zuur
#Zuur et al (2009) 'Mixed effects models and extensions in ecology with R' Springer
mod.graz1a=lme(biomass~fFENCE*fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
mod.graz1b=lme(biomass~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz1a, mod.graz1b)  #FENCE*NPK*HAB*fSITE (p=0.1077) -- because it is non-significant we should drop it

mod.graz2a=lme(biomass~fFENCE*fNPK*fHAB + fFENCE*fNPK*fSITE + fFENCE*fHAB*fSITE + fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
mod.graz2b=lme(biomass~fFENCE*fNPK+fFENCE*fHAB+fNPK*fHAB+
                 fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2b)  #fFENCE*NPK*HAB (p=0.978) <- we drop this first (least significant)
mod.graz2c=lme(biomass~fFENCE*fNPK*fHAB+ fFENCE*fNPK+fFENCE*fSITE+fNPK*fSITE+
                 fFENCE*fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2c)  #fFENCE*fNPK*fSITE (p=0.2188) 
mod.graz2d=lme(biomass~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+
                 fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2d)  #fFENCE*fHAB*fSITE (p=0.0968) -- because it is non-significant we should drop it
mod.graz2e=lme(biomass~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+
                 fNPK*fHAB+fNPK*fSITE+fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2e)  #fNPK*fHAB*fSITE (p=0.6316) -- because it is non-significant we should drop it

#cannot drop yet any of the two way interactions because they are included in other 3-way interactions
mod.graz3a=lme(biomass~fFENCE*fNPK+fFENCE*fHAB+fNPK*fHAB+
                 fFENCE*fNPK*fSITE+ fFENCE*fHAB*fSITE+ fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz3b=lme(biomass~fFENCE*fNPK+fFENCE*fHAB+fNPK*fHAB+
                 fFENCE*fNPK+fFENCE*fSITE+fNPK*fSITE+ fFENCE*fHAB*fSITE+ fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz3a, mod.graz3b)  #fFENCE*fNPK*fSITE (p=0.2188)
mod.graz3c=lme(biomass~fFENCE*fNPK+fFENCE*fHAB+fNPK*fHAB+
                 fFENCE*fNPK*fSITE+ fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+ fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz3a, mod.graz3c)  #fFENCE*fHAB*fSITE (p=0.0968)
mod.graz3d=lme(biomass~fFENCE*fNPK+fFENCE*fHAB+fNPK*fHAB+
                 fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+ fNPK*fHAB+fNPK*fSITE+fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz3a, mod.graz3d)  #fNPK*fHAB*fSITE (p=0.6316) <- we drop this first (least significant)

#we can still only drop the three-way interactions
mod.graz4a=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fNPK*fSITE+ fFENCE*fHAB*fSITE+ fNPK*fSITE + fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz4b=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fNPK+fFENCE*fSITE+fNPK*fSITE+ fFENCE*fHAB*fSITE+ fNPK*fSITE + fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz4a, mod.graz4b)  #fFENCE*fNPK*fSITE (p=0.2203) <- we drop this first (least significant)
mod.graz4c=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fNPK*fSITE+ fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+ fNPK*fSITE + fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz4a, mod.graz4c)  #fFENCE*fHAB*fSITE (p=0.0978) 

#here we can drop the three way interaction and all the two way interactions that are not included in the three way one
mod.graz5a=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK*fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz5b=lme(biomass~fFENCE+fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK*fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5b)  #fFENCE*fNPK (p=0.0118) 
mod.graz5c=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK+fHAB+
                 fFENCE*fSITE+ fNPK*fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5c)  #fNPK*fHAB (p=0.0841) 
mod.graz5d=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5d)  #fNPK*fSITE (p=0.1437) <- we drop this first (least significant)
mod.graz5e=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK*fSITE+ fFENCE*fHAB+ fFENCE*fSITE+ fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5e)  #fFENCE*fHAB*fSITE (p=0.1047) 


mod.graz6a=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz6b=lme(biomass~fFENCE+fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz6a, mod.graz6b)  #fFENCE*fNPK (p=0.0142) 
mod.graz6c=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK+fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fFENCE*fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz6a, mod.graz6c)  #fNPK*fHAB (p=0.0932) 
mod.graz6d=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fFENCE*fHAB+ fFENCE*fSITE+ fHAB*fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz6a, mod.graz6d)  #fFENCE*fHAB*fSITE (p=0.1149) <- we drop this first (least significant)



mod.graz7a=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz7b=lme(biomass~fFENCE+fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz7a, mod.graz7b)  #fFENCE*fNPK (p=0.0176) 
mod.graz7c=lme(biomass~fFENCE*fNPK+ fFENCE+fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz7a, mod.graz7c)  #fFENCE*fHAB (p=0.8948) <- we drop this first (least significant)
mod.graz7d=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK+fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz7a, mod.graz7d)  #fNPK*fHAB (p=0.1044) 
mod.graz7e=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE+fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz7a, mod.graz7e)  #fFENCE*fSITE (p=0.8594)
mod.graz7f=lme(biomass~fFENCE*fNPK+ fFENCE*fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB+fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz7a, mod.graz7f)  #fHAB*fSITE (p=0.0056) 


mod.graz8a=lme(biomass~fFENCE*fNPK+ fFENCE +fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz8b=lme(biomass~fFENCE+fNPK+ fFENCE +fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz8a, mod.graz8b)  #fFENCE*fNPK (p=0.0176) 
mod.graz8c=lme(biomass~fFENCE*fNPK+ fFENCE +fHAB+ fNPK+fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz8a, mod.graz8c)  #fNPK*fHAB (p=0.1045) 
mod.graz8d=lme(biomass~fFENCE*fNPK+ fFENCE +fHAB+ fNPK*fHAB+
                 fFENCE+fSITE+ fNPK+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz8a, mod.graz8d)  #fFENCE*fSITE (p=0.8594) <- we drop this first (least significant)
mod.graz8e=lme(biomass~fFENCE*fNPK+ fFENCE +fHAB+ fNPK*fHAB+
                 fFENCE*fSITE+ fNPK+ fSITE+ fHAB+fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz8a, mod.graz8e)  #fHAB*fSITE (p=0.0056) 


mod.graz9a=lme(biomass~fFENCE*fNPK+ fHAB+ fNPK*fHAB+
                 fFENCE+ fSITE+ fNPK+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz9b=lme(biomass~fFENCE+fNPK+ fHAB+ fNPK*fHAB+
                 fFENCE+ fSITE+ fNPK+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz9a, mod.graz9b)  #fFENCE*fNPK (p=0.0176) 
mod.graz9c=lme(biomass~fFENCE*fNPK+ fHAB+ fNPK+fHAB+
                 fFENCE+ fSITE+ fNPK+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz9a, mod.graz9c)  #fNPK*fHAB (p=0.1046) <- we drop this first (least significant)
mod.graz9d=lme(biomass~fFENCE*fNPK+ fHAB+ fNPK*fHAB+
                 fFENCE+ fSITE+ fNPK+ fHAB+fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz9a, mod.graz9d)  #fHAB*fSITE (p=0.0056) 

mod.graz10a=lme(biomass~fFENCE*fNPK+ fHAB+ fNPK+
                  fFENCE+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz10b=lme(biomass~fFENCE+fNPK+ fHAB+ fNPK+
                  fFENCE+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
dd<-anova(mod.graz10a, mod.graz10b)  #fFENCE*fNPK (p=0.0218) 

mod.graz10c=lme(biomass~fFENCE*fNPK+ fHAB+ fNPK+
                  fFENCE+ fSITE+ fHAB+fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
dd<-anova(mod.graz10a, mod.graz10c)  #fHAB*fSITE (p=0.0056) 
  #fFENCE*fNPK (p=0.0218) 

#FINAL MODEL
final.mod.graz=lme(biomass~fFENCE*fNPK+ fHAB+ fNPK+
                     fFENCE+ fSITE+ fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
summary(final.mod.graz)
dd<-summary(final.mod.graz)$tTable


####FIGURE 1 -- plant biomass 2017
par(mfrow=c(1,1))
#barplot for biomass 2017
#we want to make a bar graph for Fence*NPK and another for Hab*Site
##### FENCE*NPK
non.fertilized <- subset(grazing, grazing$fNPK=="C")
means.non.fertilized <- tapply(non.fertilized$biomass,
                               non.fertilized$fFENCE,mean)
fertilized <- subset(grazing, grazing$fNPK=="NPK")
means.fertilized <- tapply(fertilized$biomass,
                           fertilized$fFENCE,mean)
heights <- matrix(c(means.non.fertilized,means.fertilized),2,2,byrow=TRUE)

bp <- barplot(heights, las=1, ylim=c(0,650), beside=TRUE, xlab='grazing exclusion',
              names.arg=c("non-fenced","fenced"),
              ylab="total above-ground biomass (g/m^2)",
              cex.names=1.2,
              cex.lab=1.2,
              col=c("Green2","DarkGreen"))

sds.non.fertilized <- tapply(non.fertilized$biomass,
                             non.fertilized$fFENCE,sd)
sds.fertilized <- tapply(fertilized$biomass,
                         fertilized$fFENCE,sd)
length(fertilized$biomass)
length(non.fertilized$biomass)
#24 is the number of observations we have in both fertilized and non-fertilized (length tells you the number of obsv)
se.mat <- matrix(c(sds.non.fertilized,sds.fertilized),2,2,byrow=TRUE)*1.96/sqrt(24)

error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
  if(length(x)!=length(y)|length(y)!=length(lower)|length(lower)!=length(upper))
    stop("vectors must be same length")
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}  

error.bar(bp,heights,se.mat) 

legend(1,600,c("non-fertilized","fertilized"), bty="n", cex=1.5,
       fill=c("Green2","DarkGreen"))

###assessing the differences between the melur and heath bare ground bars 
grazing$BIO2 <- paste(grazing$fFENCE,grazing$fNPK, sep="")
graz.bio = lme(biomass~BIO2, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio)
#releveling the variable to compare a new baseline
grazing$BIO20npk<-factor(grazing$BIO2, levels=c("0NPK", "0C", "1NPK", "1C"))
graz.bio.0npk = lme(biomass~BIO20npk, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio.0npk)

grazing$BIO21c<-factor(grazing$BIO2, levels=c("1C","0NPK","0C","1NPK"))
graz.bio.1c = lme(biomass~BIO21c, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio.1c)

#write the letters showing which bars are similar and different
text(1.5,300,"a",cex=1.5)
text(2.5,400,"a",cex=1.5)
text(4.5,400,"a",cex=1.5)
text(5.5,550,"b",cex=1.5)




##### FIGURE 2 HAB*SITE

melur <- subset(grazing, grazing$fHAB=="melur")
means.melur <- tapply(melur$biomass, melur$fSITE,mean)
# we are calculating the means for melur in both sites (Aud and Theista)
heath <- subset(grazing, grazing$fHAB=="heath")
means.heath <- tapply(heath$biomass,
                      heath$fSITE,mean)
heights <- matrix(c(means.melur,means.heath),2,2,byrow=TRUE)

bp <- barplot(heights, las=1, ylim=c(0,600), beside=TRUE, xlab="site",
              names.arg=c("Audkuluheidi","Theistareykir"),
              ylab="plant biomass (g/m^2)",
              col=c("Purple","Blue"))

sds.melur <- tapply(melur$biomass,
                    melur$fSITE,sd)
sds.heath <- tapply(heath$biomass,
                    heath$fSITE,sd)
length(melur$biomass)
length(heath$biomass)
#24 is the number of observations we have in both melur and heath (length tells you the number of obsv)
se.mat <- matrix(c(sds.melur,sds.heath),2,2,byrow=TRUE)*1.96/sqrt(24)

error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
  if(length(x)!=length(y)|length(y)!=length(lower)|length(lower)!=length(upper))
    stop("vectors must be same length")
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}  

error.bar(bp,heights,se.mat) 



legend(1,500,c("melur","heath"), bty="n",
       fill=c("Purple","Blue"))


###assessing the differences between Hab*Site
grazing$BIO3 <- paste(grazing$fHAB,grazing$fSITE, sep="")
graz.bio = lme(biomass~BIO3, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio)
#releveling the variable to compare a new baseline
grazing$BIO30npk<-factor(grazing$BIO3, levels=c("heathTheistareykir", "melurAudkuluheidi", "melurTheistareykir", "heathAudkuluheidi"))
graz.bio.heathTheista = lme(biomass~BIO30npk, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio.heathTheista)

grazing$BIO31c<-factor(grazing$BIO3, levels=c("melurAudkuluheidi","melurTheistareykir","heathAudkuluheidi","heathTheistareykir"))
graz.bio.melurAud = lme(biomass~BIO31c, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio.melurAud)

text(1.5,250,"a",cex=1.5)
text(2.5,300,"a",cex=1.5)
text(4.5,300,"a",cex=1.5)
text(5.5,600,"b",cex=1.5)


###leave this just in case
#AUDKULUHEIDI
audk <- subset(grazing, grazing$fSITE=="Audkuluheidi")  #this dataset should have 24 observations 
audk.graz=lme(biomass~fFENCE+fNPK+fHAB, random=~1|fBLOCK, method="REML", data=audk)
summary(audk.graz)
dd<-summary(audk.graz)$tTable

#in Audkuluheidi fertilisers and fences increased biomass and there was a difference between heath and melur
#THEISTAREYKIR
theis <- subset(grazing, grazing$fSITE=="Theistareykir")  #this dataset should have 24 observations 
theis.graz=lme(biomass~fFENCE+fNPK+fHAB, random=~1|fBLOCK, method="REML", data=theis)
summary(theis.graz)
dd<-summary(theis.graz)$tTable

#in Theistareykir fences did not have an effect, while fertilizers had an effect and the differences 
#between habitats were much stronger
#(same results hold if we remove the plots fertilised by farmers TM3 and TM2C/CF)






#################################################################
#BARE GROUND
grazing$Bare<-grazing$bare
mean(grazing$Bare); sd(grazing$Bare)

#just to see what the variable looks like (this affects the choice of model we make)
par(mfrow=c(1,3)) #just telling that we want three charts arranged in a row
boxplot(grazing$Bare)
hist(grazing$Bare, breaks=25, main= "Bare ground (%)")
dotchart(grazing$Bare)

#how is bare ground related to other variables? This doesn't consider location of treatment
par(mfrow=c(1,1)) #we just want one chart now
x<-boxplot(grazing$Bare~grazing$fTTM)
x$stats
#we have there the median values of biomass in each treatment 
#(but this is not separating the different habitats)
par(mfrow=c(1,1)) 
y<-boxplot(grazing$Bare~grazing$fHAB)
y$stats
#there are definitively some difference between heath and melur :) but not really b/w aud and theista
z<-boxplot(grazing$Bare~grazing$fSITE)
z$stats
#to see the treatments in the different site/habitats we can use
library(lattice)
bwplot(grazing$Bare~grazing$fTTM|grazing$fSITE2)
#broader ranges in Theista but similar

##let's build our model :)
#we want to see the effects of treatment and habitat on bare ground, 
#and we need to tell the model that our plots are grouped in blocks
#we will do that with a linear mixed effect model (LMM)
library(nlme)
mod.graz=lme(Bare~fFENCE*fNPK*fHAB*fSITE, random=~1|fBLOCK, method="REML", data=grazing)
summary(mod.graz)
#but this summary is messy, as there are many levels in the combinations of habitat*treatment 
#it is easier if we get just one value of significance for the whole interaction term
#to do that, we can compare the models with and without the interaction term
#to do this comparison we have to use method="ML" in the formula (just believe me ;P)
#if you want to read more about models and model selection, there is a very good book by Alain Zuur
#Zuur et al (2009) 'Mixed effects models and extensions in ecology with R' Springer
mod.graz1a=lme(Bare~fFENCE*fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
mod.graz1b=lme(Bare~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz1a, mod.graz1b)  #FENCE*NPK*HAB*fSITE (p=0.5442) -- because notsignificant we should drop it in the model

mod.graz2a=lme(Bare~fFENCE*fNPK*fHAB + fFENCE*fNPK*fSITE + fFENCE*fHAB*fSITE + fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
mod.graz2b=lme(Bare~fFENCE*fNPK+fFENCE*fHAB+fNPK*fHAB+
                 fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2b)  #FENCE*NPK*HAB (p=0.9929) <- we drop this first (least significant)
mod.graz2c=lme(Bare~fFENCE*fNPK*fHAB+ fFENCE*fNPK+fFENCE*fSITE+fNPK*fSITE+
                 fFENCE*fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2c)  #fFENCE*fNPK*fSITE (p=0.3798) 
mod.graz2d=lme(Bare~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+
                 fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2d)  #fFENCE*fHAB*fSITE (p=0.6263) -- because it not-significant we drop it
mod.graz2e=lme(Bare~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+fFENCE*fHAB*fSITE+
                 fNPK*fHAB+fNPK*fSITE+fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz2a, mod.graz2e)  #fNPK*fHAB*fSITE (p=0.0159) <-- significant means we cant drop

#cannot drop yet any of the two way interactions because they are included in other 3-way interactions
mod.graz3a=lme(Bare~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+
                 fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz3b=lme(Bare~fFENCE*fNPK+fHAB*fFENCE+fNPK*fHAB+fFENCE*fNPK*fSITE+
                 fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)
anova(mod.graz3a, mod.graz3b)  #fFENCE*fNPK*fHAB (p=0.9929) <-- drop this one
mod.graz3c=lme(Bare~fFENCE*fNPK*fHAB+fFENCE*fNPK+fSITE*fHAB+fFENCE*fSITE+
                 fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz3a, mod.graz3c)  #fFENCE*fNPK*fSITE (p=0.3813)
mod.graz3d=lme(Bare~fFENCE*fNPK*fHAB+fFENCE*fNPK*fSITE+
                 fFENCE*fHAB+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB+fSITE*fNPK+fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz3a, mod.graz3d)  #fNPK*fHAB*fSITE (p=0.0162) 

#we can still only drop the three-way interactions
mod.graz4a=lme(Bare~fFENCE*fNPK+fNPK*fHAB+fFENCE*fNPK*fSITE+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing) 
mod.graz4b=lme(Bare~fFENCE*fNPK+fNPK*fHAB+fFENCE*fNPK+fFENCE*fSITE+fNPK*fSITE+
                 fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz4a, mod.graz4b)  #fFENCE*fNPK*fSITE (p=0.3825) <- we drop this 
mod.graz4c=lme(Bare~fFENCE*fNPK+fNPK*fHAB+fFENCE*fNPK*fSITE+fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB+fSITE*fNPK+fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz4a, mod.graz4c)  #fNPK*fHAB*fSITE (p=0.0165) 

#here we can drop the three way interaction and all the two way interactions that are not included in the three way one
mod.graz5a=lme(Bare~fNPK*fHAB+fFENCE*fNPK+fNPK*fSITE+
                 fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)    
mod.graz5b=lme(Bare~fNPK*fHAB+fFENCE+fNPK+fNPK*fSITE+
                 fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5b)  #fFENCE*fNPK (p=0.7913) < -- drop the two way
mod.graz5c=lme(Bare~fNPK*fHAB+fFENCE*fNPK+fNPK*fSITE+
                 fFENCE+fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5c)  #fFENCE*fSITE (p=0.6599) 
mod.graz5d=lme(Bare~fNPK*fHAB+fFENCE*fNPK+fNPK*fSITE+
                 fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB+fSITE*fNPK+fSITE*fHAB, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz5a, mod.graz5d)  #fNPK*fHAB*fSITE (p=0.0175) <-- don't drop the 3 way 


mod.graz6a=lme(Bare~fNPK*fHAB+fFENCE+fNPK+fNPK*fSITE+
                 fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
mod.graz6b=lme(Bare~fNPK*fHAB+fFENCE+fNPK+fNPK*fSITE+
                 fFENCE+fSITE+fHAB*fSITE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz6a, mod.graz6b)  #fFENCE*fSITE (p=0.6602) <--- drop this 
mod.graz6c=lme(Bare~fNPK*fHAB+fFENCE+fNPK+fNPK*fSITE+
                 fFENCE*fSITE+fHAB*fSITE+fNPK*fHAB+fSITE*fNPK+fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  
anova(mod.graz6a, mod.graz6c)  #fNPK*fHAB*SITE (p=0.0176) 



#FINAL MODEL 
mod.graz7a=lme(Bare~fFENCE+fNPK*fHAB*fSITE, random=~1|fBLOCK, method="ML", data=grazing)  


## now we have to seperate by site and do a plot for NPK + HAB 

#AUDKULUHEIDI
audk <- subset(grazing, grazing$fSITE=="Audkuluheidi")  #this dataset should have 24 observations 
audk.graz=lme(Bare~fNPK*fHAB+fFENCE, random=~1|fBLOCK, method="REML", data=audk)
summary(audk.graz)
audk.graz1a=lme(Bare~fNPK*fHAB+fFENCE, random=~1|fBLOCK, method="ML", data=audk)
audk.graz1b=lme(Bare~fNPK+fHAB+fFENCE, random=~1|fBLOCK, method="ML", data=audk)
anova(audk.graz1a, audk.graz1b)         #fNPK*fHAB (p=<0.0001) -- because it is significant we should keep it in the model

#THEISTAREYKIR
theist <- subset(grazing, grazing$fSITE=="Theistareykir")  #this dataset should have 24 observations 
theist.graz=lme(Bare~fNPK*fHAB+fFENCE, random=~1|fBLOCK, method="REML", data=theist)
summary(theist.graz)
theist.graz1a=lme(Bare~fNPK*fHAB+fFENCE, random=~1|fBLOCK, method="ML", data=theist)
theist.graz1b=lme(Bare~fNPK+fHAB+fFENCE, random=~1|fBLOCK, method="ML", data=theist)
anova(theist.graz1a, theist.graz1b)         #fNPK*fHAB (p=<0.6625) -- because it is not significant 
summary(theist.graz1b)
#this means that we only have to plot the interaction for Audkuluheidi because this cannot be directly interpreted 



audk.m <- subset(audk, audk$fHAB=="melur") 
#there is an effect of fence and fert in the melur of audkuluheidi so we can't drop that
#we use the variable fTTM (and re-leveled variations) to see what levels are different
# making a baseline with C CF and NPK
audk.m.ttm = lme(Bare~fTTM, random=~1|fBLOCK, method="REML", data=audk.m)
summary(audk.m.ttm)$tTable
audk.m.ttmcf = lme(Bare~fTTMcf, random=~1|fBLOCK, method="REML", data=audk.m)
summary(audk.m.ttmcf)$tTable
audk.m.ttmnpk = lme(Bare~fTTMnpk, random=~1|fBLOCK, method="REML", data=audk.m)
summary(audk.m.ttmnpk)$tTable
#not a huge effect of fence but there was an effect of fertilised and non-fertilised plots in the melur in Audkuluheidi


####FIGURE 3 -- bare ground in audkuluheidi 2017
par(mar=c(8,4,2,3))
#barplot for bare ground 2017 only in aud
#we want to make a bar graph for NPK*HAB
##### NPK * HAB
audk <- subset(grazing, grazing$fSITE=="Audkuluheidi") 
non.fertilized.bg <- subset(audk, audk$fNPK=="C")
means.non.fertilized.bg <- tapply(non.fertilized.bg$Bare,
                                  non.fertilized.bg$fHAB,mean)
fertilized.bg <- subset(audk, audk$fNPK=="NPK")
means.fertilized.bg <- tapply(fertilized.bg$Bare,
                              fertilized.bg$fHAB,mean)
heights <- matrix(c(means.non.fertilized.bg,means.fertilized.bg),2,2,byrow=TRUE)

bp <- barplot(heights, las=1, ylim=c(0,90), beside=TRUE, xlab='habitat',
              names.arg=c("heath","melur"),
              ylab="percent cover of bare ground",
              cex.lab=1.5,
              cex.names=1.5,
              col=c("Goldenrod3","DarkOrange"))
title("Figure 2 Effect of fertilizer on percent cover of bare ground within the habitats 
of heath (90% vegetated) and melur (<1% vegetated). Data from Audkuluheidi site 
      (i.e. outside VAZ).",line= -25, cex.lab= 0.3)

sds.non.fertilized.bg <- tapply(non.fertilized.bg$Bare,
                                non.fertilized.bg$fHAB,sd)
sds.fertilized.bg <- tapply(fertilized.bg$Bare,
                            fertilized.bg$fHAB,sd)
length(fertilized.bg$Bare)
length(non.fertilized.bg$Bare)
#24 is the number of observations we have in both fertilized and non-fertilized (length tells you the number of obsv)
se.mat <- matrix(c(sds.non.fertilized.bg,sds.fertilized.bg),2,2,byrow=TRUE)*1.96/sqrt(24)

error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
  if(length(x)!=length(y)|length(y)!=length(lower)|length(lower)!=length(upper))
    stop("vectors must be same length")
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}  

error.bar(bp,heights,se.mat) 

legend(1.3,64,c("non-fertilized","fertilized"), bty="n", cex=1.5,
       fill=c("Goldenrod3","DarkOrange"))
###assessing the differences between the melur and heath bare ground bars 
audk$fHAB2 <- paste(audk$fHAB,audk$fNPK, sep="")
audk.hab = lme(Bare~fHAB2, random=~1|fBLOCK, method="REML", data=audk)
summary(audk.hab)$tTable

#releveling the variable to compare a new baseline
audk$fHAB2hnpk<-factor(audk$fHAB2, levels=c("heathNPK", "melurC", "melurNPK", "heathC"))
audk.hab.hnpk = lme(Bare~fHAB2hnpk, random=~1|fBLOCK, method="REML", data=audk)
summary(audk.hab.hnpk)$tTable

audk$fHAB2mc<-factor(audk$fHAB2, levels=c("melurC","heathNPK", "melurNPK", "heathC"))
audk.hab.mc = lme(Bare~fHAB2mc, random=~1|fBLOCK, method="REML", data=audk)
summary(audk.hab.mc)$tTable

#write the letters showing which bars are similar and different
text(1.5,16,"a",cex=1.5)
text(2.5,13,"a",cex=1.5)
text(4.5,80,"b",cex=1.5)
text(5.5,35,"c",cex=1.5)


#### Figure for Theista showing no differences 

theist <- subset(grazing, grazing$fSITE=="Theistareykir") 
non.fertilized.bg <- subset(theist, theist$fNPK=="C")
means.non.fertilized.bg <- tapply(non.fertilized.bg$Bare,
                                  non.fertilized.bg$fHAB,mean)
fertilized.bg <- subset(theist, theist$fNPK=="NPK")
means.fertilized.bg <- tapply(fertilized.bg$Bare,
                              fertilized.bg$fHAB,mean)
heights <- matrix(c(means.non.fertilized.bg,means.fertilized.bg),2,2,byrow=TRUE)

bp <- barplot(heights, las=1, ylim=c(0,90), beside=TRUE, xlab='habitat',
              names.arg=c("heath","melur"),
              ylab="percent cover of bare ground",
              cex.lab=1.5,
              cex.names=1.5,
              col=c("Goldenrod3","DarkOrange"))


sds.non.fertilized.bg <- tapply(non.fertilized.bg$Bare,
                                non.fertilized.bg$fHAB,sd)
sds.fertilized.bg <- tapply(fertilized.bg$Bare,
                            fertilized.bg$fHAB,sd)
length(fertilized.bg$Bare)
length(non.fertilized.bg$Bare)
#24 is the number of observations we have in both fertilized and non-fertilized (length tells you the number of obsv)
se.mat <- matrix(c(sds.non.fertilized.bg,sds.fertilized.bg),2,2,byrow=TRUE)*1.96/sqrt(24)

error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
  if(length(x)!=length(y)|length(y)!=length(lower)|length(lower)!=length(upper))
    stop("vectors must be same length")
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}  

error.bar(bp,heights,se.mat) 

legend(1.3,64,c("non-fertilized","fertilized"), bty="n", cex=1.5,
       fill=c("Goldenrod3","DarkOrange"))

###assessing the differences between the melur and heath bare ground bars 
theist$fHAB2 <- paste(theist$fHAB,theist$fNPK, sep="")
theist.hab = lme(Bare~fHAB2, random=~1|fBLOCK, method="REML", data=theist)
summary(theist.hab)$tTable

#releveling the variable to compare a new baseline
theist$fHAB2hnpk<-factor(theist$fHAB2, levels=c("heathNPK", "melurC", "melurNPK", "heathC"))
theist.hab.hnpk = lme(Bare~fHAB2hnpk, random=~1|fBLOCK, method="REML", data=theist)
summary(theist.hab.hnpk)$tTable

theist$fHAB2mc<-factor(theist$fHAB2, levels=c("melurC","heathNPK", "melurNPK", "heathC"))
theist.hab.mc = lme(Bare~fHAB2mc, random=~1|fBLOCK, method="REML", data=theist)
summary(theist.hab.mc)$tTable

#write the letters showing which bars are similar and different
text(1.5,9,"a",cex=1.5)
text(2.5,9,"a",cex=1.5)
text(4.5,83,"b",cex=1.5)
text(5.5,83,"b",cex=1.5)

#run biomass.bare.2017
unfenced <- subset(grazing, grazing$fFENCE=="0")
means.unfenced <- tapply(unfenced$biomass,
                         unfenced$fNPK,mean)
fenced <- subset(grazing, grazing$fFENCE=="1")
means.fenced <- tapply(fenced$biomass,
                       fenced$fNPK,mean)
heights <- matrix(c(means.unfenced,means.fenced),2,2,byrow=TRUE)

par(mar = c(4, 5, 3, 4), family="Times New Roman")
bp <- barplot(heights, las=1, ylim=c(0,650), beside=TRUE, xlab='treatment',
              names.arg=c("non-fertilised","fertilised"),
              ylab= expression("total above-ground biomass (g m"^paste("-2")*~')'),
              cex.names=1.2,
              cex.lab=1.2,
              col=c("Green2","DarkGreen"),
              family="Times New Roman")

sds.unfenced <- tapply(unfenced$biomass,
                       unfenced$fNPK,sd)
sds.fenced <- tapply(fenced$biomass,
                     fenced$fNPK,sd)
length(fenced$biomass)
length(unfenced$biomass)
#24 is the number of observations we have in both fenced and non-fenced (length tells you the number of obsv)
se.mat <- matrix(c(sds.unfenced,sds.fenced),2,2,byrow=TRUE)*1.96/sqrt(24)

error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
  if(length(x)!=length(y)|length(y)!=length(lower)|length(lower)!=length(upper))
    stop("vectors must be same length")
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}  

error.bar(bp,heights,se.mat) 

legend(1,600,c("unfenced","fenced"), bty="n", cex=1.3,
       fill=c("Green2","DarkGreen"))

###assessing the differences between the melur and heath bare ground bars 
grazing$BIO2 <- paste(grazing$fFENCE,grazing$fNPK, sep="")
graz.bio = lme(biomass~BIO2, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio)
#releveling the variable to compare a new baseline
grazing$BIO20npk<-factor(grazing$BIO2, levels=c("0NPK", "0C", "1NPK", "1C"))
graz.bio.0npk = lme(biomass~BIO20npk, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio.0npk)

grazing$BIO21c<-factor(grazing$BIO2, levels=c("1C","0NPK","0C","1NPK"))
graz.bio.1c = lme(biomass~BIO21c, random=~1|fBLOCK, method="REML", data=grazing)
summary(graz.bio.1c)

#write the letters showing which bars are similar and different
text(5.0,550,"*",cex=1.5)