#Msweli et al_ms R1_20220510


setwd("C:/Users/samsp/Documents/...")

library(readxl)
FueltraitsDataset <- read_excel("CompleteFlammability datasets_Msweli_ST - Copy.xlsx", sheet = "FuelTraits_Dataset")
View(FueltraitsDataset)

library(Matrix)
library(lme4)

library(car)
library(effects)

library(ggplot2)

library(knitr)
library(rmarkdown)

library(patchwork)


##STEPWISE MULTIREG MODELS  

#1. Run stepwise Multiple Linear Regression -lm() Function - model for each flammability variable

maximumtemperatureRegM <- lm(scale(AverageMaximumTemperature) ~ 
                               scale(AverageFineFuels) + scale(AverageCoarseFuels) + scale(AverageDeadFuel)+
                               scale(AverageFBP) + scale(AverageDryBiomass) + scale(AverageFuelMoisture), 
                             data = FueltraitsDataset)
summary(maximumtemperatureRegM)
step(maximumtemperatureRegM, direction = "backward")



COBRegM <- lm(scale(AverageCompletenessOfBurn) ~ 
                scale(AverageFineFuels) + scale(AverageCoarseFuels) + scale(AverageDeadFuel)+
                scale(AverageFBP) + scale(AverageDryBiomass) + scale(AverageFuelMoisture), 
              data = FueltraitsDataset)
summary(COBRegM)
step(COBRegM, direction = "backward")



TimeToIgnitionRegM <- lm(scale(AverageTimeToIgnition) ~ 
                           scale(AverageFineFuels) + scale(AverageCoarseFuels) + scale(AverageDeadFuel)+
                           scale(AverageFBP) + scale(AverageDryBiomass) + scale(AverageFuelMoisture), 
                         data = FueltraitsDataset)
summary(TimeToIgnitionRegM)
step(TimeToIgnitionRegM, direction = "backward")




#1.2 run the MultiReg models with the lowest AIC 

MaximumtemperatureMod <- lm(formula = scale(AverageMaximumTemperature) ~ scale(AverageCoarseFuels) +
                              scale(AverageDryBiomass) + scale(AverageFuelMoisture), 
                            data = FueltraitsDataset)
summary(MaximumtemperatureMod)



COBMod <- lm(formula = scale(AverageCompletenessOfBurn) ~ scale(AverageFineFuels) +
               scale(AverageDryBiomass) + scale(AverageFuelMoisture), 
             data = FueltraitsDataset)
summary(COBMod)



TimeToIgnitionMod <- lm(formula = scale(AverageTimeToIgnition) ~ scale(AverageFineFuels) + scale(AverageFuelMoisture), 
                        data = FueltraitsDataset)
summary(TimeToIgnitionMod)



#Figure 1: Relationships between flammability measures (Maximum temperature, Completeness of burn, and Time to ignition) and fuel traits
#scatter plots of all models non-significant , without a line
####MAXIMUMTEMPERATURE

#ggarrange, plot_grid()
#library(cowplot)
#plot_grid(scatter_plotA, scatter_plotB, scatter_plotC)

scatter_plotA <- ggplot(FueltraitsDataset, aes(AverageMaximumTemperature, AverageFineFuels))
scatter_plotA + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = "Fine fuels (proportion)") + 
  theme_light() + 
  theme(legend.position = "none")



scatter_plotD <- ggplot(FueltraitsDataset, aes(AverageMaximumTemperature, AverageCoarseFuels))
scatter_plotD + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = "Coarse fuels (proportion)") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotG <- ggplot(FueltraitsDataset, aes(AverageMaximumTemperature, AverageDeadFuel))
scatter_plotG + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = "Dead fuels (proportion)") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotJ <- ggplot(FueltraitsDataset, aes(AverageMaximumTemperature, AverageFBP))
scatter_plotJ + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = "Fuel bed porosity (ratio)") + 
  theme_light()+
  theme(legend.position = "none")


#MAXIMUMTEMPERATURE with line

scatter_plotM <- ggplot(FueltraitsDataset, aes(AverageMaximumTemperature, AverageDryBiomass))
scatter_plotM + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = "Dry biomass (g)") + 
  geom_smooth(method = lm, color = 'black') + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotP <- ggplot(FueltraitsDataset, aes(AverageMaximumTemperature, AverageFuelMoisture))
scatter_plotP + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = expression("Maximum temperature ( "*degree*C*")"), y = "Fuel moisture (%)") + 
  geom_smooth(method = lm, color = 'black') + 
  theme_light() +
  theme(legend.position = "none")




