#########################################################
##how do species vary in their responses to fire? #######
##can we predict their responses using traits?###########
##do any show responses contingent on burn history?######
#########################################################
#scale_color_manual(values = c("skyblue4", "gold", "chocolate1", "brown", "gray28"))+
#scale_fill_manual(values = alpha(c("skyblue4", "gold", "chocolate1", "brown", "gray28"), 0.5))+
  
source("C:/Users/simha/Dropbox (Duke Bio_Ea)/Anita/Research_projects/R_code/ftbr_2020-03-18_read-in-data.R")
library(gridExtra)
library(interactions)
####################
#defining lambda####
####################
#take rows where there are burns
#add column of cover for previous year
data3 = data.table(gather(data1, species, cover, ACAGRA:ZIGGLA))
#data3[cover2 := lapply(data3[Year == Year-1 & UID == UID & species == species, cover])]
data4 = data3[!is.na(BurnDate),]
for (i in 1:dim(data4)[1]){
  yr = data4[i, Year-1]
  uid = data4[i, UID]
  sp = data4[i, species]
  if(!isEmpty(data3[Year == yr & UID == uid & species == sp, cover])){
    data4[i, cover2 := (data3[Year == yr & UID == uid & species == sp, cover])]
  }
}

#response ratio as change in cover over a year where there was fire
#not just where a species was present but where it was at least minorly abundant...
data5=data4[cover+cover2>5,] 

data5[,firedelt :=cover-cover2]
#anything less than 1 for cover/cover2 means negative change
#anything less than zero is negative change for log
data5[,firelam := log((cover+1)/(cover2+1))]
data5 = data5[firelam < Inf & firelam > -Inf,]
#data5[,numburns := as.integer(numburns)]

######################################################
########trait basis of lambda variation###############
######################################################
#take data with at least 3 points - 129 spp
tab = as.data.frame(table(data5$species))
data5_5 = data5[species %in% tab$Var1[tab$Freq>2]]
data6 = data5_5 %>% group_by(species) %>% summarize(lambda = mean(firelam), sdlam = sd(firelam), n = n())
#data6 = data5_5 %>% group_by(species) %>% summarize(lambda = mean(lam), sdlam = sd(lam), n = n())
#write.csv(data6, paste0(Sys.Date(), "_species-lambda.csv"))

###visualize species variation in log response-ratio
#png(paste0(Sys.Date(), "_ftbr-lambda-hist.png"), height = 4, width = 5, units = "in", res = 300)
a = ggplot(data = data6, aes(x = lambda))+ theme_classic()+
  theme(panel.border = element_rect(color = "black", fill = NA))+
  xlab("Mean LRR")+ylab("Frequency")+
  theme(plot.title = element_text(size = 16))+
  theme(axis.title = element_text(size = 14))+
  theme(axis.text = element_text(size = 12))+
  ggtitle("(A)")+
  geom_histogram(color = "black", fill = "gray60", binwidth = 0.2)
#dev.off()


traits = data.table(read.xlsx("CombinedTraits_SummaryStats.xlsx"))
data7 = merge(data6, traits, by.x = "species", by.y = "SpecCode", all.x = TRUE, all.y = FALSE)
species = as.data.table(read.xlsx("2019-09-13_ftbr_species_cleaned.xlsx"))
cols = c("species_code", "woodiness", "growth_form")
species = species[, ..cols]
data7 = merge(data7, species, by.x = "species", by.y = "species_code", all.x = T, all.y = F)
setDT(data7)
cols = c("species", "lambda", "sdlam", "n", "woodiness", "growth_form", "Avg_LDMC", "Avg_LeafArea",
               "Avg_SLA", "Avg_height_cm", "SD_LDMC", "SD_SLA", "SD_height_cm", "SD_LeafArea")
data7 = data7[, ..cols]
data7 = na.omit(data7)
#add CV variables
data7[,CV_SLA := SD_SLA/Avg_SLA]
data7[,CV_LDMC := SD_LDMC/Avg_LDMC]
data7$woodiness = as.factor(data7$woodiness)
#t.test(data7$Avg_height_cm[data7$woodiness == 0], data7$Avg_height_cm[data7$woodiness == 1]) #woody species are taller

