# MAXIMUM CARBON ASSIMILATION MODEL for UNDERSTORY SHRUBS ####

# Written by Erin O'Connell 2020 (oconn877@d.umn.edu)
# Citation: O'Connell, E. and Savage, J. (2020) Extended leaf phenology has limited benefits for invasive species growing at northern latitudes. Biological Invasions. In press.

# PROGRAM DESCRIPTION:
#   This program models maximum seasonal carbon assimilation for 8 species (6 plants/species) growing in a 50-year-old mixed forest at Bagley Nature Area in Duluth, MN.
#   Species included in the model are Berberis thunbergii (BETH), Corylus cornuta (COCO), Cornus sericea (COSE), Frangula alnus (FRAL), Lonicera x. bella (LOBE), Rhamnus cathartica (RHCA), Rubus ideaus (RUID), and Viburnum lentago (VILE).
#   Maximum carbon assimilation is modeled from understory light availability on sunny days, carbon dioxide assimilation rates, and leaf area over time.
#   See publication for more details.

# CITATIONS:
#   Lines 532-535 and 598-601 are from: Heberling JM, Fridley JD (2013) Resource-use strategies of native and invasive plants in Eastern North American forests. New Phytol 200(2):523-533. doi: 10.1111/nph.12388
#   and based on: Marshall B, Biscoe PV (1980) A model for C3 leaves describing the dependence of net photosynthesis on irradiance. J Exp Bot 31(1):29-39. doi: 10.1093/jxb/31.1.29

# NOTES:
#   If running this program in RStudio, be sure the "Plots" window (lower right frame) is expanded to avoid "figure margins too large" error.
#   In RStudio, use alt+o to collapse code and alt+shift+o to expand code.

# ================== SET - UP =========================================================.====
#(1) LOAD PACKAGES --------------------------------------------------------------------.----
#library(readr)
library(ggplot2)
library(broom) #for the "tidy" function to return SS read-outs
library(dplyr) #for pivot tables and other dataframe organization
library(tidyr) #for gathering columns
library(insol) # for calculating daylength
library(DescTools) #for calculating area under the curve and converting to seconds
library(directlabels) #for labelling lines in ggplot
library(gridExtra) #for arranging plots
library(ggsignif) #for adding sginificance in ggplot
library(pracma) # for calculating nthroot
library(car) #for multiple linear regression
library(lmSupport) # for correlation tests

#(2) LOAD DATA ------------------------------------------------------------------------.----
#       Set working directory ####
setwd() #Add your directory here

#       General Plant Traits  ####
#Including: individual plant numbers included in the model, plant height, stem count, leaf count, average leaf area, and total leaf area
Individuals <- read.csv("Carbon_gain_input.csv")

#       Leaf Exapansion ####
Leaf.Area <- read.csv("Spring_leaf_expansion_2017.csv")

#       Fall Leaf Phenology ####
Fall.Pheno <-  read.csv("Fall_leaf_phenology_2017.csv")

#       Photosynthesis ####
Photo_All <- read.csv("Photosynthetic_light_response_curves.csv") # Raw LiCOR 6400 data
Photo_JUN <- Photo_All[Photo_All$Month == "JUN",] # Young Leaves
Photo_AUG <- Photo_All[Photo_All$Month == "AUG",] # Mature LEaves

#       Light Calibration ####
calib_all <- read.csv("Understory_light_monthly.csv")
calib_all$Date <- as.Date(calib_all$Date, "%m/%d/%Y") 
              
# Build data frame for light calibration regression results
Season <- c("Spring", "Summer", "Summer", "Fall")
Light_Calib_Results <- as.data.frame(Season)
Light_Calib_Results$Month <- c("MAY", "JUL", "SEP", "OCT")
Light_Calib_Results$Int <- ""
Light_Calib_Results$Slope <- ""

#       Diurnal Light ####
Diurn.All <- read.csv("Understory_light_daily.csv")
Diurn <- subset.data.frame(Diurn.All, Weather == "sunny") # Subset for sunny days only 
Diurn$Date <- as.Date(Diurn$Date, "%m/%d/%Y")

#       Day Length ####
#R package (insol) citation: Corripio, J. G.: 2003, Vectorial algebra algorithms for calculating terrain parameters from DEMs and the position of the sun for solar radiation modelling in mountainous terrain, International Journal of Geographical Information Science 17(1), 1-23.
jd2017=JD(seq(ISOdate(2018,1,1),ISOdate(2018,12,31),by='day'))
Day.Length <- as.data.frame(daylength(46.7867, -92.1005,jd2017,-6))
Day.Length$Date_J <- seq(1,365,by=1)
Day.Length$Date <- as.Date(Day.Length$Date_J, origin = as.Date("2018-01-01"))

#       Species List ####
# Create dataframe of species in the model
Species <- c("BETH", "COCO", "COSE", "FRAL", "LOBE", "RHCA", "RUID", "VILE")
status <- c()
All_Species <- as.data.frame(Species, status)
for(i in 1:nrow(All_Species)) {
  if(All_Species$Species[i] == "BETH"|All_Species$Species[i] == "LOBE"|All_Species$Species[i] == "FRAL"|All_Species$Species[i] == "RHCA"){
    All_Species$status[i] <- "Invasive"
  }else {
    All_Species$status[i] <- "Native"
  } 
}

#(3) FUNCTIONS --------------------------------------------------------------- --------.----
#       Return values from fitted sigmoidal curves ####
interpolateX <- function(asym, xmid, scal, yval){
  xval <- (-scal)*log((asym/yval)-1)+xmid
  return(xval)
}

interpolateY <- function (asym, xmid, scal, xval){
  yval <- asym/(1+exp((xmid-xval)/scal))
  yval.adj <- yval/asym *100
  return(yval.adj)
}

interpolateX.JMP <- function(asym, xmid, scal, yval){
  xval <- (log((asym/yval)-1)/-scal)+xmid
  return(xval)
}

interpolateY.JMP <- function (asym, xmid, scal, xval){
  yval <- asym/(1+exp((-scal*(xval-xmid))))
  yval.adj <- yval/asym *100
  return(yval.adj)
}

#       Return values from fitted asymptotic curves ####
interpolateX.asymp <- function(asym, lrc, r0, yval){
  xval <- (log(-log((yval-asym)/(r0-asym))))/lrc
  return(xval)
}
# #####
# ================== EXPANDING LEAF AREA ==============================================.====
#(4) SPRING LEAF EXPANSION ------------------------------------------------------------.----
#       Format one measurement per row ####
Leaf.Area.adj <- gather(Leaf.Area, Leaf.No, Area, 
                        Leaf_1, Leaf_2, Leaf_3, Leaf_4, Leaf_5, Leaf_6, Leaf_7, Leaf_8, Leaf_9, Leaf_10, Leaf_11, Leaf_12, Leaf_13, Leaf_14, Leaf_15)

#       View curves ####
par(mfrow=c(3,6))
for (i in Individuals$Plant_ID){
  plot(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID == i] ~
         Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID == i],
       xlab = "Julian Date", ylab = "Area (cm2)", main = i)
}

#       Fit Sigmoidal Curves ####
# Loop for all plants
par(mfrow=c(4,5))
for (i in Individuals$Plant_ID[Individuals$Plant_ID !=13 &
                               Individuals$Plant_ID != 35 &
                               Individuals$Plant_ID != 51 &
                               Individuals$Plant_ID != 58 &
                               Individuals$Plant_ID != 73 &
                               Individuals$Plant_ID != 74 &
                               Individuals$Plant_ID != 85]) {
  #Plot leaf area for each individual
  plot(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID==i] ~
         Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID==i],
       xlab = "Julian Date", ylab = "Leaf Area (cm2)", 
       main = paste(i))
  
  #Fit Sigmoidal Curve
  Area.adj <- Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID == i]
  Date.adj <- Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID == i]
  fit1 <- nls(Area.adj ~ SSlogis(Date.adj, Asym, xmid, scal), 
              data = Leaf.Area.adj)
  lines(seq(0, 200, 1), 
        predict(fit1, newdata = data.frame(Date.adj = seq(0, 200, 1))))
  
  #Return values
  results1 <- data.frame(tidy(fit1))
  Individuals$Asym.EXPAN[Individuals$Plant_ID==i] <- results1$estimate[results1$term=="Asym"]
  Individuals$Xmid.EXPAN[Individuals$Plant_ID==i] <- results1$estimate[results1$term=="xmid"]
  Individuals$Scal.EXPAN[Individuals$Plant_ID==i] <- results1$estimate[results1$term=="scal"]
}

#       Add missing values from JMP ####
Individuals$Asym.EXPAN[Individuals$Plant_ID == 13] <- 1.5556369
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 13] <- 142.12722
Individuals$Scal.EXPAN[Individuals$Plant_ID == 13] <- 0.1044352

Individuals$Asym.EXPAN[Individuals$Plant_ID == 35] <- 0.9998121
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 35] <- 140.43884
Individuals$Scal.EXPAN[Individuals$Plant_ID == 35] <- 0.1616791

Individuals$Asym.EXPAN[Individuals$Plant_ID == 51] <- 35.04168
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 51] <- 148.445
Individuals$Scal.EXPAN[Individuals$Plant_ID == 51] <- 0.1580033

Individuals$Asym.EXPAN[Individuals$Plant_ID == 58] <- 55.750177
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 58] <- 157.93622
Individuals$Scal.EXPAN[Individuals$Plant_ID == 58] <- 1.5121374
Individuals$Scal.EXPAN[Individuals$Plant_ID == 58] <- 0.5 #estimated to more realistically fit data

Individuals$Asym.EXPAN[Individuals$Plant_ID == 73] <- 42.778604
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 73] <- 155.74798
Individuals$Scal.EXPAN[Individuals$Plant_ID == 73] <- 0.5 #estimated to more realistically fit data

Individuals$Asym.EXPAN[Individuals$Plant_ID == 74] <- 38.376347
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 74] <- 151.43508
Individuals$Scal.EXPAN[Individuals$Plant_ID == 74] <- 2.1036005
Individuals$Scal.EXPAN[Individuals$Plant_ID == 74] <- 0.5 #estimated to more realistically fit data

# Convert to numeric
Individuals$Asym.EXPAN <- as.numeric(Individuals$Asym.EXPAN)
Individuals$Xmid.EXPAN <- as.numeric(Individuals$Xmid.EXPAN)
Individuals$Scal.EXPAN <- as.numeric(Individuals$Scal.EXPAN)

#Adjust for outlier
Individuals$Asym.EXPAN[Individuals$Plant_ID == 4] <- Individuals$Ave_leaf_area[Individuals$Plant_ID==4]
Individuals$Xmid.EXPAN[Individuals$Plant_ID == 4] <- ave(Individuals$Xmid.EXPAN[Individuals$Species=="LOBE"])[1]

#       Check fits of JMP fit data and adjusted fits ####
# LOBE 4
par(mfrow = c(1,1))
Asym.a <- Individuals$Asym.EXPAN[Individuals$Plant_ID==4][1] 
xmid.a <- Individuals$Xmid.EXPAN[Individuals$Plant_ID==4][1]
scal.a <- Individuals$Scal.EXPAN[Individuals$Plant_ID==4][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((xmid-x)/scal))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=125, to=190,  xlab = "Julian Date", ylab = "Leaf Area (cm2)",
      ylim=c(0,30))
points(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID==4] ~ Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID==4],
     xlab = "Julian Date", ylab = "Leaf Area (cm2)")

# COSE 73
Asym.a <- Individuals$Asym.EXPAN[Individuals$Plant_ID==73][1] 
xmid.a <- Individuals$Xmid.EXPAN[Individuals$Plant_ID==73][1]
scal.a <- Individuals$Scal.EXPAN[Individuals$Plant_ID==73][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((-scal*(x-xmid))))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=125, to=190,  xlab = "Julian Date", ylab = "Leaf Area (cm2)",
      ylim=c(0,70))
points(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID==73] ~ Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID==73],
     xlab = "Julian Date", ylab = "Leaf Area (cm2)")

# RUID 74
Asym.a <- Individuals$Asym.EXPAN[Individuals$Plant_ID==74][1] 
xmid.a <- Individuals$Xmid.EXPAN[Individuals$Plant_ID==74][1]
scal.a <- Individuals$Scal.EXPAN[Individuals$Plant_ID==74][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((-scal*(x-xmid))))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=125, to=190,  xlab = "Julian Date", ylab = "Leaf Area (cm2)",
      ylim=c(0,60))
points(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID==74] ~ Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID==74],
     xlab = "Julian Date", ylab = "Leaf Area (cm2)")

# COCO 58
Asym.a <- Individuals$Asym.EXPAN[Individuals$Plant_ID==58][1] 
xmid.a <- Individuals$Xmid.EXPAN[Individuals$Plant_ID==58][1]
scal.a <- Individuals$Scal.EXPAN[Individuals$Plant_ID==58][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((-scal*(x-xmid))))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=125, to=190,  xlab = "Julian Date", ylab = "Leaf Area (cm2)",
      ylim=c(0,80))
points(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID==58] ~ Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID==58],
     xlab = "Julian Date", ylab = "Leaf Area (cm2)")

# VILE 51
Asym.a <- Individuals$Asym.EXPAN[Individuals$Plant_ID==51][1] 
xmid.a <- Individuals$Xmid.EXPAN[Individuals$Plant_ID==51][1]
scal.a <- Individuals$Scal.EXPAN[Individuals$Plant_ID==51][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((-scal*(x-xmid))))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=125, to=190,  xlab = "Julian Date", ylab = "Leaf Area (cm2)",
      ylim=c(0,50))
points(Leaf.Area.adj$Area[Leaf.Area.adj$Plant_ID==51] ~ Leaf.Area.adj$Date_J[Leaf.Area.adj$Plant_ID==51],
     xlab = "Julian Date", ylab = "Leaf Area (cm2)")

