chapter 4 version 10 NRCS graphs code

1 Set up

1.1 1) Load packages

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(tidyverse)
library(janitor)
library(corrplot)
library(lavaan)       #
library(semPlot)
library(plyr)


#mapping
# library(ggmap)
# library(maps)
# library(mapdata)    
# library(ggspatial)
# library(gridExtra)

1.2 2) Import data

#import data
soilhealth <- read_excel("7-30-21 part 11 nrcs soil health final data.xlsx") 

#change column names
soilhealth <- soilhealth %>% 
  clean_names()
#str(soilhealth)
#View(soilhealth)

#convert variables to factors
soilhealth$nhorizon <- as.factor(soilhealth$nhorizon)
soilhealth$blk <- as.factor(soilhealth$blk)
soilhealth$horizon <- as.factor(soilhealth$horizon)
soilhealth$thorizon <- as.factor(soilhealth$thorizon)
soilhealth$rlanduse <- as.factor(soilhealth$rlanduse)
soilhealth$landuse <- as.factor(soilhealth$landuse)
soilhealth$prec <- as.factor(soilhealth$prec)
soilhealth$ag <- as.factor(soilhealth$ag)
soilhealth$tillage <- as.factor(soilhealth$tillage)

#convert variables to numeric
soilhealth$bdepth=as.numeric(soilhealth$bdepth)
soilhealth$precip <- as.numeric(soilhealth$precip)
soilhealth$t_np <- as.numeric(soilhealth$t_np)
soilhealth$to_cp <- as.numeric(soilhealth$to_cp)
soilhealth$ca <- as.numeric(soilhealth$ca)
soilhealth$cu <- as.numeric(soilhealth$cu)
soilhealth$mg <- as.numeric(soilhealth$mg)
soilhealth$mn <- as.numeric(soilhealth$mn)
soilhealth$na <- as.numeric(soilhealth$na)
soilhealth$p <- as.numeric(soilhealth$p)
soilhealth$p_h <- as.numeric(soilhealth$p_h)
soilhealth$k <- as.numeric(soilhealth$k)
soilhealth$zn <- as.numeric(soilhealth$zn)
soilhealth$fe <- as.numeric(soilhealth$fe)
soilhealth$cec <- as.numeric(soilhealth$cec)
soilhealth$socstock <- as.numeric(soilhealth$socstock)
soilhealth$tnstock <- as.numeric(soilhealth$tnstock)
soilhealth$mb <- as.numeric(soilhealth$mb)
soilhealth$gramp <- as.numeric(soilhealth$gramp)
soilhealth$gramn <- as.numeric(soilhealth$gramn)
soilhealth$actin <- as.numeric(soilhealth$actin)
soilhealth$amf <- as.numeric(soilhealth$amf)
soilhealth$fungi <- as.numeric(soilhealth$fungi)
soilhealth$fb <- as.numeric(soilhealth$fb)


#vartable(soilhealth)

1.3 3) Data Wrangling

# removes ottawa
soilhealth1 <- soilhealth %>%
  filter(location!="Ottawa") 

#str(soilhealth1)

soilhealth1 <- as.data.frame(soilhealth1)

# keep nutrients

nut <- soilhealth1 %>%
 dplyr::filter(location!="Ottawa") %>%
  dplyr::select(location, treatment, bdepth, horizon, blk, replication, clay, nh4, no3, t_np, to_cp, ca, cu, mg, mn, na, p, p_h, k, zn, fe, gwc, cec, socstock, tnstock) %>%
  dplyr::mutate(cn=(to_cp/t_np), ca_g=(ca/1000))
#str(nut)
#summary(nut)
#vartable(nut)


# keeps aggregate data
agg <- soilhealth1 %>%
  dplyr::select(location, treatment, bdepth, horizon, blk, replication, x20wsa2000, x20wsa250, x20wsa53, x20wsa20, x20mwd, x5wsa2000, x5wsa250, x5wsa53, x5wsa20, x5mwd, nagg, bd)
#str(agg)
#summary(agg)
#plot(agg)
#vartable(agg)


# Keep respiration data, enzymes, proteins
# normalize the enzyme distributions

#vartable(agg)

epr <- soilhealth1 %>%
  dplyr::select(location, treatment, bdepth, horizon, blk, replication, ac, glucosidase, glucosaminidase, acidphosphotase, alkphosphatase, arylsulfatase, phosphodiesterase, protein, respiration, mb, gramp, gramn, actin, amf, fungi, fb, pgrmp, pgrmn, pactin, pamf, pfungi) %>%
  dplyr::mutate(lnac=log(ac), lngluc=log(glucosidase), lnglucosam=log(glucosaminidase), lnacidp=log(acidphosphotase), lnalkp=log(alkphosphatase), lnary=log(arylsulfatase), lnpho=log(phosphodiesterase), lnprotein=log(protein), lngluc_glucosam= (lngluc/lnglucosam), lngluc_acidp= (lngluc/lnacidp), lngluc_alkp=(lngluc/lnalkp), lngluc_pho=(lngluc/lnpho), lngluc_ary=(lngluc/lnary))
#str(epr)
#vartable(epr)


ratios <- soilhealth1 %>%
  dplyr::select(location, treatment, bdepth, horizon, blk, replication, respiration, to_cg, t_nm, protein) %>%
  dplyr::mutate(t_ng=(t_nm/1000), ptn=(protein/t_ng), rec=(respiration/to_cg))
#vartable(ratios)
#To deal with egative natural log values from enzymes, I added 1.1 to each variable.

epr2 <- soilhealth1 %>%
  dplyr::select(location, treatment, bdepth, horizon, blk, replication, ac, glucosidase, glucosaminidase, acidphosphotase, alkphosphatase, arylsulfatase, phosphodiesterase, protein, respiration) %>%
  mutate(lnac=log(ac), lngluc=log(glucosidase*100), lnglucosam=log(glucosaminidase*100), lnacidp=log(acidphosphotase*100), lnalkp=log(alkphosphatase*100), lnary=log(arylsulfatase*100), lnpho=log(phosphodiesterase*100), lnprotein=log(protein), lngluc_glucosam= (lngluc/lnglucosam), lngluc_acidp= (lngluc/lnacidp), lngluc_alkp=(lngluc/lnalkp), lngluc_pho=(lngluc/lnpho), lngluc_ary=(lngluc/lnary))



#summary(epr1)

1.4 4) Theme James

library(grid)
library(ggthemes)

#fcol <- c( "AG" = "black", "EA" = "grey40", "NP" = "grey90", "IR"="blue")
#kstate purple - 512888
fcol <- c( "AG" = "#FF3333", "EA" = "#FF9933", "NP" = "#512888", "IR"="#0033FF")

