######################################################################## # Analysis for effect of grazing #chapter 2 thesis # 22-Mar-2018 ######################################################################## rm(list=ls()) #first thing is to specify your working directory setwd("~/Desktop/R/chapter 1") #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("chapter2_completedata.xls") #sites in rows, species in columns, no blank cells #the file has 84 observations of 76 variables plus 1 for bare #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:86)] #this dataset should have 48 observations of 9 variables plus one for bare #Run to remove the accidentally Fertilised melur #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 Fertiliser and grazing (4 treatments # per block: C control unfenced, CF control fenced, # NPK Fertilised unfenced, NPKF Fertilised fenced) #"plot": specific ID code for each plot #"NPK": Fertilisation treatment, either Fertilised (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 grazing$pair<-paste(grazing$block, grazing$NPK, sep="") fenced<-subset(grazing, grazing$fence=="1") unfenced<-subset(grazing, grazing$fence=="0") effG.unstd <- data.frame(habitat=unfenced$habitat, site=unfenced$site, NPK=unfenced$NPK, pair=unfenced$pair, biomass=unfenced$biomass-(fenced$biomass), bare=unfenced$bare-(fenced$bare), graminoid=unfenced$graminoid-(fenced$graminoid), grass=unfenced$grass-(fenced$grass), sedge.rush=unfenced$sedges+unfenced$rush-(fenced$sedges+fenced$rush), decid.shrub=unfenced$`decid.shrub`-(fenced$`decid.shrub`), shrub=unfenced$shrub-(fenced$shrub), herb=unfenced$herb-(fenced$herb)) ##standardized effect size values effG.unordered <- data.frame(habitat=effG.unstd$habitat, site=effG.unstd$site, NPK=effG.unstd$NPK, pair=effG.unstd$pair, biomass=(effG.unstd$biomass-mean(effG.unstd$biomass))/sd(effG.unstd$biomass), bare=(effG.unstd$bare-mean(effG.unstd$bare))/sd(effG.unstd$bare), graminoid=(effG.unstd$graminoid-mean(effG.unstd$graminoid))/sd(effG.unstd$graminoid), herb=(effG.unstd$herb-mean(effG.unstd$herb))/sd(effG.unstd$herb), shrub=(effG.unstd$shrub-mean(effG.unstd$shrub))/sd(effG.unstd$shrub), grass=(effG.unstd$grass-mean(effG.unstd$grass))/sd(effG.unstd$grass), sedge.rush=(effG.unstd$sedge.rush-mean(effG.unstd$sedge.rush))/sd(effG.unstd$sedge.rush), decid.shrub=(effG.unstd$decid.shrub-mean(effG.unstd$decid.shrub))/sd(effG.unstd$decid.shrub) ) ### making effG.unordered into alphabetisized order by pair effG.ordered<-effG.unordered[order(effG.unordered$pair),] ###adding lichen and moss dataset and reordering the pair so we can merge the datasets setwd("~/Desktop/ICB format") library("readxl") lichen.moss <- read_excel("Nutnet2017_icb_copy.xls", sheet="moss.lichen") names(lichen.moss) graz.lichen.moss <- subset(lichen.moss, lichen.moss$chapter=="grazing")[,c(1:9)] names(graz.lichen.moss) ##sum of cover together by each plot t<-aggregate(cover~plot+site+taxa+habitat+block+treatment+fence, graz.lichen.moss, sum) graz.ag<-as.data.frame(t) lichen.moss <- graz.ag lichen.moss$pair<-paste(lichen.moss$block, lichen.moss$treatment, sep="") fenced.l.m<-subset(lichen.moss, lichen.moss$fence=="1") unfenced.l.m<-subset(lichen.moss, lichen.moss$fence=="0") #split into mosses and lichen graz.moss<-subset(lichen.moss,lichen.moss$taxa=="moss") fenced.moss<-subset(graz.moss, graz.moss$fence=="1") unfenced.moss<-subset(graz.moss, graz.moss$fence=="0") graz.lichen<-subset(lichen.moss, lichen.moss$taxa=="lichen") fenced.lichen<-subset(graz.lichen, graz.lichen$fence=="1") unfenced.lichen<-subset(graz.lichen, graz.lichen$fence=="0") effG.cover.unstd <- data.frame(habitat=unfenced.moss$habitat, site=unfenced.moss$site, NPK=unfenced.moss$treatment, pair=unfenced.moss$pair, lichen=unfenced.lichen$cover-(fenced.lichen$cover), moss=unfenced.moss$cover-(fenced.moss$cover), total.moss.lichen=(unfenced.moss$cover+unfenced.lichen$cover)-(fenced.moss$cover+fenced.lichen$cover)) effG.cover <- data.frame(habitat=effG.cover.unstd$habitat, site=effG.cover.unstd$site, NPK=effG.cover.unstd$NPK, pair=effG.cover.unstd$pair, lichen=(effG.cover.unstd$lichen-mean(effG.cover.unstd$lichen))/sd(effG.cover.unstd$lichen), moss=(effG.cover.unstd$moss-mean(effG.cover.unstd$moss))/sd(effG.cover.unstd$moss), total.moss.lichen=(effG.cover.unstd$total.moss.lichen-mean(effG.cover.unstd$total.moss.lichen))/sd(effG.cover.unstd$total.moss.lichen) ) View(effG.cover) effG.ordered.cover<-effG.cover[order(effG.cover$pair),] ##add lichen and moss to dataset effG.ordered$moss<-paste(effG.ordered.cover$moss) effG.ordered$lichen<-paste(effG.ordered.cover$lichen) effG.ordered$total.moss.lichen<-paste(effG.ordered.cover$total.moss.lichen) ##final dataset name names(effG.ordered) ##this dataset has the effects of grazing on all variables effG<-effG.ordered ### move to end of code eventually... ### measuring the overall effect of grazing on biomass, bare, graminoid, shrub, herb ### to do this we can compare the fenced and unfenced plots for all variables without a significant effect ## first take the fenced and unfenced from dataset : "grazing" unfenced<-subset(grazing,grazing$fence=="0") fenced<-subset(grazing, grazing$fence=="1") t.test(unfenced$biomass, fenced$biomass, paired=TRUE) t.test(unfenced$graminoid, fenced$graminoid, paired=TRUE) t.test(unfenced$bare, fenced$bare, paired=TRUE) t.test(unfenced$shrub, fenced$shrub, paired=TRUE) t.test(unfenced$herb, fenced$herb, paired=TRUE) t.test(unfenced.l.m$cover, fenced.l.m$cover, paired=TRUE) ##combining this effect with the effect of Fertiliser for biomass and graminoids that had a significant interaction fert.unfenced<-subset(unfenced, unfenced$NPK=="NPK") fert.fenced<-subset(fenced, fenced$NPK=="NPK") unfert.unfenced<-subset(unfenced, unfenced$NPK=="C") unfert.fenced<-subset(fenced, fenced$NPK=="C") t.test(fert.fenced$biomass, unfert.fenced$biomass, paired=TRUE) t.test(unfert.unfenced$biomass, unfert.fenced$biomass, paired=TRUE) t.test(unfert.unfenced$bare, unfert.fenced$bare, paired=TRUE) t.test(fert.fenced$graminoid, fert.fenced$graminoid, paired=TRUE) t.test(unfert.unfenced$graminoid, unfert.fenced$graminoid, paired=TRUE) ##effect of Fertiliser t.test(fert.fenced$biomass, unfert.fenced$biomass) t.test(fert.fenced$graminoid, unfert.fenced$graminoid) t.test(fert.fenced$shrub, unfert.fenced$shrub) t.test(fert.fenced$herb, unfert.fenced$herb) t.test(fert.fenced$total.moss.lichen, unfert.fenced$total.moss.lichen) ##subset unfenced and untreated plots to give site description without grazing unfenced.control<-subset(unfenced, unfenced$NPK=="C") ##test bg in heath and melur heath.control<-subset(unfenced.control, unfenced.control$habitat=="heath") melur.control<-subset(unfenced.control, unfenced.control$habitat=="melur") t.test(heath.control$bare, melur.control$bare, paired=TRUE) ##test bg in aud and theist aud.control<-subset(unfenced.control, unfenced.control$site=="Audkuluheidi") theist.control<-subset(unfenced.control, unfenced.control$site=="Theistareykir") t.test(aud.control$bare, theist.control$bare, paired=TRUE) ##biomass in heath and melur t.test(heath.control$biomass, melur.control$biomass, paired=TRUE) ##biomass in aud and theist t.test(aud.control$biomass, theist.control$biomass, paired=TRUE) aud.heath<-subset(aud.control, aud.control$habitat=="heath") aud.melur<-subset(aud.control, aud.control$habitat=="melur") theist.heath<-subset(theist.control, theist.control$habitat=="heath") theist.melur<-subset(theist.control, theist.control$habitat=="melur") mean(aud.heath$biomass) mean(aud.melur$biomass) mean(theist.heath$biomass) mean(theist.melur$biomass) sum(melur.control$graminoid)/sum(melur.control$biomass) sum(heath.control$shrub)/sum(heath.control$biomass) ##subset ## testing the effect of grazing on biomass model.bio<-lm(biomass~NPK*habitat*site, data=effG) summary(model.bio) model.bio.1a<-lm(biomass~NPK*habitat*site, data=effG) model.bio.1b<-lm(biomass~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.bio.1a, model.bio.1b) #NPK*habitat*site (p=0.1918) model.bio.2a<-lm(biomass~NPK*habitat+NPK*site+habitat*site, data=effG) model.bio.2b<-lm(biomass~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.bio.2a, model.bio.2b) #NPK*habitat (p=0.9817) <- drop model.bio.2c<-lm(biomass~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.bio.2a, model.bio.2c) #NPK*site (p=0.3103) model.bio.2d<-lm(biomass~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.bio.2a, model.bio.2d) #habitat*site (p=0.1723) model.bio.3a<-lm(biomass~NPK*site+habitat*site, data=effG) model.bio.3b<-lm(biomass~NPK+site+habitat*site, data=effG) anova(model.bio.3a, model.bio.3b) #NPK*site (p=0.2961) <- drop model.bio.3c<-lm(biomass~NPK*site+habitat+site, data=effG) anova(model.bio.3a, model.bio.3c) #habitat*site (p=0.1599) model.bio.4a<-lm(biomass~NPK+site+habitat*site, data=effG) model.bio.4b<-lm(biomass~NPK+site+habitat+site, data=effG) anova(model.bio.4a, model.bio.4b) #habitat*site (p=0.1606) model.bio.5a<-lm(biomass~NPK+site+habitat, data=effG) summary.bio<-summary(model.bio.5a) model.bio.5b<-lm(biomass~site+habitat, data=effG) anova(model.bio.5a,model.bio.5b) model.bio.5c<-lm(biomass~NPK+habitat, data=effG) anova(model.bio.5a,model.bio.5c) model.bio.5d<-lm(biomass~NPK+site, data=effG) anova(model.bio.5a,model.bio.5d) ##plot for total biomass library(ggplot2) summary.lm.fit <- summary.bio coef.matrix <- summary.lm.fit$coefficients ## coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 ##multiply std.error by 1.96 lm.df <- data.frame(Effect.total.biomass= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.total.biomass)) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) plot.bio<-g #testing the effect of grazing on bare-ground model.bare<-lm(bare~NPK*habitat*site, data=effG) summary(model.bare) model.bare.1a<-lm(bare~NPK*habitat*site, data=effG) model.bare.1b<-lm(bare~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.bare.1a, model.bare.1b) #NPK*habitat*site (p=0.7196) model.bare.2a<-lm(bare~NPK*habitat+NPK*site+habitat*site, data=effG) model.bare.2b<-lm(bare~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.bare.2a, model.bare.2b) #NPK*habitat (p=0.4924) model.bare.2c<-lm(bare~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.bare.2a, model.bare.2c) #NPK*site (p=0.3925) model.bare.2d<-lm(bare~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.bare.2a, model.bare.2d) #habitat*site (p=0.8722) <- drop model.bare.3a<-lm(bare~NPK*habitat+NPK*site+habitat+site, data=effG) model.bare.3b<-lm(bare~NPK*habitat+NPK+site+habitat+site, data=effG) anova(model.bare.3a, model.bare.3b) #NPK*site (p=0.3789) model.bare.3c<-lm(bare~NPK+habitat+NPK*site+habitat+site, data=effG) anova(model.bare.3a, model.bare.3c) #NPK*habitat (p=0.4799) <- drop model.bare.4a<-lm(bare~NPK+habitat+NPK*site+habitat+site, data=effG) model.bare.4b<-lm(bare~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.bare.4a, model.bare.4b) #habitat*site (p=0.3723) model.bare.5a<-lm(bare~NPK+habitat+site, data=effG) model.bare.5b<-lm(bare~habitat+site, data=effG) anova(model.bare.5a, model.bare.5b) model.bare.5c<-lm(bare~NPK+site, data=effG) anova(model.bare.5a, model.bare.5c) model.bare.5d<-lm(bare~NPK+habitat, data=effG) anova(model.bare.5a, model.bare.5d) summary.bare<-summary(model.bare.5a) summary.lm.fit <- summary.bare coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.bare)) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) plot.bare<-g #testing the effect of grazing on graminoids model.gram<-lm(graminoid~NPK*habitat*site, data=effG) summary(model.gram) model.gram.1a<-lm(graminoid~NPK*habitat*site, data=effG) model.gram.1b<-lm(graminoid~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.gram.1a, model.gram.1b) #NPK*habitat*site (p=0.3597) model.gram.2a<-lm(graminoid~NPK*habitat+NPK*site+habitat*site, data=effG) model.gram.2b<-lm(graminoid~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.gram.2a, model.gram.2b) #NPK*habitat (p=0.9918) <- drop model.gram.2c<-lm(graminoid~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.gram.2a, model.gram.2c) #NPK*site (p=0.3919) model.gram.2d<-lm(graminoid~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.gram.2a, model.gram.2d) #habitat*site (p=0.6549) model.gram.3a<-lm(graminoid~NPK+habitat+NPK*site+habitat*site, data=effG) model.gram.3b<-lm(graminoid~NPK+habitat+NPK+site+habitat*site, data=effG) anova(model.gram.3a, model.gram.3b) #NPK*site (p=0.378) model.gram.3c<-lm(graminoid~NPK+habitat+NPK*site+habitat+site, data=effG) anova(model.gram.3a, model.gram.3c) #habitat*site (p=0.6453) <- drop model.gram.4a<-lm(graminoid~NPK+habitat+NPK*site+habitat+site, data=effG) model.gram.4b<-lm(graminoid~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.gram.4a, model.gram.4b) #NPK*site (p=0.3675) model.gram.5a<-lm(graminoid~NPK+habitat+site, data=effG) model.gram.5b<-lm(graminoid~habitat+site, data=effG) anova(model.gram.5a,model.gram.5b) model.gram.5a<-lm(graminoid~NPK+habitat+site, data=effG) model.gram.5c<-lm(graminoid~NPK+site, data=effG) anova(model.gram.5a,model.gram.5c) model.gram.5a<-lm(graminoid~NPK+habitat+site, data=effG) model.gram.5d<-lm(graminoid~NPK+habitat, data=effG) anova(model.gram.5a,model.gram.5d) summary.gram<-summary(model.gram.5a) ## NPK 0.0417 t- value = -2.177 summary.lm.fit <- summary.gram coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.gram= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.gram)) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) plot.gram<-g ##by grasses #testing the effect of grazing on grasses that are more palatable than sedges and rushes model.grass<-lm(grass~NPK*habitat*site, data=effG) summary(model.grass) model.grass.1a<-lm(grass~NPK*habitat*site, data=effG) model.grass.1b<-lm(grass~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.grass.1a, model.grass.1b) #NPK*habitat*site (p=0.335) model.grass.2a<-lm(grass~NPK*habitat+NPK*site+habitat*site, data=effG) model.grass.2b<-lm(grass~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.grass.2a, model.grass.2b) #NPK*habitat (p=0.89) <- drop model.grass.2c<-lm(grass~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.grass.2a, model.grass.2c) #NPK*site (p=0.42) model.grass.2d<-lm(grass~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.grass.2a, model.grass.2d) #habitat*site (p=0.81) model.grass.3a<-lm(grass~NPK+habitat+NPK*site+habitat*site, data=effG) model.grass.3b<-lm(grass~NPK+habitat+NPK+site+habitat*site, data=effG) anova(model.grass.3a, model.grass.3b) #NPK*site (p=0.42) model.grass.3c<-lm(grass~NPK+habitat+NPK*site+habitat+site, data=effG) anova(model.grass.3a, model.grass.3c) #habitat*site (p=0.81) <- drop model.grass.4a<-lm(grass~NPK+habitat+NPK*site+habitat+site, data=effG) model.grass.4b<-lm(grass~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.grass.4a, model.grass.4b) #NPK*site (p=0.39) model.grass.5a<-lm(grass~NPK+habitat+NPK+site, data=effG) summary.grass<-summary(model.grass.5a) ##significant p = 0.0347 NPK ##sedged and rushes model.sedge.rush<-lm(sedge.rush~NPK*habitat*site, data=effG) summary(model.sedge.rush) model.sedge.rush.1a<-lm(sedge.rush~NPK*habitat*site, data=effG) model.sedge.rush.1b<-lm(sedge.rush~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.sedge.rush.1a, model.sedge.rush.1b) #NPK*habitat*site (p=0.335) model.sedge.rush.2a<-lm(sedge.rush~NPK*habitat+NPK*site+habitat*site, data=effG) model.sedge.rush.2b<-lm(sedge.rush~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.sedge.rush.2a, model.sedge.rush.2b) #NPK*habitat (p=0.5) model.sedge.rush.2c<-lm(sedge.rush~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.sedge.rush.2a, model.sedge.rush.2c) #NPK*site (p=0.8)<- drop model.sedge.rush.2d<-lm(sedge.rush~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.sedge.rush.2a, model.sedge.rush.2d) #habitat*site (p=0.1) model.sedge.rush.3a<-lm(sedge.rush~NPK*habitat+NPK+site+habitat*site, data=effG) model.sedge.rush.3b<-lm(sedge.rush~NPK*habitat+NPK+site+habitat+site, data=effG) anova(model.sedge.rush.3a, model.sedge.rush.3b) #habitat*site (p=0.09) model.sedge.rush.3c<-lm(sedge.rush~NPK+habitat+NPK+site+habitat*site, data=effG) anova(model.sedge.rush.3a, model.sedge.rush.3c) #habitat*NPK (p=0.55) <- drop model.sedge.rush.4a<-lm(sedge.rush~NPK+habitat+NPK+site+habitat*site, data=effG) model.sedge.rush.4b<-lm(sedge.rush~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.sedge.rush.4a, model.sedge.rush.4b) #hab*site (p=0.09) model.sedge.rush.5a<-lm(sedge.rush~NPK+habitat+NPK+site, data=effG) summary.sedge.rush<-summary(model.sedge.rush.5a) ##sedges and rushes had no significant effect of grazing on NPK p=0.886 #testing the effect of grazing on herbs (note we changed the name to forbs in the manuscript) model.herb<-lm(herb~NPK*habitat*site, data=effG) summary(model.herb) model.herb.1a<-lm(herb~NPK*habitat*site, data=effG) model.herb.1b<-lm(herb~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.herb.1a, model.herb.1b) #NPK*habitat*site (p=0.1367) model.herb.2a<-lm(herb~NPK*habitat+NPK*site+habitat*site, data=effG) model.herb.2b<-lm(herb~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.herb.2a, model.herb.2b) #NPK*habitat (p=0.151) model.herb.2c<-lm(herb~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.herb.2a, model.herb.2c) #NPK*site (p=0.3653) <- drop model.herb.2d<-lm(herb~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.herb.2a, model.herb.2d) #habitat*site (p=0.2375) model.herb.3a<-lm(herb~NPK*habitat+NPK+site+habitat*site, data=effG) model.herb.3b<-lm(herb~NPK+habitat+NPK+site+habitat*site, data=effG) anova(model.herb.3a, model.herb.3b) #NPK*habitat (p=0.1486) model.herb.3c<-lm(herb~NPK*habitat+NPK+site+habitat+site, data=effG) anova(model.herb.3a, model.herb.3c) #habitat*site (p=0.2349) <- drop model.herb.4a<-lm(herb~NPK*habitat+NPK+site+habitat+site, data=effG) model.herb.4b<-lm(herb~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.herb.4a, model.herb.4b) #NPK*habitat (p=0.1528) model.herb.5a<-lm(herb~NPK+habitat+site, data=effG) model.herb.5b<-lm(herb~habitat+site, data=effG) anova(model.herb.5a, model.herb.5b) model.herb.5c<-lm(herb~NPK+site, data=effG) anova(model.herb.5a, model.herb.5c) model.herb.5d<-lm(herb~habitat+NPK, data=effG) anova(model.herb.5a, model.herb.5d) summary.herb<-summary(model.herb.5a) #########################PLOT FOR HERBS######################### summary.lm.fit <- summary.herb coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) plot.herb<- g ## nothing significant for herbs #testing the effect of grazing on shrubs model.shrub<-lm(shrub~NPK*habitat*site, data=effG) summary(model.shrub) model.shrub.1a<-lm(shrub~NPK*habitat*site, data=effG) model.shrub.1b<-lm(shrub~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.shrub.1a, model.shrub.1b) #NPK*habitat*site (p=0.4397) model.shrub.2a<-lm(shrub~NPK*habitat+NPK*site+habitat*site, data=effG) model.shrub.2b<-lm(shrub~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.shrub.2a, model.shrub.2b) #NPK*habitat (p=0.7487) model.shrub.2c<-lm(shrub~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.shrub.2a, model.shrub.2c) #NPK*site (p=0.8646) <- drop model.shrub.2d<-lm(shrub~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.shrub.2a, model.shrub.2d) #habitat*site (p=0.2025) model.shrub.3a<-lm(shrub~NPK*habitat+NPK+site+habitat*site, data=effG) model.shrub.3b<-lm(shrub~NPK+habitat+NPK+site+habitat*site, data=effG) anova(model.shrub.3a, model.shrub.3b) #NPK*habitat (p=0.7417) <- drop model.shrub.3c<-lm(shrub~NPK*habitat+NPK+site+habitat+site, data=effG) anova(model.shrub.3a, model.shrub.3c) #habitat*site (p=0.1898) model.shrub.4a<-lm(shrub~NPK+habitat+NPK+site+habitat*site, data=effG) model.shrub.4b<-lm(shrub~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.shrub.4a, model.shrub.4b) #habitat*site (p=0.1789) model.shrub.5a<-lm(shrub~NPK+habitat+site, data=effG) model.shrub.5b<-lm(shrub~habitat+site, data=effG) anova(model.shrub.5a,model.shrub.5b) model.shrub.5c<-lm(shrub~habitat+NPK, data=effG) anova(model.shrub.5a,model.shrub.5c) model.shrub.5c<-lm(shrub~NPK+site, data=effG) anova(model.shrub.5a,model.shrub.5c) summary.shrub<-summary(model.shrub.5a) #########################PLOT FOR SHRUBS######################### summary.lm.fit <- summary.shrub coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.shrub= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.shrub)) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + geom_hline(yintercept= 0)+ theme(legend.position="none",axis.title.x=element_blank()) g plot.shrub<-g ## nothing significant for shrubs ### decidious shrubs model.decid.shrub<-lm(decid.shrub~NPK*habitat*site, data=effG) summary(model.decid.shrub) model.decid.shrub.1a<-lm(decid.shrub~NPK*habitat*site, data=effG) model.decid.shrub.1b<-lm(decid.shrub~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.decid.shrub.1a, model.decid.shrub.1b) #NPK*habitat*site (p=0.4397) model.decid.shrub.2a<-lm(decid.shrub~NPK*habitat+NPK*site+habitat*site, data=effG) model.decid.shrub.2b<-lm(decid.shrub~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.decid.shrub.2a, model.decid.shrub.2b) #NPK*habitat (p=0.3602) model.decid.shrub.2c<-lm(decid.shrub~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.decid.shrub.2a, model.decid.shrub.2c) #NPK*site (p=0.7743) <- drop model.decid.shrub.2d<-lm(decid.shrub~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.decid.shrub.2a, model.decid.shrub.2d) #habitat*site (p=0.5636) model.decid.shrub.3a<-lm(decid.shrub~NPK*habitat+NPK+site+habitat*site, data=effG) model.decid.shrub.3b<-lm(decid.shrub~NPK+habitat+NPK+site+habitat*site, data=effG) anova(model.decid.shrub.3a, model.decid.shrub.3b) #NPK*habitat (p=0.3472) model.decid.shrub.3c<-lm(decid.shrub~NPK*habitat+NPK+site+habitat+site, data=effG) anova(model.decid.shrub.3a, model.decid.shrub.3c) #habitat*site (p=0.553) <- drop model.decid.shrub.4a<-lm(decid.shrub~NPK+habitat+NPK+site+habitat*site, data=effG) model.decid.shrub.4b<-lm(decid.shrub~NPK+habitat+NPK+site+habitat+site, data=effG) anova(model.decid.shrub.4a, model.decid.shrub.4b) #habitat*site (p=0.5519) model.decid.shrub.5a<-lm(decid.shrub~NPK+habitat+NPK+site, data=effG) summary.decid.shrub<-summary(model.decid.shrub.5a) ##deciduous shrubs also have no effect of grazing #testing the effect of grazing on lichen, moss, and lichen and moss together #start by looking at lichens model.bio<-lm(lichen~NPK*habitat*site, data=effG) summary(model.bio) model.bio.1a<-lm(lichen~NPK*habitat*site, data=effG) model.bio.1b<-lm(lichen~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.bio.1a, model.bio.1b) #NPK*habitat*site (p=0.05884) close but not significant model.bio.2a<-lm(lichen~NPK*habitat+NPK*site+habitat*site, data=effG) model.bio.2b<-lm(lichen~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.bio.2a, model.bio.2b) #NPK*habitat (p=0.5252) <- drop model.bio.2c<-lm(lichen~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.bio.2a, model.bio.2c) #NPK*site (p=0.07893) model.bio.2d<-lm(lichen~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.bio.2a, model.bio.2d) #habitat*site (p=0.07761) model.bio.3a<-lm(lichen~NPK*site+habitat*site, data=effG) model.bio.3b<-lm(lichen~NPK+site+habitat*site, data=effG) anova(model.bio.3a, model.bio.3b) #NPK*site (p=0.0735) <- drop model.bio.3c<-lm(lichen~NPK*site+habitat+site, data=effG) anova(model.bio.3a, model.bio.3c) #habitat*site (p=0.07231) model.bio.4a<-lm(lichen~NPK+site+habitat*site, data=effG) model.bio.4b<-lm(lichen~NPK+site+habitat+site, data=effG) anova(model.bio.4a, model.bio.4b) #habitat*site (p=0.08937) model.bio.5a<-lm(lichen~NPK+site+habitat, data=effG) summary.bio<-summary(model.bio.5a) ##nothing significant for lichen ## try and explore the nearly signficant 3way interaction with lichen p = 0.058 heath.lichen<-subset(effG, effG$habitat=="heath") melur.lichen<-subset(effG,effG$habitat=="melur") model.lichen.2a<-lm(lichen~NPK*site, data=heath.lichen) model.lichen.2b<-lm(lichen~NPK+site, data=heath.lichen) anova(model.lichen.2a, model.lichen.2b) model.lichen.2a<-lm(lichen~NPK*site, data=melur.lichen) model.lichen.2b<-lm(lichen~NPK+site, data=melur.lichen) anova(model.lichen.2a, model.lichen.2b) aud.lichen<-subset(effG, effG$site=="Audkuluheidi") theist.lichen<-subset(effG,effG$site=="Theistareykir") model.lichen.2a<-lm(lichen~NPK*habitat, data=aud.lichen) model.lichen.2b<-lm(lichen~NPK+habitat, data=aud.lichen) anova(model.lichen.2a, model.lichen.2b) model.lichen.2a<-lm(lichen~NPK*habitat+NPK+habitat, data=theist.lichen) model.lichen.2b<-lm(lichen~NPK+habitat, data=theist.lichen) anova(model.lichen.2a, model.lichen.2b) #(NPK*habitat p=0.01404) summary.lichen<-summary(model.lichen.2a) ##now look at mosses model.bio<-lm(moss~NPK*habitat*site, data=effG) summary(model.bio) model.bio.1a<-lm(moss~NPK*habitat*site, data=effG) model.bio.1b<-lm(moss~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.bio.1a, model.bio.1b) #NPK*habitat*site (p=0.2003) model.bio.2a<-lm(moss~NPK*habitat+NPK*site+habitat*site, data=effG) model.bio.2b<-lm(moss~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.bio.2a, model.bio.2b) #NPK*habitat (p=0.9276) model.bio.2c<-lm(moss~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.bio.2a, model.bio.2c) #NPK*site (p=0.6225) model.bio.2d<-lm(moss~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.bio.2a, model.bio.2d) #habitat*site (p=0.977)<- drop model.bio.3a<-lm(moss~NPK*habitat+NPK*site, data=effG) model.bio.3b<-lm(moss~NPK*habitat+NPK+site, data=effG) anova(model.bio.3a, model.bio.3b) #NPK*site (p=0.6122) model.bio.3c<-lm(moss~NPK+habitat+NPK*site, data=effG) anova(model.bio.3a, model.bio.3c) #habitat*NPK (p=0.9254) <- drop model.bio.4a<-lm(moss~NPK+habitat+NPK*site, data=effG) model.bio.4b<-lm(moss~NPK+habitat+NPK+site, data=effG) anova(model.bio.4a, model.bio.4b) #NPK*site (p=0.6023) model.bio.5a<-lm(moss~NPK+site+habitat, data=effG) summary.bio<-summary(model.bio.5a) ###nothing significant ## testing the effect of grazing on moss and lichen biomass combined model.total.m.l<-lm(total.moss.lichen~NPK*habitat*site, data=effG) summary(model.total.m.l) model.total.m.l.1a<-lm(total.moss.lichen~NPK*habitat*site, data=effG) model.total.m.l.1b<-lm(total.moss.lichen~NPK*habitat+NPK*site+habitat*site, data=effG) anova(model.total.m.l.1a, model.total.m.l.1b) #NPK*habitat*site (p=0.09954) model.total.m.l.2a<-lm(total.moss.lichen~NPK*habitat+NPK*site+habitat*site, data=effG) model.total.m.l.2b<-lm(total.moss.lichen~NPK+habitat+NPK*site+habitat*site, data=effG) anova(model.total.m.l.2a, model.total.m.l.2b) #NPK*habitat (p=0.9339) model.total.m.l.2c<-lm(total.moss.lichen~NPK*habitat+NPK+site+habitat*site, data=effG) anova(model.total.m.l.2a, model.total.m.l.2c) #NPK*site (p=0.8786)<- drop model.total.m.l.2d<-lm(total.moss.lichen~NPK*habitat+NPK*site+habitat+site, data=effG) anova(model.total.m.l.2a, model.total.m.l.2d) #habitat*site (p=0.723) model.total.m.l.3a<-lm(total.moss.lichen~NPK*habitat+habitat*site, data=effG) model.total.m.l.3b<-lm(total.moss.lichen~NPK+habitat+habitat*site, data=effG) anova(model.total.m.l.3a, model.total.m.l.3b) #NPK+habitat (p=0.932) <- drop model.total.m.l.3c<-lm(total.moss.lichen~NPK*habitat+habitat+site, data=effG) anova(model.total.m.l.3a, model.total.m.l.3c) #habitat*site (p=0.9226) model.total.m.l.4a<-lm(total.moss.lichen~NPK+habitat+habitat*site, data=effG) model.total.m.l.4b<-lm(total.moss.lichen~NPK+habitat+habitat+site, data=effG) anova(model.total.m.l.4a, model.total.m.l.4b) #habitat*site (p=0.9204) model.total.m.l.5a<-lm(total.moss.lichen~NPK+site+habitat, data=effG) model.total.m.l.5b<-lm(total.moss.lichen~site+habitat, data=effG) model.total.m.l.5c<-lm(total.moss.lichen~NPK+habitat, data=effG) model.total.m.l.5d<-lm(total.moss.lichen~NPK+site, data=effG) anova(model.total.m.l.5a,model.total.m.l.5b) #NPK anova(model.total.m.l.5a,model.total.m.l.5c) #site anova(model.total.m.l.5a,model.total.m.l.5d) #habitat summary.m.l<-summary(model.total.m.l.5a) ##nothing significant for lichen and moss combined ########################PLOT MOSS LICHEN############################ summary.lm.fit <- summary.m.l coef.matrix <- summary.m.l$coefficients coef.matrix[2:4,"Std. Error"]<-(coef.matrix[2:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.moss.lichen)) g <- g + geom_hline(yintercept= 0)+ theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) plot.m.l<-g ###############################Figure 2 ################################ #run this code for figure 2 rm(list=ls()) #first thing is to specify your working directory setwd("~/Desktop/R/chapter 1") #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("chapter2_completedata.xls") #sites in rows, species in columns, no blank cells #the file has 84 observations of 76 variables plus 1 for bare #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:86)] #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 Fertiliser and grazing (4 treatments # per block: C control unfenced, CF control fenced, # NPK Fertilised unfenced, NPKF Fertilised fenced) #"plot": specific ID code for each plot #"NPK": Fertilisation treatment, either Fertilised (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 grazing$pair<-paste(grazing$block, grazing$NPK, sep="") fenced<-subset(grazing, grazing$fence=="1") unfenced<-subset(grazing, grazing$fence=="0") effG.unstd <- data.frame(habitat=unfenced$habitat, site=unfenced$site, NPK=unfenced$NPK, pair=unfenced$pair, biomass=unfenced$biomass-(fenced$biomass), bare=unfenced$bare-(fenced$bare), graminoid=unfenced$graminoid-(fenced$graminoid), grass=unfenced$grass-(fenced$grass), sedge.rush=unfenced$sedges+unfenced$rush-(fenced$sedges+fenced$rush), decid.shrub=unfenced$`decid.shrub`-(fenced$`decid.shrub`), shrub=unfenced$shrub-(fenced$shrub), herb=unfenced$herb-(fenced$herb)) ##standardized effect size values effG.unordered <- data.frame(habitat=effG.unstd$habitat, site=effG.unstd$site, NPK=effG.unstd$NPK, pair=effG.unstd$pair, biomass=(effG.unstd$biomass-mean(effG.unstd$biomass))/sd(effG.unstd$biomass), bare=(effG.unstd$bare-mean(effG.unstd$bare))/sd(effG.unstd$bare), graminoid=(effG.unstd$graminoid-mean(effG.unstd$graminoid))/sd(effG.unstd$graminoid), herb=(effG.unstd$herb-mean(effG.unstd$herb))/sd(effG.unstd$herb), shrub=(effG.unstd$shrub-mean(effG.unstd$shrub))/sd(effG.unstd$shrub), grass=(effG.unstd$grass-mean(effG.unstd$grass))/sd(effG.unstd$grass), sedge.rush=(effG.unstd$sedge.rush-mean(effG.unstd$sedge.rush))/sd(effG.unstd$sedge.rush), decid.shrub=(effG.unstd$decid.shrub-mean(effG.unstd$decid.shrub))/sd(effG.unstd$decid.shrub) ) ### making effG.unordered into alphabetisized order by pair effG.ordered<-effG.unordered[order(effG.unordered$pair),] ###adding lichen and moss dataset and reordering the pair so we can merge the datasets setwd("~/Desktop/ICB format") library("readxl") lichen.moss <- read_excel("Nutnet2017_icb_copy.xls", sheet="moss.lichen") names(lichen.moss) graz.lichen.moss <- subset(lichen.moss, lichen.moss$chapter=="grazing")[,c(1:9)] names(graz.lichen.moss) ##sum of cover together by each plot t<-aggregate(cover~plot+site+taxa+habitat+block+treatment+fence, graz.lichen.moss, sum) graz.ag<-as.data.frame(t) lichen.moss <- graz.ag lichen.moss$pair<-paste(lichen.moss$block, lichen.moss$treatment, sep="") fenced.l.m<-subset(lichen.moss, lichen.moss$fence=="1") unfenced.l.m<-subset(lichen.moss, lichen.moss$fence=="0") #split into mosses and lichen graz.moss<-subset(lichen.moss,lichen.moss$taxa=="moss") fenced.moss<-subset(graz.moss, graz.moss$fence=="1") unfenced.moss<-subset(graz.moss, graz.moss$fence=="0") graz.lichen<-subset(lichen.moss, lichen.moss$taxa=="lichen") fenced.lichen<-subset(graz.lichen, graz.lichen$fence=="1") unfenced.lichen<-subset(graz.lichen, graz.lichen$fence=="0") effG.cover.unstd <- data.frame(habitat=unfenced.moss$habitat, site=unfenced.moss$site, NPK=unfenced.moss$treatment, pair=unfenced.moss$pair, lichen=unfenced.lichen$cover-(fenced.lichen$cover), moss=unfenced.moss$cover-(fenced.moss$cover), total.moss.lichen=(unfenced.moss$cover+unfenced.lichen$cover)-(fenced.moss$cover+fenced.lichen$cover)) effG.cover <- data.frame(habitat=effG.cover.unstd$habitat, site=effG.cover.unstd$site, NPK=effG.cover.unstd$NPK, pair=effG.cover.unstd$pair, lichen=(effG.cover.unstd$lichen-mean(effG.cover.unstd$lichen))/sd(effG.cover.unstd$lichen), moss=(effG.cover.unstd$moss-mean(effG.cover.unstd$moss))/sd(effG.cover.unstd$moss), total.moss.lichen=(effG.cover.unstd$total.moss.lichen-mean(effG.cover.unstd$total.moss.lichen))/sd(effG.cover.unstd$total.moss.lichen) ) effG.ordered.cover<-effG.cover[order(effG.cover$pair),] ##add lichen and moss to dataset effG.ordered$moss<-paste(effG.ordered.cover$moss) effG.ordered$lichen<-paste(effG.ordered.cover$lichen) effG.ordered$total.moss.lichen<-paste(effG.ordered.cover$total.moss.lichen) ##final dataset name names(effG.ordered) ##this dataset has the effects of grazing on all variables effG<-effG.ordered ##plot for biomass model.bio.5a<-lm(biomass~NPK+site+habitat, data=effG) summary.bio<-summary(model.bio.5a) ##plot for total biomass library(ggplot2) summary.lm.fit <- summary.bio coef.matrix <- summary.lm.fit$coefficients ## coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 ##multiply std.error by 1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + coord_cartesian( ylim=(c(-2,2))) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + ylab("Effect of grazing") g <- g + ggtitle("total biomass") g <- g + theme(axis.title.y = element_text(size = rel(0.7), angle = 90)) plot.bio<-g #plot for bare model.bare.5a<-lm(bare~NPK+habitat+NPK+site, data=effG) summary.bare<-summary(model.bare.5a) summary.lm.fit <- summary.bare coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + coord_cartesian( ylim=(c(-2,2))) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + ylab("Effect of grazing") g <- g + ggtitle("bare ground") g <- g + theme(axis.title.y = element_text(size = rel(0.7), angle = 90)) plot.bare<-g ##plot for graminoids model.gram.5a<-lm(graminoid~NPK+habitat+NPK+site, data=effG) summary.gram<-summary(model.gram.5a) ## NPK 0.0417 t- value = -2.177 summary.lm.fit <- summary.gram coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + coord_cartesian( ylim=(c(-2,2))) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + ylab("Effect of grazing") g <- g + ggtitle("graminoids") g <- g + theme(axis.title.y = element_text(size = rel(0.7), angle = 90)) plot.gram<-g #plot herb (changed to forb) model.herb.5a<-lm(herb~NPK+habitat+NPK+site, data=effG) summary.herb<-summary(model.herb.5a) summary.lm.fit <- summary.herb coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + coord_cartesian( ylim=(c(-2,2))) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + ylab("Effect of grazing") g <- g + ggtitle("forbs") g <- g + theme(axis.title.y = element_text(size = rel(0.7), angle = 90)) plot.herb<- g #plot shrub model.shrub.5a<-lm(shrub~NPK+habitat+NPK+site, data=effG) summary.shrub<-summary(model.shrub.5a) summary.lm.fit <- summary.shrub coef.matrix <- summary.lm.fit$coefficients coef.matrix[1:4,"Std. Error"]<-(coef.matrix[1:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + coord_cartesian( ylim=(c(-2,2))) g <- g + geom_hline(yintercept= 0)+theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + ylab("Effect of grazing") g <- g + ggtitle("shrubs") g <- g + theme(axis.title.y = element_text(size = rel(0.7), angle = 90)) g plot.shrub<-g ##plot lichen moss model.total.m.l.5a<-lm(total.moss.lichen~NPK+site+habitat, data=effG) summary.m.l<-summary(model.total.m.l.5a) ##nothing significant for lichen and moss combined ########################PLOT MOSS LICHEN############################ summary.lm.fit <- summary.m.l coef.matrix <- summary.m.l$coefficients coef.matrix[2:4,"Std. Error"]<-(coef.matrix[2:4,"Std. Error"])*1.96 lm.df <- data.frame(Effect.of.grazing= c(coef.matrix[2:4, "Estimate"], sd.lower = coef.matrix[2:4, "Estimate"] - coef.matrix[2:4, "Std. Error"], sd.upper = coef.matrix[2:4, "Estimate"] + coef.matrix[2:4, "Std. Error"]), Variable = rep(c("Fertiliser", "Habitat", "Site"), 3), type = c(rep("estimate", 3), rep("sd", 3), rep("sd", 3))) g <- ggplot(lm.df, aes(x = Variable, y = Effect.of.grazing)) g <- g + coord_cartesian( ylim=(c(-2,2))) g <- g + geom_hline(yintercept= 0)+ theme(legend.position="none",axis.title.x=element_blank()) g <- g + geom_point(aes(shape = type, color = Variable), size = 3) + scale_shape_manual(values = c(16, 95)) #+ scale_shape_manual g <- g + geom_path(aes(color = Variable)) g <- g + ylab("Effect of grazing") g <- g + theme(axis.title.y = element_text(size = rel(0.7), angle = 90)) g <- g + ggtitle("moss and lichen") plot.m.l<-g library(gridExtra) library(cowplot) mygrid<-plot_grid(plot.bare, plot.bio, plot.gram, plot.shrub, plot.herb, plot.m.l, labels="auto", ncol=2) grid.arrange(mygrid,bottom= "Effect of grazing on a) aboveground biomass, b) bareground, c) graminoids, d) shrubs, e) forbs, f) moss and lichen in plots treated with Fertiliser, between habitat (heath and melur), and across sites (in/outside the volcanic active zone)")