#(5) CALCULATE DATE of 95% EXPANSION --------------------------------------------------.----
#       Calulation ####
Individuals$Full.Expan <- ""
for(i in 1:nrow(Individuals)){
  asym.I <- Individuals$Asym.EXPAN[i]
  xmid.I <- Individuals$Xmid.EXPAN[i]
  scal.I <- Individuals$Scal.EXPAN[i]
  if(Individuals$Plant_ID[i] == 13 | Individuals$Plant_ID[i] == 35 |
     Individuals$Plant_ID[i] == 51 | Individuals$Plant_ID[i] == 58 |
     Individuals$Plant_ID[i] == 73 | Individuals$Plant_ID[i] == 74){
    Individuals$Full.Expan[i] <- interpolateX.JMP(asym.I, xmid.I, scal.I, asym.I*0.95)
  } else {
    Individuals$Full.Expan[i] <- interpolateX(asym.I, xmid.I, scal.I, asym.I*0.95)
    Individuals$Start.Expan[i] <- interpolateX(asym.I, xmid.I, scal.I, asym.I*0.05)
  }
}
Individuals$Full.Expan <- as.numeric(Individuals$Full.Expan)

#       Figure (one outlier needs to be addressed) ####
ggplot(Individuals,aes(Species, Full.Expan, fill = Status, color = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = "Date of 95% Leaf Expansion", x = "Species") +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.85,.85),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  coord_flip()

#(6) WEEKLY PERCENT EXPANSION ---------------------------------------------------------.----
#       Create dataframes ####
# Create temporary dataframe
Temp <- data.frame(Plant_ID = "", Date_J = "", Area.Perc = "")
# Create Expansion dataframe
Expansion <- data.frame(Plant_ID = "", Date_J = "", Area.Perc = "")
Expansion <- Expansion[-1,]

#       Build observation lists ####
for (i in Individuals$Plant_ID){
  t <- ceiling(min(Individuals$Start.Expan))
  Area.Perc <- 0
  asym.I <- Individuals$Asym.EXPAN[Individuals$Plant_ID == i]
  xmid.I <- Individuals$Xmid.EXPAN[Individuals$Plant_ID == i]
  scal.I <- Individuals$Scal.EXPAN[Individuals$Plant_ID == i]
  
  while(t < ceiling(Individuals$Full.Expan[Individuals$Plant_ID == i])){
    Temp$Plant_ID <- i
    Temp$Date_J <- t
    if(i == 13 | i == 35 | i == 51 | i == 58 | i == 73 | i == 74){
      Area.Perc <- interpolateY.JMP(asym.I, xmid.I, scal.I, t)
    } else {
      Area.Perc <- interpolateY(asym.I, xmid.I, scal.I, t)
    }
    Area.Perc <- as.numeric(Area.Perc)
    Temp$Area.Perc <- Area.Perc
    Expansion <- rbind(Expansion, Temp)
    t <- t + 1
  }
}
Expansion$Area.Perc <- as.numeric(Expansion$Area.Perc)

# ####
# ================== SENESCING LEAF AREA ==============================================.====
#(7) FALL LEAF SENESCENCE CURVES -------------------------------------------------------.----
#       Build Sigmoidal curve for all plants ####
#Avoid zeros
Fall.Pheno$PLS.adj <- Fall.Pheno$Percent_leaf_senescence + 0.001

#Loop for all plants
par(mfrow=c(4,5))
for (i in Individuals$Plant_ID) {
  #Plot each plant
  plot(Fall.Pheno$Percent_leaf_senescence[Fall.Pheno$Plant_ID==i] ~ 
         Fall.Pheno$Date_J[Fall.Pheno$Plant_ID==i],
       xlab = "Julian Date", ylab = "% Clr Chge of Ttl", 
       main = i)
  
  #Fit sigmoid curve to each plant
  PLS.adj2 <- Fall.Pheno$PLS.adj[Fall.Pheno$Plant_ID==i]
  Date_J.adj <- Fall.Pheno$Date_J[Fall.Pheno$Plant_ID==i]
  fit2 <- nls(PLS.adj2 ~ SSlogis(Date_J.adj, Asym, xmid, scal), 
              data = Fall.Pheno)
  lines(seq(220, 350, 1), 
        predict(fit2, newdata = data.frame(Date_J.adj = seq(220, 350, 1))))
  
  #Return 50% leaf drop (curve midpoint)
  results1 <- data.frame(tidy(fit2))
  asym2 <- results1$estimate[results1$term=="Asym"]
  xmid2 <- results1$estimate[results1$term=="xmid"]
  scal2 <- results1$estimate[results1$term=="scal"]
  
  Individuals$Asym.SENES[Individuals$Plant_ID==i] <- results1$estimate[results1$term=="Asym"]
  Individuals$Xmid.SENES[Individuals$Plant_ID==i] <- results1$estimate[results1$term=="xmid"]
  Individuals$Scal.SENES[Individuals$Plant_ID==i] <- results1$estimate[results1$term=="scal"]
}

#Convert to numbers
Individuals$Asym.SENES <- as.numeric(Individuals$Asym.SENES)
Individuals$Xmid.SENES <- as.numeric(Individuals$Xmid.SENES)
Individuals$Scal.SENES <- as.numeric(Individuals$Scal.SENES)

#Correct for values over 100%
for(i in Individuals$Asym.SENES){
  if(i > 100){
    Individuals$Asym.SENES[Individuals$Asym.SENES==i] <- 100
  } else if (i <= 100){
    Individuals$Asym.SENES[Individuals$Asym.SENES==i] <- i
  }
}

#       Adjust for outliers due to bad fits ####
# LOBE 4
Individuals$Scal.SENES[Individuals$Plant_ID == 16] <- 15

par(mfrow = c(1,1))
Asym.a <- Individuals$Asym.SENES[Individuals$Plant_ID==16][1] 
xmid.a <- Individuals$Xmid.SENES[Individuals$Plant_ID==16][1]
scal.a <- Individuals$Scal.SENES[Individuals$Plant_ID==16][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((xmid-x)/scal))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=240, to=350,  xlab = "Julian Date", ylab = "% Clr Chge of Ttl")
points(Fall.Pheno$Percent_leaf_senescence[Fall.Pheno$Plant_ID==16] ~ Fall.Pheno$Date_J[Fall.Pheno$Plant_ID==16],
     xlab = "Julian Date", ylab = "% Clr Chge of Ttl")

#RUID 84
Individuals$Xmid.SENES[Individuals$Plant_ID == 84] <- mean(Individuals$Xmid.SENES[Individuals$Species == "RUID"])
Individuals$Scal.SENES[Individuals$Plant_ID == 84] <- 15

Asym.a <- Individuals$Asym.SENES[Individuals$Plant_ID==84][1] 
xmid.a <- Individuals$Xmid.SENES[Individuals$Plant_ID==84][1]
scal.a <- Individuals$Scal.SENES[Individuals$Plant_ID==84][1]
eq = function(x, Asym, xmid, scal){Asym/(1+exp((xmid-x)/scal))}
curve(eq(x, Asym.a, xmid.a, scal.a), from=240, to=350,  xlab = "Julian Date", ylab = "% Clr Chge of Ttl")
points(Fall.Pheno$Percent_leaf_senescence[Fall.Pheno$Plant_ID==84] ~ Fall.Pheno$Date_J[Fall.Pheno$Plant_ID==84],
     xlab = "Julian Date", ylab = "% Clr Chge of Ttl")

#(8) CALCULATE DATE of 5% SENESCENCE (START SENES) ------------------------------------.----
Individuals$Start.Senes <- ""
for(i in 1:nrow(Individuals)){
    asym.II <- Individuals$Asym.SENES[i]
    xmid.II <- Individuals$Xmid.SENES[i]
    scal.II <- Individuals$Scal.SENES[i]
    Individuals$Start.Senes[i] <- interpolateX(asym.II, xmid.II, scal.II, asym.II*0.05)
}
Individuals$Start.Senes <- as.numeric(Individuals$Start.Senes)

#       Figure ####  
ggplot(Individuals,aes(Species, Start.Senes, fill = Status, color=Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = "Date of 5% Leaf Senescence", x = "Species") +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.15,.85),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  coord_flip()

#(9) CALCULATE DATE of 95% SENESCENCE (FIN SENES)--------------------------------------.----
Individuals$Fin.Senes <- ""
for(i in 1:nrow(Individuals)){
    asym.II <- Individuals$Asym.SENES[i]
    xmid.II <- Individuals$Xmid.SENES[i]
    scal.II <- Individuals$Scal.SENES[i]
    Individuals$Fin.Senes[i] <- interpolateX(asym.II, xmid.II, scal.II, asym.II*0.95)
}
Individuals$Fin.Senes <- as.numeric(Individuals$Fin.Senes)

#       Figure ####  
ggplot(Individuals,aes(Species, Fin.Senes, fill = Status, color=Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = "Date of 95% Leaf Senescence", x = "Species") +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.15,.85),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  coord_flip()

#(10) WEEKLY PERCENT SENESCENCE -------------------------------------------------------.----
#       Create dataframes ####
# Create temporary dataframe
Temp <- data.frame(Plant_ID = "", Date_J = "", Senes = "")
# Create Expansion dataframe
Senescence <- data.frame(Plant_ID = "", Date_J = "", Senes = "")
Senescence <- Senescence[-1,]

#       Build observation lists ####
for (i in Individuals$Plant_ID){
  t <- ceiling(Individuals$Start.Senes[Individuals$Plant_ID == i])
  asym.I <- Individuals$Asym.SENES[Individuals$Plant_ID == i]
  xmid.I <- Individuals$Xmid.SENES[Individuals$Plant_ID == i]
  scal.I <- Individuals$Scal.SENES[Individuals$Plant_ID == i]
  senes <- interpolateY(asym.I, xmid.I, scal.I, t)
  
  while(t <= ceiling(max(Individuals$Fin.Senes))){ ##THIS LINE CHANGES WHETHER THE MODELS ENDS AT 5% FOR EACH PLANT OR THE LAST PLANT
    Temp$Plant_ID <- i
    Temp$Date_J <- t
    senes <-  interpolateY(asym.I, xmid.I, scal.I, t)
    senes <- as.numeric(senes)
    Temp$Senes <- senes
    Senescence <- rbind(Senescence, Temp)
    t <- t + 1
  }
}
Senescence$Senes <- as.numeric(Senescence$Senes)
Senescence$Date_J <- as.numeric(Senescence$Date_J)
Senescence$Area.Perc <- 100 - Senescence$Senes

# ####
# ================== SUMMER LEAF AREA =================================================.====
#(11) SUMMER AREA OVER TIME -----------------------------------------------------------.----
#       Create dataframes ####
Temp <- data.frame(Plant_ID = "", Date_J = "", Area.Perc = "")
FullLeafArea <- data.frame(Plant_ID = "", Date_J = "", Area.Perc = "")
FullLeafArea <- FullLeafArea[-1,]

#       Build observation lists ####
for (i in Individuals$Plant_ID[Individuals$Plant_ID != 85]){
  t <- ceiling(Individuals$Full.Expan[Individuals$Plant_ID == i])
  while(t < ceiling(Individuals$Start.Senes[Individuals$Plant_ID == i])){
    Temp$Plant_ID <- i
    Temp$Date_J <- t
    Temp$Area.Perc <- 100
    FullLeafArea <- rbind(FullLeafArea, Temp)
    t <- t + 1
  }
}

# ####
# ================== TOTAL LEAF AREA ==================================================.====
#(12) Calc total area (in m2) ---------------------------------------------------------.----
# Leaf count times average leaf area 
# convert to m2 -> divide by 1000 cm2 per meter 
Individuals$Total_leaf_area <- Individuals$Leaf_count * Individuals$Ave_leaf_area / 1000

#(13) SeasonSeriesArea ----------------------------------------------------------------.----
#       Merge with Individuals ####
Expansion$Leaf.Age <- "YOUNG"
Expansion <- merge.data.frame(Expansion, Individuals, by = "Plant_ID", all.x =  TRUE)

FullLeafArea$Leaf.Age <- "MATURE" 
FullLeafArea <- merge.data.frame(FullLeafArea, Individuals, by = "Plant_ID", all.x = TRUE)

Senescence$Senes <- NULL
Senescence$Leaf.Age <- "SENES"
Senescence <- merge.data.frame(Senescence, Individuals, by = "Plant_ID", all.x = TRUE)

#       Stack all instances of leaf area ####
SeasonSeriesArea <- rbind(Expansion, FullLeafArea, Senescence)

# ####
# ================== PHOTOSYNTHETIC LIGHT RESPONSE CURVES =============================.====
#(14) FIT CURVES ----------------------------------------------------------------------.----
# Nonlinear least squares regression (non-rectangular hyperbola) from Jason Heberling 
#       Individuals in June ####
#Create new variables
Individuals$Amax.JUN <- ""      #Am
Individuals$QuantYld.JUN <- ""  #AQY
Individuals$Resp.JUN <- ""      #Rd
Individuals$Theta.JUN <- ""     #theta