data7[woodiness == 0, woodiness := "Not Woody"]
data7[woodiness == 1, woodiness := "Woody"]
t.test(data7$lambda[data7$woodiness == "Woody"], data7$lambda[data7$woodiness == "Not Woody"])

##visualize effect of woodiness on LRR
#png(paste0(Sys.Date(), "_ftbr-lambda-woodiness.png"), height = 6, width = 8, units = "in", res = 300)
b = ggplot(data = data7, aes(x = woodiness, y = lambda))+   theme_classic()+
  theme(panel.border = element_rect(color = "black", fill = NA))+
  xlab("Woodiness")+ylab("LRR")+
  theme(plot.title = element_text(size = 16))+
  theme(axis.title = element_text(size = 14))+
  theme(axis.text = element_text(size = 12))+
  #scale_color_manual(values = c("darkseagreen4", "chocolate4"))+
  geom_boxplot(width = 0.4)+ #position = position_nudge(x=-0.1))+
  geom_point(data = data7, aes(x = woodiness, y = lambda), position = position_jitter(height = 0, width =0.07), shape = 16, size = 1.8, alpha = 0.4)+
  #stat_summary(fun = "mean", na.rm = T, size = 1, geom="line", aes(group = species, color = woodiness))+
  theme(legend.position="none")+ggtitle("(B)")
#dev.off()

#i think SLA, height, and woodiness will be imp traits here
#Avg_SLA, CV_SLA, Avg_height, woodiness
mod1 = lm(lambda ~ woodiness*CV_SLA+Avg_height_cm, data = data7)
summary(mod1)
AIC(mod1) #111
mod2 = lm(lambda ~ woodiness*CV_SLA,  data = data7)
summary(mod2)
AIC(mod2) #111
mod3 = lm(lambda ~ woodiness*CV_SLA*Avg_height_cm,  data = data7)
summary(mod3)
AIC(mod3) #109
mod4 = lm(lambda ~ woodiness+CV_SLA*Avg_height_cm,  data = data7)
summary(mod4)
AIC(mod4) #104.9
mod5 = lm(lambda ~ woodiness*Avg_SLA*CV_SLA*Avg_height_cm,  data = data7)
summary(mod5)
AIC(mod5) #107
mod6 = lm(lambda ~ woodiness*Avg_SLA*CV_SLA,  data = data7)
summary(mod6)
AIC(mod6) #107

#interaction with height and CV_SLA is sig
car::Anova(mod4, type = 3)

#visualize interaction between CV_SLA and Avg height
c = interact_plot(mod4, CV_SLA, modx = Avg_height_cm, 
              plot.points = TRUE, point.alpha = 0.5, point.size = 1.5,
              interval = TRUE, legend.main = "Avg. Height", 
              x.label = "CV SLA", y.label = "LRR",
              colors = c("darkred", "gray70", "gray10"))+
  theme_classic()+ggtitle("(C)")+  theme(plot.title = element_text(size = 16))+
  theme(panel.border = element_rect(color = "black", fill = NA))+
  theme(axis.title = element_text(size = 14))+
  theme(legend.key.size = unit(0.5, "line"))+
  theme(legend.title = element_text(size = 12), 
        legend.text = element_text(size = 10))+
  theme(axis.text = element_text(size = 12)) +
  scale_x_continuous(expand = c(0,0))

####################################################
##which LRR slopes are statistically significant? ##
####################################################
tab = as.data.frame(table(data5$species))
tab = tab[order(-tab$Freq),]
#species with at least, like, 5 data points
spp = tab$Var1[tab$Freq > 4] #113 spp total
data5_1 = data5[species %in% spp, ]
data5_2 = data5_1 %>% group_by(species) %>% summarize(yint = lm(firelam ~ FRI)$coefficients[1],
                                                      slope = lm(firelam ~ FRI)$coefficients[2],
                                                      mean = mean(firelam),
                                                      pval = summary(lm(firelam ~ FRI))$coefficients[,4][2],
                                                      n = length(firelam))