pfcol <- c( "AG" = "#FF3333", "EA" = "#FF9933", "NP" = "#512888")


pd <- position_dodge(0.1)

theme_James <- function(base_size=14, base_family="TT Times New Roman") {
  (theme_foundation(base_size=base_size, base_family=base_family)+
      theme_bw()+ 
      theme(panel.background = element_rect(colour = NA),
            plot.background = element_rect(colour = NA),
            axis.title = element_text(color="black",size=rel(1.2)),
            axis.text = element_text(color="black", size = 12),
            legend.key = element_rect(colour = NA),
            legend.spacing = unit(0, "cm"),
            legend.text = element_text(size=12),
            legend.title = element_blank(),
            panel.grid = element_blank(),
            plot.title = element_text(color="Black",size = rel(1.5),face = "bold",hjust = 0.5),
            strip.text = element_text(color="Black",size = rel(1),face="bold")
      ))
}

2 Extra

2.1 Correlation Matrix

library(ggcorrplot)

# Compute a correlation matrix

names(soilhealth)

##  [1] "sample_name"       "location"          "treatment"        
##  [4] "depth"             "replication"       "bdepth"           
##  [7] "nhorizon"          "blk"               "horizon"          
## [10] "thorizon"          "rlanduse"          "landuse"          
## [13] "elevation"         "clay"              "temp"             
## [16] "precip"            "prec"              "ag"               
## [19] "tillage"           "nh4"               "no3"              
## [22] "t_np"              "t_nm"              "to_cp"            
## [25] "to_cg"             "ca"                "cu"               
## [28] "mg"                "mn"                "na"               
## [31] "p"                 "p_h"               "carbonate"        
## [34] "k"                 "zn"                "fe"               
## [37] "cec"               "gwc"               "odad"             
## [40] "bd"                "ac"                "x20wsa2000"       
## [43] "x20wsa250"         "x20wsa53"          "x20wsa20"         
## [46] "x20mwd"            "x5wsa2000"         "x5wsa250"         
## [49] "x5wsa53"           "x5wsa20"           "x5mwd"            
## [52] "nagg"              "respiration"       "mb"               
## [55] "gramp"             "gramn"             "actin"            
## [58] "amf"               "fungi"             "fb"               
## [61] "glucosidase"       "glucosaminidase"   "acidphosphotase"  
## [64] "alkphosphatase"    "arylsulfatase"     "phosphodiesterase"
## [67] "protein"           "socstock"          "tnstock"          
## [70] "pgrmp"             "pgrmn"             "pactin"           
## [73] "pamf"              "pfungi"

soilmatrix <- soilhealth %>%
  filter(location!="Ottawa", treatment!="IR") %>% 
  dplyr::select(-sample_name, -blk, -bdepth, -elevation, -horizon, -prec, -carbonate, - precip, -nhorizon, -odad, -gwc, -temp, -to_cg, -t_nm, -x5wsa53, -x5wsa20, -x20wsa53, -x20wsa20, -nh4, -no3, -x20wsa2000, -x20wsa250, -x5mwd, -nagg, -x5wsa250,-socstock, -tnstock, -clay, -x5wsa2000, -ca, -cu, -mg, -mn, -na, -zn, -fe, -k) %>%
  dplyr::select(where(is.numeric)) %>% drop_na()

soilmatrix <- na.omit(soilmatrix)

names(soilmatrix)

##  [1] "t_np"              "to_cp"             "p"                
##  [4] "p_h"               "cec"               "bd"               
##  [7] "ac"                "x20mwd"            "respiration"      
## [10] "mb"                "gramp"             "gramn"            
## [13] "actin"             "amf"               "fungi"            
## [16] "fb"                "glucosidase"       "glucosaminidase"  
## [19] "acidphosphotase"   "alkphosphatase"    "arylsulfatase"    
## [22] "phosphodiesterase" "protein"           "pgrmp"            
## [25] "pgrmn"             "pactin"            "pamf"             
## [28] "pfungi"

colnames(soilmatrix)[1] <- "TN"
colnames(soilmatrix)[2] <- "SOC"
colnames(soilmatrix)[3] <- "P"
colnames(soilmatrix)[4] <- "pH"
colnames(soilmatrix)[5] <- "CEC"
colnames(soilmatrix)[6] <- "bd"
colnames(soilmatrix)[7] <- "POXC"
colnames(soilmatrix)[8] <- "MWD"
colnames(soilmatrix)[9] <- "Resp"
colnames(soilmatrix)[10] <- "MB"
colnames(soilmatrix)[11] <- "GmP"
colnames(soilmatrix)[12] <- "GmN"
colnames(soilmatrix)[13] <- "Act"
colnames(soilmatrix)[14] <- "AMF"
colnames(soilmatrix)[15] <- "Fun"
colnames(soilmatrix)[16] <- "FB"
colnames(soilmatrix)[17] <- "bG"
colnames(soilmatrix)[18] <- "NAG"
colnames(soilmatrix)[19] <- "AP"
colnames(soilmatrix)[20] <- "ALK"
colnames(soilmatrix)[21] <- "ARY"
colnames(soilmatrix)[22] <- "PHO"
colnames(soilmatrix)[23] <- "Pro"

soilmatrix <- soilmatrix %>% 
  dplyr::select(-CEC)

corr <- round(cor(soilmatrix, method = "pearson"), 2)

#head(corr[,])

# Compute a matrix of correlation p-values
p.mat <- cor_pmat(soilmatrix , conf.level = 0.95)

#head(p.matr[,])

#colnames(corr) <- c("Bulk Density", "Potassium", "Phosphorus", "Manganese", "pH", "Calcium", "Clay", "5 min WSA 250", "20 min MWD", "20 min WSA 2000", "NRCS", "5 min MWD", "5 min WSA 2000", "Iron", "Sodium", "Copper", "Magnesium", "20 min WSA 250", "Acid Phosphatase", "Glucosaminidase", "Total Carbon", "Total Nitrogen", "Zinc", "Respiration", "Phosphodiesterase", "Arylsulfatase", "Alkaline Phosphatase", "Protein", "Glucosidase", "Active Carbon")

#rownames(corr) <- c("Bulk Density", "Potassium", "Phosphorus", "Manganese", "pH", "Calcium", "Clay", "5 min WSA 250", "20 min MWD", "20 min WSA 2000", "NRCS", "5 min MWD", "5 min WSA 2000", "Iron", "Sodium", "Copper", "Magnesium", "20 min WSA 250", "Acid Phosphatase", "Glucosaminidase", "Total Carbon", "Total Nitrogen", "Zinc", "Respiration", "Phosphodiesterase", "Arylsulfatase", "Alkaline Phosphatase", "Protein", "Glucosidase", "Active Carbon")