#Loop for all plants
par(mfrow=c(4,6), mar=c(4,4,1,1))
for (i in Individuals$Plant_ID) {
  #Plot each plant
  plot(Photo_JUN$Photo[Photo_JUN$Plant_ID == i] ~ Photo_JUN$PARi[Photo_JUN$Plant_ID == i],
       xlab = "Light (PPFD)", ylab = "A (umol co2/m2/s)")
  title(paste(Individuals$Species[Individuals$Plant_ID ==i], i), adj = .7, line = -5)
  
  #Fit monomolecular curve to each plant
  Photo.adj <- Photo_JUN$Photo[Photo_JUN$Plant_ID == i]
  PAR.adj <- Photo_JUN$PARi[Photo_JUN$Plant_ID == i]
  
  fit1 <- nls(Photo.adj ~ (1/(2*theta))*(AQY*PAR.adj+Am-sqrt((AQY*PAR.adj+Am)^2-4*AQY*theta*Am*PAR.adj))-Rd,start=list(Am=(max(Photo.adj)-min(Photo.adj)),AQY=0.05,Rd=-min(Photo.adj),theta=1)) 
  curve((1/(2*summary(fit1)$coef[4,1]))*(summary(fit1)$coef[2,1]*x+summary(fit1)$coef[1,1]-sqrt((summary(fit1)$coef[2,1]*x+summary(fit1)$coef[1,1])^2-4*summary(fit1)$coef[2,1]*summary(fit1)$coef[4,1]*summary(fit1)$coef[1,1]*x))-summary(fit1)$coef[3,1],lwd=2,col="blue",add=T)
      #fit1 code from: Heberling JM, Fridley JD (2013) Resource-use strategies of native and invasive plants in Eastern North American forests. New Phytol 200(2):523-533. doi: 10.1111/nph.12388
      #and based on: Marshall B, Biscoe PV (1980) A model for C3 leaves describing the dependence of net photosynthesis on irradiance. J Exp Bot 31(1):29-39. doi: 10.1093/jxb/31.1.29

  #Return values
  results1 <- data.frame(tidy(fit1))
  Amax <- results1$estimate[results1$term=="Am"]
  QuantYld <- results1$estimate[results1$term=="AQY"]
  r0 <- results1$estimate[results1$term=="Rd"]
  Theta <- results1$estimate[results1$term=="theta"]
  
  Individuals$Amax.JUN[Individuals$Plant_ID==i] <- Amax
  Individuals$QuantYld.JUN[Individuals$Plant_ID==i] <- QuantYld
  Individuals$Resp.JUN[Individuals$Plant_ID==i] <- r0
  Individuals$Theta.JUN[Individuals$Plant_ID==i] <- Theta
}
results1 <- data.frame(tidy(fit1))

#Convert to numbers
Individuals$Amax.JUN <- as.numeric(Individuals$Amax.JUN)
Individuals$QuantYld.JUN  <- as.numeric( Individuals$QuantYld.JUN)
Individuals$Resp.JUN <- as.numeric(Individuals$Resp.JUN)
Individuals$Theta.JUN <- as.numeric(Individuals$Theta.JUN)

# Status Figure
ggplot(Individuals,aes(Status, Amax.JUN, fill=Status, color=Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = "June Amax (umol co2/s/m2)", x = "JUNE") +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL,  guide=FALSE, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.5,.9),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#Species Figure
ggplot(Individuals,aes(Species, Amax.JUN, fill=Status, color=Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = " June Amax (umol co2/s/m2)", x = "Species") +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL,  guide=FALSE, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.2,.9),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#       Individuals in August ####
#Create new variables
Individuals$Amax.AUG <- ""      #Am
Individuals$QuantYld.AUG <- ""  #AQY
Individuals$Resp.AUG <- ""      #Rd
Individuals$Theta.AUG <- ""     #theta

#Loop for all plants
par(mfrow=c(4,6), mar=c(4,4,1,1))
for (i in Individuals$Plant_ID) {
  #Plot each plant
  plot(Photo_AUG$Photo[Photo_AUG$Plant_ID == i] ~ Photo_AUG$PARi[Photo_AUG$Plant_ID == i],
       xlab = "Light (PPFD)", ylab = "A (umol co2/m2/s)")
  title(paste(Individuals$Species[Individuals$Plant_ID ==i], i), adj = .7, line = -5)
  
  #Fit monomolecular curve to each plant
  Photo.adj <- Photo_AUG$Photo[Photo_AUG$Plant_ID == i]
  PAR.adj <- Photo_AUG$PARi[Photo_AUG$Plant_ID == i]
  
  fit1 <- nls(Photo.adj ~ (1/(2*theta))*(AQY*PAR.adj+Am-sqrt((AQY*PAR.adj+Am)^2-4*AQY*theta*Am*PAR.adj))-Rd,start=list(Am=(max(Photo.adj)-min(Photo.adj)),AQY=0.05,Rd=-min(Photo.adj),theta=1)) 
  curve((1/(2*summary(fit1)$coef[4,1]))*(summary(fit1)$coef[2,1]*x+summary(fit1)$coef[1,1]-sqrt((summary(fit1)$coef[2,1]*x+summary(fit1)$coef[1,1])^2-4*summary(fit1)$coef[2,1]*summary(fit1)$coef[4,1]*summary(fit1)$coef[1,1]*x))-summary(fit1)$coef[3,1],lwd=2,col="blue",add=T)
       #fit1 code from: Heberling JM, Fridley JD (2013) Resource-use strategies of native and invasive plants in Eastern North American forests. New Phytol 200(2):523-533. doi: 10.1111/nph.12388
       #and based on: Marshall B, Biscoe PV (1980) A model for C3 leaves describing the dependence of net photosynthesis on irradiance. J Exp Bot 31(1):29-39. doi: 10.1093/jxb/31.1.29
  
  #Return values
  results1 <- data.frame(tidy(fit1))
  Amax <- results1$estimate[results1$term=="Am"]
  QuantYld <- results1$estimate[results1$term=="AQY"]
  r0 <- results1$estimate[results1$term=="Rd"]
  Theta <- results1$estimate[results1$term=="theta"]
  
  Individuals$Amax.AUG[Individuals$Plant_ID==i] <- Amax
  Individuals$QuantYld.AUG[Individuals$Plant_ID==i] <- QuantYld
  Individuals$Resp.AUG[Individuals$Plant_ID==i] <- r0
  Individuals$Theta.AUG[Individuals$Plant_ID==i] <- Theta
}
results1 <- data.frame(tidy(fit1))

#Convert to numbers
Individuals$Amax.AUG <- as.numeric(Individuals$Amax.AUG)
Individuals$QuantYld.AUG  <- as.numeric( Individuals$QuantYld.AUG)
Individuals$Resp.AUG <- as.numeric(Individuals$Resp.AUG)
Individuals$Theta.AUG <- as.numeric(Individuals$Theta.AUG)

# Status Figure
ggplot(Individuals,aes(Status, Amax.AUG, fill=Status, color=Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = "AUG Amax (umol co2/s/m2)", x = "AUGUST") +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL,  guide=FALSE, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.5,.9),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#Species Figure
ggplot(Individuals,aes(Species, Amax.AUG, fill=Status, color=Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = " AUG Amax (umol co2/s/m2)", x = NULL) +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL,  guide=FALSE, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.75,.9),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#Check Distributions of Amax
par(mfrow = c(1,2))
hist(Individuals$Amax.JUN, main = "Young Leaves", xlab = "Amax (umol C m-2 s-1)")
hist(Individuals$Amax.AUG, main = "Mature Leaves", xlab = "Amax (umol C m-2 s-1)")

# ####
# ================== UNDERSTORY LIGHT ENVIRONMENT =====================================.====
#(15) CALIBRATION CURVES for Spot vs. LUX ---------------------------------------------.----
#       Combine seasons ####
calib_all$Season2 <- ""
for(i in 1:nrow(calib_all)){
  if (calib_all$Season[i] == "SummerI" | calib_all$Season[i] == "SummerII"){
    calib_all$Season2[i] <- "Summer"
  } else if (calib_all$Season[i] == "Spring"){
    calib_all$Season2[i] <- "Spring"
  } else if (calib_all$Season[i] == "Fall"){
    calib_all$Season2[i] <- "Fall"
  }
}

#       View data as scatterplot ####
# Untransformed
ggplot(calib_all, aes(colour = calib_all$Season2, calib_all$LUX_plant, calib_all$PPFD_plant)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "PPFD (umol photons m-2 s-1)", x = "LUX (lumens m-2)") +
  scale_colour_manual(name = NULL, values = c("orange", "forest green", "black")) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

# Logged
ggplot(calib_all, aes(colour = calib_all$Season2, log(calib_all$LUX_plant, base=10), 
                      log(calib_all$PPFD_plant, base=10))) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "PPFD (umol photons m-2 s-1)", x = "LUX (lumens m-2)") +
  scale_colour_manual(name = NULL, values = c("orange", "forest green", "black")) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

#       Fit Regessions ####
# Spring
par(mfrow=c(1,2))
hist(sqrt(calib_all$PPFD_plant[calib_all$Season2 == "Spring"]), main="Spring PPFD", xlab=NULL)
hist(sqrt(calib_all$LUX_plant[calib_all$Season2 == "Spring"]), main="Spring LUX", xlab=NULL)

par(mfrow=c(2,2))
Season.Calib <- lm(log(calib_all$PPFD_plant[calib_all$Season2 == "Spring"], base = 10)~log(calib_all$LUX_plant[calib_all$Season2 == "Spring"], base =10))
plot(Season.Calib)
summary(Season.Calib)
Calib.R <- data.frame(tidy(Season.Calib))
Light_Calib_Results$Int[Light_Calib_Results$Season == "Spring"] <- Calib.R[1,2]
Light_Calib_Results$Slope[Light_Calib_Results$Season == "Spring"] <- Calib.R[2,2]

# Summer
par(mfrow=c(1,2))
hist(log(calib_all$PPFD_plant[calib_all$Season2 == "Summer"], base=10), main = "Summer PPFD", xlab=NULL)
hist(log(calib_all$LUX_plant[calib_all$Season2 == "Summer"], base=10), main = "Summer LUX", xlab=NULL)

par(mfrow=c(2,2))
Season.Calib <- lm(log(calib_all$PPFD_plant[calib_all$Season2 == "Summer"], base=10)~log(calib_all$LUX_plant[calib_all$Season2 == "Summer"], base=10))
plot(Season.Calib)
summary(Season.Calib)
Calib.R <- data.frame(tidy(Season.Calib))
Light_Calib_Results$Int[Light_Calib_Results$Season == "Summer"] <- Calib.R[1,2]
Light_Calib_Results$Slope[Light_Calib_Results$Season == "Summer"] <- Calib.R[2,2]

# Fall 
par(mfrow=c(1,2))
hist(log(calib_all$PPFD_plant[calib_all$Season2 == "Fall"], base=10), main = "Fall PPFD", xlab=NULL)
hist(log(calib_all$LUX_plant[calib_all$Season2 == "Fall"], base=10), main = "Fall LUX", xlab=NULL)

par(mfrow=c(2,2))
Season.Calib <- lm(log(calib_all$PPFD_plant[calib_all$Season == "Fall"], base = 10)~log(calib_all$LUX_plant[calib_all$Season == "Fall"], base=10))
plot(Season.Calib)
summary(Season.Calib)
Calib.R <- data.frame(tidy(Season.Calib))
Light_Calib_Results$Int[Light_Calib_Results$Season == "Fall"] <- Calib.R[1,2]
Light_Calib_Results$Slope[Light_Calib_Results$Season == "Fall"] <- Calib.R[2,2]

# Convert to numbers
Light_Calib_Results$Int <- as.numeric(Light_Calib_Results$Int)
Light_Calib_Results$Slope <- as.numeric(Light_Calib_Results$Slope)

#(16) CONVERT DIURNAL LIGHT from LUX to PAR -------------------------------------------.----
#       Diurnal Meas ####
# Remove 0 light
Diurn2 <- subset.data.frame(Diurn, Diurn$LUX !=0)

# Tranform LUX (log)
Diurn2$LUX.log <- log(Diurn2$LUX, base = 10)

# Convert LUX to PAR (for diurnal curves)
Diurn2$PPFD.log <- "" 
for (i in 1:nrow(Diurn2)){
  m <- Light_Calib_Results$Slope[Light_Calib_Results$Month == Diurn2$Month[i]][1]
  b <- Light_Calib_Results$Int[Light_Calib_Results$Month == Diurn2$Month[i]][1]
  l <- Diurn2$LUX.log[i][1]
  Diurn2$PPFD.log[i] <- m*l + b 
}

# Transform PAR (un-log)
Diurn2$PPFD.log <- as.numeric(Diurn2$PPFD.log)
Diurn2$PPFD <- 10^Diurn2$PPFD.log

#       Spot Meas ####
# Tranform LUX (log)
calib_all$LUX_full_sun_log <- log(calib_all$LUX_full_sun, base = 10)

# Convert LUX to PAR (for spot meas)
calib_all$PPFD_full_sun_log <- "" 
for (i in 1:nrow(calib_all)){
  m <- Light_Calib_Results$Slope[Light_Calib_Results$Season == calib_all$Season2[i]][1]
  b <- Light_Calib_Results$Int[Light_Calib_Results$Season == calib_all$Season2[i]][1]
  l <- calib_all$LUX_full_sun_log[i][1]
  calib_all$PPFD_full_sun_log[i] <- m*l + b 
}

# Transform PAR (un-log)
calib_all$PPFD_full_sun_log <- as.numeric(calib_all$PPFD_full_sun_log)
calib_all$PPFD_full_sun <- 10^calib_all$PPFD_full_sun_log

#Calulate % PPFD
calib_all$PPFD_percent <- calib_all$PPFD_plant / calib_all$PPFD_full_sun *100 

# ================= UNDERSTORY LIGHT SUMMARY ==========================================.====
#(17) SUMMARIZE DAILY LIGHT (area under diurnal curve)---------------------------------.----
#       Preparation ####
# Convert to seconds
Diurn2$Time.s <- HmsToSec(Diurn2$Time)

# Plot
par(mfrow=(c(4,6)))
for (i in unique(Diurn2$Logger_No)){
  for (j in unique(Diurn2$Date)){
    plot(Diurn2$Time.s[Diurn2$Date==j & Diurn2$Logger_No == i],
         Diurn2$PPFD[Diurn2$Date==j & Diurn2$Logger_No == i],
         main = paste(i,", ", j),
         xlab = "Time (sec)", ylab = "PPFD",
         type = "l", col="blue")
  }
}

#       Calculate Area under the Curve ####
# AUC Calulation by HOBO and Date
DailyLight <- Diurn2 %>%
  group_by(Logger_No, Date) %>%
  summarize(PPFD = AUC(Time.s, PPFD, method = "spline"),
            LUX = AUC(Time.s, LUX, method = "spline"))