setDT(data5_2) #58 species
#png(paste0(Sys.Date(), "_lambda-FRI-statsig-spp.png"), width = 6, height = 6, units="in", res=300)
#ggplot(data = data5[species %in% data5_2[pval<0.05, species]], aes(y = firelam, x = FRI))+
#  facet_wrap(~species, ncol = 3)+geom_hline(yintercept=0, color="gray40", linetype = "longdash")+
#  geom_point()+geom_smooth(method="lm")+theme_bw()+
#  xlab("FRI")+ylab("Lambda")
#dev.off()

d = ggplot(aes(x = FRI, y = firelam), data = data5_1[species == "ARISTR",])+
  theme_classic()+
  theme(panel.border = element_rect(color = "black", fill = NA))+
  theme(plot.title = element_text(size = 16))+
  theme(axis.title = element_text(size = 14))+
  theme(axis.text = element_text(size = 12))+
  geom_smooth(color = "gray30", method = "lm", na.rm = T, fullrange = T, size = 1.2)+
  geom_point(color = "gray30", alpha = 0.5, size = 1.5)+
  xlab("Fire Return Interval")+ylab("A. stricta LRR")+ggtitle("(D)")+
  scale_x_continuous(expand = c(0,0))

#lay <- rbind(c(1, 1, 1, NA, 2, 2,2),
#             c(3,3,3,3,4,4,4))
#png(paste0(Sys.Date(), "_Fig5.png"), width = 7, height = 5.5, units="in", res=300)
#grid.arrange(a, b, c, d, layout_matrix = lay)
#dev.off()
#####################################################
####generate species table of most abundant species##
#####################################################
#data5_FRI = setDT(data5_2[pval < 0.051])
#data5_FRI = dplyr::select(data5_FRI, species, slope)
#names(data5_FRI) = c("species", "FRI sensitivity")
#
#ab = data3 %>% group_by(species) %>% summarize(abundance = mean(cover[cover > 0], na.rm = T),
#                                               incidence = sum(cover>0, na.rm = T)/dim(data1)[1])
#
#traits = data.table(read.xlsx("CombinedTraits_SummaryStats.xlsx"))
#traits = dplyr::select(traits, "SpecCode", "Avg_SLA", "Avg_height_cm", "SD_SLA")
#
#species = as.data.table(read.xlsx("2019-09-13_ftbr_species_cleaned.xlsx"))
#species = dplyr::select(species, species_code, species_name, growth_form, woodiness)
#x = dplyr::select(data5_2, species, mean)
#names(x) = c("species", "Mean LRR")

#splist = merge(species, traits, by.y = "SpecCode", by.x = "species_code", all.x = T, all.y = F)
#splist = merge(ab, splist, by.x = "species", by.y = "species_code")
#setDT(splist)
#splist[,CV_SLA := SD_SLA/Avg_SLA]
#splist = merge(splist, data5_FRI, by="species", all = T)
#splist = merge(splist, x, by = "species", all.x = T, all.y = F)
#splist = dplyr::select(splist, species_name, growth_form, abundance, incidence, woodiness, Avg_height_cm, Avg_SLA, CV_SLA,
#                "FRI sensitivity", "Mean LRR")
#names(splist) = c("Species", "Growth Form", "Abundance", "Incidence", "Woodiness", "Average Height", "Average SLA", "CV SLA",
#                 "FRI sensitivity", "Mean LRR")

#splist = splist[order(-Abundance),]
#splist = splist[!is.na(splist$`Average Height`)]
#splist = splist[!is.na(splist$`CV SLA`)]
#splist = splist[!is.na(splist$`Mean LRR`)]
#
#write.xlsx(splist, paste0(Sys.Date(), "_Supplement_Species-List.xlsx"))
#with FRI sensitivity, #burns sensitivity, woodiness, height, CV_SLA