WLC Data analysis - Chemical RScript written by Erinne Stirling 2018 R version 3.5.1 (2018-07-02) -- "Feather Spray" Copyright (C) 2018 The R Foundation for Statistical Computing Platform: x86_64-apple-darwin15.6.0 (64-bit) #Create the environment library(readr) library(dplyr) library(ggplot2) library(agricolae) ##########----------########## C:N RATIO ##########----------########## #Call the data WLC<-read_csv("2018 WLC chemical data.csv")%>%#Read the WLC file transmute(TIME, REALTIME, CHAR, ID, CN=TC/TN, LITTER=ifelse(LITTER=="BLK", "Blank", ifelse(LITTER=="UL", "Unscorched", "Scorched"))) #Test for interactions anova(lm(CN~CHAR*LITTER, data=WLC)) #------- WLC0<-WLC%>% filter(TIME==0) anova(lm(CN~CHAR*LITTER, data=WLC0))#1 way: Char HSD.test(aov(CN~CHAR, data=WLC0), c("CHAR"), console=TRUE, alpha=0.05) #------- WLC4<-WLC%>% filter(TIME==4) anova(lm(CN~CHAR*LITTER, data=WLC4))# Litter:char HSD.test(aov(CN~CHAR*LITTER, data=WLC4), c("CHAR","LITTER"), console=TRUE, alpha=0.05) #------- WLC5<-WLC%>% filter(TIME==5) anova(lm(CN~CHAR*LITTER, data=WLC4))#2way: Litter:char HSD.test(aov(CN~CHAR*LITTER, data=WLC5), c("CHAR","LITTER"), console=TRUE, alpha=0.05) #------ WLC6<-WLC%>% filter(TIME==6) anova(lm(CN~CHAR*LITTER, data=WLC6)) HSD.test(aov(CN~CHAR*LITTER, data=WLC6), c("CHAR","LITTER"), console=TRUE, alpha=0.05) #------Adjust data for plotting CN_sd<-aggregate(WLC$CN, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=sd) CN_mean<-aggregate(WLC$CN, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=mean) CN_mean<-CN_mean%>% mutate(sd=CN_sd$x, Time=c("00","00","00","00","00","00", "02","02","02","02","02","02", "07","07","07","07","07","07", "13","13","13","13","13","13", "65","65","65","65","65","65", "80","80","80","80","80","80", "94","94","94","94","94","94"))%>% mutate(X=x, ymax=X+sd,ymin=X-sd, Litter=as.factor(Litter))%>% mutate(ymax=ifelse(x0,ymin==0,ymin),ymin)) levels(CN_mean$Litter)[c(1,3,2)] CN_mean<-CN_mean[order(CN_mean$Litter),] CN_mean<-CN_mean[order(CN_mean$Char),] ggplot(CN_mean, aes(CN_mean$Time, CN_mean$X)) + facet_grid(Char~Litter) + geom_bar(stat="identity") + guides(fill="none") + geom_errorbar(aes(ymax=CN_mean$ymax, ymin=CN_mean$ymin), width=0.2) + xlab("Time (days)") + ylab("C:N ratio") + scale_y_continuous(limits=c(0,50), breaks=c(0,10,20,30,40,50)) + scale_x_discrete(labels=c("0","65","80","94")) theme(panel.background=element_rect(fill="white"), panel.border = element_rect(linetype = "solid", fill = NA), strip.text.y=element_text(angle=0, size=12), legend.title = element_text(size=10), legend.text=element_text(size=8), legend.position=c(0.86,0.1)) ##########----------########## Mineral Nitrogen ##########----------########## #Call the data WLC<-read_csv("2018 WLC chemical data.csv")%>%#Read the WLC file transmute(TIME, REALTIME, CHAR, ID, AN=NO3+NH4, LITTER=ifelse(LITTER=="BLK", "Blank", ifelse(LITTER=="UL", "Unscorched", "Scorched"))) #Test for interactions anova(lm(AN~CHAR*LITTER, data=WLC)) WLC0<-WLC%>% filter(TIME==0) anova(lm(AN~CHAR*LITTER, data=WLC0))# HSD.test(aov(AN~LITTER, data=WLC0), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC1<-WLC%>% filter(TIME==1) anova(lm(AN~CHAR*LITTER, data=WLC1)) HSD.test(aov(AN~LITTER, data=WLC1), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC2<-WLC%>% filter(TIME==2) anova(lm(AN~CHAR*LITTER, data=WLC2)) #------ WLC3<-WLC%>% filter(TIME==3) anova(lm(AN~CHAR*LITTER, data=WLC3)) #------ WLC4<-WLC%>% filter(TIME==4) anova(lm(AN~CHAR*LITTER, data=WLC4)) HSD.test(aov(AN~LITTER, data=WLC4), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC5<-WLC%>% filter(TIME==5) anova(lm(AN~CHAR*LITTER, data=WLC5)) HSD.test(aov(AN~LITTER, data=WLC5), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC6<-WLC%>% filter(TIME==6) anova(lm(AN~CHAR*LITTER, data=WLC6)) HSD.test(aov(AN~LITTER, data=WLC6), c("LITTER"), console=TRUE, alpha=0.05) #------ #Adjust data for plotting AN_mean<-aggregate(WLC$NH4, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=mean) AN_mean$Type<-"NH4" AN_mean$Type<-factor(AN_mean$Type) NH4_sd<-aggregate(WLC$NH4, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=sd) colnames(NH4_sd)<-c("Litter","Char","Time","sd") AN_mean$sd<-NH4_sd$sd NO3<-aggregate(WLC$NO3, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=mean) NO3$Type<-"NO3" NO3$Type<-factor(NO3$Type) NO3_sd<-aggregate(WLC$NO3, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=sd) colnames(NO3_sd)<-c("Litter","Char","Time","sd") NO3$sd<-NO3_sd$sd AN_mean$joinx<-AN_mean$x+NO3$x NO3$joinx<-AN_mean$x+NO3$x AN_mean<-rbind(AN_mean,NO3) AN_mean$Time<-c("00","00","00","00","00","00", "02","02","02","02","02","02", "07","07","07","07","07","07", "13","13","13","13","13","13", "65","65","65","65","65","65", "80","80","80","80","80","80", "94","94","94","94","94","94", "00","00","00","00","00","00", "02","02","02","02","02","02", "07","07","07","07","07","07", "13","13","13","13","13","13", "65","65","65","65","65","65", "80","80","80","80","80","80", "94","94","94","94","94","94") AN_mean$ymax<-ifelse(AN_mean$Type=="NH4",AN_mean$X+AN_mean$sd,AN_mean$joinX+AN_mean$sd) AN_mean$ymin<-ifelse(AN_mean$Type=="NH4",AN_mean$X-AN_mean$sd,AN_mean$joinX-AN_mean$sd) AN_mean$Litter<-factor(AN_mean$Litter,levels(AN_mean$Litter)[c(1,3,2)]) AN_mean$ymax<-ifelse(AN_mean$x0,AN_mean$ymin==0,AN_mean$ymin),AN_mean$ymin) AN_mean<-AN_mean[order(AN_mean$Time),] AN_mean<-AN_mean[order(AN_mean$Litter),] AN_mean<-AN_mean[order(AN_mean$Char),] #Plotting ggplot(AN_mean, aes(AN_mean$Time, AN_mean$X)) + facet_grid(Char~Litter) + geom_bar(stat="identity") + scale_fill_manual(values=c("red","light blue","maroon","blue")) + geom_errorbar(aes(ymax=AN_mean$ymax, ymin=AN_mean$ymin), width=0.2) + xlab("Time (days)") + ylab(expression(paste("AN (",mu,"gN/gSoil)"))) + scale_y_continuous(limits=c(0,13),breaks=c(0,4,8,12),labels=c(0,4,8,12)) + scale_x_discrete(labels=c("0","2","7","13","65","80","94")) theme(panel.background=element_rect(fill="white"), panel.border = element_rect(linetype = "solid", fill = NA), strip.text.y=element_text(angle=0, size=12), legend.title = element_text(size=10), legend.text=element_text(size=8), legend.position=c(0.86,0.1)) ##########----------########## Potentially Mineralisable Nitrogen ##########----------########## WLC<-read_csv("2018 WLC chemical data.csv")%>%#Read the WLC file transmute(TIME, REALTIME, CHAR, ID, PMN, LITTER=ifelse(LITTER=="BLK", "Blank", ifelse(LITTER=="UL", "Unscorched", "Scorched"))) #Test for interactions anova(lm(PMN~CHAR*LITTER, data=WLC)) #------ WLC0<-WLC%>% filter(TIME==0) anova(lm(PMN~CHAR*LITTER, data=WLC0)) HSD.test(aov(PMN~LITTER, data=WLC0), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC1<-WLC%>% filter(TIME==1) anova(lm(PMN~CHAR*LITTER, data=WLC1)) HSD.test(aov(PMN~LITTER, data=WLC1), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC2<-WLC%>% filter(TIME==2) anova(lm(PMN~CHAR*LITTER, data=WLC2)) HSD.test(aov(PMN~LITTER, data=WLC2), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC3<-WLC%>% filter(TIME==3) anova(lm(PMN~CHAR*LITTER, data=WLC3)) HSD.test(aov(PMN~LITTER, data=WLC3), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC4<-WLC%>% filter(TIME==4) anova(lm(PMN~CHAR*LITTER, data=WLC4)) HSD.test(aov(PMN~CHAR*LITTER, data=WLC4), c("CHAR","LITTER"), console=TRUE, alpha=0.05) #------ WLC5<-WLC%>% filter(TIME==5) anova(lm(PMN~CHAR*LITTER, data=WLC5)) HSD.test(aov(PMN~LITTER, data=WLC5), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC6<-WLC%>% filter(TIME==6) anova(lm(PMN~CHAR*LITTER, data=WLC6)) HSD.test(aov(PMN~CHAR*LITTER, data=WLC6), c("CHAR","LITTER"), console=TRUE, alpha=0.05) #Adjust data for plotting PMN_sd<-aggregate(WLC$PMN, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=sd) PMN_mean<-aggregate(WLC$PMN, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=mean) PMN_mean<-PMN_mean%>% mutate(sd=PMN_sd$x, Time=c("00","00","00","00","00","00", "02","02","02","02","02","02", "07","07","07","07","07","07", "13","13","13","13","13","13", "65","65","65","65","65","65", "80","80","80","80","80","80", "94","94","94","94","94","94"))%>% mutate(X=x, ymax=X+sd,ymin=X-sd, Litter=as.factor(Litter))%>% mutate(ymax=ifelse(x0,ymin==0,ymin),ymin)) levels(PMN_mean$Litter)[c(1,3,2)] PMN_mean<-PMN_mean[order(PMN_mean$Litter),] PMN_mean<-PMN_mean[order(PMN_mean$Char),] #Plotting ggplot(PMN_mean, aes(PMN_mean$Time, PMN_mean$X)) + facet_grid(Char~Litter) + geom_bar(stat="identity") + scale_fill_manual(values=c("grey50","grey75")) + geom_errorbar(aes(ymax=PMN_mean$ymax, ymin=PMN_mean$ymin), width=0.2) + xlab("Time (days)") + ylab(expression(paste("PMN (",mu,"gN/gSoil)"))) + scale_y_continuous(limits=c(0,10.7),breaks=c(0,2,4,6,8,10),labels=c(0,2,4,6,8,10)) + scale_x_discrete(labels=c("0","2","7","13","65","80","94")) + theme(panel.background=element_rect(fill="white"), panel.border = element_rect(linetype = "solid", fill = NA), strip.text.y=element_text(angle=0, size=12), legend.title = element_text(size=10), legend.text=element_text(size=8), legend.position=c(0.86,0.1)) ##########----------########## Microbial Biomass Nitrogen ##########----------########## #Call the data WLC<-read_csv("2018 WLC chemical data.csv")%>%#Read the WLC file transmute(TIME, REALTIME, CHAR, ID, MBN=MBN/0.54, LITTER=ifelse(LITTER=="BLK", "Blank", ifelse(LITTER=="UL", "Unscorched", "Scorched"))) #Statistics #Test for interactions anova(lm(MBN~CHAR*LITTER, data=WLC)) WLC0<-WLC%>% filter(TIME==0) anova(lm(MBN~CHAR*LITTER, data=WLC0)) HSD.test(aov(MBN~CHAR*LITTER, data=WLC0), c("LITTER","CHAR"), console=TRUE, alpha=0.05) #------ WLC1<-WLC%>% filter(TIME==1) anova(lm(MBN~CHAR*LITTER, data=WLC1)) HSD.test(aov(MBN~CHAR*LITTER, data=WLC1), c("LITTER","CHAR"), console=TRUE, alpha=0.05) #------ WLC2<-WLC%>% filter(TIME==2) anova(lm(MBN~CHAR*LITTER, data=WLC2)) HSD.test(aov(MBN~LITTER, data=WLC1), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC3<-WLC%>% filter(TIME==3) anova(lm(MBN~CHAR*LITTER, data=WLC3)) HSD.test(aov(MBN~CHAR*LITTER, data=WLC3), c("LITTER","CHAR"), console=TRUE, alpha=0.05) #------ WLC4<-WLC%>% filter(TIME==4) anova(lm(MBN~CHAR*LITTER, data=WLC4)) HSD.test(aov(MBN~LITTER, data=WLC4), c("LITTER"), console=TRUE, alpha=0.05) #------ WLC5<-WLC%>% filter(TIME==5) anova(lm(MBN~CHAR*LITTER, data=WLC5)) WLC6<-WLC%>% filter(TIME==6) anova(lm(MBN~CHAR*LITTER, data=WLC6)) HSD.test(aov(MBN~CHAR*LITTER, data=WLC6), c("LITTER","CHAR"), console=TRUE, alpha=0.05) #Adjust data for plotting MBN_sd<-aggregate(WLC$MBN, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=sd) MBN_mean<-aggregate(WLC$MBN, list(Litter = WLC$LITTER, Char = WLC$CHAR, Time=WLC$TIME), FUN=mean) MBN_mean<-MBN_mean%>% mutate(sd=MBN_sd$x, Time=c("00","00","00","00","00","00", "02","02","02","02","02","02", "07","07","07","07","07","07", "13","13","13","13","13","13", "65","65","65","65","65","65", "80","80","80","80","80","80", "94","94","94","94","94","94"))%>% mutate(X=x, ymax=X+sd,ymin=X-sd, Litter=as.factor(Litter))%>% mutate(ymax=ifelse(x0,ymin==0,ymin),ymin)) levels(MBN_mean$Litter)[c(1,3,2)] MBN_mean<-MBN_mean[order(MBN_mean$Litter),] MBN_mean<-MBN_mean[order(MBN_mean$Char),] #Plotting ggplot(MBN_mean, aes(MBN_mean$Time, MBN_mean$X)) + facet_grid(Char~Litter) + geom_bar(stat="identity") + geom_errorbar(aes(ymax=MBN_mean$ymax, ymin=MBN_mean$ymin), width=0.2) + xlab("Time (days)") + ylab(expression(paste("MBN (",mu,"gN/gSoil)"))) + scale_y_continuous(limits=c(0,17),breaks=c(0,4,8,12,16),labels=c(0,4,8,12,16)) + scale_x_discrete(labels=c("0","2","7","13","65","80","94")) + theme(panel.background=element_rect(fill="white"), panel.border = element_rect(linetype = "solid", fill = NA), strip.text.y=element_text(angle=0, size=12), legend.title = element_text(size=10), legend.text=element_text(size=8), legend.position=c(0.86,0.9))