#COMPLETENESSOFBURN with line

scatter_plotB <- ggplot(FueltraitsDataset, aes(AverageCompletenessOfBurn, AverageFineFuels))
scatter_plotB + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  geom_smooth(method = lm, color = 'black') + 
  theme_light() +
  theme(legend.position = "none")


#####COMPLETENESSOFBURN without line

scatter_plotE <- ggplot(FueltraitsDataset, aes(AverageCompletenessOfBurn, AverageCoarseFuels))
scatter_plotE + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light()+
  theme(legend.position = "none")



scatter_plotH <- ggplot(FueltraitsDataset, aes(AverageCompletenessOfBurn, AverageDeadFuel))
scatter_plotH + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " " ,y = " ") + 
  theme_light()+
  theme(legend.position = "none")



scatter_plotK <- ggplot(FueltraitsDataset, aes(AverageCompletenessOfBurn, AverageFBP))
scatter_plotK + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotN <- ggplot(FueltraitsDataset, aes(AverageCompletenessOfBurn, AverageDryBiomass))
scatter_plotN + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light() +
  theme(legend.position = "none")


scatter_plotQ <- ggplot(FueltraitsDataset, aes(AverageCompletenessOfBurn, AverageFuelMoisture))
scatter_plotQ + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = "Completeness of burn (proportion)", y = " ") + 
  theme_light() +
  theme(legend.position = "none")






#####TIMETOIGNITION without line

scatter_plotC <- ggplot(FueltraitsDataset, aes(AverageTimeToIgnition, AverageFineFuels))
scatter_plotC + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotF <- ggplot(FueltraitsDataset, aes(AverageTimeToIgnition, AverageCoarseFuels))
scatter_plotF + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light()+
  theme(legend.position = "none")


scatter_plotI <- ggplot(FueltraitsDataset, aes(AverageTimeToIgnition, AverageDeadFuel))
scatter_plotI + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotL <- ggplot(FueltraitsDataset, aes(AverageTimeToIgnition, AverageFBP))
scatter_plotL + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotO <- ggplot(FueltraitsDataset, aes(AverageTimeToIgnition, AverageDryBiomass))
scatter_plotO + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = " ", y = " ") + 
  theme_light() +
  theme(legend.position = "none")



scatter_plotR <- ggplot(FueltraitsDataset, aes(AverageTimeToIgnition, AverageFuelMoisture))
scatter_plotR + geom_point(aes(shape = Vegetationgroup, color = Vegetationgroup), size = 3) +
  scale_shape_manual(values = c(17, 15, 5)) +
  scale_color_manual(values = c('black', 'grey', 'black')) +
  labs(x = "Time to ignition (seconds)", y = " ") + 
  theme_light() +
  theme(legend.position = "none")




##Figure 2: BOXPLOTS per vegetation group (IAPs, Fynbos and Thicket) (EXPORT jpeg in 300x or 400x)

boxplotA <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageMaximumTemperature)) +
  geom_boxplot() + ylab(expression("Maximum temperature ( "*degree*C*")")) + 
  theme_light() 
plot(boxplotA)


boxplotB <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageCompletenessOfBurn)) +
  geom_boxplot() + ylab("Completeness of burn (proportion)") + 
  theme_light() 
plot(boxplotB)


boxplotC <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageTimeToIgnition)) +
  geom_boxplot() + ylab("Time to ignition (seconds)") + 
  theme_light() 
plot(boxplotC)


boxplotD <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageFineFuels)) +
  geom_boxplot() + ylab("Fine fuels (proportion)") + 
  theme_light() 
plot(boxplotD)


boxplotE <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageCoarseFuels)) +
  geom_boxplot() + ylab("Coarse fuels (proportion)") + 
  theme_light() 
plot(boxplotE)


boxplotF <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageDeadFuel)) +
  geom_boxplot() + ylab("Dead fuels (proportion)") + 
  theme_light() 
plot(boxplotF)


boxplotG <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageFBP)) +
  geom_boxplot() + ylab("Fuel bed porosity (ratio)") + 
  theme_light() 
plot(boxplotG)


boxplotH <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageDryBiomass)) +
  geom_boxplot() + ylab("Dry biomass (g)") + 
  theme_light() 
plot(boxplotH)


boxplotI <- ggplot(data = FueltraitsDataset, aes(x = Vegetationgroup, y = AverageFuelMoisture)) +
  geom_boxplot() + ylab("Fuel moisture (%)") + 
  theme_light() 
plot(boxplotI)


##################################The End########################################################