#load packages
library(ggplot2)
library(dplyr)

#First set directory to where your files are located

#Aboveground Productivity##############################################
#read csv
bio.productivity<-read.csv('BIO II biomass.csv')
str(bio.productivity)#check structure of dataframe

#Cleaning up data
#Use subset to exclude some '0 species' plots, create new df
bio.productivity <- subset(bio.productivity, NumSp > 0) 

sum(is.na(bio.productivity$Biomass..g.m2.)) #count number of rows with NA -- there are 2
bio.productivity<- bio.productivity[!is.na(bio.productivity$Biomass..g.m2.), ] #remove the rows with NA, create new df
sum(is.na(bio.productivity$Biomass..g.m2.)) #now there are no NAs, 


#Check the Years of data and the treatments and numbers of plots
length(unique(bio.productivity$Year))#22 years of data
unique(bio.productivity$Year)#lists the years, final year is 2017
length(unique(bio.productivity$NumSp))#confirming there are 5 species treatments
bio.productivity$Plot<-as.factor(bio.productivity$Plot)#make Plot a factor
length(unique(bio.productivity$Plot))#168 plots
unique(bio.productivity$NumSp)#unique treatments are 1,2,4,8,16 spp


bio.productivity$Plot<-as.factor(bio.productivity$Plot)#make Plot a factor

#Example 1 -- Subset one year only, create new df
first.year <- subset(bio.productivity, Year == 1996) 
str(first.year)
first.year$Plot<- as.factor(first.year$Plot)

#plot graph of biomass vs. species number
ggplot(first.year, aes(x=NumSp, y=Biomass..g.m2.))+
  geom_point() + 
  geom_smooth(se=FALSE)+
  ylab(Productivity~(g/m2))+
  xlab(Number~of~species)+
  theme_bw()


#Example 2 -- Use dplyr to summarize data  
bio.prod.summary<- 
  bio.productivity %>% 
  group_by(NumSp, Plot, Year)%>%
  summarise(mean=mean(Biomass..g.m2.), max=max(Biomass..g.m2.), min=min(Biomass..g.m2.))

str(bio.prod.summary)
bio.prod.summary<-as.data.frame(bio.prod.summary)
  

#plot mean response of biomass for each treatment over time
ggplot(bio.prod.summary, aes(x=Year, y=mean, group=NumSp, col=NumSp))+
  geom_point()+
  geom_smooth(se=FALSE)+
  ylab(Productivity~(g/m2)) +
  theme_bw()

#you can also use geom_ribbon to plot the max and min of each



#Soil Carbon#############################################################################
bio.carbon<-read.csv('BIO II soil carbon.csv')
str(bio.carbon)#check structure of df
#Use subset to exclude some '0 species' plots, create new df
bio.carbon <- subset(bio.carbon, NumSp > 0) 

sum(is.na(bio.carbon$pTotalCs)) #count number of rows with NA -- there are 7
bio.carbon<- bio.carbon[!is.na(bio.carbon$pTotalCs), ] #remove the rows with NA, create new df
sum(is.na(bio.carbon$pTotalCs)) #now there are no NAs
bio.carbon$Plot<-as.factor(bio.3$Plot)#make Plot a factor

#Check the Years of data
length(unique(bio.carbon$Year))#7 years of data
unique(bio.carbon$Year)#lists the years, final year is 2015
length(unique(bio.carbon$NumSp))#confirming there are 5 species treatments



#Example 1 -- Subset last year only, create new df
last.year.carbon <- subset(bio.carbon, Year == 2015) 
str(last.year.carbon)
last.year.carbon$Plot<- as.factor(last.year.carbon$Plot)

#plot graph of biomass vs. species number
ggplot(last.year.carbon, aes(x=NumSp, y=pTotalCs))+
  geom_point() + 
  geom_smooth(se=FALSE)

#Example 2 -- Use dplyr to summarize data  
bio.carbon.summary<- 
  bio.carbon %>% 
  group_by(NumSp, Plot, Year)%>%
  summarise(mean=mean(pTotalCs), max=max(pTotalCs), min=min(pTotalCs))

str(bio.carbon.summary)
bio.carbon.summary<-as.data.frame(bio.carbon.summary)


#plot mean response of biomass for each treatment over time
ggplot(bio.carbon.summary, aes(x=Year, y=mean, group=NumSp, col=NumSp))+
  geom_point()+
  geom_smooth(se=FALSE)

#you can also use geom_ribbon to plot the max and min of each