# Types of correlogram layout
# Get the lower triangle type=lower
## using hierarchical clustering for hc.order
#ggcorrplot(corr, hc.order = TRUE, type = "lower")
#outline white
#ggcorrplot(corr, hc.order = TRUE, type = "lower", outline.col = "white")
#Add correlation labels
#ggcorrplot(corr, hc.order = TRUE, type = "lower", outline.col = "white", lab=TRUE)
#Remove non-significance
#ggcorrplot(corr, hc.order = TRUE, type = "lower", p.mat=p.mat, outline.col = "white", insig = "blank", lab=TRUE)

ggcorrplot(corr, hc.order = TRUE, type = "lower", p.mat=p.mat, outline.col = "white", insig = "blank", lab=TRUE, lab_size = 2)

corpsoil <- ggcorrplot(corr, hc.order = TRUE, type = "lower", p.mat=p.mat, outline.col = "white", insig = "blank", lab=TRUE, lab_size = 2, title = "Spearman Correlation Matrix p > 0.05") 

ggsave("corplot.png", corpsoil, width=15, height=10)


#Correlation with p values
  
#p.mat <-  soilmatrix %>% dplyr::select(where(is.numeric)) %>% cor_pmat()

#soilmatrix %>% dplyr::select(where(is.numeric)) %>% cor() %>% ggcorrplot(hc.order = TRUE, type = "lower", p.mat = p.mat, insig = "blank", lab= TRUE)



ggcorrplot(corr, sig.level = 0.05, method = "circle", hc.order = F, outline.color = "white", 
    type = "lower", lab = T, lab_size = 3, colors = c("tomato2", "white", "springgreen3"), 
    p.mat = p.mat, insig = "blank", ggtheme = theme_bw) + theme(panel.grid = element_blank(), 
    legend.position = c(0.1, 0.8), axis.text = element_text(face = "bold"))

ggcorrplot(corr, sig.level = 0.05, hc.order = F, type = "lower", lab=T, lab_size= 3,   colors = c("tomato2", "white", "springgreen3"), p.mat=p.mat, outline.col = "white", insig = "blank") + theme(panel.grid = element_blank(),
    legend.position = c(0.1, 0.8), axis.text = element_text(face = "bold"))

corsoil <- ggcorrplot(corr, tl.cex = 20, sig.level = 0.05, hc.order = F, type = "lower", lab=T, lab_size= 4.5,   colors = c("tomato2", "white", "springgreen3"), p.mat=p.mat, outline.col = "white", insig = "blank") + theme(panel.grid = element_blank(),
    legend.position = c(0.1, 0.8), axis.text = element_text(face = "bold"))

ggsave("ncorsoilplot.png", corsoil, width=10, height=10)

3 Biological Graphs

3.1 B-Glucosidase

#mean, standard deviation, and std error
bgluc<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(glucosidase)), mean=mean(glucosidase, na.rm=T), sd = sd(glucosidase, na.rm=T), se = sd / sqrt(N)) 


bgluc$location_f =factor(bgluc$location, levels=c('Tribune', 'Hays', 'Manhattan'))
bgluc$depth =factor(bgluc$bdepth, levels=c(0,20,40,60,80,100))


ggplot(bgluc, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(beta~-Glucosidase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(bgluc, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(beta~-Glucosidase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

#+   ggsave("gluc1.png",  height=6, width=9)

##N-acetyl-β-D-glucosaminidase

#mean, standard deviation, and std error
bglum<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(glucosaminidase)), mean=mean(glucosaminidase, na.rm=T), sd = sd(glucosaminidase, na.rm=T), se = sd / sqrt(N)) 


bglum$location_f =factor(bglum$location, levels=c('Tribune', 'Hays', 'Manhattan'))
bglum$depth =factor(bglum$bdepth, levels=c(0,20,40,60,80,100))


ggplot(bglum, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,25) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(N~-acetyl~-beta~-Glucosaminidase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(bglum, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,25) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(N~-acetyl~-beta~-Glucosaminidase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+ggsave("glucosam1.png",  height=6, width=9)

3.2 Acid Phosphatase

#mean, standard deviation, and std error
acpho<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(acidphosphotase
)), mean=mean(acidphosphotase, na.rm=T), sd = sd(acidphosphotase, na.rm=T), se = sd / sqrt(N)) 


acpho$location_f =factor(acpho$location, levels=c('Tribune', 'Hays', 'Manhattan'))
acpho$depth =factor(acpho$bdepth, levels=c(0,20,40,60,80,100))


ggplot(acpho, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,400) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Acid~ Phosphotase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(acpho, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,400) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Acid~ Phosphatase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+  ggsave("acpho1.png",  height=6, width=9)

3.3 Alkaline Phosphatase

#mean, standard deviation, and std error
akpho<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(alkphosphatase
)), mean=mean(alkphosphatase, na.rm=T), sd = sd(alkphosphatase, na.rm=T), se = sd / sqrt(N)) 


akpho$location_f =factor(akpho$location, levels=c('Tribune', 'Hays', 'Manhattan'))
akpho$depth =factor(akpho$bdepth, levels=c(0,20,40,60,80,100))


ggplot(akpho, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,100) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Alkaline~ Phosphotase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(akpho, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,100) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Alkaline~ Phosphotase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

#+  ggsave("akpho1.png",  height=6, width=9)

3.4 Arylsulfatase

#mean, standard deviation, and std error
ary<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(arylsulfatase
)), mean=mean(arylsulfatase, na.rm=T), sd = sd(arylsulfatase, na.rm=T), se = sd / sqrt(N)) 


ary$location_f =factor(ary$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ary$depth =factor(ary$bdepth, levels=c(0,20,40,60,80,100))


ggplot(ary, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Arylsulfatase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(ary, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Arylsulfatase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+ggsave("ary1.png",  height=6, width=9)

3.5 Phosphodiesterase

#mean, standard deviation, and std error
phos<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(phosphodiesterase
)), mean=mean(phosphodiesterase, na.rm=T), sd = sd(phosphodiesterase, na.rm=T), se = sd / sqrt(N)) 


phos$location_f =factor(phos$location, levels=c('Tribune', 'Hays', 'Manhattan'))
phos$depth =factor(phos$bdepth, levels=c(0,20,40,60,80,100))


ggplot(phos, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,80) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Phosphodiesterase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(phos, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,80) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Phosphodiesterase ~(mg ~kg^{-1} ~hr^{-1}))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+  ggsave("phos1.png",  height=6, width=9)

3.6 Protein Content