# Add Months back in
Diurn.Month <- subset.data.frame(Diurn2, select=c(Date, Month))
Diurn.Month <- Diurn.Month[!duplicated(Diurn.Month), ] # remove repeated rows
DailyLight <- merge.data.frame(DailyLight, Diurn.Month, by = "Date", all.y = FALSE)
DailyLight$Logger_No <- as.character(DailyLight$Logger_No)

# Convert to mols
DailyLight$PPFD.umol <- DailyLight$PPFD
DailyLight$PPFD <- DailyLight$PPFD/1000000

#       All Bin Average Light ####
DailyLight.bin.ave <- DailyLight %>%
  group_by(Date) %>%
  summarise(PPFD = mean (PPFD),
            LUX = mean (LUX))

# Add Months back in
DailyLight.bin.ave <- merge.data.frame(DailyLight.bin.ave, Diurn.Month, by = "Date", all.y = FALSE)

#       Summarize ####
ave(DailyLight.bin.ave$PPFD[DailyLight.bin.ave$Month=="MAY"])[1]
ave(DailyLight.bin.ave$PPFD[DailyLight.bin.ave$Month=="JUL"|DailyLight.bin.ave$Month=="SEP"])[1]
ave(DailyLight.bin.ave$PPFD[DailyLight.bin.ave$Month=="OCT"])[1]

min(DailyLight$PPFD[DailyLight$Month=="MAY"])[1]
min(DailyLight$PPFD[DailyLight$Month=="JUL"|DailyLight$Month=="SEP"])[1]
min(DailyLight$PPFD[DailyLight$Month=="OCT"])[1]

max(DailyLight$PPFD[DailyLight$Month=="MAY"])[1]
max(DailyLight$PPFD[DailyLight$Month=="JUL"|DailyLight$Month=="SEP"])[1]
max(DailyLight$PPFD[DailyLight$Month=="OCT"])[1]

#       Visualize ####
# Canopy context (from the prior year)
Can.Light <-data.frame(xmin=as.Date(c(120,161,269), origin = as.Date("2018-01-01")), 
                       xmax=as.Date(c(161,269,310), origin = as.Date("2018-01-01")),
                       Can.Open=c("Open", "Closed", "Open"))

# All Bins
DailyPPFD.plot <- ggplot() +
  geom_rect(data = Can.Light, aes(xmin = xmin, ymin=-Inf, xmax = xmax, ymax=Inf, fill = Can.Open)) +
  geom_point(data = DailyLight, aes(Date, PPFD, color = Logger_No), size = 4) + 
  #geom_line() +
  labs(x = NULL, y = "Daily PPFD (mol photons day-1)") +
  theme_classic() +
  scale_color_manual(name = "Light Bin", values = c("#33FF00", "993333", "#FF9933", "#003399", "#9933CC", "#334455")) +
  scale_fill_manual(name = "Canopy", values = c("Open" = "white", "Closed" = "#e6e6e6")) +
  theme(text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  scale_x_date(date_breaks = '1 month', date_labels = '%b') 
DailyPPFD.plot

# Averages
ave.DailyPPFD.plot <- ggplot(DailyLight.bin.ave, aes(Date, PPFD)) +
  geom_point(size = 4) + 
  labs(x = NULL, y = "Average Total Daily Light \n(mol photons day-1)") +
  theme_classic() +
  theme(text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))
ave.DailyPPFD.plot

#       Average Total Daily Light ####
#Add light bins to main data
for (i in 1:nrow(Individuals)){
  Individuals$Bin.Spr[i]    <- mean(calib_all$Bin[calib_all$Season=="Spring"   & calib_all$Plant_ID==Individuals$Plant_ID[i]], na.rm = T)
  Individuals$Bin.Sum.I[i]  <- mean(calib_all$Bin[calib_all$Season=="SummerI"  & calib_all$Plant_ID==Individuals$Plant_ID[i]], na.rm = T)
  Individuals$Bin.Sum.II[i] <- mean(calib_all$Bin[calib_all$Season=="SummerII" & calib_all$Plant_ID==Individuals$Plant_ID[i]], na.rm = T)
  Individuals$Bin.Fall[i]   <- mean(calib_all$Bin[calib_all$Season=="Fall"     & calib_all$Plant_ID==Individuals$Plant_ID[i]], na.rm = T)
}

Individuals$Bin.Fall[Individuals$Plant_ID == 74] <- 1

#Average per month
DailyPPFD <- DailyLight %>%
  group_by(Logger_No, Month) %>%
  summarize(PPFD = mean(PPFD))
DailyPPFD.bin <- spread(DailyPPFD, Month, PPFD)

#MAY
DailyPPFD.MAY <- subset.data.frame(DailyPPFD.bin, select = c(Logger_No, MAY))
colnames(DailyPPFD.MAY)[1] <- "Bin.Spr"
colnames(DailyPPFD.MAY)[2] <- "PPFD.MAY"
Individuals <- merge.data.frame(Individuals, DailyPPFD.MAY, by = "Bin.Spr", all.x=TRUE)

#JUL
DailyPPFD.JUL <- subset.data.frame(DailyPPFD.bin, select = c(Logger_No, JUL))
colnames(DailyPPFD.JUL)[1] <- "Bin.Sum.I"
colnames(DailyPPFD.JUL)[2] <- "PPFD.JUL"
Individuals <- merge.data.frame(Individuals, DailyPPFD.JUL, by = "Bin.Sum.I", all.x=TRUE)

#SEP
DailyPPFD.SEP <- subset.data.frame(DailyPPFD.bin, select = c(Logger_No, SEP))
colnames(DailyPPFD.SEP)[1] <- "Bin.Sum.II"
colnames(DailyPPFD.SEP)[2] <- "PPFD.SEP"
Individuals <- merge.data.frame(Individuals, DailyPPFD.SEP, by = "Bin.Sum.II", all.x=TRUE)

#OCT
DailyPPFD.OCT <- subset.data.frame(DailyPPFD.bin, select = c(Logger_No, OCT))
colnames(DailyPPFD.OCT)[1] <- "Bin.Fall"
colnames(DailyPPFD.OCT)[2] <- "PPFD.OCT"
Individuals <- merge.data.frame(Individuals, DailyPPFD.OCT, by = "Bin.Fall", all.x=TRUE)

#(18) COMPARE WITH DAYLENGTH ----------------------------------------------------------.----
#       Add daylength to AUC data ####
DailyLight <- merge.data.frame(DailyLight, Day.Length, by = "Date", all.x =TRUE)
DailyLight.bin.ave <- merge.data.frame(DailyLight.bin.ave, Day.Length, by = "Date", all.x = TRUE)

#       Visualize ####
# Bin PPFD total
daylen.plot <- ggplot(DailyLight, aes(DailyLight$daylen, DailyLight$PPFD)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "Light (mol photons day-1)", x = "Day Length (minutes)") +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

# Bin PPFD by month
daylen.bin.mon.plot <-ggplot(DailyLight, aes(color = DailyLight$Month, DailyLight$daylen, DailyLight$PPFD)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "Light (mol photons day-1)", x = "Day Length (minutes)") +
  scale_colour_manual(name = NULL, values = c("forest green", "blue", "orange", "black")) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

# Ave PPFD by Month
daylen.ave.mon.plot <- ggplot(DailyLight.bin.ave, aes(color = DailyLight.bin.ave$Month, DailyLight.bin.ave$daylen, DailyLight.bin.ave$PPFD)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "Light (mol photons day-1)", x = "Day Length (minutes)") +
  scale_colour_manual(name = NULL, values = c("forest green", "blue", "orange", "black")) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

# Bin PPFD by bin
daylen.bin.plot <- ggplot(DailyLight, aes(color = DailyLight$Logger_No, DailyLight$daylen, DailyLight$PPFD)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "Light (mol photons day-1)", x = "Day Length (minutes)") +
  scale_colour_manual(name = NULL, values = blues9) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

grid.arrange(daylen.plot, daylen.bin.mon.plot, daylen.bin.plot, daylen.ave.mon.plot,
             ncol = 2, nrow = 2)

#(19) DAILY CO2 SEASONAL REGRESSIONS (SeasonSeriesLight)-------------------------------.----
#       Create Data frames ####
# Create temporary dataframe
Temp.L <- data.frame(Date_J = "", PPFD = "")
# Create Expansion dataframe
SeasonSeriesLight <- data.frame(Date_J = "", PPFD = "")
SeasonSeriesLight <- SeasonSeriesLight[-1,]

#       Fit polynomial regressions and build data frame ####
x.data = DailyLight.bin.ave$Date_J
y.data = DailyLight.bin.ave$PPFD
par(mfrow=c(1,1))
plot(x.data, y.data, xlab = "Date", ylab = "PPFD")
Light.Mod1 <- lm(y.data ~ poly(x.data, 5))
predict.Aday <- predict(Light.Mod1, data.frame(x=x.data))
lines(x.data, predict.Aday, col='green',lwd=2)
  
Date_J.temp <- 100
while(Date_J.temp <= 344){
  Temp.L$Date_J <- Date_J.temp
  A <- predict(Light.Mod1, data.frame(x.data=c(Date_J.temp)))
  A <- as.numeric(A)
  Temp.L$PPFD <- A
  SeasonSeriesLight <- rbind(SeasonSeriesLight, Temp.L)
  Date_J.temp <- Date_J.temp + 1
}
SeasonSeriesLight$PPFD <- as.numeric(SeasonSeriesLight$PPFD)
SeasonSeriesLight$Date_J <- as.numeric(SeasonSeriesLight$Date_J)

#Correct for unrealistically high and low light levels
SeasonSeriesLight$PPFD[SeasonSeriesLight$PPFD > max(DailyLight.bin.ave$PPFD)] <- max(DailyLight.bin.ave$PPFD)
SeasonSeriesLight$PPFD[SeasonSeriesLight$PPFD < min(DailyLight.bin.ave$PPFD)] <- min(DailyLight.bin.ave$PPFD)

#(20) SEASONAL % LIGHT (area under the seasonal curve) --------------------------------.----
#       Define Seasons ####
#Average spring and fall phenology
Bud.Burst <- 110 #Average Julian date of all plants observed in 2017
Leaf.Drop <- 291 #Average Julian date of all plants observed in 2017 (75% leaf senescence)

#Subset daily A data
SeasonSeriesLight <- SeasonSeriesLight[SeasonSeriesLight$Date_J >= Bud.Burst & 
                                         SeasonSeriesLight$Date_J <= Leaf.Drop,]

Spring.Series.L <- SeasonSeriesLight[SeasonSeriesLight$Date_J <= 161,] 
Sum.Series.L <- SeasonSeriesLight[SeasonSeriesLight$Date_J > 161 & SeasonSeriesLight$Date_J < 269,] 
Fall.Series.L <- SeasonSeriesLight[SeasonSeriesLight$Date_J > 269,] 

#       Spring ####
par(mfrow=c(1,3))
plot(Spring.Series.L$Date_J, Spring.Series.L$PPFD,
       xlab = "Date", ylab = "PPFD",
       type = "l", col="blue")
SPR.L <- AUC(Spring.Series.L$Date_J, Spring.Series.L$PPFD, method = "spline")

#       Summer ####
plot(Sum.Series.L$Date_J, Sum.Series.L$PPFD,
     xlab = "Date", ylab = "PPFD",
     type = "l", col="blue")
SUM.L <- AUC(Sum.Series.L$Date_J, Sum.Series.L$PPFD, method = "spline")

#       Fall ####
plot(Fall.Series.L$Date_J, Fall.Series.L$PPFD,
     xlab = "Date", ylab = "PPFD",
     type = "l", col="blue")
FALL.L <- AUC(Fall.Series.L$Date_J, Fall.Series.L$PPFD, method = "spline")

#       Total Percents ####
Tot.L <- SPR.L + SUM.L + FALL.L
(SPR.L/Tot.L)*100
(SUM.L/Tot.L)*100
(FALL.L/Tot.L)*100

#(21) SUMMARY DATA --------------------------------------------------------------------.----
#       % PPFD and % LUX ####
#Spring
calib_all.spr <- subset.data.frame(calib_all, calib_all$Season=="Spring", select = c(Plant_ID, PPFD_percent, LUX_percent))
calib_all.spr <- calib_all.spr %>%
  group_by(Plant_ID) %>%
  summarize(PPFD_percent = mean(PPFD_percent, na.rm=TRUE),
            LUX_percent = mean(LUX_percent, na.rm=TRUE))
colnames(calib_all.spr)[2] <- "PPFD_percent.spr"
colnames(calib_all.spr)[3] <- "LUX_percent.spr"
Individuals <- merge.data.frame(Individuals, calib_all.spr, by = "Plant_ID",  all.x = TRUE, all.y = FALSE)

ave(calib_all$PPFD_percent[calib_all$Season=="Spring"])[1]
min(calib_all$PPFD_percent[calib_all$Season=="Spring"])[1]
max(calib_all$PPFD_percent[calib_all$Season=="Spring"])[1]

#Summer I
calib_all.sumI <- subset.data.frame(calib_all, calib_all$Season=="SummerI", select = c(Plant_ID, PPFD_percent, LUX_percent))
calib_all.sumI <- calib_all.sumI %>%
  group_by(Plant_ID) %>%
  summarize(PPFD_percent = mean(PPFD_percent, na.rm=TRUE),
            LUX_percent = mean(LUX_percent, na.rm=TRUE))
colnames(calib_all.sumI)[2] <- "PPFD_percent.sumI"
colnames(calib_all.sumI)[3] <- "LUX_percent.sumI"
Individuals <-merge.data.frame(Individuals, calib_all.sumI, by = "Plant_ID",   all.x = TRUE, all.y = FALSE)

ave(calib_all$PPFD_percent[calib_all$Season=="SummerI"])[1]
min(calib_all$PPFD_percent[calib_all$Season=="SummerI"])[1]
max(calib_all$PPFD_percent[calib_all$Season=="SummerI"])[1]

