exec/SurfaceArea_roots.R

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

#######calculate root surface area
library(ggplot2)
library(gridExtra)
library(cowplot)
library(lme4)
library(lmerTest)
library(lubridate)
library(emmeans)
setwd("D:/R_Workspace/SDEF/SDEF.analysis/data")
load("SDEF_roots_pre_aggregation.rda")

data_roots <- SDEF_roots_pre_aggregation
data_roots$SA <- ((pi*((data_roots$diameter)^2)*data_roots$length))
data_roots$month <- month(as.Date(data_roots$date, format = "%m/%d/%Y"))

mean_SA <- aggregate(
  x = data_roots[["SA"]], by = list(
    data_roots[["site"]],
    data_roots[["month"]],
    data_roots[["depth_bin"]],
    data_roots[["tube"]]

  ), FUN = function(x) { mean(x, na.rm = TRUE) })
colnames(mean_SA) <- c("site", "month", "depth_bin", "mean_SA")


# Find s.d. of surface area, binned by site/date/depth_bin

sd_SA <- aggregate(
  x = data_roots[["SA"]], by = list(
    data_roots[["site"]],
    data_roots[["month"]],
    data_roots[["depth_bin"]]
  ), FUN = function(x) { sd(x, na.rm = TRUE) })
colnames(sd_SA) <- c("site", "month", "depth_bin", "sd_SA")
# Find sample size ('n') of root surface area, binned by site/date/depth_bin
n_SA <- aggregate(
  x = data_roots[["SA"]], by = list(
    data_roots[["site"]],
    data_roots[["month"]],
    data_roots[["depth_bin"]]
  ), FUN = function(x) { length(x) })
colnames(n_SA) <- c("site", "month", "depth_bin", "n_SA")

root_SA <- cbind(mean_SA, sd_SA$sd_SA, n_SA$n_SA)
colnames(root_SA) <- c("site", "month", "depth_bin", "mean_SA", "sd_SA", "n_SA")

root_SA$se_SA <- ((root_SA$sd_SA)/sqrt(root_SA$n_SA))
root_SA$minse <- abs((root_SA$mean_SA -root_SA$se_SA))
root_SA$maxse <- (root_SA$mean_SA +root_SA$se_SA)
rm(mean_SA, sd_SA, n_SA, SDEF_roots_pre_aggregation)
# Section making months
Dec <- subset(root_SA, root_SA$month == 12)
Dec$Date <- rep(as.Date("12/01/2015", format = "%m/%d/%Y"), times = 7)
Jan <- subset(root_SA, root_SA$month == 1)
Jan$Date <- rep(as.Date("01/01/2016", format = "%m/%d/%Y"), times = 23)
Feb <- subset(root_SA, root_SA$month == 2)
Feb$Date <- rep(as.Date("02/01/2016", format = "%m/%d/%Y"), times = 15)
Mar <- subset(root_SA, root_SA$month == 3)
Mar$Date <- rep(as.Date("03/01/2016", format = "%m/%d/%Y"), times = 8)
May <- subset(root_SA, root_SA$month == 5)
May$Date <- rep(as.Date("05/01/2016", format = "%m/%d/%Y"), times = 16)
root_SA <- rbind(Dec, Jan, Feb, Mar, May)
rm(Dec, Jan, Feb, Mar, May)

### add in native v invasive
ad <- subset(root_SA, root_SA$site == "adenostoma")
ad$veg_type <- rep("Native", times = 37)
g <- subset(root_SA, root_SA$site == "grass")
g$veg_type <- rep("Invasive", times = 36)
alldates <- rbind(ad, g)
rm(ad, g)

### Subset by depth
zt10 <- subset(alldates, alldates$depth_bin == 1)
tent20 <- subset(alldates, alldates$depth_bin == 2)
twent30 <- subset(alldates, alldates$depth_bin == 3)
thirt40 <- subset(alldates, alldates$depth_bin == 4)

#######COLOR###########
##########Make the same figure in COLOR
ad10 <- subset(zt10, zt10$site == "adenostoma")
g10 <- subset(zt10, zt10$site == "grass")

p1 <- ggplot(data = ad10, aes(y = mean_SA, x = Date)) +
  scale_x_date(labels = date_format("%b-%y"), date_breaks = "2 month")+
  geom_point(data = ad10, aes(colour = "Native", size = 1), show.legend = F)+
  geom_errorbar(data = ad10, aes(x = Date, ymin=minse, ymax=maxse, width = 10,  colour = "Native") ) +
  ylab(expression(paste("Root biomass per soil volume (g   ", dm^-3, ")"))) +
  xlab("") +
  ggtitle("0-10 cm") +
  geom_point(data = g10,position = position_nudge(x = 10),  aes(colour = "Invasive", size = 1), show.legend = F)+
  geom_errorbar(data = g10, position = position_nudge(x = 10), aes(x = Date, ymin=minse, ymax=maxse, width = 15, colour = "Invasive") ) +
  scale_colour_manual("", breaks = c("Native", "Invasive"), values=c("#E58601", "#E2D200"), guide = F) +
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5),
        text = element_text(size=18))