#mean, standard deviation, and std error
protein<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(protein
)), mean=mean(protein, na.rm=T), sd = sd(protein, na.rm=T), se = sd / sqrt(N)) 


protein$location_f =factor(protein$location, levels=c('Tribune', 'Hays', 'Manhattan'))
protein$depth =factor(protein$bdepth, levels=c(0,20,40,60,80,100))


ggplot(protein, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,15) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Protein ~Content ~(g~protein ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(protein, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,15) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Protein ~Content ~(g~protein ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+
 # ggsave("protein1.png",  height=6, width=9)

3.7 Permanganate Oxidizable Organic Carbon

#mean, standard deviation, and std error
ac<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(ac
)), mean=mean(ac, na.rm=T), sd = sd(ac, na.rm=T), se = sd / sqrt(N)) 


ac$location_f =factor(ac$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ac$depth =factor(ac$bdepth, levels=c(0,20,40,60,80,100))


ggplot(ac, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,100) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Active ~Carbon~(mg ~Oxidizable ~C ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.2,0.3))

ggplot(ac, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,100) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Active ~Carbon ~(mg ~Oxidizable ~C ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.2,0.3))

#+  ggsave("ac1.png",  height=6, width=9)

3.8 Soil Respiration

#mean, standard deviation, and std error
respiration<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(respiration
)), mean=mean(respiration, na.rm=T), sd = sd(respiration, na.rm=T), se = sd / sqrt(N)) 


respiration$location_f =factor(respiration$location, levels=c('Tribune', 'Hays', 'Manhattan'))
respiration$depth =factor(respiration$bdepth, levels=c(0,20,40,60,80,100))


ggplot(respiration, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Respiration ~(mg ~CO[2]~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(respiration, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Respiration ~(mg ~CO[2]~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

#+  ggsave("respiration1.png",  height=6, width=9)

3.9 Microbial Biomass

#mean, standard deviation, and std error
mb<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(mb
)), mean=mean(mb, na.rm=T), sd = sd(mb, na.rm=T), se = sd / sqrt(N)) 


mb$location_f =factor(mb$location, levels=c('Tribune', 'Hays', 'Manhattan'))
mb$depth =factor(mb$bdepth, levels=c(0,20,40,60,80,100))

mbir <- mb %>%
  filter(treatment!="IR")

ggplot(mbir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Microbial ~Biomass ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=pfcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(mb, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Microbial ~Biomass ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+  ggsave("microbialbiomass1.png",  height=6, width=9)

3.10 Gram Negative Bacteria

#mean, standard deviation, and std error
gramn<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(gramn)), mean=mean(gramn, na.rm=T), sd = sd(gramn, na.rm=T), se = sd / sqrt(N)) 


gramn$location_f =factor(gramn$location, levels=c('Tribune', 'Hays', 'Manhattan'))
gramn$depth =factor(gramn$bdepth, levels=c(0,20,40,60,80,100))

gramnir <- gramn %>%
  filter(treatment!="IR")


ggplot(gramnir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Gram~"(-)" ~Bacteria ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(gramn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Gram~"(-)" ~Bacteria ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

#+  ggsave("gramneg1.png",  height=6, width=9)

3.11 Gram Positive Bacteria

#mean, standard deviation, and std error
gramp<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(gramp)), mean=mean(gramp, na.rm=T), sd = sd(gramp, na.rm=T), se = sd / sqrt(N)) 


gramp$location_f =factor(gramp$location, levels=c('Tribune', 'Hays', 'Manhattan'))
gramp$depth =factor(gramp$bdepth, levels=c(0,20,40,60,80,100))

grampir <- gramp %>%
  filter(treatment!="IR")


ggplot(grampir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Gram~("+") ~Bacteria ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(gramp, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Gram~("+") ~Bacteria ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

# +  ggsave("grampos1.png",  height=6, width=9)

3.12 Actinomycetes

#mean, standard deviation, and std error
actin<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(actin)), mean=mean(actin, na.rm=T), sd = sd(actin, na.rm=T), se = sd / sqrt(N)) 


actin$location_f =factor(actin$location, levels=c('Tribune', 'Hays', 'Manhattan'))
actin$depth =factor(actin$bdepth, levels=c(0,20,40,60,80,100))

actinir <- actin %>%
  filter(treatment!="IR")

ggplot(actinir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Actinomycetes ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(actin, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Actinomycetes ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

#+ ggsave("actino1.png",  height=6, width=9)

3.13 Arbuscular Mycorrhizal Fungi

#mean, standard deviation, and std error
amf<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(amf
)), mean=mean(amf, na.rm=T), sd = sd(amf, na.rm=T), se = sd / sqrt(N)) 


amf$location_f =factor(amf$location, levels=c('Tribune', 'Hays', 'Manhattan'))
amf$depth =factor(amf$bdepth, levels=c(0,20,40,60,80,100))

amfir <- amf %>%
  filter(treatment!="IR")


ggplot(amfir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~AMF ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(respiration, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~AMF ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+  ggsave("amf1.png",  height=6, width=9)

3.14 Saprophytic Fungi

#mean, standard deviation, and std error
fungi<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(fungi)), mean=mean(fungi, na.rm=T), sd = sd(fungi, na.rm=T), se = sd / sqrt(N)) 


fungi$location_f =factor(fungi$location, levels=c('Tribune', 'Hays', 'Manhattan'))
fungi$depth =factor(fungi$bdepth, levels=c(0,20,40,60,80,100))

fungiir <- fungi %>%
  filter(treatment!="IR")


ggplot(fungiir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Saprophytic ~Fungi ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(fungi, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Saprophytic ~Fungi ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3)) 

# +  ggsave("sapfungi1.png",  height=6, width=9)

3.15 Fungal to Bacteria ratio

#mean, standard deviation, and std error
fb<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(fb
)), mean=mean(fb, na.rm=T), sd = sd(fb, na.rm=T), se = sd / sqrt(N)) 


fb$location_f =factor(fb$location, levels=c('Tribune', 'Hays', 'Manhattan'))
fb$depth =factor(fb$bdepth, levels=c(0,20,40,60,80,100))

fbir <- fb %>%
  filter(treatment!="IR")


ggplot(fbir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Fungi ~to ~Bacteria ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.60,0.54))

ggplot(fb, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Fungi ~to ~Bacteria ~(nmol ~PLFA ~g^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+  ggsave("fb1.png",  height=6, width=9)

3.16 Composition Gram Negative Bacteria

#mean, standard deviation, and std error
pgrmn<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(pgrmn)), mean=mean(pgrmn, na.rm=T), sd = sd(pgrmn, na.rm=T), se = sd / sqrt(N)) 