#Summer II
calib_all.sumII <- subset.data.frame(calib_all, calib_all$Season=="SummerII", select = c(Plant_ID, PPFD_percent, LUX_percent))
calib_all.sumII <- calib_all.sumII %>%
  group_by(Plant_ID) %>%
  summarize(PPFD_percent = mean(PPFD_percent, na.rm=TRUE),
            LUX_percent = mean(LUX_percent, na.rm=TRUE))
colnames(calib_all.sumII)[2] <- "PPFD_percent.sumII"
colnames(calib_all.sumII)[3] <- "LUX_percent.sumII"
Individuals <-merge.data.frame(Individuals, calib_all.sumII, by = "Plant_ID",   all.x = TRUE, all.y = FALSE)

ave(calib_all$PPFD_percent[calib_all$Season=="SummerII"])[1]
min(calib_all$PPFD_percent[calib_all$Season=="SummerII"])[1]
max(calib_all$PPFD_percent[calib_all$Season=="SummerII"])[1]

ave(calib_all$PPFD_percent[calib_all$Season=="SummerI" | calib_all$Season=="Summer II"])[1]
min(calib_all$PPFD_percent[calib_all$Season=="SummerI" | calib_all$Season=="Summer II"])[1]
max(calib_all$PPFD_percent[calib_all$Season=="SummerI" | calib_all$Season=="Summer II"])[1]
SD(calib_all$PPFD_percent[calib_all$Season=="SummerI" | calib_all$Season=="Summer II"])[1]

#Fall
calib_all.fall <- subset.data.frame(calib_all, calib_all$Season=="Fall", select = c(Plant_ID, PPFD_percent, LUX_percent))
calib_all.fall <- calib_all.fall %>%
  group_by(Plant_ID) %>%
  summarize(PPFD_percent = mean(PPFD_percent, na.rm=TRUE),
            LUX_percent = mean(LUX_percent, na.rm=TRUE))
colnames(calib_all.fall)[2] <- "PPFD_percent.fall"
colnames(calib_all.fall)[3] <- "LUX_percent.fall"
Individuals <-merge.data.frame(Individuals, calib_all.fall, by = "Plant_ID",  all.x = TRUE, all.y = FALSE)

ave(calib_all$PPFD_percent[calib_all$Season=="Fall"])[1]
min(calib_all$PPFD_percent[calib_all$Season=="Fall"])[1]
max(calib_all$PPFD_percent[calib_all$Season=="Fall"])[1]

#       Spot PPFD ####
ave(calib_all$PPFD_plant[calib_all$Season2=="Spring"])[1]
min(calib_all$PPFD_plant[calib_all$Season2=="Spring"])[1]
max(calib_all$PPFD_plant[calib_all$Season2=="Spring"])[1]

ave(calib_all$PPFD_plant[calib_all$Season2=="Summer"])[1]
min(calib_all$PPFD_plant[calib_all$Season2=="Summer"])[1]
max(calib_all$PPFD_plant[calib_all$Season2=="Summer"])[1]

ave(calib_all$PPFD_plant[calib_all$Season2=="Fall"])[1]
min(calib_all$PPFD_plant[calib_all$Season2=="Fall"])[1]
max(calib_all$PPFD_plant[calib_all$Season2=="Fall"])[1]

# ####
# ================== DAILY CO2 ASSIMILATION ===========================================.====
#(22) PAIR LIGHT DATA WITH PHOTOSYNTHESIS PARAMETERS ----------------------------------.----
#Combine photosynthesis data for each indivual with diurnal light data

# Select necessary columns
Diurn3 <- subset.data.frame(Diurn2, select = c(Logger_No, Month:Time.s))

#       Break down by monthly diurnal measurements ####       
# Spring
Diurn_MAY_conv <- subset.data.frame(Diurn3, Diurn3$Month == "MAY")
Diurn_MAY_conv$Bin.Spr <- Diurn_MAY_conv$Logger_No #Spring Logger_No is same as Bin.No.Spr, but name needs to match "Individuals" dataset
Diurn_MAY_Ind <- merge.data.frame(x = Individuals, y = Diurn_MAY_conv, by = "Bin.Spr", y.all = TRUE)

# Summer I
Diurn_JUL_conv <- subset.data.frame(Diurn3, Diurn3$Month == "JUL")
Diurn_JUL_conv$Bin.Sum.I <- Diurn_JUL_conv$Logger_No 
Diurn_JUL_Ind <- merge.data.frame(x = Individuals, y = Diurn_JUL_conv, by = "Bin.Sum.I")

# Summer II
Diurn_SEP_conv <- subset.data.frame(Diurn3, Diurn3$Month == "SEP")
Diurn_SEP_conv$Bin.Sum.II <- Diurn_SEP_conv$Logger_No 
Diurn_SEP_Ind <- merge.data.frame(x = Individuals, y = Diurn_SEP_conv, by = "Bin.Sum.II")

# Fall
Diurn_OCT_conv <- subset.data.frame(Diurn3, Diurn3$Month == "OCT")
Diurn_OCT_conv$Bin.Fall <- Diurn_OCT_conv$Logger_No 
Diurn_OCT_Ind <- merge.data.frame(x = Individuals, y = Diurn_OCT_conv, by = "Bin.Fall")

#       Reassemble ####
Diurn.Ind <- rbind.data.frame(Diurn_MAY_Ind, Diurn_JUL_Ind, Diurn_SEP_Ind, Diurn_OCT_Ind)
Diurn.Ind <- subset.data.frame(Diurn.Ind, select = c(Plant_ID:Ave_leaf_area,Amax.JUN:Theta.AUG,
                                                     Month:Time.s))

#(23) CONVERT PPFD TO A ---------------------------------------------------------------.----
#       Calculate net carbon gain for young leaves (takes ~2 min to run) ####
Diurn.Ind$A.sec.JUN <- "" 
for (i in 1:nrow(Diurn.Ind)) {
  Am <- Diurn.Ind$Amax.JUN[i]
  AQY <- Diurn.Ind$QuantYld.JUN[i]
  Rd <- Diurn.Ind$Resp.JUN[i]
  theta <- Diurn.Ind$Theta.JUN[i]
  PARlrc <- Diurn.Ind$PPFD[i]
  Diurn.Ind$A.sec.JUN[i] <- (1/(2*theta))*(AQY*PARlrc+Am-sqrt((AQY*PARlrc+Am)^2-4*AQY*theta*Am*PARlrc))-Rd
}
Diurn.Ind$A.sec.JUN <- as.numeric(Diurn.Ind$A.sec.JUN)

# Remove respiration counfounding factor
Diurn.Ind$A.sec.JUN.adj <- Diurn.Ind$A.sec.JUN + Diurn.Ind$Resp.JUN

#       Calculate net carbon gain for mature leaves (takes ~2 min to run) ####
Diurn.Ind$A.sec.AUG <- "" 
for (i in 1:nrow(Diurn.Ind)) {
  Am <- Diurn.Ind$Amax.AUG[i]
  AQY <- Diurn.Ind$QuantYld.AUG[i]
  Rd <- Diurn.Ind$Resp.AUG[i]
  theta <- Diurn.Ind$Theta.AUG[i]
  PARlrc <- Diurn.Ind$PPFD[i]
  Diurn.Ind$A.sec.AUG[i] <- (1/(2*theta))*(AQY*PARlrc+Am-sqrt((AQY*PARlrc+Am)^2-4*AQY*theta*Am*PARlrc))-Rd
}
Diurn.Ind$A.sec.AUG <- as.numeric(Diurn.Ind$A.sec.AUG)

# Remove respiration counfounding factor
Diurn.Ind$A.sec.AUG.adj <- Diurn.Ind$A.sec.AUG + Diurn.Ind$Resp.AUG

#(24) CALCULATE TOTAL A FOR EACH DAY (area under curve) -------------------------------.----
Diurn.Ind$Date_Time <- as.POSIXct(paste(Diurn.Ind$Date, Diurn.Ind$Time), format="%Y-%m-%d %H:%M:%S")
Diurn.Ind <- Diurn.Ind[order(Diurn.Ind$Date_Time),]

#       Plot ####
Diurn.Ind <- Diurn.Ind[with(Diurn.Ind, order(Plant_ID, Date, Time.s)),] #order rows

par(mfrow=(c(5,6)))
for (i in unique(Diurn.Ind$Plant_ID)){
  for (j in unique(Diurn.Ind$Date)){
    plot(Diurn.Ind$Time.s[Diurn.Ind$Date==j & Diurn.Ind$Plant_ID == i],
         Diurn.Ind$A.sec.JUN.adj[Diurn.Ind$Date==j & Diurn.Ind$Plant_ID == i],
         main = paste(i,", ", j),
         xlab = "Time (sec)", ylab = "A sec (umol C m-2 s-1)",
         type = "l", col="blue")
  }
}

#       Calculate area under the curve (Total Aday for young leaves) ####
DailyA.date <- Diurn.Ind %>%
  group_by(Plant_ID, Date) %>%
  summarize(Aday.YOUNG = AUC(Time.s, A.sec.JUN.adj, method = "spline"),
            Aday.MATURE = AUC(Time.s, A.sec.AUG.adj, method = "spline"))

#       Convert from umol to mol ####
DailyA.date$Aday.YOUNG <- DailyA.date$Aday.YOUNG / 1000000
DailyA.date$Aday.MATURE <- DailyA.date$Aday.MATURE / 1000000

#       Add months back in ####
Diurn.Month <- subset.data.frame(Diurn.Ind, select=c(Date, Month))
Diurn.Month <- Diurn.Month[!duplicated(Diurn.Month), ] # remove repeated rows
DailyA.date <- merge.data.frame(DailyA.date, Diurn.Month, by = "Date", all.y = FALSE)

#(25) COMPARE DAILY A to DAYLENGTH ----------------------------------------------------.----
#       Add daylength to AUC data ####
DailyA.date <- merge.data.frame(DailyA.date, Day.Length, by = "Date", all.x =TRUE)

#       Visualize ####
# Daily A by month
ggplot(DailyA.date, aes(color = DailyA.date$Month, DailyA.date$daylen, DailyA.date$Aday.YOUNG)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "Daily A (mol C s-1 m-2)", x = "Day Length (minutes)") +
  scale_colour_manual(name = NULL, values = c("forest green", "blue", "orange", "black")) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

# Daily A by month
ggplot(DailyA.date, aes(color = DailyA.date$Month, DailyA.date$Date, DailyA.date$Aday.YOUNG)) +
  geom_point(size = 4) +
  theme_classic() +
  labs(y = "Daily A (mol C s-1 m-2)", x = NULL) +
  scale_colour_manual(name = NULL, values = c("forest green", "blue", "orange", "black")) +
  stat_quantile(method = "lm", formula=y~x, size=2) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18),
        panel.border = element_rect(colour = "black", fill=NA, size=.5))

#(26) DAILY CO2 SEASONAL REGRESSIONS (SeasonSeriesA)-----------------------------------.----
#       Stack young and mature daily A ####
DailyA.L <- gather(DailyA.date, Age, Aday, Aday.YOUNG, Aday.MATURE)
for (i in DailyA.L$Age){
  if(i == "Aday.YOUNG"){
    DailyA.L$Age[DailyA.L$Age == i] <- "YOUNG"
  }else if(i == "Aday.MATURE"){
    DailyA.L$Age[DailyA.L$Age == i] <- "MATURE"
  }
}
#       Create Data frames ####
# Create temporary dataframe
Temp <- data.frame(Plant_ID = "", Date_J = "", Age = "", Aday.lf = "")
# Create Expansion dataframe
SeasonSeriesA <- data.frame(Plant_ID = "", Date_J = "", Age = "", Aday.lf = "")
SeasonSeriesA <- SeasonSeriesA[-1,]
#Create preliminary dataset
SSL.pre <- DailyA.L[0,]

#       Fit polynomial regressions and build data frame ####
# Subset DailyA to include proper leaf ages
for (i in 1:nrow(DailyA.L)){
  expan.date <- Individuals$Full.Expan[Individuals$Plant_ID == DailyA.L$Plant_ID[i]]
  if(DailyA.L$Date_J[i] <= expan.date & DailyA.L$Age[i] == "YOUNG" |
     DailyA.L$Date_J[i] > expan.date & DailyA.L$Age[i] == "MATURE"){
    SSL.pre <- rbind(SSL.pre, DailyA.L[i,])
  }
}
SSL.pre$Color[SSL.pre$Age== "MATURE"] <- "black"
SSL.pre$Color[SSL.pre$Age=="YOUNG"] <- "blue"

# Loop that fits 4 parameter polynomial curves to DailyA
# and then creates a new data frame with daily predicted DailyA values
# blue open circles represent young leaves and black represent mature leaves
par(mfrow=c(4,6))
for (i in unique(SSL.pre$Plant_ID)){
    y.data <- SSL.pre$Aday[SSL.pre$Plant_ID == i]
    x.data <- SSL.pre$Date_J[SSL.pre$Plant_ID == i]
    colorchoice <- SSL.pre$Color[SSL.pre$Plant_ID == i]
    plot(y.data ~ x.data, col=colorchoice, xlab = "Date", ylab = "Aday", main = i)
    Light.Mod <- lm(y.data ~ poly(x.data, 5))
    predict.Aday <- predict(Light.Mod,data.frame(x=x.data))
    lines(x.data, predict.Aday, col='green',lwd=2)
    
    Date_J.temp <- 100
    while(Date_J.temp <= 344){
      Temp$Plant_ID <- i
      Temp$Date_J <- Date_J.temp
      if(Date_J.temp <= Individuals$Full.Expan[Individuals$Plant_ID==i]){
        Temp$Age <- "YOUNG"
      } else if (Date_J.temp > Individuals$Full.Expan[Individuals$Plant_ID==i]){
        Temp$Age <- "MATURE"
      }
      A <- predict(Light.Mod,data.frame(x.data=Date_J.temp))
      A <- as.numeric(A)
      Temp$Aday.lf <- A
      SeasonSeriesA <- rbind(SeasonSeriesA, Temp)
      Date_J.temp <- Date_J.temp + 1
    }
}
SeasonSeriesA$Aday.lf <- as.numeric(SeasonSeriesA$Aday.lf)
SeasonSeriesA$Date_J <- as.numeric(SeasonSeriesA$Date_J)