p1

# 10 -20
ad20 <- subset(tent20, tent20$site == "adenostoma")
g20 <- subset(tent20, tent20$site == "grass")

p2 <- ggplot(data = ad20, aes(y = mean_biodensity_dm, x = Date)) +
  scale_x_date(labels = date_format("%b-%y"), date_breaks = "2 month")+
  geom_point(data = ad20, aes(colour = "Native", size = 1), show.legend = F)+
  geom_errorbar(data = ad20, aes(x = Date, ymin=minse, ymax=maxse, width = 10,  colour = "Native") ) +
  xlab("") +
  ggtitle("10-20 cm") +
  scale_y_continuous(limits = c(0, .75))+
  geom_point(data = g20,position = position_nudge(x = 10),  aes(colour = "Invasive", size = 1), show.legend = F)+
  geom_errorbar(data = g20, position = position_nudge(x = 10), aes(x = Date, ymin=minse, ymax=maxse, width = 15, colour = "Invasive") ) +
  scale_colour_manual("", breaks = c("Native", "Invasive"), values=c("#E58601", "#E2D200"), guide = F) +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        plot.title = element_text(hjust = 0.5),
        text = element_text(size=18))
p2

### 20-30
ad30 <- subset(twent30, twent30$site == "adenostoma")
g30 <- subset(twent30, twent30$site == "grass")

p3 <- ggplot(data = ad30, aes(y = mean_biodensity_dm, x = Date)) +
  scale_x_date(labels = date_format("%b-%y"), date_breaks = "2 month")+
  geom_point(data = ad30, aes(colour = "Native", size = 1), show.legend = F)+
  geom_errorbar(data = ad30, aes(x = Date, ymin=minse, ymax=maxse, width = 10,  colour = "Native") ) +
  xlab("") +
  ggtitle("20-30 cm") +
  scale_y_continuous(limits = c(0, .75))+
  geom_point(data = g30,position = position_nudge(x = 10),  aes(colour = "Invasive", size = 1), show.legend = F)+
  geom_errorbar(data = g30, position = position_nudge(x = 10), aes(x = Date, ymin=minse, ymax=maxse, width = 15, colour = "Invasive") ) +
  scale_colour_manual("", breaks = c("Native", "Invasive"), values=c("#E58601", "#E2D200"), guide = F) +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        plot.title = element_text(hjust = 0.5),
        text = element_text(size=18))
p3

# 30-40
ad40 <- subset(thirt40, thirt40$site == "adenostoma")
g40 <- subset(thirt40, thirt40$site == "grass")

p4 <- ggplot(data = ad40, aes(y = mean_biodensity_dm, x = Date)) +
  scale_x_date(labels = date_format("%b-%y"), date_breaks = "2 month")+
  geom_point(data = ad40, aes(colour = "Native", size = 1), show.legend = F)+
  geom_errorbar(data = ad40, aes(x = Date, ymin=minse, ymax=maxse, width = 10,  colour = "Native") ) +
  xlab("") +
  ggtitle("30-40 cm") +
  scale_y_continuous(limits = c(0, .75))+
  geom_point(data = g40,position = position_nudge(x = 10),  aes(colour = "Invasive", size = 1), show.legend = F)+
  geom_errorbar(data = g40, position = position_nudge(x = 10), aes(x = Date, ymin=minse, ymax=maxse, width = 15, colour = "Invasive") ) +
  scale_colour_manual("", breaks = c("Native", "Invasive"), values=c("#E58601", "#E2D200"), guide = F) +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        plot.title = element_text(hjust = 0.5),
        text = element_text(size=18))
p4

## Total

adtot <- subset(Totals, Totals$veg_type == "Native")
gtot <- subset(Totals, Totals$veg_type == "Invasive")

p5 <- ggplot(data = adtot, aes(y = biomass, x = Date)) +
  scale_x_date(labels = date_format("%b-%y"), date_breaks = "2 month")+
  geom_point(data = adtot, aes(colour = "Native", size = 1), show.legend = F)+
  geom_errorbar(data = adtot, aes(x = Date, ymin=minse, ymax=maxse, width = 10,  colour = "Native") ) +
  xlab("") +
  ggtitle("Total") +
  scale_y_continuous(limits = c(0, .75))+
  geom_point(data = gtot,position = position_nudge(x = 10),  aes(colour = "Invasive", size = 1), show.legend = F)+
  geom_errorbar(data = gtot, position = position_nudge(x = 10), aes(x = Date, ymin=minse, ymax=maxse, width = 15, colour = "Invasive") ) +
  scale_colour_manual("", breaks = c("Native", "Invasive"), values=c("#E58601", "#E2D200"), guide =F) +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        plot.title = element_text(hjust = 0.5),
        text = element_text(size=18),
        legend.position = "bottom")
p5

grid.arrange(p1, p2, p3, p4, p5, ncol=5)