pgrmn$location_f =factor(pgrmn$location, levels=c('Tribune', 'Hays', 'Manhattan'))
pgrmn$depth =factor(pgrmn$bdepth, levels=c(0,20,40,60,80,100))

pgrmnir <- pgrmn %>%
  filter(treatment!="IR")


ggplot(pgrmnir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
    ylab(expression(Gram~"(-)" ~Bacteria ~Composition ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#ggsave("pgramneg1.png",  height=6, width=9)

3.17 Composition Gram Positive Bacteria

#mean, standard deviation, and std error
pgrmp<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(pgrmp)), mean=mean(pgrmp, na.rm=T), sd = sd(pgrmp, na.rm=T), se = sd / sqrt(N)) 


pgrmp$location_f =factor(pgrmp$location, levels=c('Tribune', 'Hays', 'Manhattan'))
pgrmp$depth =factor(pgrmp$bdepth, levels=c(0,20,40,60,80,100))

pgrmpir <- pgrmp %>%
  filter(treatment!="IR")


ggplot(pgrmpir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Gram~("+") ~Bacteria ~Composition ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.77,0.3))

#  ggsave("pgrampos1.png",  height=6, width=9)

3.18 Composition Actinomycetes

#mean, standard deviation, and std error
pactin<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(pactin)), mean=mean(pactin, na.rm=T), sd = sd(pactin, na.rm=T), se = sd / sqrt(N)) 


pactin$location_f =factor(pactin$location, levels=c('Tribune', 'Hays', 'Manhattan'))
pactin$depth =factor(pactin$bdepth, levels=c(0,20,40,60,80,100))

pactinir <- pactin %>%
  filter(treatment!="IR")

ggplot(pactinir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Actinomycetes ~Composition ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#ggsave("pactino1.png",  height=6, width=9)

3.19 Composition Arbuscular Mycorrhizal Fungi

#mean, standard deviation, and std error
pamf<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(pamf
)), mean=mean(pamf, na.rm=T), sd = sd(pamf, na.rm=T), se = sd / sqrt(N)) 

pamf$location_f =factor(pamf$location, levels=c('Tribune', 'Hays', 'Manhattan'))
pamf$depth =factor(pamf$bdepth, levels=c(0,20,40,60,80,100))

pamfir <- pamf %>%
  filter(treatment!="IR")


ggplot(pamfir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Arbuscular ~Mycorrhizal ~Fungi ~Composition ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.57,0.7))

# ggsave("pamf1.png",  height=6, width=9)

3.20 Composition Saprophytic Fungi

#mean, standard deviation, and std error
pfungi<- ddply(epr, c("location","treatment","bdepth"), mutate, N=sum(!is.na(pfungi)), mean=mean(pfungi, na.rm=T), sd = sd(pfungi, na.rm=T), se = sd / sqrt(N)) 


pfungi$location_f =factor(pfungi$location, levels=c('Tribune', 'Hays', 'Manhattan'))
pfungi$depth =factor(pfungi$bdepth, levels=c(0,20,40,60,80,100))

pfungiir <- pfungi %>%
  filter(treatment!="IR")


ggplot(pfungiir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(~Saprophytic ~Fungi ~Composition ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.7))

# ggsave("psapfungi1.png",  height=6, width=9)

4 Chemical Graphs

4.1 Total Nitrogen

#mean, standard deviation, and std error
tnp<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(t_np)), mean=mean(t_np, na.rm=T), sd = sd(t_np, na.rm=T), se = sd / sqrt(N)) 


tnp$location_f =factor(tnp$location, levels=c('Tribune', 'Hays', 'Manhattan'))
tnp$depth =factor(tnp$bdepth, levels=c(0,20,40,60,80,100))


ggplot(tnp, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Total~ Nitrogen ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

ggplot(tnp, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Total~ Nitrogen ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

#+ ggsave("tnp1.png",  height=6, width=9)

4.2 Total Nitrogen stock

#mean, standard deviation, and std error
tnstock<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(tnstock)), mean=mean(tnstock, na.rm=T), sd = sd(tnstock, na.rm=T), se = sd / sqrt(N)) 


tnstock$location_f =factor(tnstock$location, levels=c('Tribune', 'Hays', 'Manhattan'))
tnstock$depth =factor(tnstock$bdepth, levels=c(0,20,40,60,80,100))


tnstockir <- tnstock %>%
  filter(treatment!="IR")

ggplot(tnstockir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Total~ Nitrogen ~Stocks ~(Mg ~N ~ha^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.75,0.5))

ggplot(tnstockir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Total~ Nitrogen ~Stocks ~(Mg ~N ~ha^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

#+  ggsave("tnstock1.png",  height=6, width=9)

4.3 TN stocks graph

#sum all depths per treatment for carbon stocks

tnstock_bar<- ddply(tnstock, c("location","treatment", "bdepth"), mutate,N=sum(!is.na(tnstock)), mean=mean(tnstock, na.rm=T), sd = sd(tnstock, na.rm=T), se = sd / sqrt(N)) 
 
 tnstock_bar=tnstock_bar[seq(1, nrow(tnstock_bar), 4), ]

 tnstock_bar1 <- tnstock_bar %>%
    dplyr::select(location, treatment, replication, bdepth, socstock, mean,  sd, se )

 
tnstock_bar11<- ddply(tnstock_bar1, c("location","treatment"), mutate,N=sum(!is.na(mean)), sum=sum(mean, na.rm=T)) 

 tnstock_barr11 <- tnstock_bar11[c(1,8,15,21,28,35,42,48,61),]

tnstock_barsum<- tnstock_barr11%>%
    dplyr::select(location, treatment, sum)
tnstock_barsum

##     location treatment       sum
## 1       Hays        AG  9.605797
## 8       Hays        EA 13.934172
## 15      Hays        NP 14.907049
## 21 Manhattan        AG 15.655787
## 28 Manhattan        EA 14.194287
## 35 Manhattan        NP 16.675074
## 42   Tribune        AG 10.723844
## 48   Tribune        EA 11.075330
## 61   Tribune        NP 12.955694

sfcol <- c( "AG" = "#FF3333", "EA" = "#FF9933", "NP" = "#512888")

 #colbio <- c( "AG" = "#ff6500", "EA" = "deepskyblue1", "NP" = "green4")

tnstock_barsum$location_f =factor(tnstock_barsum$location, levels=c('Tribune', 'Hays', 'Manhattan'))

 ggplot(data=tnstock_barsum, aes(x=treatment, y=sum, fill = treatment)) + 
  geom_bar(position=position_dodge(), stat="identity", colour = "black") +
  scale_y_continuous(limits=c(0,25)) +  
  facet_wrap(facets=vars(location_f), strip.position="bottom")  +
  xlab("")+
  scale_fill_manual(values = sfcol) + 
     ggtitle("Total TN Stocks")+ 
  theme_James() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        strip.background=element_rect(size=0.5, colour = "black"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.position= c(0.1,0.86)) +
  ylab(expression(Total ~TN ~Stocks ~(Mg ~N ~ha^{-1} ~soil))) +
  guides(fill = guide_legend(override.aes = aes(label="")))+
 geom_label(aes(label=round(sum,0), y = sum+1),
            label.padding = unit(.3,"lines"), fill=NA, show.legend=NA, label.size = NA, font="bold" )

4.4 Soil Organic Carbon

#mean, standard deviation, and std error
to_cp<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(to_cp)), mean=mean(to_cp, na.rm=T), sd = sd(to_cp, na.rm=T), se = sd / sqrt(N)) 