# ####
# ================== CO2 ASSIMILATION ACCROSS THE SEASON ==============================.====
#(27) DETERMINE DAILY Aday BASED ON SEASONAL LEAF AREA AND DAYLENGTH-------------------.----
#         Add Aday based on leaf age and seasonal light ####
SeasonSeries <- SeasonSeriesArea
SeasonSeries$Aday <- ""
SeasonSeries <- subset.data.frame(SeasonSeries, Date_J >= 100 & Date_J <= 344)
for(i in 1:nrow(SeasonSeries)){
  #Assign daily leaf assimilation
  SeasonSeries$Aday[i] <- SeasonSeriesA$Aday.lf[SeasonSeriesA$Date_J == SeasonSeries$Date_J[i] & 
                                                      SeasonSeriesA$Plant_ID == SeasonSeries$Plant_ID[i]]
}
SeasonSeries$Aday <- as.numeric(SeasonSeries$Aday)

#Remove Aday values below zero
SeasonSeries <- subset.data.frame(SeasonSeries, SeasonSeries$Aday >=0 )

#         Calulate daily A with various parameters ####
# Daily assimilation ignoring leaf area (mol m-2 day-1)
SeasonSeries$Aleaf <- ((SeasonSeries$Area.Perc/100) * SeasonSeries$Aday)

# Daily assimilion for whole plant (mol plant-1 day-1)
SeasonSeries$Aplant <- (SeasonSeries$Total_leaf_area * SeasonSeries$Aleaf)

# Daily assimilion for one stem (mol stem-1 day-1)
SeasonSeries$Astem <- (SeasonSeries$Aplant/SeasonSeries$Stem_count)

#(28) CALCULATE TOTAL C ASSIMILATED PER SEASON (AUC) ----------------------------------.----
#         Denote Season ####
SeasonSeries$Season <- ""
for(i in SeasonSeries$Date_J){
  if(i <= 161){
    SeasonSeries$Season[SeasonSeries$Date_J == i] <- "SPR"
  } else if (i > 161 & i <= 269){
    SeasonSeries$Season[SeasonSeries$Date_J == i] <- "SUM"
  }else if (i > 269){
    SeasonSeries$Season[SeasonSeries$Date_J == i] <- "FALL"
  }
}

#         Spring ####
# Plot all curves by individual
par(mfrow=(c(3,4)))
for (i in unique(SeasonSeries$Plant_ID[SeasonSeries$Season == "SPR"])){
  plot(SeasonSeries$Date_J[SeasonSeries$Plant_ID == i & SeasonSeries$Season == "SPR"],
       SeasonSeries$Aplant[SeasonSeries$Plant_ID == i &SeasonSeries$Season == "SPR"],
       main = i,
       xlab = "Date", ylab = "A plant (mmol C day-1)",
       type = "l", col="blue")
}

# Calculate area under the curve (Total A.spr)
Spring.Series <- subset.data.frame(SeasonSeries,SeasonSeries$Season == "SPR")
C.Tot.Spr <- Spring.Series %>%
  group_by(Plant_ID) %>%
  summarize(Aplant.spr = AUC(Date_J, Aplant, method = "spline"),
            Aleaf.spr = AUC(Date_J, Aleaf, method = "spline"),
            Astem.spr = AUC(Date_J, Astem, method = "spline"))
# Merge with individuals
Individuals <- merge.data.frame(Individuals, C.Tot.Spr, by = "Plant_ID", all = TRUE)

#         Summer ####
# Plot all curves by individual
par(mfrow=(c(4,5)))
for (i in unique(SeasonSeries$Plant_ID[SeasonSeries$Season == "SUM"])){
  plot(SeasonSeries$Date_J[SeasonSeries$Plant_ID == i & SeasonSeries$Season == "SUM"],
       SeasonSeries$Aplant[SeasonSeries$Plant_ID == i &SeasonSeries$Season == "SUM"],
       main = i,
       xlab = "Date", ylab = "A plant (mmol C day-1)",
       type = "l", col="blue")
}

# Calculate area under the curve (Total A.sm)
Summer.Series <- subset.data.frame(SeasonSeries,SeasonSeries$Season == "SUM")
C.Tot.Sm <- Summer.Series %>%
  group_by(Plant_ID) %>%
  summarize(Aplant.sm = AUC(Date_J, Aplant, method = "spline"),
            Aleaf.sm = AUC(Date_J, Aleaf, method = "spline"),
            Astem.sm = AUC(Date_J, Astem, method = "spline"))
# Merge with individuals
Individuals <- merge.data.frame(Individuals, C.Tot.Sm, by = "Plant_ID", all = TRUE)

#         Fall ####
# Plot all curves by individual
for (i in unique(SeasonSeries$Plant_ID[SeasonSeries$Season == "FALL"])){
  plot(SeasonSeries$Date_J[SeasonSeries$Plant_ID == i & SeasonSeries$Season == "FALL"],
       SeasonSeries$Aplant[SeasonSeries$Plant_ID == i &SeasonSeries$Season == "FALL"],
       main = i,
       xlab = "Date", ylab = "A plant (mmol C day-1)",
       type = "l", col="blue")
}

# Calculate area under the curve (Total A.fall)
Fall.Series <- subset.data.frame(SeasonSeries,SeasonSeries$Season == "FALL")
C.Tot.Fall <- Fall.Series %>%
  group_by(Plant_ID) %>%
  summarize(Aplant.fall = AUC(Date_J, Aplant, method = "spline"),
            Aleaf.fall = AUC(Date_J, Aleaf, method = "spline"),
            Astem.fall = AUC(Date_J, Astem, method = "spline"))

# Merge with individuals
Individuals <- merge.data.frame(Individuals, C.Tot.Fall, by = "Plant_ID", all = TRUE)

#Some Plants senesced before canopy closed
Individuals$Aplant.fall[Individuals$Plant_ID == 57 | Individuals$Plant_ID == 58 | Individuals$Plant_ID == 83] <- 0
Individuals$Aleaf.fall[Individuals$Plant_ID == 57 | Individuals$Plant_ID == 58 | Individuals$Plant_ID == 83] <- 0
Individuals$Astem.fall[Individuals$Plant_ID == 57 | Individuals$Plant_ID == 58 | Individuals$Plant_ID == 83] <- 0

#         Total ####
#mols of C ssimilated in entire season
Individuals$Aplant.Tot <- Individuals$Aplant.spr + Individuals$Aplant.sm + Individuals$Aplant.fall

#convert to kg
Individuals$Aplant.Tot.kg <- Individuals$Aplant.Tot * 44.01 / 1000

#Percents
Individuals$Aplant.spr.perc <- Individuals$Aplant.spr/Individuals$Aplant.Tot * 100
Individuals$Aplant.sm.perc <- Individuals$Aplant.sm/Individuals$Aplant.Tot * 100
Individuals$Aplant.fall.perc <- Individuals$Aplant.fall/Individuals$Aplant.Tot * 100

Individuals$Aplant.Tot.perc <- Individuals$Aplant.spr.perc + Individuals$Aplant.sm.perc + Individuals$Aplant.fall.perc 

#(29) CLEAN UP FINAL DATA SETS --------------------------------------------------------.----
#         SeasonSeries ####
#Add Date
SeasonSeries$Date <- as.Date(SeasonSeries$Date_J, origin=as.Date("2017-01-01"))
SeasonSeries.clean <- subset.data.frame(SeasonSeries, select = c(Plant_ID,Species,Status,Date_J,Date,
                                                                 Leaf.Age,Season,Stem_count,Leaf_count,
                                                                 Ave_leaf_area,Total_leaf_area,Aday,Aleaf,Aplant,Astem))
SeasonSeries.clean <- SeasonSeries.clean[order(SeasonSeries.clean$Plant_ID, SeasonSeries.clean$Date_J),]

# ####
#(30) NOTES ---------------------------------------------------------------------------.----
# "Individuals" dataframe has summary time points and the parameters used in the model
# "SeasonSeries" and ""SeasonSeries.clean" has all data over time

