exec/Isotpe_Figures.R

In bmcnellis/SDEF.analysis: Data Analysis for 2016 San Dimas Experimental Forest Project

#
library(ggplot2)
setwd("D:/R_Workspace/SDEF/SDEF_Isotope_Data")

I <- read.csv(file = "13C_Data_Submission_2_17.csv")
head(I)
Iad <- subset(I, I$Species == "ADEFAS")
Ig <- subset(I, I$Species == "EHRCAL")
plot(Iad$D13C, type = "p", col = 4)
points(Ig$D13C, type = "p", col = 51)



iso <- ggplot(I, aes(x = Month , y = D13C)) + geom_boxplot(aes(fill = Species))
iso <- iso + theme_bw(base_family = "sans", base_size = 13)+ ylab(expression(delta^{13}*"C (‰)"))+
  scale_fill_manual(values=c("White", "Grey"))+ xlab("Time (Month)")
iso


I.lm <- lm(I$D13C ~ I$Season + I$Species)
shapiro.test(I.lm$resid)
I.anova <- anova(I.lm)
I.anova


iso <- ggplot(I, aes(x = Month , y = D15N)) + geom_boxplot(aes(fill = Species))
iso <- iso +theme_bw(base_family = "sans", base_size = 13)+ ylab(expression(delta^{15}*"N (‰)")) +
scale_fill_manual(values=c("White", "Grey"))+ xlab("Time (Month)")
iso

I.lm <- lm(I$D15N ~ I$Season + I$Species)
step(I.lm)
shapiro.test(I.lm$resid)
I.anova <- anova(I.lm)
I.anova


#### D18 O figure
se <- function(x) sqrt(var(x)/length(x))
O <- read.csv(file = "18_Oisotope.csv")
head(O)
Odry <- subset(O, O$Season == "Dry")
dry <- mean(Odry[,2])
dryse <- se(Odry[,2])
Ostem <- subset(O, O$Type == "Stem")
Owetstem <- subset(Ostem, Ostem$Season == "Wet")
wetstem <- mean(Owetstem[,2])
wetstemse <- se(Owetstem[,2])
Odrystem <- subset(Ostem, Ostem$Season == "Dry")
Owell <- subset(O, O$Season == "Well")
well <- mean(Owell[,2])
wellse <- se(Owell[,2])
rain <- subset(O, O$Season =="Rain")
rainme <- mean(rain[,2])
rainse <- se(rain[,2])


Mean18O <- as.data.frame(c(rainme, well, wetstem, dry))
colnames(Mean18O) <- "Mean"
SE <- as.data.frame(c(rainse, wellse, wetstemse, dryse))
colnames(SE)<-"SE"
O18 <- cbind(Mean18O, SE)
O18$Source <- c("Rainwater", "Well Water", "Stem (Feb)", "Stem (Jun)" )

#O <- O[c(1:20),c(1:4)]


o18 <- ggplot(O18, aes(x = Source, y = Mean)) + geom_bar(stat = "identity", fill = "grey") +
  geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width=0.2)
o18 <- o18+ theme_bw(base_family = "sans", base_size = 13)+ ylab(expression(delta^{18}*"O (‰)"))
o18

# 5/18?2018 try to make this figure into a box plot.
index <- c("Rain", "Well","Wet", "Dry")
values <- c( "Rainwater",  "Well Water", "Stem (Feb)", "Stem (Jun)")
O$Source <- values[match(O$Season, index)]

O$Source_1 <- factor(O$Source, c("Rainwater",  "Well Water", "Stem (Feb)", "Stem (Jun)"))

o18 <- ggplot(O, aes(y = X18O, x = Source_1)) + geom_boxplot(aes(fill = Source_1)) + scale_fill_manual(values = c("gray82", "black", "white", "gray40"))
o18 <- o18 + theme_bw(base_family = "sans", base_size = 13)+ ylab(expression(delta^{18}*"O (‰)")) + guides(fill = FALSE)+
  xlab("Source")  + annotate("text", x = 4, y = 5.5, label = "***")
o18

o18mod <- lm(X18O ~ Source_1, data = O)
shapiro.test(o18mod$residuals)
aov <- aov(o18mod)
TukeyHSD(x = aov, "Source_1")