to_cp$location_f =factor(to_cp$location, levels=c('Tribune', 'Hays', 'Manhattan'))
to_cp$depth =factor(to_cp$bdepth, levels=c(0,20,40,60,80,100))


ggplot(to_cp, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Soil~ Organic ~Carbon ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

ggplot(to_cp, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Soil~ Organic ~Carbon ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

#+ ggsave("tcp1.png",  height=6, width=9)

4.5 Soil Organic Carbon stocks

#mean, standard deviation, and std error
socstock<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(socstock)), mean=mean(socstock, na.rm=T), sd = sd(socstock, na.rm=T), se = sd / sqrt(N)) 


socstock$location_f =factor(socstock$location, levels=c('Tribune', 'Hays', 'Manhattan'))
socstock$depth =factor(socstock$bdepth, levels=c(0,20,40,60,80,100))

socstockir <- socstock %>%
  filter(treatment!="IR")

ggplot(socstockir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Soil~ Organic ~Carbon ~Stocks ~(Mg ~C ~ha^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.23,0.8))

ggplot(socstockir, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Soil~ Organic ~Carbon ~Stocks ~(Mg ~C ~ha^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

#+ ggsave("socstockp1.png",  height=6, width=9)

4.6 SOC stocks graph

#sum all depths per treatment for carbon stocks

socstock_bar<- ddply(socstock, c("location","treatment", "bdepth"), mutate,N=sum(!is.na(socstock)), mean=mean(socstock, na.rm=T), sd = sd(socstock, na.rm=T), se = sd / sqrt(N)) 
 
 socstock_barr=socstock_bar[seq(1, nrow(socstock_bar), 4), ]

 socstock_bar1 <- socstock_barr %>%
    dplyr::select(location, treatment, replication, bdepth, socstock, mean,  sd, se )

socstock_bar11<- ddply(socstock_bar1, c("location","treatment"), mutate,N=sum(!is.na(mean)), sum=sum(mean, na.rm=T)) 
 
 socstock_barr11 <- socstock_bar11[c(1,8,15,21,28,35,42,48,61),]

socstock_barsum<- socstock_barr11%>%
    dplyr::select(location, treatment, sum)
socstock_barsum

##     location treatment       sum
## 1       Hays        AG  80.33702
## 8       Hays        EA  98.93570
## 15      Hays        NP 134.47705
## 21 Manhattan        AG 148.80752
## 28 Manhattan        EA 149.22606
## 35 Manhattan        NP 185.81541
## 42   Tribune        AG 103.66324
## 48   Tribune        EA 111.08303
## 61   Tribune        NP 139.13003

sfcol <- c( "AG" = "#FF3333", "EA" = "#FF9933", "NP" = "#512888")

 #colbio <- c( "AG" = "#ff6500", "EA" = "deepskyblue1", "NP" = "green4")

socstock_barsum$location_f =factor(socstock_barsum$location, levels=c('Tribune', 'Hays', 'Manhattan'))

 ggplot(data=socstock_barsum, aes(x=treatment, y=sum, fill = treatment)) + 
  geom_bar(position=position_dodge(), stat="identity", colour = "black") +
  scale_y_continuous(limits=c(0,250)) +  
  facet_wrap(facets=vars(location_f), strip.position="bottom")  +
  xlab("")+
  scale_fill_manual(values = sfcol) + 
     ggtitle("Total SOC Stocks")+ 
  theme_James() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),
        strip.background=element_rect(size=0.5, colour = "black"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.position= c(0.1,0.86)) +
  ylab(expression(Total ~SOC ~Stocks ~(Mg ~C ~ha^{-1} ~soil))) +
  guides(fill = guide_legend(override.aes = aes(label="")))+
 geom_label(aes(label=round(sum,0), y = sum+15),
            label.padding = unit(.3,"lines"), fill=NA, show.legend=NA, label.size = NA, font="bold" )

4.7 Extractable Ca Content (g/kg)

#mean, standard deviation, and std error
ca_g <- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(ca_g)), mean=mean(ca_g, na.rm=T), sd = sd(ca_g, na.rm=T), se = sd / sqrt(N)) 


ca_g$location_f =factor(ca_g$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ca_g$depth =factor(ca_g$bdepth, levels=c(0,20,40,60,80,100))


ggplot(ca_g, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,6000) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
    ylab(expression(~Extractable ~Ca ~Content ~(g ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.07,0.2))

ggplot(ca_g, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,6000) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
    ylab(expression(~Extractable ~Ca ~Content ~(g ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.07,0.2))

#+ ggsave("ca1.png",  height=6, width=9)

4.8 Extractable Cu Content (mg/kg)

#mean, standard deviation, and std error
cu<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(cu)), mean=mean(cu, na.rm=T), sd = sd(cu, na.rm=T), se = sd / sqrt(N)) 


cu$location_f =factor(cu$location, levels=c('Tribune', 'Hays', 'Manhattan'))
cu$depth =factor(cu$bdepth, levels=c(0,20,40,60,80,100))


ggplot(cu, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
      ylab(expression(~Extractable ~Cu ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

ggplot(cu, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
      ylab(expression(~Extractable ~Cu ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

#+ ggsave("cu1.png",  height=6, width=9)

4.9 Extractable Mg Content (mg/kg)

#mean, standard deviation, and std error
mg<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(mg)), mean=mean(mg, na.rm=T), sd = sd(mg, na.rm=T), se = sd / sqrt(N)) 


mg$location_f =factor(mg$location, levels=c('Tribune', 'Hays', 'Manhattan'))
mg$depth =factor(mg$bdepth, levels=c(0,20,40,60,80,100))


ggplot(mg, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,1000) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Mg ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.07,0.2))

ggplot(mg, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,1000) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Mg ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.07,0.2))

#+ ggsave("mg1.png",  height=6, width=9)

4.10 Extractable Mn Content (mg/kg)

#mean, standard deviation, and std error
mn<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(mn)), mean=mean(mn, na.rm=T), sd = sd(mn, na.rm=T), se = sd / sqrt(N)) 