#------------------------ END OF CARBON ASSIMILATION MODEL-----------------------------.----
# ============================== RESULTS FIGURES ======================================.====
#(31) INDIVIDUALS TIME SERIES  --------------------------------------------------------.----
#       Aplant ####
ggplot(SeasonSeries.clean, aes(Date_J, Aplant, colour = Status, group = Plant_ID)) +
  geom_line(size = .2) +
  labs(x = "Julian Day", y = "A plant (mol C plant-1 day-1)") + 
  geom_dl(aes(label = Plant_ID), method = list(dl.combine("first.points", "last.points"), cex = 0.8)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#454f11"),labels = c("Invasive", "Native")) +
  theme_classic() + xlim(100,400)
  theme(legend.position = c(0.15,.9),text = element_text(size=14), panel.border = element_rect(colour = "black", fill=NA, size=.5))

par(mfrow=c(4,5))
for(i in unique(SeasonSeries.clean$Plant_ID)){
  plot(SeasonSeries.clean$Date_J[SeasonSeries.clean$Plant_ID == i],
       SeasonSeries.clean$Aplant[SeasonSeries.clean$Plant_ID == i],
       main = i, xlab = "Date", ylab = "A plant (mmol C day-1)",
       type = "l", col="blue")
}

#       Astem ####
SeasonSeries.clean.89 <- subset.data.frame(SeasonSeries.clean, SeasonSeries.clean$Plant_ID != 89)
ggplot(SeasonSeries.clean, aes(Date_J, Astem, colour = Status, group = Plant_ID)) +
  geom_line(size = 1) +
  labs(x = "Julian Day", y = "A stem (mol C stem-1 day-1)") +
  geom_dl(aes(label = Plant_ID), method = list(dl.combine("first.points", "last.points"), cex = 0.8)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#454f11"),labels = c("Invasive", "Native")) +
  theme_classic() + xlim(100,400) +
  theme(legend.position = c(0.15,.9),text = element_text(size=14), panel.border = element_rect(colour = "black", fill=NA, size=.5)) 

#       Aleaf ####
ggplot(SeasonSeries.clean, aes(Date_J, Aleaf, colour = Status, group = Plant_ID)) +
  geom_line(size = 1) +
  labs(x = "Julian Day", y = "A (mol C leaf-1 day-1)") +
  geom_dl(aes(label = Plant_ID), method = list(dl.combine("first.points", "last.points"), cex = 0.8)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#454f11"),labels = c("Invasive", "Native")) +
  theme_classic() + xlim(100,400) +
  theme(legend.position = c(0.15,.9),text = element_text(size=14), panel.border = element_rect(colour = "black", fill=NA, size=.5)) 

#       Aday - does not account for leaf area ####
ggplot(SeasonSeries.clean, aes(Date_J, Aday, colour = Status, group = Plant_ID)) +
  geom_line(size = 1) +
  labs(x = "Julian Day", y = "A day (mmol C m-2 day-1)") +
  geom_dl(aes(label = Plant_ID), method = list(dl.combine("first.points", "last.points"), cex = 0.8)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#454f11"),labels = c("Invasive", "Native")) +
  theme_classic() + xlim(100,400) +
  theme(legend.position = c(0.15,.9),text = element_text(size=14),panel.border = element_rect(colour = "black", fill=NA, size=.5)) 

#(32) SPECIES TIME SERIES  ------------------------------------------------------------.----
#       Species pivot table ####
SeasonSeries.Sum <- SeasonSeries.clean %>%
  group_by(Species, Date_J) %>%
  summarize(Aday = mean(Aday, na.rm = TRUE),
            Aplant = mean(Aplant, na.rm = TRUE),
            Aleaf = mean(Aleaf, na.rm = TRUE), 
            Astem = mean(Astem, na.rm = TRUE))

# Add status
SeasonSeries.Sum$Status <- ""
for(i in SeasonSeries.Sum$Species) {
  if(i == "BETH"|i == "LOBE"|i == "FRAL"|i == "RHCA"){
    SeasonSeries.Sum$Status[SeasonSeries.Sum$Species == i] <- "Invasive"
  }else if(i == "RUID"|i == "COCO"|i == "COSE"|i == "VILE"){
    SeasonSeries.Sum$Status[SeasonSeries.Sum$Species == i] <- "Native"
  } else
    return("AN ERROR MESSEGE")
}

# Add Genus
SeasonSeries.Sum$Genus <- ""
for(i in SeasonSeries.Sum$Species) {
  if(i == "BETH"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Berberis"
  }else if(i == "COCO"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Corylus"
  }else if(i == "COSE"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Cornus"
  }else if(i == "FRAL"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Frangula"
  }else if(i == "LOBE"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Lonicera"
  }else if(i == "RHCA"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Rhamnus"
  }else if(i == "RUID"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Rubus"
  }else if(i == "VILE"){
    SeasonSeries.Sum$Genus[SeasonSeries.Sum$Species == i] <- "Viburnum"
  } else
    return("AN ERROR MESSEGE")
}

#Add Date
SeasonSeries.Sum$Date <- as.Date(SeasonSeries.Sum$Date_J, origin=as.Date("2017-01-01"))

#   ****Aplant (MAIN SUMMARY FIGURE)**** ####
#Add overstory canopy context
Can.Light <-data.frame(xmin=as.Date(c(100,161,269), origin = as.Date("2017-01-01")), 
                       xmax=as.Date(c(161,269,340), origin = as.Date("2017-01-01")),
                       Can.Open=c("Open", "Closed", "Open"))
ggplot() +
  geom_rect(data = Can.Light, aes(xmin = xmin, ymin=-Inf, xmax = xmax, ymax=Inf, fill = Can.Open)) +
  geom_line(data = SeasonSeries.Sum, aes(Date, Aplant, group = Genus, colour = Status, linetype = Genus), size = 1) +
  #geom_point(data = SeasonSeries.Sum, aes(Date, Aplant, shape = Genus, group = Genus, colour = Status), size = 3) +
  labs(x = NULL, y=bquote(A['plant'] ~'(mol' ~CO[2] ~day^-1*')')) + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b', lim = c(as.Date("2017-04-10"),as.Date("2017-12-10"))) +
  scale_linetype_manual(name = "Genus", values = c("Berberis"="solid", "Frangula"="longdash","Lonicera"="dotted","Rhamnus"="dotdash",
                                                  "Corylus"="dotted","Cornus"="dotdash","Rubus"="solid","Viburnum"="longdash"),
                     guide=guide_legend(label.theme = element_text(angle = 0, face = "italic"))) +
  scale_color_manual(name = "Genus",  values = c("Invasive"="#79964c", "Native"="#133a00")) +
  scale_fill_manual(name = "Canopy", values = c("Open" = "white", "Closed" = "#e8e8e8")) +
  theme_classic() +
  theme(legend.position = "right",text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        plot.margin=unit(c(0,.05,.05,.05), "null"))

#       Aleaf ####
# Photosynethsis rates on m-2 basis, taking into account % changes in leaf area, but not the area itself (excludes leaf number and leaf size)
ggplot() +
  geom_rect(data = Can.Light, aes(xmin = xmin, ymin=-Inf, xmax = xmax, ymax=Inf, fill = Can.Open)) +
  geom_line(data = SeasonSeries.Sum, aes(Date, Aleaf, group = Species, colour = Status, linetype = Genus), size = 1) +
  labs(x = NULL, y=bquote(A['leaf'] ~'(mol' ~CO[2] ~day^-1*')')) + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b', lim = c(as.Date("2017-04-10"),as.Date("2017-12-10"))) +
  scale_linetype_manual(name = "Genus", values = c("Berberis"="solid", "Frangula"="longdash","Lonicera"="dotted","Rhamnus"="dotdash",
                                                   "Corylus"="dotted","Cornus"="dotdash","Rubus"="solid","Viburnum"="longdash"),
                        guide=guide_legend(label.theme = element_text(angle = 0, face = "italic"))) +
  scale_color_manual(name = "Genus",  values = c("Invasive"="#79964c", "Native"="#133a00")) +
  scale_fill_manual(name = "Canopy", values = c("Open" = "white", "Closed" = "#e8e8e8")) +
  theme_classic() +
  theme(legend.position = "right",text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        plot.margin=unit(c(0,.05,.05,.05), "null"))

#       Astem ####
ggplot() +
  geom_rect(data = Can.Light, aes(xmin = xmin, ymin=-Inf, xmax = xmax, ymax=Inf, fill = Can.Open)) +
  geom_line(data = SeasonSeries.Sum, aes(Date, Astem, group = Species, colour = Status, linetype = Genus), size = 1) +
  labs(x = NULL, y=bquote(A['stem'] ~'(mol' ~CO[2] ~day^-1*')')) + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b', lim = c(as.Date("2017-04-10"),as.Date("2017-12-10"))) +
  scale_linetype_manual(name = "Genus", values = c("Berberis"="solid", "Frangula"="longdash","Lonicera"="dotted","Rhamnus"="dotdash",
                                                 "Corylus"="dotted","Cornus"="dotdash","Rubus"="solid","Viburnum"="longdash"),
                   guide=guide_legend(label.theme = element_text(angle = 0, face = "italic"))) +
  scale_color_manual(name = "Genus",  values = c("Invasive"="#79964c", "Native"="#133a00")) +
  scale_fill_manual(name = "Canopy", values = c("Open" = "white", "Closed" = "#e8e8e8")) +
  theme_classic() +
  theme(legend.position = "right",text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        plot.margin=unit(c(0,.05,.05,.05), "null"))

#       Aday - does not account for leaf area ####
ggplot(SeasonSeries.Sum, aes(Date, Aday, colour = Status, group = Species)) +
  geom_line(size = 1) +
  labs(x = NULL, y = "A day (mol C m-2 day-1)") + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b', lim = c(as.Date("2017-04-10"),as.Date("2017-12-10"))) +
  geom_dl(aes(label = Species), method = list(dl.combine("first.points", "last.points"), cex = 0.8)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#454f11"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(legend.position = c(0.15,.9),text = element_text(size=14), panel.border = element_rect(colour = "black", fill=NA, size=.5))

ggplot() +
  geom_rect(data = Can.Light, aes(xmin = xmin, ymin=-Inf, xmax = xmax, ymax=Inf, fill = Can.Open)) +
  geom_line(data = SeasonSeries.Sum, aes(Date, Aday, group = Species, colour = Status, linetype = Genus), size = 1) +
  labs(x = NULL, y=bquote(A['day'] ~'(mol' ~CO[2] ~day^-1* ~m^-2*')')) + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b', lim = c(as.Date("2017-04-10"),as.Date("2017-12-10"))) +
  scale_linetype_manual(name = "Genus", values = c("Berberis"="solid", "Frangula"="longdash","Lonicera"="dotted","Rhamnus"="dotdash",
                                                   "Corylus"="dotted","Cornus"="dotdash","Rubus"="solid","Viburnum"="longdash"),
                        guide=guide_legend(label.theme = element_text(angle = 0, face = "italic"))) +
  scale_color_manual(name = "Genus",  values = c("Invasive"="#79964c", "Native"="#133a00")) +
  scale_fill_manual(name = "Canopy", values = c("Open" = "white", "Closed" = "#e8e8e8")) +
  theme_classic() +
  theme(legend.position = "right",text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        plot.margin=unit(c(0,.05,.05,.05), "null"))

#(33) TOTAL CARBON - SEASON COMPARISONS -----------------------------------------------.----
#       Total ####
C.tot.fig.sp <- ggplot(Individuals,aes(Species, Aplant.Tot.kg, fill = Status, color = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = bquote("Total Carbon Gain (kg "~CO[2] ~yr^-1*')'), x = NULL) +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() + ylim(0,400) +
  scale_color_manual(name = NULL,  guide=FALSE, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  scale_fill_manual(name = NULL, guide = FALSE, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  annotate("text", x = c(1,2,3,4,5,6,7,8), y = c(60, 390, 180, 120, 50, 320, 300, 130), 
           label = c("C", "A", "A", "AB", "BC", "A", "A", "A"), color = "blue", size = 5) +
  theme(legend.position = c(0.15,.85),text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x = element_text(angle = 45, face = "italic", hjust=1),
        plot.margin=unit(c(.05,0,.05,.05), "null"))

C.tot.fig.st <- ggplot(Individuals,aes(Status, Aplant.Tot.kg, fill=Status, color = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = .75)+
  geom_boxplot(lwd = .75) +
  labs(y = NULL, x = NULL) +
  theme_classic() + ylim(0,400) +
  scale_fill_manual(guide = FALSE, name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL,  guide=FALSE, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  geom_signif(y_position=380, xmin=1, xmax=2, annotation=c("NS"), tip_length=0, colour = "blue", textsize = 5)+
  theme(legend.position = c(0.5,.9),text = element_text(size=18),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x = element_text(angle = 45, hjust=1),
        plot.margin=unit(c(.05,.05,.08,0), "null")) 
  
grid.arrange(C.tot.fig.sp, C.tot.fig.st, widths = c(3,1), ncol = 2, nrow = 1)

#       Spring ####
ggplot(Individuals,aes(Species, Aplant.spr, fill = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25)+
  geom_boxplot() +
  labs(y = "Total C assimilated in Spring (mol)", x = NULL) +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#       Summer ####
ggplot(Individuals,aes(Species, Aplant.sm, fill = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25)+
  geom_boxplot() +
  labs(y = "Total C assimilated in Summer (mol)", x = NULL) +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#       Fall ####
ggplot(Individuals,aes(Species, Aplant.fall, fill = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25)+
  geom_boxplot() +
  labs(y = "Total C assimilated in Fall (mol)", x = NULL) +
  scale_x_discrete(name = NULL, limits=c("BETH","LOBE","RHCA","FRAL","RUID","COCO", "COSE", "VILE"),
                   label=c("BETH"="Berberis", "COCO"="Corylus", "COSE"="Cornus", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  theme(legend.position = c(0.15,.8),text = element_text(size=18), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#(34) TOTAL CARBON - SEASON PERCENTAGES -----------------------------------------------.----
#       Stacked Box plot ####
# Summarize Percents
Species.Summary <- Individuals %>%
  group_by(Species) %>%
  summarize(A.spr.perc = mean(Aplant.spr.perc, na.rm = TRUE),
            A.sum.perc = mean(Aplant.sm.perc, na.rm = TRUE),
            A.fall.perc = mean (Aplant.fall.perc, na.rm = TRUE))

# Stack Percents
Percentage <- gather(Species.Summary, Season, Percent,
                        A.spr.perc, A.sum.perc, A.fall.perc)

# Subset by status
Percentage.N <- subset.data.frame(Percentage, Species == "COCO"|Species == "COSE"|Species == "RUID"|Species == "VILE")
Percentage.I <- subset.data.frame(Percentage, Species == "BETH"|Species == "LOBE"|Species == "FRAL"|Species == "RHCA")

# Invasive
A.perc.I <- ggplot(Percentage.I, aes(x = Species, y = Percent, fill=Season)) +
  geom_bar(stat='identity') +
  labs(x="Invasive", y="Total Carbon Gain (%)") +
  scale_x_discrete(name = "Invasive", limits=c("BETH","LOBE","RHCA","FRAL"),
                   label=c("BETH"="Berberis", "FRAL"="Frangula", "LOBE"="Lonicera", "RHCA"="Rhamnus")) +
  theme_classic() +
  scale_fill_manual(name = "Season", guide = FALSE, values = c("#e8c16d", "#CADD6D", "#5F6D19"),labels = c("Fall","Spring", "Summer")) +
  theme(legend.position = "right", text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x = element_text(face = "italic"),
        plot.margin=unit(c(.05,0,.05,.05), "null")) 

# Native
A.perc.N <- ggplot(Percentage.N, aes(x = Species, y = Percent, fill=Season)) +
  geom_bar(stat='identity') +
  labs(x="Native", y=NULL) +
  scale_x_discrete(name = "Native", limits=c("RUID","COCO", "COSE", "VILE"),
                   label=c("COCO"="Corylus", "COSE"="Cornus", "RUID"="Rubus", "VILE"="Viburnum")) +
  theme_classic() +
  scale_fill_manual(name = "Season", values = c("#e8c16d", "#CADD6D", "#5F6D19"),labels = c("Fall","Spring", "Summer")) +
  theme(legend.position = "right", text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x = element_text(face = "italic"),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        plot.margin=unit(c(.05,.05,.05,0), "null")) 

grid.arrange(A.perc.I, A.perc.N, widths = c(1,1.1), ncol = 2, nrow = 1)

#       Percents by status ####
#Spring
Spring.perc.st.fig <- ggplot(Individuals,aes(Status, Aplant.spr.perc, fill=Status, color = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = 0.75)+
  geom_boxplot(lwd = 0.75) +
  labs(y = "Carbon Gain (%)", title = "Spring", x = NULL)+
  theme_classic() +  ylim(0,105) +
  scale_fill_manual(guide = FALSE, name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(guide = FALSE, name = NULL, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.5,.9),text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        plot.margin=unit(c(.05,0,.05,.15), "null"),
        plot.title = element_text(hjust = 0.5))+
  geom_signif(y_position=100, xmin=1, xmax=2, annotation=c("***"), tip_length=0, colour = "blue", textsize = 8)

#Summer
Sum.perc.st.fig <- ggplot(Individuals,aes(Status, Aplant.sm.perc, fill=Status, color = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = 0.75)+
  geom_boxplot(lwd = 0.75) +
  labs(y = NULL, title = "Summer", x = NULL) +
  theme_classic() +  ylim(0,105) +
  scale_fill_manual(guide = FALSE, name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(guide = FALSE, name = NULL, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = c(0.5,.9),text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        plot.margin=unit(c(.05,0,.05,0), "null"),
        plot.title = element_text(hjust = 0.5),
        panel.background = element_rect(fill = "#e8e8e8")) +
  geom_signif(y_position=100, xmin=1, xmax=2, annotation=c("*"), tip_length=0, colour = "blue", textsize = 8)

#Fall
Fall.perc.st.fig <- ggplot(Individuals,aes(Status, Aplant.fall.perc, fill=Status, color = Status))+
  stat_boxplot(geom ='errorbar', width = 0.25, lwd = 0.75)+
  geom_boxplot(lwd = 0.75) +
  labs(y = NULL, title = "Fall", x = NULL) +
  theme_classic() + ylim(0,105) +
  scale_fill_manual(name = NULL, values = c("#CADD6D", "#5F6D19"),labels = c("Invasive", "Native")) +
  scale_color_manual(name = NULL, values = c("Invasive" = "#79964c", "Native" = "#133a00")) +
  theme(legend.position = "right",text = element_text(size=18), 
        panel.border = element_rect(colour = "black", fill=NA, size=.5),
        axis.text.x=element_blank(), axis.ticks.x=element_blank(),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        plot.margin=unit(c(.05,.05,.05,0), "null"),
        plot.title = element_text(hjust = 0.5)) +
  geom_signif(y_position=100, xmin=1, xmax=2, annotation=c("**"), tip_length=0, colour = "blue", textsize = 8)

#Assemble
grid.arrange(Spring.perc.st.fig, Sum.perc.st.fig, Fall.perc.st.fig, 
             widths = c(1.5, 1.8, 1.9), ncol = 3, nrow = 1)
# ####
# ================== EXPLORATION OF RESULTS: PREDICTOR VARIABLES ======================.====
# Check A Total distribution ####
par(mfrow=c(1,2))
hist(Individuals$Aplant.Tot, xlab = "Total Carbon Assimilated (mol)", main = "Total A")
hist(log(Individuals$Aplant.Tot, base = 10), xlab = "log Total Carbon Assimilated (mol)", main = "log Total A")

#(35) PLANT HEIGHT --------------------------------------------------------------------.----
#       Set-up ####
Species.Summary.new <- Individuals %>%
  group_by(Species) %>%
  summarise(Height = mean(Plant_height, na.rm = TRUE),
            Aplant.Tot = mean(Aplant.Tot, na.rm = TRUE))

status <- c()
for(i in Species.Summary.new$Species) {
  if(i == "BETH"|i == "LOBE"|i == "FRAL"|i == "RHCA"){
    status <- c(status, "Invasive")
  }else if(i == "RUID"|i == "COCO"|i == "COSE"|i == "VILE"){
    status <- c(status, "Native")
  } else
    return("AN ERROR MESSEGE")
}
Species.Summary.new <- cbind(Species.Summary.new, status)

#       All Individuals ####
ggplot(Individuals,aes(Plant_height, Aplant.Tot)) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "Total Carbon (mol)", x = "Plant Height", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#       Species ####
ggplot(Species.Summary.new,aes(Height, Aplant.Tot)) +
  geom_point(aes(colour = status), size=6) +
  stat_quantile(method = "lm", formula = y ~ x, size = 2) +
  labs(y = "Total Carbon (mol)", x = "Plant Height (cm)", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14), panel.border = element_rect(colour = "black", fill=NA, size=.5))

#(36) LEAF AREA -----------------------------------------------------------------------.----
#       Leaf Count ####
par(mfrow=c(1,2))
hist(Individuals$Leaf_count, xlab = "# of Leaves", main = "Leaf Count")
hist(log(Individuals$Leaf_count, base = 10), xlab = "log # of Leaves", main = "log Leaf Count")

ggplot(Individuals,aes(Leaf_count, Aplant.Tot)) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "Total Carbon Gain (mol)", x = "Leaf Count", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

ggplot(Individuals,aes(log(Leaf_count, base = 10), log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "log Leaf Count", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#       Leaf Size ####
par(mfrow=c(1,2))
hist(Individuals$Ave_leaf_area, xlab = "Leaf Size (cm2)", main = "Leaf Size")
hist(sqrt(Individuals$Ave_leaf_area), xlab = "sqrt Leaf Size (cm2)", main = "sqrt Leaf Size")

ggplot(Individuals, aes(Ave_leaf_area, Aplant.Tot)) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "Total Carbon Gain (mol)", x = "Leaf Area", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

ggplot(Individuals, aes(sqrt(Ave_leaf_area), log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "log Leaf Area", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#       Total Area ####
par(mfrow=c(1,2))
hist(Individuals$Total_leaf_area, xlab = "Total Area (m2)", main = "Total Leaf Area")
hist(log(Individuals$Total_leaf_area, base = 10), xlab = "sqrt Total Area (m2)", main = "sqrt Total Leaf Area")

ggplot(Individuals, aes(Total_leaf_area, Aplant.Tot)) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "Total Carbon Gain (mol)", x = "Total Leaf Area", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

ggplot(Individuals, aes(log(Total_leaf_area, base = 10), log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "log Total Leaf Area", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

par(mfrow=c(1,1))
plot(log(Individuals$Aplant.Tot, base = 10) ~ log(Individuals$Total_leaf_area, base = 10),
     xlab = "log Total Leaf Area (m2)", ylab = "log Total C Gain (mol)", main = "Total Leaf Area")

#(37) PHOTOSYNTHESIS RATES ------------------------------------------------------------.----
#       Check Distributions ####
par(mfrow=c(1,2))
hist(Individuals$Amax.JUN, xlab = "Amax June (mmol C s-1 m-2)", main = "Amax JUNE")
hist(Individuals$Amax.AUG, xlab = "Amax August (mmol C s-1 m-2)", main = "Amax AUGUST")

plot(log(Individuals$Aplant.Tot, base = 10) ~ Individuals$Amax.JUN,
     xlab = "Amax June (mmol C s-1 m-2)", ylab ="log Total Carbon Gain (mol)", main = "Young Amax")
plot(log(Individuals$Aplant.Tot, base = 10) ~ Individuals$Amax.AUG,
     xlab = "Amax August (mmol C s-1 m-2)", ylab ="log Total Carbon Gain (mol)", main = "Mature Amax")

#       Young ####
ggplot(Individuals, aes(Amax.JUN, log(Aplant.Tot, base =10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "Amax June (mmol C s-1 m-2)", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#       Mature ####
ggplot(Individuals, aes(Amax.AUG, log(Aplant.Tot, base =10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "Amax August (mmol C s-1 m-2)", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#(38) UNDERSTORY LIGHT LEVELS (bins and % lux) ----------------------------------------.----
#       Check Distributions ####
# PPFD
par(mfrow=c(1,3))
hist(Individuals$PPFD.MAY, xlab = "Daily PPFD Spring", main = "Spring")
hist(log(Individuals$PPFD.JUL, base = 10), xlab = "Daily PPFD Summer", main = "Summer")
hist(Individuals$PPFD.JUL, xlab = "Daily PPFD Fall", main = "Fall")

plot(log(Individuals$Aplant.Tot, base = 10) ~ Individuals$PPFD.MAY,
     xlab = "Daily PPFD Spring", ylab ="log Total Carbon Gain (mol)", main = "Spring")
plot(log(Individuals$Aplant.Tot, base = 10) ~ log(Individuals$PPFD.JUL, base = 10),
     xlab = "log Daily PPFD Summer", ylab ="log Total Carbon Gain (mol)", main = "Summer")
plot(log(Individuals$Aplant.Tot, base = 10) ~ Individuals$PPFD.OCT,
     xlab = "Daily PPFD Fall", ylab ="log Total Carbon Gain (mol)", main = "Fall")

# % LUX
par(mfrow=c(1,4))
hist(Individuals$LUX_percent.spr, xlab = "% LUX Spring", main = "Spring")
hist(log(Individuals$LUX_percent.sumII, base = 10), xlab = "% LUX Summer", main = "Summer I")
hist(log(Individuals$LUX_percent.sumII, base = 10), xlab = "% LUX Summer", main = "Summer II")
hist(Individuals$PPFD_percent.fall, xlab = "Daily PPFD Fall", main = "Fall")

par(mfrow=c(1,4))
plot(log(Individuals$Aplant.Tot, base = 10) ~ Individuals$LUX_percent.spr,
     xlab = "% LUX Spring", ylab ="log Total Carbon Gain (mol)", main = "Spring % LUX")
plot(log(Individuals$Aplant.Tot, base = 10) ~ log(Individuals$LUX_percent.sumI, base = 10),
     xlab = "log % LUX Summer I", ylab ="log Total Carbon Gain (mol)", main = "Summer I % LUX")
plot(log(Individuals$Aplant.Tot, base = 10) ~ log(Individuals$LUX_percent.sumII, base = 10),
     xlab = "log % LUX Summer II", ylab ="log Total Carbon Gain (mol)", main = "Summer II % LUX")
plot(log(Individuals$Aplant.Tot, base = 10) ~ Individuals$LUX_percent.fall,
     xlab = "Daily PPFD Fall", ylab ="log Total Carbon Gain (mol)", main = "Fall")

#       Spring ####
#Bins
ggplot(Individuals, aes(PPFD.MAY, log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "Daily PPFD Spring", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

# % LUX
ggplot(Individuals, aes(LUX_percent.spr, log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "% LUX Spring", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#       Summer ####
#Bins
ggplot(Individuals, aes(log(PPFD.JUL, base = 10), log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "log Daily PPFD Summer", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

# % LUX
ggplot(Individuals, aes(log(LUX_percent.sumI, base = 10), log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "log % LUX Summer", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#       Fall ####
#Bins
ggplot(Individuals, aes(PPFD.OCT, log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "Daily PPFD Fall", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

# % LUX
ggplot(Individuals, aes(log(LUX_percent.fall, base = 10), log(Aplant.Tot, base = 10))) +
  geom_point(aes(shape = Species, colour = Status), size=4) +
  stat_quantile(aes(group = Species, colour = Status), method = "lm", formula = y ~ x) +
  stat_quantile(method = "lm", formula = y ~ x) +
  labs(y = "log Total Carbon Gain (mol)", x = "log % LUX Summer", shape = "Species") +
  scale_shape_manual(values = c(15,4,1,16,17,8,3,2)) +
  scale_colour_manual(name = NULL, values = c("#94a34e", "#333a0d"),labels = c("Invasive", "Native")) +
  theme_classic() +
  theme(text = element_text(size=14))

#(39) MULTIPLE LINEAR REGRESSION ------------------------------------------------------.----
#       Transform the variables ####
par(mfrow=c(1,1))
Individuals$Aplant.Tot.t <- log(Individuals$Aplant.Tot)
shapiro.test(Individuals$Aplant.Tot.t)
hist(Individuals$Aplant.Tot.t)

Individuals$Total_leaf_area.t <- log(Individuals$Total_leaf_area, base=10)
shapiro.test(Individuals$Total_leaf_area.t)
hist(Individuals$Total_leaf_area.t)

Individuals$Plant_height.t <- sqrt(Individuals$Plant_height)
shapiro.test(Individuals$Plant_height.t)
hist(Individuals$Plant_height.t)

Individuals$Stem_count.t <- nthroot(Individuals$Stem_count, 4)
shapiro.test(Individuals$Stem_count.t)
hist(Individuals$Stem_count.t)
  #Not transformable

shapiro.test(Individuals$Amax.JUN)
hist(Individuals$Amax.JUN)

Individuals$Amax.AUG.t <- log(Individuals$Amax.AUG)
shapiro.test(Individuals$Amax.AUG.t)
hist(Individuals$Amax.AUG.t)

Individuals$PPFD.MAY.t <- nthroot(Individuals$PPFD.MAY,2)
shapiro.test(Individuals$PPFD.MAY.t)
hist(Individuals$PPFD.MAY.t)
  #Not transpformable

Individuals$PPFD.JUL.t <- nthroot(Individuals$PPFD.JUL,4)
shapiro.test(Individuals$PPFD.JUL.t)
hist(Individuals$PPFD.JUL.t)
  #Not transpformable

Individuals$PPFD.SEP.t <- nthroot(Individuals$PPFD.SEP,4)
shapiro.test(Individuals$PPFD.SEP.t)
hist(Individuals$PPFD.SEP.t)
  #Not transpformable

#       Pick out important variables ####
mod.all <- lm(Individuals$Aplant.Tot.t ~ Individuals$Total_leaf_area.t + Individuals$Stem_count + Individuals$Plant_height +
                Individuals$Amax.JUN + Individuals$Amax.AUG.t + 
                Individuals$PPFD.MAY + Individuals$PPFD.JUL +Individuals$PPFD.SEP + Individuals$PPFD.OCT)
avPlots(mod.all)

#       Test for correlations with response ####
panel.cor <- function(x,y, digits=2, cex.cor) 
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0,1,0,1))
  r <- abs(cor(x,y, method = "pearson"))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  test <- cor.test(x,y)
  p <- test$p.value
  p <- format(p, digits=digits)[1]
  text(0.5, 0.5, paste("r=", txt, "\np=", p))
}

pairs(Individuals$Aplant.Tot.t ~ Individuals$Total_leaf_area.t + Individuals$Stem_count + Individuals$Plant_height +
        Individuals$Amax.JUN + Individuals$Amax.AUG.t + 
        Individuals$PPFD.MAY + Individuals$PPFD.JUL +Individuals$PPFD.SEP + Individuals$PPFD.OCT,
      lower.panel = panel.smooth, upper.panel = panel.cor)
      #Leaf area, stem count, and height corrleated with A.Tot
      #But Amax AUG is not !?

#       Test for correlations within predicors ####
#(can't have predictor variables correlated with each other)
pairs(~ Individuals$Total_leaf_area.t + Individuals$Stem_count + Individuals$Plant_height +
        Individuals$Amax.JUN + Individuals$Amax.AUG.t + 
        Individuals$PPFD.MAY + Individuals$PPFD.JUL +Individuals$PPFD.SEP + Individuals$PPFD.OCT,
      lower.panel = panel.smooth, upper.panel = panel.cor)
      #Amax in JUN and AUG are correlated
      #leaf area and stem count are correlated
      #leaf area and height are correlated
      #Amax JUN and PPFD JUl are correlated

#       VIF's to double check colinearity ####
#       Create models ####
# Model 1 - everything
mod1 <- lm(Individuals$Aplant.Tot.t ~ Individuals$Total_leaf_area.t + Individuals$Stem_count + Individuals$Plant_height.t +
             Individuals$Amax.JUN + Individuals$Amax.AUG.t + 
             Individuals$PPFD.MAY + Individuals$PPFD.JUL +Individuals$PPFD.SEP + Individuals$PPFD.OCT)
avPlots(mod1) 
par(mfrow = c(2,2))
plot(mod1)
modelEffectSizes(mod1)
summary(mod1)

# Model 2 - most significant of each caterogry
mod2 <- lm(Individuals$Aplant.Tot.t ~ Individuals$Total_leaf_area.t + Individuals$Amax.AUG.t + 
             Individuals$PPFD.OCT)
avPlots(mod2) 
par(mfrow = c(2,2))
plot(mod2)
modelEffectSizes(mod2)
summary(mod2)

# Model 3 - leaf.area and Amax
mod3 <- lm(Individuals$Aplant.Tot.t ~ Individuals$Total_leaf_area.t + Individuals$Amax.AUG.t)
avPlots(mod3) 
par(mfrow = c(2,2))
plot(mod3)
modelEffectSizes(mod3)
summary(mod3)

# Model 3 - only significantly correlated with predictor
mod4 <- lm(Individuals$Aplant.Tot.t ~ Individuals$Total_leaf_area.t + Individuals$Stem_count.t +
             Individuals$Plant_height.t)
avPlots(mod4) 
par(mfrow = c(2,2))
plot(mod4)
modelEffectSizes(mod4)
summary(mod4)

#       Comparing Models ####
AIC(mod1, mod2, mod3, mod4)