############################################################### ##do communities respond differently to subsequent burns####### ##based on the recency of the last burn?####################### ############################################################### #scale_color_manual(values = c("skyblue4", "gold", "chocolate1", "brown", "gray28"))+ # scale_fill_manual(values = alpha(c("skyblue4", "gold", "chocolate1", "brown", "gray28"), 0.5))+ ##########load data############### source("../R_code/ftbr_2020-03-18_read-in-data.R") #we only have cover data beginning in 2011 burns = burns[BurnYear > 2010] data5 = data1 data5$Burns = data5$numburns #what to do w plots burned twice between cover samples? remove those years x = as.data.frame(table(data1[,c("UID", "Year")])) x=x[x$Freq == 2,] data5[, tag := "tag"] for (i in 1:dim(x)[1]){ data5[UID == x$UID[i] & Year == x$Year[i], tag := NA_character_] } data5 = data5[!is.na(tag)] ####################################################### ##########change in species richness################### ####################################################### ####delta-richness over a burn data5[,deltaS := as.numeric(NA_integer_)] for (i in 1:dim(burns)[1]){ plot = burns[i, Plot] year = burns[i, BurnYear] if (!isEmpty(data5[UID == plot & Year == year, richness]) & !isEmpty(data5[UID == plot & Year == year-1, richness])){ data5[UID == plot & Year == year, deltaS := (data5[UID == plot & Year == year, richness] - data5[UID == plot & Year == year-1, richness])] } } data6 = data5[!is.na(deltaS),] data6 = data6[!is.na(FRI),] data6$Burns = as.factor(data6$Burns) #assessing fits fit_noint = lm(deltaS ~ FRI+DSB, data = data6) fit_FRI = lm(deltaS ~ FRI, data = data6) fit_all = lm(deltaS ~ FRI*DSB, data = data6) fit_dsb = lm(deltaS ~ DSB, data = data6) AIC(fit_all) #969 AIC(fit_noint) #994 AIC(fit_FRI) #1003 AIC(fit_dsb) #1001 summary(fit_all) car::Anova(fit_all, type = 3) #Call: #lm(formula = deltaS ~ DSB * FRI, data = data6) # #Residuals: # Min 1Q Median 3Q Max #-20.0552 -2.2606 0.2411 2.5033 13.7974 # #Coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) 2.717e+00 1.287e+00 2.111 0.0363 * # DSB -4.662e-02 1.517e-02 -3.073 0.0025 ** # FRI -9.965e-03 1.586e-03 -6.282 3.02e-09 *** # DSB:FRI 1.172e-04 2.194e-05 5.341 3.13e-07 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # #Residual standard error: 4.563 on 160 degrees of freedom #Multiple R-squared: 0.2645, Adjusted R-squared: 0.2507 #F-statistic: 19.18 on 3 and 160 DF, p-value: 1.127e-10 # ###visualize #png(paste0(Sys.Date(), "_lm-richness.png"), width = 8, height = 6, units = "in", res = 300) richness_lm = interact_plot(fit_all, FRI, modx = DSB, plot.points = TRUE, point.alpha = 0.6, jitter = c(20, 2), interval = TRUE, legend.main = "DSB", x.label = "Fire Return Interval", y.label = "Change in Richness", colors = c("darkred", "gray70", "gray10"))+ theme_classic()+ggtitle("(A)")+ theme(plot.title = element_text(size = 20))+ theme(panel.border = element_rect(color = "black", fill = NA))+ theme(axis.title = element_text(size = 14))+ theme(legend.title = element_text(size = 14), legend.text = element_text(size = 12))+ theme(axis.text = element_text(size = 12))+ scale_x_continuous(breaks = c(seq(0,2500,500)), expand = c(0,0))#+ #geom_hline(yintercept = 0, color = "black", size = 0.5) #dev.off() #Sanova = aov(deltaS ~ as.factor(Burns), data = data6) #summary(Sanova) #TukeyHSD(Sanova, conf.level=.95) ###################################################### ##########change in % veg cover ###################### ###################################################### ####delta-pctcover over a burn data5[,deltaPC := as.numeric(NA_integer_)] for (i in 1:dim(burns)[1]){ plot = burns[i, Plot] year = burns[i, BurnYear] if (!isEmpty(data5[UID == plot & Year == year, pctcover]) & !isEmpty(data5[UID == plot & Year == year-1, pctcover])){ data5[UID == plot & Year == year, deltaPC := (data5[UID == plot & Year == year, pctcover] - data5[UID == plot & Year == year-1, pctcover])] } } data8 = data5[!is.na(deltaPC),] data8 = data8[!is.na(FRI),] data8$Burns = as.factor(data8$Burns) fit_noint_PC = lm(deltaPC ~ DSB+FRI, data = data8) fit_FRI_PC = lm(deltaPC ~ FRI, data = data8) fit_dsb_PC = lm(deltaPC ~ DSB, data = data8) fit_all_PC = lm(deltaPC ~ FRI*DSB, data = data8) AIC(fit_all_PC) #2114 AIC(fit_noint_PC) #2126 AIC(fit_FRI_PC) #2135 AIC(fit_dsb_PC) #2139 summary(fit_all_PC) #adj R^2 = 0.21 car::Anova(fit_all_PC, type = 3) # #Call: # lm(formula = deltaPC ~ DSB * FRI, data = data8) # #Residuals: # Min 1Q Median 3Q Max #-487.80 -68.41 16.33 84.32 372.25 # #Coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) 74.8643662 46.5906674 1.607 0.110232 #DSB -1.0943116 0.5408822 -2.023 0.044865 * # FRI -0.3274240 0.0619790 -5.283 4.5e-07 *** # DSB:FRI 0.0031784 0.0008287 3.835 0.000185 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # #Residual standard error: 151.5 on 147 degrees of freedom #Multiple R-squared: 0.2207, Adjusted R-squared: 0.2048 #F-statistic: 13.88 on 3 and 147 DF, p-value: 5.119e-08 #visualize PC_lm = interact_plot(fit_all_PC, FRI, modx = DSB, plot.points = TRUE, point.alpha = 0.6, jitter = c(10, 10), interval = TRUE, legend.main = "DSB", x.label = "Fire Return Interval", y.label = "Change in % Veg Cover", colors = c("darkred", "gray70", "gray10"))+ theme_classic()+ggtitle("(B)")+ theme(plot.title = element_text(size = 20))+ theme(panel.border = element_rect(color = "black", fill = NA))+ theme(axis.title = element_text(size = 14))+ theme(axis.text = element_text(size = 12))+ theme(legend.position = "none")+ scale_x_continuous(breaks = c(seq(0,2500,500)), expand = c(0,0))#+ #geom_hline(yintercept = 0, color = "black", size = 0.5) #lay = rbind(c(1, 1, 1,1,2,2, 2), c(1, 1, 1,1, 2,2, 2)) #png(paste0(Sys.Date(), "_Fig4.png"), width = 7, height = 3.5, units="in", res=300) #grid.arrange(richness_lm, PC_lm, layout_matrix = lay) #dev.off()