mn$location_f =factor(mn$location, levels=c('Tribune', 'Hays', 'Manhattan'))
mn$depth =factor(mn$bdepth, levels=c(0,20,40,60,80,100))


ggplot(mn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Mn ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

ggplot(mn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Mn ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

#+ ggsave("mn1.png",  height=6, width=9)

4.11 Extractable Na Content (mg/kg)

#mean, standard deviation, and std error
na<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(na)), mean=mean(na, na.rm=T), sd = sd(na, na.rm=T), se = sd / sqrt(N)) 


na$location_f =factor(na$location, levels=c('Tribune', 'Hays', 'Manhattan'))
na$depth =factor(na$bdepth, levels=c(0,20,40,60,80,100))


ggplot(na, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,1200) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Na ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.8))

ggplot(na, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,1200) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Na ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.8))

#+ ggsave("na1.png",  height=6, width=9)

4.12 Extractable P Content (mg/kg)

#mean, standard deviation, and std error
p<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(p)), mean=mean(p, na.rm=T), sd = sd(p, na.rm=T), se = sd / sqrt(N)) 


p$location_f =factor(p$location, levels=c('Tribune', 'Hays', 'Manhattan'))
p$depth =factor(p$bdepth, levels=c(0,20,40,60,80,100))


ggplot(p, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,180) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~P ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

ggplot(p, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,180) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~P ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

#+ ggsave("p1.png",  height=6, width=9)

4.13 Extractable K Content (mg/kg)

#mean, standard deviation, and std error
k<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(k)), mean=mean(k, na.rm=T), sd = sd(k, na.rm=T), se = sd / sqrt(N)) 


k$location_f =factor(k$location, levels=c('Tribune', 'Hays', 'Manhattan'))
k$depth =factor(k$bdepth, levels=c(0,20,40,60,80,100))


ggplot(k, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,1200) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~K ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

ggplot(k, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,1200) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~K ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

#+ ggsave("k1.png",  height=6, width=9)

4.14 Extractable Zn Content (mg/kg)

#mean, standard deviation, and std error
zn<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(zn)), mean=mean(zn, na.rm=T), sd = sd(zn, na.rm=T), se = sd / sqrt(N)) 


zn$location_f =factor(zn$location, levels=c('Tribune', 'Hays', 'Manhattan'))
zn$depth =factor(zn$bdepth, levels=c(0,20,40,60,80,100))


ggplot(zn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Zn ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

ggplot(zn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Zn ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.9,0.2))

#+ ggsave("zn1.png",  height=6, width=9)

4.15 Extractable Fe Content (mg/kg)

#mean, standard deviation, and std error
fe<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(fe)), mean=mean(fe, na.rm=T), sd = sd(fe, na.rm=T), se = sd / sqrt(N)) 


fe$location_f =factor(fe$location, levels=c('Tribune', 'Hays', 'Manhattan'))
fe$depth =factor(fe$bdepth, levels=c(0,20,40,60,80,100))


ggplot(fe, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,100) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Fe ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.2,0.2))

ggplot(fe, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,100) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
        ylab(expression(~Extractable ~Fe ~Content ~(mg ~kg^{-1} ~soil))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.2,0.2))

#+ ggsave("fe1.png",  height=6, width=9)

4.16 pH

#mean, standard deviation, and std error
p_h<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(p_h)), mean=mean(p_h, na.rm=T), sd = sd(p_h, na.rm=T), se = sd / sqrt(N)) 


p_h$location_f =factor(p_h$location, levels=c('Tribune', 'Hays', 'Manhattan'))
p_h$depth =factor(p_h$bdepth, levels=c(0,20,40,60,80,100))


ggplot(p_h, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)",
       y="pH") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.1,0.2))

ggplot(p_h, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)",
       y="pH") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.1,0.2))

#+ ggsave("ph1.png",  height=6, width=9)

5 Physical graphs

5.1 Clay

#mean, standard deviation, and std error
clay<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(clay)), mean=mean(clay, na.rm=T), sd = sd(clay, na.rm=T), se = sd / sqrt(N)) 


clay$location_f =factor(clay$location, levels=c('Tribune', 'Hays', 'Manhattan'))
clay$depth =factor(clay$bdepth, levels=c(0,20,40,60,80,100))


ggplot(clay, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Clay ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.8))

ggplot(clay, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)") +
  ylab(expression(Clay ~("%"))) +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.8))

#+ ggsave("clay1.png",  height=6, width=9)


ggplot(clay, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,40) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)")+
  ylab(expression(Clay ~("%"))) +
scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.55,0.19))

5.2 20 minute MWD

#mean, standard deviation, and std error
x20mwd<- ddply(agg, c("location","treatment","bdepth"), mutate, N=sum(!is.na(x20mwd)), mean=mean(x20mwd, na.rm=T), sd = sd(x20mwd, na.rm=T), se = sd / sqrt(N)) 


x20mwd$location_f =factor(x20mwd$location, levels=c('Tribune', 'Hays', 'Manhattan'))
x20mwd$depth =factor(x20mwd$bdepth, levels=c(0,20,40,60,80,100))


twe<- ggplot(x20mwd, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)",
       y="20 Minute Mean Weight Diameter (mm)") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.93,0.2))




twe<- ggplot(x20mwd, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)",
       y="20 Minute Mean Weight Diameter (mm)") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.93,0.2))
#+ ggsave("20mwd1.png",  height=6, width=9)

ggplot(x20mwd, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)",
       y="Mean Weight Diameter (mm)") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.55,0.19))

5.3 Bulk Density

#mean, standard deviation, and std error
bd<- ddply(agg, c("location","treatment","bdepth"), mutate, N=sum(!is.na(bd)), mean=mean(bd, na.rm=T), sd = sd(bd, na.rm=T), se = sd / sqrt(N)) 

bd$location_f =factor(bd$location, levels=c('Tribune', 'Hays', 'Manhattan'))
bd$depth =factor(bd$bdepth, levels=c(0,20,40,60,80,100))


ggplot(bd, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)")+
ylab(expression(Bulk ~Density ~(g ~cm^{-3} ))) +  
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.1,0.2))

ggplot(bd, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,4) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)")+
ylab(expression(Bulk ~Density ~(g ~cm^{-3} ))) +  
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.1,0.2))

#+ ggsave("bd1.png",  height=6, width=9)

6 Ratio Graphs

6.1 Soil Organic Carbon to Nitrogen

