### Manoela Machado ##### ## 17/01/2023 ## # Thermal study and flammability experiment #### require(nlme) require(mgcv) require(ggplot2) require(tidyverse) require(gridExtra) require(effects) # ─── Understory air temperature --------------------------------------------- ## replacing log 0.1 for log of 1. New figures. ## effect of logging and its size on the thermal environment. rm(list = ls()) setwd('I:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/AirTemp_IButtons/') load('Analysis_EdgeEffect_iButton_H3b_RealDistances.RData') # no zeros for log: IB2017HeatRealDistMax$newStudyPoint[IB2017HeatRealDistMax$newStudyPoint%in%c(0)] <- 1 sort(unique(IB2017HeatRealDistMax$newStudyPoint)) # data head(IB2017HeatRealDistMax) unique(IB2017HeatRealDistMax$GapID) summary(IB2017HeatRealDistMax) ## revisions in April 2025 # range of gap size sort(unique(IB2017HeatRealDistMax$Area)) # 135.3650 to 1172.2326 ## Create size cat of gaps ## didn't solve anything. sort(unique(IB2017HeatRealDistMax$Area)) IB2017HeatRealDistMax %>% ggplot(.) + geom_histogram(aes(x=Area)) gArea <- IB2017HeatRealDistMax %>% ggplot(., aes(y=MaxTemp, x=Area, colour=newStudyPoint)) + geom_point() + geom_smooth() #### IB Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } ## air temp ~ distance from logging gap ## Q: Do distance from gap and gap area affect understorey air temperature? ## simple (+ Area) mod.1 <- lme(MaxTemp ~ log(newStudyPoint) + Area, random = list(~1|local, ~1|plotID), # random date is terrible data = IB2017HeatRealDistMax, method = 'ML') res(mod.1) # 30259.82 summary(mod.1) # area is not sig # interaction with area mod.2 <- lme(MaxTemp ~ log(newStudyPoint)*Area, random = list(~1|local, ~1|plotID), data = IB2017HeatRealDistMax, method = 'ML') res(mod.2) # 30256.86 summary(mod.2) # Area and the interaction with Area are sig # comparing the models anova(mod.2, mod.1) # sig, we go for complete # mod.3 <- lme(MaxTemp ~ log(newStudyPoint) + SizeCat, # random = list(~1|local, ~1|plotID), # # data = IB2017HeatRealDistMax, # method = 'ML') # # res(mod.3) # 30262.71 # summary(mod.3) # Size is not sig # # ## # mod.4 <- lme(MaxTemp ~ log(newStudyPoint)*SizeCat, # random = list(~1|local, ~1|plotID), # # data = IB2017HeatRealDistMax, # method = 'ML') # # res(mod.4) # 30264.7 # summary(mod.4) # Size and the interaction with SizeCat is not sig # # # comparing the models # anova(mod.3, mod.4) # not sig, go for simple (without interaction with GapSize) # ## Predicts with Area ## ## 30/10 -- predict from new model (log(1) and new figs with quantiles) # figure with normal x newdatIB <- expand.grid(newStudyPoint = seq(min(IB2017HeatRealDistMax$newStudyPoint), max(IB2017HeatRealDistMax$newStudyPoint), length.out = 200), Area = seq(min(IB2017HeatRealDistMax$Area), max(IB2017HeatRealDistMax$Area))) head(newdatIB) newdatIB$predsmod.2 <- predict(mod.2, newdatIB, level=0) head(newdatIB) #### Plot Predicts #### gIB <- newdatIB %>% ggplot(.) + geom_point(data=IB2017HeatRealDistMax, aes(y=MaxTemp, x=newStudyPoint), alpha=.5, colour='grey') + geom_point(data=newdatIB, aes(y=predsmod.2, x=newStudyPoint, colour=Area)) + scale_colour_gradient(low = "black", high = "grey50", name='Gap area (m2)') + labs(y='Understory ambient air (?C)', x='Distance from logging gap center (m)') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gIB ggsave(gIB, file='../Figs/FigAirTemp.png', width = 8, height = 6) #### figure log x #### # mod.2 <- lme(MaxTemp ~ log(newStudyPoint)*Area, data = IB2017HeatRealDistMax) interaction1 <- effect("log(newStudyPoint)*Area", mod.2, xlevels=3) plot(interaction1, multiline=TRUE) interaction2_IB <- data.frame(interaction1) interaction2_IB$Area <- as.factor(interaction2_IB$Area) newyNMDS <- predict(mod.2, type="response",level=0, se.fit = TRUE) glimpse(newyNMDS) gIB_log <- ggplot(IB2017HeatRealDistMax) + geom_point(aes(x=log(newStudyPoint), y=MaxTemp), alpha=.5, colour='grey') + geom_line(data=interaction2_IB, aes(x=log(newStudyPoint), y=fit, colour=Area), linetype="solid", linewidth=0.8, alpha=0.7) + scale_colour_manual(values = c('black','black','black'), name='Gap area (m2)') + labs(y='Understory ambient air (?C)', x='Log transformed distance from \nlogging gap center (m)') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gIB_log # ggsave(gIB_log, file='../Figs/FigAirTemp_log.png', width = 8, height = 6) #### save IB RData #### save.image('AirTemp_IB_analysis_results_Oct2023.RData') # ─── Surface temperature --------------------------------------------- rm(list = ls()) setwd('H:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/SurfaceTemp_FLIR/') load("datFLIR_FinalForAnalysis.RData") glimpse(datFLIRMax) glimpse(datFLIR) # Log(x) # # no zeros: datFLIRMax$newStudyPoint[datFLIRMax$newStudyPoint%in%c(0)] <- 1 sort(unique(datFLIRMax$newStudyPoint)) # ## Create size cat of gaps # sort(unique(datFLIRMax$Area)) # datFLIRMax %>% # ggplot(.) + # geom_histogram(aes(x=Area)) #+ # # geom_point(aes(y=Area,x=SizeCat)) # # datFLIRMax <- datFLIRMax %>% # mutate(SizeCat = ifelse(Area>=1000, 'Large', 'Small')) # # datFLIRMax %>% ggplot(.) + # geom_point(aes(y=MaxUpper.tail, x=Area, colour=SizeCat)) # gDist <- datFLIRMax %>% ggplot(.,aes(y=MaxUpper.tail, x=newStudyPoint)) + # geom_point(aes(colour=SizeCat)) + # geom_smooth(aes(colour=SizeCat)) + # labs(x='Distance from gap center') gArea <- datFLIRMax %>% ggplot(.,aes(y=MaxUpper.tail, x=Area, colour=newStudyPoint)) + geom_point() + geom_smooth() # gridExtra::grid.arrange(gDist, gArea, ncol=1) # #### FLIR Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } ## ground max temp ~ distance from logging gap ## Q: Do distance from gap and gap area affect ground surface temperature? ## simple mod.1 <- lme(MaxUpper.tail ~ log(newStudyPoint) + Area, random = list(~1|GapID, ~1|PlotID), data = datFLIRMax, method = 'ML') res(mod.1) summary(mod.1) # area is not sig # interaction with area mod.2 <- lme(MaxUpper.tail ~ log(newStudyPoint)*Area, random = list(~1|GapID, ~1|PlotID), data = datFLIRMax, method = 'ML') res(mod.2) summary(mod.2) # the interaction is sig # comparing models anova(mod.2, mod.1) # sig, go for complete # # size cat instead of area # mod.3 <- lme(MaxUpper.tail ~ log(newStudyPoint) + SizeCat, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.3) # 2838.576 # summary(mod.3) # size cat is not sig # # ## interaction with Size Cat # mod.4 <- lme(MaxUpper.tail ~ log(newStudyPoint)*SizeCat, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.4) # 2839.91 # summary(mod.4) # size cat and the interaction are not sig # # # comparing the models # anova(mod.3, mod.4) # not sig, go for simple (without the interaction) # figure with normal x newdat_FLIR_MAX <- expand.grid(newStudyPoint = seq(min(datFLIRMax$newStudyPoint), max(datFLIRMax$newStudyPoint), length.out = 200), Area = seq(min(datFLIRMax$Area), max(datFLIRMax$Area), length.out = 200)) head(newdat_FLIR_MAX) newdat_FLIR_MAX$predsmod.2_MAX <- predict(mod.2, newdat_FLIR_MAX, level=0) #### Plot Predicts #### gFLIR_Max <- newdat_FLIR_MAX %>% ggplot(.) + geom_point(data=datFLIRMax, aes(y=MaxUpper.tail, x=newStudyPoint), alpha=.5, colour='grey') + geom_point(aes(y=predsmod.2_MAX, x=newStudyPoint, colour=Area)) + scale_colour_gradient(low = "black", high = "grey50", name='Gap area (m2)') + labs(y='Ground surface temperature (?C)', x='Distance from logging gap center (m)') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gFLIR_Max ggsave(gFLIR_Max, file='../Figs/SurfaceTemp_Max.png', width = 8, height = 6) #### figure log x #### # mod.2 <- lme(MaxTemp ~ log(newStudyPoint)*Area, data = IB2017HeatRealDistMax) interaction1 <- effect("log(newStudyPoint)*Area", mod.2, xlevels=3) plot(interaction1, multiline=TRUE) interaction2_FLIR <- data.frame(interaction1) interaction2_FLIR$Area <- as.factor(interaction2_FLIR$Area) newyNMDS <- predict(mod.2, type="response",level=0, se.fit = TRUE) head(newyNMDS) gFLIR_Max_log <- ggplot(datFLIRMax) + geom_point(aes(x=log(newStudyPoint), y=MaxUpper.tail), alpha=.5, colour='grey') + geom_line(data=interaction2_FLIR, aes(x=log(newStudyPoint), y=fit, colour=Area), linetype="solid", size=0.8, alpha=0.7) + scale_colour_manual(values = c('black','black','black'), name='Gap area (m2)') + labs(y='Understory ambient air (?C)', x='Log transformed distance from \nlogging gap center (m)') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gFLIR_Max_log ggsave(gFLIR_Max_log, file='../Figs/SurfaceTemp_Max_log.png', width = 8, height = 6) #### save FLIR max RData #### save.image('FLIRmax_analysis_results_Oct2023.RData') # ─── Air Humidity --------------------------------------------- rm(list = ls()) setwd('H:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/AirHumidity/') load("datFLIR_FinalForAnalysis.RData") head(datFLIRMax) # no zeros: datFLIRMax$newStudyPoint[datFLIRMax$newStudyPoint%in%c(0)] <- 1 sort(unique(datFLIRMax$newStudyPoint)) # ## Create size cat of gaps # sort(unique(datFLIRMax$Area)) # datFLIRMax %>% # ggplot(.) + # geom_histogram(aes(x=Area)) #+ # # geom_point(aes(y=Area,x=SizeCat)) # # datFLIRMax <- datFLIRMax %>% # mutate(SizeCat = ifelse(Area>=1000, 'Large', 'Small')) # # datFLIRMax %>% ggplot(.) + # geom_point(aes(y=MinLower.tail, x=Area, colour=SizeCat)) # gDist <- datFLIRMax %>% ggplot(.,aes(y=MeanHumidity, x=newStudyPoint)) + # geom_point(aes(colour=SizeCat)) + # geom_smooth(aes(colour=SizeCat)) + # geom_smooth() + # labs(x='Distance from gap center') gArea <- datFLIRMax %>% ggplot(.,aes(y=MeanHumidity, x=Area, colour=newStudyPoint)) + geom_point() + geom_smooth() # gridExtra::grid.arrange(gDist, gArea, ncol=1) ## so, no one is log# #### Air Hum Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } ## ground max temp ~ distance from logging gap ## Q: Do distance from gap and gap area affect ground surface temperature? ## simple mod.1 <- lme(MeanHumidity ~ (newStudyPoint) + Area, random = list(~1|GapID, ~1|PlotID), data = datFLIRMax, method = 'ML') res(mod.1) # decent summary(mod.1) # area is not sig # # log ## # mod.1Log <- lme(MeanHumidity ~ log(newStudyPoint) + Area, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.1Log) # res are worse # summary(mod.1Log) # # AIC(mod.1, mod.1Log) # freaking same # # # comparing models # anova(mod.1Log, mod.1) # # # # interaction with area # mod.2 <- lme(MeanHumidity ~ (newStudyPoint)*Area, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.2) # summary(mod.2) # the interaction is not sig # # # log # mod.2Log <- lme(MeanHumidity ~ log(newStudyPoint)*Area, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.2Log) # summary(mod.2) # the interaction is not sig # # ## comparing the interaction term # anova(mod.2, mod.1) # not sig, go for simple mod 1 # # # ## Null model # mod.1 <- lme(MeanHumidity ~ (newStudyPoint) + Area, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # mod.1.NULL <- lme(MeanHumidity ~ 1, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # anova(mod.1, mod.1.NULL) ## they are different, so i reject the null # # # log # mod.1Log <- lme(MeanHumidity ~ log(newStudyPoint) + Area, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # mod.1Log.NULL <- lme(MeanHumidity ~ 1, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # anova(mod.1Log, mod.1Log.NULL) # they are diff, so i reject the null # # # # ## using Size Cat instead of Area (numeric) # mod.3 <- lme(MeanHumidity ~ (newStudyPoint)+SizeCat, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.3) # 3190.164 # summary(mod.3) # Size Cat is not sig # # ## interaction with Size Cat # mod.4 <- lme(MeanHumidity ~ (newStudyPoint)*SizeCat, # random = list(~1|GapID, ~1|PlotID), # data = datFLIRMax, # method = 'ML') # # res(mod.4) # 3192.155 # summary(mod.4) # Size Cat is not sig, interaction neither # # ## compraring models # anova(mod.4, mod.3) ## not sig, go for simple (without the interaction) # # summary(mod.1) newdat_Hum <- expand.grid(newStudyPoint = seq(min(datFLIRMax$newStudyPoint), max(datFLIRMax$newStudyPoint), length.out = 300), Area = max(datFLIRMax$Area)) head(newdat_Hum) newdat_Hum$predsmod.1 <- predict(mod.1, newdat_Hum, level=0) # install.packages("remotes") # remotes::install_github("bsurial/bernr") CI <- bernr::bolker_ci(mod.1, newdat_Hum, pred_int = FALSE, conf_level = 0.95) head(CI) #### Plot Predicts #### gFLIRhum <- datFLIRMax %>% ggplot(.) + geom_point(data=datFLIRMax, aes(y=MeanHumidity, x=newStudyPoint), alpha = 0.5, colour='grey') + geom_point(data = newdat_Hum, aes(y=predsmod.1, x=newStudyPoint), size=2) + geom_point(data = CI, aes(y=predsmod.1, x=newStudyPoint), size=2) + geom_ribbon(data=CI, aes(y=predsmod.1, x=newStudyPoint, ymin=ci_l, ymax=ci_h), alpha=0.2) + scale_colour_gradient(low = "black", high = "grey50", name='Gap area (m2)') + labs(y='Air humidity (%)', x='Distance from logging gap center (m)') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.2)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gFLIRhum ggsave(gFLIRhum, file='../Figs/AirHumidity.png', width = 8, height = 6) #### save AirHum RData #### save.image('AirHumid_analysis_results_Oct2023.RData') # ─── Fuel moisture --------------------------------------------- rm(list = ls()) setwd('H:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/MoistureContent') load("datMCnoNA.RData") glimpse(datMCnoNA) range(datMCnoNA$days_since_rain) datMCnoNA$days_since_rain[datMCnoNA$days_since_rain%in%c(0.1)] <- 1 datMCnoNA$newStudyPoint[datMCnoNA$newStudyPoint%in%c(0.1)] <- 1 gDays <- datMCnoNA %>% filter(days_since_rain < 40) %>% ggplot(.,aes(y=moisture_content, x=days_since_rain)) + geom_point(aes(colour=Area)) + geom_smooth() gDist <- datMCnoNA %>% ggplot(.,aes(y=moisture_content, x=newStudyPoint, colour=Area)) + geom_point() + geom_smooth() + labs(x='Distance from gap center') gArea <- datMCnoNA %>% ggplot(.,aes(y=moisture_content, x=Area, colour=days_since_rain)) + geom_point() + geom_smooth() ## stop the experiment at 40 days after rain because lack of data for all the range of gap sizes range(datMCnoNA$days_since_rain) range(datMCnoNA$newStudyPoint) datMCnoNA <- datMCnoNA %>% filter(days_since_rain <= 40) # ## So: x=days since rain is log ## x= distance from gap centrer is not log ! #### Fuel Moisture Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } ## rate of moisture content loss = MC ~ days since rain ## Q: Do distance from gap and gap area affect leaf litter desiccation? ## stop the experiment at 40 days after rain because lack of data for all the range of gap sizes datMCnoNA <- datMCnoNA %>% filter(days_since_rain <= 40) ## mod.1 <- lme(moisture_content ~ log(days_since_rain) + newStudyPoint + Area, random = list(~1|DateG, ~1|GapID), data = datMCnoNA, method = 'ML') res(mod.1) ## ok summary(mod.1) # area is not sig # ## interaction log(days since rain) * distance # mod.2 <- lme(moisture_content ~ log(days_since_rain)*newStudyPoint + Area, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA, # method = 'ML') # res(mod.2) ## ok # summary(mod.2) ## inter is not sig! # # ## comparing mod 1 and mod 2 # anova(mod.2, mod.1) # not sig , go for simple # # # ## interaction log(days since rain) * area # mod.3 <- lme(moisture_content ~ log(days_since_rain)*Area + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA, # method = 'ML') # res(mod.3) ## ok # summary(mod.3) ## inter is not sig! # # # comparing mod 1 and mod 3 # anova(mod.3, mod.1) # not sig, go for simple mod.1 # # ## two way interaction # mod.4 <- lme(moisture_content ~ log(days_since_rain)*Area + log(days_since_rain)*newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA, # method = 'ML') # res(mod.4) ## ok # summary(mod.4) ## inter with distance is not sig. inter with are is sig. # # ## # anova(mod.3, mod.4) # not sig, go for simple mod.3 # # # ## REML # mod.3REML <- lme(moisture_content ~ log(days_since_rain)*Area + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA) # # ## Size Cat instead of Area # mod.5 <- lme(moisture_content ~ log(days_since_rain) + SizeCat + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA, # method = 'ML') # res(mod.5) ## ok # summary(mod.5) ## all sig # # # interaction # mod.6 <- lme(moisture_content ~ log(days_since_rain)*SizeCat + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA, # method = 'ML') # res(mod.6) ## ok # summary(mod.6) ## size cat and interaction are not sig # # ## comparing models # anova(mod.5, mod.6) # not sig, go for simple mod.5 # # # two way interaction # mod.7 <- lme(moisture_content ~ log(days_since_rain)*SizeCat + log(days_since_rain)*newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datMCnoNA, # method = 'ML') # res(mod.7) ## ok # summary(mod.7) ## size cat and interaction are not sig # # anova(mod.5, mod.7) # not sig, go for simple mod 5 ## ## Predicts with Area ## ## New dat for Minimum distance (centre) newdat_MC_MinDist <- expand.grid(days_since_rain = seq(min(datMCnoNA$days_since_rain), max(datMCnoNA$days_since_rain), length.out = 500), Area = max(datMCnoNA$Area), newStudyPoint = min(datMCnoNA$newStudyPoint)) # = centre range(newdat_MC_MinDist$days_since_rain) newdat_MC_MinDist$predsmod.1 <- predict(mod.1, newdat_MC_MinDist, level=0) head(newdat_MC_MinDist) CI_min <- bernr::bolker_ci(mod.1, newdat_MC_MinDist, pred_int = FALSE, conf_level = 0.95) head(CI_min) ## New dat for mean distance newdat_MC_MeanDist <- expand.grid(days_since_rain = seq(min(datMCnoNA$days_since_rain), max(datMCnoNA$days_since_rain), length.out = 500), Area = max(datMCnoNA$Area), newStudyPoint = mean(datMCnoNA$newStudyPoint)) # newdat_MC_MeanDist$predsmod.1 <- predict(mod.1, newdat_MC_MeanDist, level=0) head(newdat_MC_MeanDist) CI_mean <- bernr::bolker_ci(mod.1, newdat_MC_MeanDist, pred_int = FALSE, conf_level = 0.95) head(CI_mean) ## New dat for max dsitance newdat_MC_MaxDist <- expand.grid(days_since_rain = seq(min(datMCnoNA$days_since_rain), max(datMCnoNA$days_since_rain), length.out = 500), Area = max(datMCnoNA$Area), newStudyPoint = max(datMCnoNA$newStudyPoint)) # newdat_MC_MaxDist$predsmod.1 <- predict(mod.1, newdat_MC_MaxDist, level=0) head(newdat_MC_MaxDist) CI_max <- bernr::bolker_ci(mod.1, newdat_MC_MaxDist, pred_int = FALSE, conf_level = 0.95) head(CI_max) ## dat_CI <- CI_min %>% rbind(CI_mean, CI_max) head(dat_CI) dat_CI$newStudyPoint <- as.factor(dat_CI$newStudyPoint) unique(dat_CI$newStudyPoint) dat_CI$Distance <- NA dat_CI$Distance[dat_CI$newStudyPoint%in%c(1)] <- 'Centre' dat_CI$Distance[dat_CI$newStudyPoint%in%c(26.6228784648188)] <- 'Mean distance' dat_CI$Distance[dat_CI$newStudyPoint%in%c(75.25 )] <- 'Maximum distance' head(dat_CI) unique(dat_CI$Distance) dat_CI$Distance <- factor(dat_CI$Distance, c('Centre', 'Mean distance', 'Maximum distance')) head(dat_CI) # summary(mod.1) #### Plot predicts ##### gMC <- datMCnoNA %>% # filter(days_since_rain <40 ) %>% ggplot(.) + geom_point(aes(y=moisture_content, x=days_since_rain), alpha=.5, colour='grey') + # distances geom_line(data= newdat_MC_MinDist, aes(y=predsmod.1, x=days_since_rain), linetype='dashed', colour='black', alpha=.7) + geom_line(data= newdat_MC_MeanDist, aes(y=predsmod.1, x=days_since_rain), linetype='dashed', colour='red', alpha=.7) + geom_line(data= newdat_MC_MaxDist, aes(y=predsmod.1, x=days_since_rain), linetype='dashed', colour='blue', alpha=.7) + geom_point(aes(y=predsmod.1, x=days_since_rain), colour='black', alpha=.7) + labs(y='Leaf litter moisture content (%)', x='Days after rain') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + #scale_color_continuous(name='Gap area (m2)') + #coord_cartesian(xlim=c(0,3)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gMC ggsave(gMC, file='../Figs/MC.png', width = 8, height = 5) #### figure with log x #### head(dat) gMC_log <- ggplot(datMCnoNA) + geom_point(aes(x=log(days_since_rain), y=moisture_content), alpha=.5, colour='grey') + geom_line(data=dat_CI, aes(x=log(days_since_rain), y=predsmod.1, linetype=Distance), size=1) + geom_ribbon(data=dat_CI, aes(y=predsmod.1, x=log(days_since_rain), ymin=ci_l, ymax=ci_h, fill=Distance), alpha=0.2) + labs(y='Leaf litter moisture content (%)', x='Ln(number of days after rain)') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gMC_log ggsave(gMC_log, file='../Figs/MC_log.png', width = 8, height = 6) #### save Fuel Moisture Content RData #### save.image('MC_analysis_results_Oct2023.RData') # ─── Probability of ignition --------------------------------------------- # Proportion of ignition (5 attempts per sample) vs Moisture Content setwd('H:/Meu Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/Probabilty_Ignition/') load('dat_flamm_PropIG.RData') glimpse(dat_flamm_PropIG) class(dat_flamm_PropIG$PropStatusIG) # num class(dat_flamm_PropIG$moisture_content) # num dat_flamm_PropIG$Plot <- factor(dat_flamm_PropIG$Plot, levels=c("centre","edge","10m","20m","50m","PF", ordered=TRUE)) dat_flamm_PropIGnoNA <- dat_flamm_PropIG %>% filter(!is.na(dat_flamm_PropIG$PropStatusIG) & !is.na(dat_flamm_PropIG$moisture_content)) # 617 obs datThresIG <- dat_flamm_PropIGnoNA head(datThresIG) # plotting ggplot(datThresIG) + geom_point(aes(y=PropStatusIG, x=moisture_content, colour=Plot), size=3, alpha=0.6) + geom_smooth(aes(y=PropStatusIG, x=moisture_content), se = F) + labs(x="Leaf-litter moisture content", y="Proportion of ignition") #geom_vline(xintercept = 17, linetype="dashed", color = "red") #### Ig prop Analysis ##### #### glm #### # # Proportional data 0-1 -> GLM (p.254 Zuur) # # my samples are independent? each value of MC is from each individual sample, from a plot nested in a GapID # # maybe GLMM? (chap 13) # sort(unique(datThresIG$PropStatusIG)) # 0 1 # # # mod1 <- glm(cbind(PropStatusIG*5,5-PropStatusIG*5) ~ moisture_content, # adjusting the success / failure # data = datThresIG, family = binomial) # # summary(mod1) # qqnorm(residuals(mod1)) # # # Check prediction from the model: page 257 # # the thick line in the middle shows the estimated probabilities for a range of length values for the female data # MyData <- data.frame(moisture_content = seq(from = min(datThresIG$moisture_content), # to = max(datThresIG$moisture_content), by = 1)) # P1 <- predict(mod1, newdata = MyData, type = "link", se = TRUE) # plot(MyData$moisture_content, exp(P1$fit) / (1+exp(P1$fit)), # type='l', ylim=c(0,1), xlab="MC", ylab="Prob IG") # lines(MyData$moisture_content, exp(P1$fit+1.96*P1$se.fit)/ # (1+exp(P1$fit+1.96*P1$se.fit)), lty=2) # lines(MyData$moisture_content, exp(P1$fit-1.96*P1$se.fit)/ # (1+exp(P1$fit-1.96*P1$se.fit)), lty=2) # points(datThresIG$moisture_content, datThresIG$PropStatusIG) # #### mass #### # ### adding random part # require(MASS) # mod2 <- glmmPQL(cbind(PropStatusIG*5,5-PropStatusIG*5) ~ moisture_content, # random = ~1 | GapID, family = binomial, # data = datThresIG) # # summary(mod2) # # Fixed effects: cbind(PropStatusIG * 5, 5 - PropStatusIG * 5) ~ moisture_content # # Value Std.Error DF t-value p-value # # (Intercept) -4.25557307639 0.2779839693184 598 -15.3086995873 0 # # moisture_content -0.12519323553 0.0103182913704 598 12.1331362956 0 # # ## plotting # g <- (4.2555730) - 0.1251932 * datThresIG$moisture_content # p.averageGapID <- exp(g) / (1 + exp(g)) # I1 <- order(datThresIG$moisture_content) # avoid spaghetti plot # plot(datThresIG$moisture_content, datThresIG$PropStatusIG, # ylab = "Probability of NOT Ignition (1=NOT ignited)", xlab = "Leaf-litter moisture content") # lines(datThresIG$moisture_content[I1], p.averageGapID[I1], lwd=3) # p.Upp<-exp(g+1.96*0.694417213018)/(1+exp(g+1.96*0.694417213018)) # StdDev: 0.694417213018 (from summary) # p.Low<-exp(g-1.96*0.694417213018)/(1+exp(g-1.96*0.694417213018)) # StdDev: 0.694417213018 (from summary) # lines(datThresIG$moisture_content[I1], p.Upp[I1]) # lines(datThresIG$moisture_content[I1], p.Low[I1]) # #### glmer lme4 #### ## same, but using lme4: require(lme4) modProbIg1 <- glmer(cbind(PropStatusIG*5,5-PropStatusIG*5) ~ moisture_content + (1 | GapID), family = binomial, data = datThresIG) AIC(modProbIg1) # 1555 qqnorm(residuals(modProbIg1)) summary(modProbIg1) # # Fixed effects: # Estimate Std. Error z value Pr(>|z|) # (Intercept) -4.36134357751 0.26310358654 -16.57653 < 2.22e-16 *** # moisture_content 0.12502677542 0.00827133754 15.11567 < 2.22e-16 *** #### plotting modProbIg1 #### # g <- (4.36134869299) + -0.12502691220 * datThresIG$moisture_content # from summary # p.averageGapID <- exp(g) / (1 + exp(g)) # I1 <- order(datThresIG$moisture_content) # avoid spaghetti plot # plot(datThresIG$moisture_content, datThresIG$PropStatusIG, # ylab = "Probability of NOT Ignition (1=NOT ignited)", xlab = "Leaf-litter moisture content") # lines(datThresIG$moisture_content[I1], p.averageGapID[I1], lwd=3) # p.Upp<-exp(g+1.96*0.784803165)/(1+exp(g+1.96*0.784803165)) # StdDev: from summary # p.Low<-exp(g-1.96*0.784803165)/(1+exp(g-1.96*0.784803165)) # StdDev: from summary # lines(datThresIG$moisture_content[I1], p.Upp[I1]) # lines(datThresIG$moisture_content[I1], p.Low[I1]) # #### random str ##### # add nested Study Plot modProbIg2 <- glmer(cbind(PropStatusIG*5,5-PropStatusIG*5) ~ moisture_content + (1 | GapID / Plot), family = binomial, data = datThresIG) AIC(modProbIg2) # 1441 better par(mfrow=c(2,2)) plot(fitted(modProbIg2), residuals(modProbIg2)) # qqnorm(residuals(modProbIg2)) qqline(residuals(modProbIg2)) acf(residuals(modProbIg2)) # summary(modProbIg2) # Fixed effects: # Estimate Std. Error z value Pr(>|z|) # (Intercept) -4.70270065214 0.27939386552 -16.8318 < 2.22e-16 *** # moisture_content 0.13699766477 0.00933366356 14.6778 < 2.22e-16 *** #### plotting modProbIg2 ##### dev.off() g <- (4.70269522197) - 0.13699786734 * datThresIG$moisture_content # from summary p.averageGapID <- exp(g) / (1 + exp(g)) I1 <- order(datThresIG$moisture_content) # avoid spaghetti plot plot(datThresIG$moisture_content, datThresIG$PropStatusIG, ylab = "Probability of Ignition (1= ignited)", xlab = "Leaf-litter moisture content") lines(datThresIG$moisture_content[I1], p.averageGapID[I1], lwd=3) p.Upp<-exp(g+1.96*0.709129068)/(1+exp(g+1.96*0.709129068)) # StdDev: from summary p.Low<-exp(g-1.96*0.709129068)/(1+exp(g-1.96*0.709129068)) # StdDev: from summary lines(datThresIG$moisture_content[I1], p.Upp[I1]) lines(datThresIG$moisture_content[I1], p.Low[I1]) ### getting more results mod 4 # my data with only one gapID: MyData <- data.frame(moisture_content = seq(from = min(datThresIG$moisture_content), to = max(datThresIG$moisture_content), by = 1), GapID=datThresIG$GapID[1]) # my data with a fixed value of moisture content: MyData <- data.frame(moisture_content = 12) # , -> Holdsworth and Uhl 1997 Ecol Appl, Threshold = 12% MyData <- data.frame(moisture_content = 15) # , -> Uhl and Kauffman 1990 Ecology, Threshold = 15% MyData <- data.frame(moisture_content = 50) # MC = 50 %, 0.1045 chance of burning predict(modProbIg2, newdata = MyData, re.form = NA, type='response') ### for moisture content of 12%, prob of ignition: 0.955157539234 !!!!!!!!!!! ### for moisture content of 15%, prob of ignition: 0.933870770186 !!!!!!!!!!! # or: ### MyData <- data.frame(moisture_content = seq(from = min(datThresIG$moisture_content), to = max(datThresIG$moisture_content), by = 0.5), GapID='PF', Plot='PF') predict(modProbIg2, newdata = MyData, type='response') # for the overall relationship -> to get the value of MC that gives 0.5 prob of not igniting MyData <- data.frame(moisture_content = seq(from = min(datThresIG$moisture_content), to = max(datThresIG$moisture_content), by = 0.1)) MyData$predictedData <- predict(modProbIg2, newdata = MyData, re.form = NA, type='response') MyData %>% filter(predictedData >= 0.49 & predictedData < 0.51) # for 50% change of ignition, moisture content = 34.27 lines(predict(modProbIg2, newdata = MyData, re.form = NA, type='response') ~ MyData$moisture_content, type='p') # overal prediction curve ## summary(modProbIg2) prob_fct <- function(x) boot::inv.logit(-0.13699786734*x + 4.70269522197) # from modProbIg2 halfThreshold <- function(x) 0.5 plot(prob_fct, 0,50, add=T,col="darkgreen") # demand equation. eqx <- function(x) halfThreshold(x) - prob_fct(x) EQ <- uniroot(eqx,c(0,300)) # Price where supp & demand intersect (Price at equilibrium) EQ$root # 34.3267517061 exactly! ## alternative newdata <- with(modProbIg2, expand.grid(moisture_content = seq(min(datThresIG$moisture_content), to = max(datThresIG$moisture_content), len = 100))) Xmat <- model.matrix(~ moisture_content, newdata) fixest <- fixef(modProbIg2) fit <- as.vector(fixest %*% t(Xmat)) SE <- sqrt(diag(Xmat %*% vcov(modProbIg2) %*% t(Xmat))) q <- qt(0.975, df = df.residual(modProbIg2)) linkinv <- binomial()$linkinv newdata <- cbind(newdata, fit = linkinv(fit), lower = linkinv(fit - q * SE), upper = linkinv(fit + q * SE)) ggplot(newdata, aes(y=fit, x=moisture_content)) + geom_line() + geom_ribbon(aes(ymin=lower, ymax=upper), color=NA, alpha=0.3)+ geom_point(data = datThresIG, aes(y=PropStatusIG), size=3, alpha=0.6) + # , colour=Plot geom_segment(x = -1, y = 0.50, xend = 34.326, yend = 0.50, linetype="dashed", color = "grey10", size=.5) + # x = 34.326 for y = 0.50 geom_segment(x = 34.326, y = 0.50, xend = 34.326, yend = -1, linetype="dashed", color = "grey10", size=.5) + geom_segment(x = -1, y = 0.955, xend = 12, yend = 0.955, linetype="dashed", color = "grey10", size=.5) + # x = 12 for y = 0.955 geom_segment(x = 12, y = 0.955, xend = 12, yend = -1, linetype="dashed", color = "grey10", size=.5) newdataProbIG <- newdata #### save Probability of Ignition RData ##### save.image(file='Prob_Ignition_vs_MoistureContent_Jan2023.RData') # ─── Ignitability --------------------------------------------- rm(list = ls()) setwd('H:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/Ignitability/') # data load("datDTnoNA.RData") glimpse(datDTnoNA) range(datDTnoNA$days_since_rain) datDTnoNA$days_since_rain[datDTnoNA$days_since_rain%in%c(0.1)] <- 1 datDTnoNA$newStudyPoint[datDTnoNA$newStudyPoint%in%c(0.1)] <- 1 # ## create a separate dataset with Gap and Area only # unique(datDTnoNA$GapID) # datArea <- datDTnoNA %>% summarise(GapID= unique(GapID), # save(datArea, file='../Combustability/datArea.RData') # ## Create size cat of gaps ## didnt solve anything, not using it. # sort(unique(datDTnoNA$Area)) # datDTnoNA <- datDTnoNA %>% # mutate(SizeCat = ifelse(Area>=1000, 'Large', 'Small')) # # datDTnoNA %>% # ggplot(.) + # #geom_histogram(aes(x=Area)) + # geom_point(aes(y=Area,x=SizeCat)) ## graphs explor ## # gDays <- datDTnoNA %>% ggplot(.,aes(y=DTPosix, x=days_since_rain)) + # geom_point(aes(colour=Area)) + # geom_smooth() # datDTnoNA %>% # # filter(GapID != '104117' & days_since_rain!=46 & DTPosix < "2019-06-08 00:00:02 GMT") %>% # ggplot(.,aes(y=DTPosix, x=days_since_rain)) + # geom_point(aes(shape=SizeCat, colour=GapID)) + # geom_smooth(aes(colour=SizeCat)) # # gDist <- datDTnoNA %>% ggplot(.,aes(y=DTPosix, x=newStudyPoint)) + # geom_point(aes(colour=Area)) + # geom_smooth() + labs(x='Distance from gap center') # # gArea <- datDTnoNA %>% ggplot(.,aes(y=DTPosix, x=Area)) + # geom_point(aes(colour=days_since_rain)) + # geom_smooth() # # gridExtra::grid.arrange(gDays, gDist, gArea, ncol=1) #### stop the experiment at 40 days after rain because lack of data for all the range of gap sizes range(datDTnoNA$days_since_rain) # 1 - 72 datDTnoNA <- datDTnoNA %>% filter(days_since_rain <= 40) ## Create a column for flame duration in seconds glimpse(datDTnoNA) datDTnoNA$DTSecsWithDate <- lubridate::time_length(datDTnoNA$DTnum, unit='second') g <- min(datDTnoNA$DTSecsWithDate) as.POSIXct(as.numeric(g+1e-6), origin = "1970-01-01", tz="GMT") # zero seconds # Remove the min to get rid of all the date part datDTnoNA <- datDTnoNA %>% mutate(DTSecs = DTSecsWithDate - g) glimpse(datDTnoNA) # ## So: x=days since rain is log ## x= distance from gap centre is not log ! #### Ignitability Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } ## delay time to ignition ~ days since rain ## Q: Do distance from gap and gap area affect the delay time to ignition? range(datDTnoNA$days_since_rain) ## range(datDTnoNA$newStudyPoint) mod.1 <- lme(DTSecs ~ log(days_since_rain) + newStudyPoint + Area, random = list(~1|DateG, ~1|GapID), #random = ~1|GapID, data = datDTnoNA, method = 'ML', control=lmeControl(msMaxIter = 30, tolerance = 1e-04, niterEM = 25, msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) res(mod.1) ## 664 (without date as random 670) summary(mod.1) # distance is not sig. are is barely sig. ## interaction log(days since rain) * distance ## # mod.2 <- lme(DTSecs ~ log(days_since_rain)*newStudyPoint + Area, # random = list(~1|DateG, ~1|GapID), # data = datDTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # # res(mod.2) ## 666 (without date as random 672) # summary(mod.2) ## inter is not sig! # # ## comparing mod 1 and mod 2 # anova(mod.2, mod.1) # not sig , go for simple mod.1 # # # ## interaction log(days since rain) * area # mod.3 <- lme(DTSecs ~ log(days_since_rain)*Area + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datDTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.3) ## 666 # summary(mod.3) ## inter is not sig # # # comparing mod 1 and mod 3 # anova(mod.3, mod.1) # not sig, go for simple mod 1 # # # ## Size Cat instead of Area # mod.4 <- lme(DTSecs ~ log(days_since_rain) + SizeCat + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datDTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.4) ## 650 # summary(mod.4) ## Size cat is sig. Distance from center is not sig. # # # interaction with size cat # mod.5 <- lme(DTSecs ~ log(days_since_rain)*SizeCat + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datDTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.5) ## 652 # summary(mod.5) ## Size cat is sig. Distance from center is not sig. # # # comparing models # anova(mod.5, mod.4) # not sig, go for simple mod 4 (without interaction with Size Cat) # ## ## Predicts with Area ## range(datDTnoNA$newStudyPoint) newdat_DT <- expand.grid(days_since_rain = seq(min(datDTnoNA$days_since_rain), max(datDTnoNA$days_since_rain), length.out = 1000), Area = max(datDTnoNA$Area), newStudyPoint = min(datDTnoNA$newStudyPoint)) # = centre head(newdat_DT) newdat_DT$predsmod.1 <- predict(mod.1, newdat_DT, level=0) head(newdat_DT) ## transform DT num into Posix newdat_DT$predsmod.1Posix <- as.POSIXct(as.numeric(newdat_DT$predsmod.1+1e-6), origin = "1970-01-01", tz="GMT") CI_DT <- bernr::bolker_ci(mod.1, newdat_DT, pred_int = FALSE, conf_level = 0.95) CI_DT$predsmod.1Posix <- as.POSIXct(as.numeric(CI_DT$predsmod.1+1e-6), origin = "1970-01-01", tz="GMT") head(CI_DT) # #### Plot predicts ##### gDT <- newdat_DT %>% ggplot(.) + geom_point(data=datDTnoNA, aes(y=DTSecs, x=days_since_rain), colour='grey', alpha=.5) + geom_point(aes(y=predsmod.1Posix, x=days_since_rain), colour='black', alpha=.7) + labs(y='Delay time to ignition (sec)', x='Days after rain') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + scale_color_continuous(name='Gap area (m2)') + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gDT ## why not only gap center and larger gap to predicts? ggsave(gDT, file='../Figs/gDT.png', width = 8, height = 5) #### figure with log x #### gDT_log <- ggplot(datDTnoNA) + geom_point(aes(x=log(days_since_rain), y=DTSecs), alpha=.5, colour='grey') + geom_line(data=CI_DT, aes(x=log(days_since_rain), y=predsmod.1Posix), linetype="solid", size=1, alpha=0.7, col='black') + geom_ribbon(data=CI_DT, aes(x=log(days_since_rain), y=predsmod.1Posix, ymin=ci_l, ymax=ci_h), alpha=0.2) + labs(y='Delay time to ignition (sec)', x='Log transformed number of days after rain') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gDT_log ggsave(gDT_log, file='../Figs/DT_log.png', width = 8, height = 6) #### save Ignitability RData #### save.image('DT_analysis_results_Oct2023.RData') # ─── Combustability --------------------------------------------- rm(list = ls()) setwd('H:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/Combustability') # OVEN: weight before - weight after = evaporated # BURN: weight before - weight after = weight loss during the burn (wldb) # wldb - evaporated = mass lost # if positive, material was lost; if negative, only water evaporated. ## flame is always positive. IF "mass loss during burn" is positive, material was lost. IF negative, only water was burnt. ## Combustibility is mass loss / flame duration, so it shouldn't be negative, right? Negative were forced to zero (zero mass was consumed) # flame data from the fire experiment load("dat_flammability_Jan2023.RData") glimpse(dat_flammnoNA) range(dat_flammnoNA$newStudyPoint) dat_flammnoNA$days_since_rain_burn[dat_flammnoNA$days_since_rain_burn%in%c(0.1)] <- 1 dat_flammnoNA$newStudyPoint[dat_flammnoNA$newStudyPoint%in%c(0.1)] <- 1 datCombust <- dat_flammnoNA %>% filter(days_since_rain_burn <= 40) glimpse(datCombust) range(datCombust$newStudyPoint) rm(dat_flammnoNA) datCombust %>% # filter(mass_lost_burn>=0) %>% filter(days_since_rain_burn <= 40) %>% ggplot(.) + geom_point(aes(y=combust_SumFlame , x=days_since_rain_burn, colour=newStudyPoint)) + geom_smooth(aes(y=combust_SumFlame , x=days_since_rain_burn)) + labs(y='Combustion rate (g/s)', x='Days after rain') #### Comb Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } head(datCombust) ## mod.1 <- lme(combust_SumFlame ~ days_since_rain_burn + newStudyPoint + Area, random = list(~1|date_collectedG, ~1|GapID), data = datCombust, method = 'ML', control=lmeControl(msMaxIter = 30, tolerance = 1e-04, niterEM = 25, msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) res(mod.1) ## summary(mod.1) # ## null model mod.1.NULL <- lme(combust_SumFlame ~ 1, random = list(~1|date_collectedG, ~1|GapID), data = datCombust, method = 'ML', control=lmeControl(msMaxIter = 30, tolerance = 1e-04, niterEM = 25, msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) res(mod.1.NULL) summary(mod.1.NULL) # anova(mod.1, mod.1.NULL) ## they are not different, so i can't reject the null # ## interaction log(days since rain) * distance # mod.2 <- lme(combust_SumFlame ~ days_since_rain_burn*newStudyPoint + Area, # random = list(~1|date_collectedG, ~1|GapID), # data = datCombust, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.2) ## # summary(mod.2) # # # ## comparing mod 1 and mod 2 # anova(mod.2, mod.1) # not sig, stick with mod.1 # # # ## interaction log(days since rain) * area # mod.3 <- lme(combust_SumFlame ~ days_since_rain_burn*Area + newStudyPoint, # random = list(~1|date_collectedG, ~1|GapID), # data = datCombust, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.3) ## # summary(mod.3) ## none sig # # # comparing mod 1 and mod 3 # anova(mod.3, mod.1) # not sig, mod.1 # # # ## ## Predicts ## range(datCombust$newStudyPoint) newdat_FT <- expand.grid(days_since_rain_burn = seq(min(datCombust$days_since_rain_burn), max(datCombust$days_since_rain_burn), length.out = 1000), Area = max(datCombust$Area), newStudyPoint = min(datCombust$newStudyPoint)) # = centre head(newdat_FT) newdat_FT$predsmod.1.NULL <- predict(mod.1.NULL, newdat_FT, level=0) head(newdat_FT) #### Plot predicts ##### gComb <- newdat_FT %>% ggplot(.) + geom_point(data=datCombust, aes(y=combust_SumFlame, x=days_since_rain_burn),colour='grey') + geom_line(aes(y=predsmod.1.NULL, x=days_since_rain_burn), size=1.5, col='black') + labs(y='Combustion rate (g/sec)', x='Days after rain') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + scale_color_continuous(name='Gap area (m2)') + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gComb ggsave(gComb, file='../Figs/gComb.png', width = 8, height = 5) #### save Combustability RData #### save.image(file='FT_analysis_results_Oct2023.RData') # ─── Sustainability --------------------------------------------- rm(list = ls()) setwd('H:/My Drive/PhD_Sheffield/Research_thesis/3rd chapter Fire experim/AfterThesis_forPub/Analysis_Done_in_2023/Sustainability') # data load("datBTnoNA.RData") glimpse(datBTnoNA) summary(datBTnoNA) range(datBTnoNA$days_since_rain) datBTnoNA$days_since_rain[datBTnoNA$days_since_rain%in%c(0.1)] <- 1 datBTnoNA$newStudyPoint[datBTnoNA$newStudyPoint%in%c(0.1)] <- 1 #### stop the experiment at 40 days after rain because lack of data for all the range of gap sizes range(datBTnoNA$days_since_rain) # 1 - 72 datBTnoNA <- datBTnoNA %>% filter(days_since_rain <= 40) range(datBTnoNA$newStudyPoint) ## graphs explor ## range(datBTnoNA$BTPosix) gDays <- datBTnoNA %>% ggplot(.,aes(y=BTPosix, x=days_since_rain)) + geom_point(aes(colour=Area)) + geom_smooth() gDist <- datBTnoNA %>% ggplot(.,aes(y=BTPosix, x=newStudyPoint)) + geom_point(aes(colour=Area)) + geom_smooth() + labs(x='Distance from gap center') gArea <- datBTnoNA %>% ggplot(.) + geom_point(aes(y=BTPosix, x=Area, colour=days_since_rain)) + geom_smooth(aes(y=BTPosix, x=Area)) # gridExtra::grid.arrange(gSizeCat, gDays, gDist, gArea, ncol=1) ### Create a column for flame duration in seconds glimpse(datBTnoNA) datBTnoNA$BTSecsWithDate <- lubridate::time_length(datBTnoNA$BTnum, unit='second') # g <- min(datBTnoNA$BTSecsWithDate) asa <- 1559952000 ## zero seconds as.POSIXct(as.numeric(asa+1e-6), origin = "1970-01-01", tz="GMT") # zero seconds # Remove the min to get rid of all the date part datBTnoNA <- datBTnoNA %>% mutate(BTSecs = BTSecsWithDate - asa) glimpse(datBTnoNA) ## So: x= is ~ days since rain log ? #### Sustainability Analysis #### res <- function(mod){ a <- par(mfrow=c(2,2)) plot(fitted(mod), residuals(mod)) qqnorm(residuals(mod)) qqline(residuals(mod)) acf(residuals(mod)) pacf(residuals(mod)) b <- AIC(mod) print(b) } ## delay time to ignition ~ days since rain ## Q: Do distance from gap and gap area affect the delay time to ignition? ## mod.1 <- lme(BTSecs ~ log(days_since_rain) + newStudyPoint + Area, random = list(~1|DateG, ~1|GapID), data = datBTnoNA, method = 'ML', control=lmeControl(msMaxIter = 30, tolerance = 1e-04, niterEM = 25, msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) res(mod.1) ## summary(mod.1) # distance and area are not sig. # ## interaction log(days since rain) * distance # mod.2 <- lme(BTSecs ~ log(days_since_rain)*newStudyPoint + Area, # random = list(~1|DateG, ~1|GapID), # data = datBTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # # res(mod.2) ## # summary(mod.2) ## inter is not sig! # # ## comparing mod 1 and mod 2 # anova(mod.2, mod.1) # not sig , go for simple mod.1 # # # ## interaction log(days since rain) * area # mod.3 <- lme(BTSecs ~ log(days_since_rain)*Area + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datBTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.3) ## # summary(mod.3) ## inter is not sig # # # comparing mod 1 and mod 3 # anova(mod.3, mod.1) # not sig, go for simple mod 1 # # ## two way interaction # mod.4 <- lme(BTSecs ~ log(days_since_rain)*Area + log(days_since_rain)*newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datBTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.4) ## 6924.709 # summary(mod.4) ## nothing is sig, except days since rain # # # comparing mod 1 and mod 3 # anova(mod.4, mod.1) # not sig, go for simple mod 1 # # # ## using categorical Size Cat instead of numeric Area ## # mod.5 <- lme(BTSecs ~ log(days_since_rain) + newStudyPoint + SizeCat, # random = list(~1|DateG, ~1|GapID), # data = datBTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.5) ## 6913.625 # summary(mod.5) # distance is not sig. Gap size cat and days are sig. # # ## interaction with gap size # mod.6 <- lme(BTSecs ~ log(days_since_rain)*SizeCat + newStudyPoint, # random = list(~1|DateG, ~1|GapID), # data = datBTnoNA, # method = 'ML', # control=lmeControl(msMaxIter = 30, # tolerance = 1e-04, niterEM = 25, # msMaxEval = 50, msTol = 1e-04,msVerbose = TRUE, opt= "optim")) # res(mod.6) ## 6913.699 # summary(mod.6) # nothing is sig. geez # # # anova(mod.6, mod.5) # not sig, keep simple ## Predicts with Area ## newdat_BT <- expand.grid(days_since_rain = seq(min(datBTnoNA$days_since_rain), max(datBTnoNA$days_since_rain), length.out = 1000), Area = max(datBTnoNA$Area), newStudyPoint = min(datBTnoNA$newStudyPoint)) # = centre head(newdat_BT) newdat_BT$predsmod.1 <- predict(mod.1, newdat_BT, level=0) CI_BT <- bernr::bolker_ci(mod.1, newdat_BT, pred_int = FALSE, conf_level = 0.95) CI_BT$predsmod.1Posix <- as.POSIXct(as.numeric(CI_BT$predsmod.1+1e-6), origin = "1970-01-01", tz="GMT") head(CI_BT) # ## transform BT num into Posix # newdat_BT$predsmod.1Posix <- as.POSIXct(as.numeric(newdat_BT$predsmod.1+1e-6), origin = "1970-01-01", tz="GMT") # head(newdat_BT) #### Plot predicts ##### # glimpse(datBTnoNA) # gBT <- newdat_BT %>% ggplot(.) + # geom_point(data=datBTnoNA, aes(y=BTSecs, x=days_since_rain), col='grey') + # geom_point(aes(y=predsmod.1Posix, x=days_since_rain)) + # scale_y_continuous(breaks=seq(0,2000,50)) + # scale_y_time() + # breaks = c('300','600','900') # labs(y='Burning time', x='Days after rain') + # theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + # scale_color_continuous(name='Gap area (m2)') + # theme(legend.box.background = element_rect(fill = "white", colour = "black"), # legend.box.margin = margin(4, 4, 4, 4), # legend.text = element_text(size=10), # legend.title = element_text(size=10)) # gBT # ggsave(gBT, file='../Figs/gBT.png', width = 8, height = 5) #### figure with log x gBT_log <- ggplot(datBTnoNA) + geom_point(aes(x=log(days_since_rain), y=BTSecs), alpha=.5, colour='grey') + geom_line(data=CI_BT, aes(x=log(days_since_rain), y=predsmod.1Posix), linetype="solid", size=1, alpha=0.7, col='black') + geom_ribbon(data=CI_BT, aes(x=log(days_since_rain), y=predsmod.1Posix, ymin=ci_l, ymax=ci_h), alpha=0.2) + # geom_line(data=newdat_BT, aes(x=log(days_since_rain), y=predsmod.1Posix), linetype="solid", alpha=0.7, col='black') + labs(y='Delay time to ignition (sec)', x='Log transformed number of days after rain') + theme_bw(base_size = 18) + theme(legend.position = c(0.8,0.8)) + theme(legend.box.background = element_rect(fill = "white", colour = "black"), legend.box.margin = margin(4, 4, 4, 4), legend.text = element_text(size=10), legend.title = element_text(size=10)) gBT_log #### save Sustainability RData #### save.image(file='BT_analysis_results_Oct2023.RData')