#mean, standard deviation, and std error
cn<- ddply(nut, c("location","treatment","bdepth"), mutate, N=sum(!is.na(cn)), mean=mean(cn, na.rm=T), sd = sd(cn, na.rm=T), se = sd / sqrt(N)) 


cn$location_f =factor(cn$location, levels=c('Tribune', 'Hays', 'Manhattan'))
cn$depth =factor(cn$bdepth, levels=c(0,20,40,60,80,100))


ggplot(cn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="Soil Organic Carbon to Total Nitrogen Ratio") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

ggplot(cn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,0.5) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="C to N ratio") +
  scale_color_manual(values=fcol) +
  theme_James() +
  theme(legend.position = c(0.92,0.5))

#+ ggsave("cn1.png",  height=6, width=9)

6.2 Protein to Total Nitrogen

#mean, standard deviation, and std error
ptn<- ddply(ratios, c("location","treatment","bdepth"), mutate, N=sum(!is.na(ptn)), mean=mean(ptn, na.rm=T), sd = sd(ptn, na.rm=T), se = sd / sqrt(N)) 


ptn$location_f =factor(ptn$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ptn$depth =factor(ptn$bdepth, levels=c(0,20,40,60,80,100))


ggplot(ptn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="Protein to Total Nitrogen ratio")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(ptn, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="Protein to Total Nitrogen ratio")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+ ggsave("ptn1.png",  height=6, width=9)

6.3 Soil Respiration to Soil Organic Carbon

#mean, standard deviation, and std error
rec<- ddply(ratios, c("location","treatment","bdepth"), mutate, N=sum(!is.na(rec)), mean=mean(rec, na.rm=T), sd = sd(rec, na.rm=T), se = sd / sqrt(N)) 


rec$location_f =factor(rec$location, levels=c('Tribune', 'Hays', 'Manhattan'))
rec$depth =factor(rec$bdepth, levels=c(0,20,40,60,80,100))


ggplot(rec, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="Respiration to Soil Organic Carbon ratio")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

ggplot(rec, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="Respiration to Soil Organic Carbon ratio")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.9,0.3))

#+ ggsave("rec1.png",  height=6, width=9)

6.4 Glucosidase to Glucosaminidase ratio

#mean, standard deviation, and std error
lncne2<- ddply(epr2, c("location","treatment","bdepth"), mutate, N=sum(!is.na(lngluc_glucosam   )), mean=mean(lngluc_glucosam-1.1   , na.rm=T), sd = sd(lngluc_glucosam   , na.rm=T), se = sd / sqrt(N)) 


lncne2$location_f =factor(lncne2$location, levels=c('Tribune', 'Hays', 'Manhattan'))
lncne2$depth =factor(lncne2$bdepth, levels=c(0,20,40,60,80,100))


ggplot(lncne2, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(NAG)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

ggplot(lncne2, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(NAG)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

#+ ggsave("incne1.png",  height=6, width=9)

6.5 Glucosidase to Acid Phosphotase ratio

#mean, standard deviation, and std error
lncape<- ddply(epr2, c("location","treatment","bdepth"), mutate, N=sum(!is.na(lngluc_acidp)), mean=mean(lngluc_acidp , na.rm=T), sd = sd(lngluc_acidp , na.rm=T), se = sd / sqrt(N)) 


lncape$location_f =factor(lncape$location, levels=c('Tribune', 'Hays', 'Manhattan'))
lncape$depth =factor(lncape$bdepth, levels=c(0,20,40,60,80,100))


ggplot(lncape, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(AP)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

ggplot(lncape, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(AP)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

#+ ggsave("incape1.png",  height=6, width=9)

6.6 Glucosidase to Alkaline Phosphatase ratio

#mean, standard deviation, and std error
lncakpe<- ddply(epr2, c("location","treatment","bdepth"), mutate, N=sum(!is.na(lngluc_alkp )), mean=mean(lngluc_alkp   , na.rm=T), sd = sd(lngluc_alkp   , na.rm=T), se = sd / sqrt(N)) 


lncakpe$location_f =factor(lncakpe$location, levels=c('Tribune', 'Hays', 'Manhattan'))
lncakpe$depth =factor(lncakpe$bdepth, levels=c(0,20,40,60,80,100))


ggplot(lncakpe, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(ALP)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

ggplot(lncakpe, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(ALP)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

#+ ggsave("incakpe1.png",  height=6, width=9)

6.7 Glucosidase to Arylsulfatase ratio

#mean, standard deviation, and std error
lncarle<- ddply(epr2, c("location","treatment","bdepth"), mutate, N=sum(!is.na(lngluc_ary)), mean=mean(lngluc_ary-1.1   , na.rm=T), sd = sd(lngluc_ary   , na.rm=T), se = sd / sqrt(N)) 


lncarle$location_f =factor(lncarle$location, levels=c('Tribune', 'Hays', 'Manhattan'))
lncarle$depth =factor(lncarle$bdepth, levels=c(0,20,40,60,80,100))


ggplot(lncarle, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(ARY)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

ggplot(lncarle, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(ARY)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

#+ ggsave("incarle1.png",  height=6, width=9)

6.8 Glucosidase to Phosphodiesterase ratio

#mean, standard deviation, and std error
lncppe<- ddply(epr2, c("location","treatment","bdepth"), mutate, N=sum(!is.na(lngluc_pho   )), mean=mean(lngluc_pho , na.rm=T), sd = sd(lngluc_pho   , na.rm=T), se = sd / sqrt(N)) 


lncppe$location_f =factor(lncppe$location, levels=c('Tribune', 'Hays', 'Manhattan'))
lncppe$depth =factor(lncppe$bdepth, levels=c(0,20,40,60,80,100))


ggplot(lncppe, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(PHO)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

ggplot(lncppe, aes(y=mean, x=bdepth, color=treatment)) +
  geom_line(position=pd, size=1.3, aes(color=treatment),se=F) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), color="black",width=3,size=0.9,position = position_dodge(.1)) + 
  geom_point( pch=21, size=3, color="black") +
  geom_point(aes(color=treatment), size=2.5)+
  ylim(0,70) +
  scale_x_reverse() +
  facet_wrap(facets=vars(location_f), strip.position="bottom") +
  coord_flip() +
  scale_x_continuous(trans = "reverse", breaks = seq(0,100,10)) +
  scale_y_continuous(position = "right")+
  labs(x="Depth (cm)", y="ln(bG) : ln(PHO)")+
  scale_color_manual(values=fcol) +
  theme_James()+
  theme(legend.position = c(0.93,0.3))+
  geom_hline(yintercept=1,linetype="dotted")

#+ ggsave("incppe1.png",  height=6, width=9)