Effect of seasonality of fire and duration of rest from grazing on mountain big sage plant communities
DAM
2022-04-19
Fire effects on community composition
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Deferment | 1 | 0.05 | 0.01 | 0.78 | 0.47 |
| BurnSeason | 2 | 1.64 | 0.22 | 14.04 | 0.01 |
| TSF | 1 | 1.07 | 0.14 | 18.28 | 0.01 |
| Deferment:BurnSeason | 2 | 0.25 | 0.03 | 2.17 | 0.04 |
| Deferment:TSF | 1 | 0.01 | 0 | 0.1 | 0.99 |
| BurnSeason:TSF | 2 | 0.38 | 0.05 | 3.21 | 0.01 |
| Deferment:BurnSeason:TSF | 2 | -0.01 | 0 | -0.08 | 1 |
| Residual | 69 | 4.04 | 0.54 | NA | NA |
| Total | 80 | 7.43 | 1 | NA | NA |
Take-aways:
- No apparent effect of grazing deferment. Dropping from subsequent plots and analysis
- Significant effect of Time since fire and interaction between Burn season and TSF. Subsequent analyses parse community comp differences among fire treatments within each year separately.
Ordination results of species-level communities in response to two seasons of fire treatment vs. unburned controls presented as dissimilarity among treatment groups in the year prior to fire treatments, one year after fire, and four years after fire. Blue ‘+’ signs in each panel of the top row correspond to species scores as in bottom panel.
| year | Pairwise comparison | Fstat | P |
|---|---|---|---|
| 2008 | FB_vs_NB | 3.422 | 0.01 |
| 2008 | FB_vs_SB | 2.628 | 0.023 |
| 2008 | SB_vs_NB | 0.7826 | 0.625 |
| 2009 | FB_vs_NB | 14.71 | 0.001 |
| 2009 | FB_vs_SB | 2.446 | 0.026 |
| 2009 | SB_vs_NB | 8.375 | 0.001 |
| 2012 | FB_vs_NB | 22.27 | 0.001 |
| 2012 | FB_vs_SB | 3.708 | 0.001 |
| 2012 | SB_vs_NB | 8.614 | 0.001 |
Fire effects by functional groups
Regression results
Whether grazers were returned to burned pastures in the season after fire or deferred for 1 or 2 years did not affect most plant functional group responses either 1 yr or 4 yr after burning. The deferment term was only statistically significant for litter cover 4 yrs after fire, with lower litter accumulation where grazing had not been deferred. In no case, however, did mean litter cover fall below 75%.
Few functional groups showed statistically-significant differences in cover 1 yr or 4 yr after fire when compared to pre-burn levels:
One year after fire, AGI was lower in fall burn plots; ASH was higher in both spring and fall burn plots; FN increased in all plots, Litter declined in fall burn plots; PGI increased in spring burn plots; PGN decreased in fall burn plots; and PSN declined in both spring and fall burn plots.
Four years after fire, Ash cover remained elevated under both burn seasons; FN declined in unburned and spring burn plots, but not fall burn plots; Litter was lower in all plots; PGI increased in spring burn plots; PGI increased in spring and fall burn plots; PSN declined in both spring and fall burn plots; and PSSN declined without fire but remained constant in spring and fall burned plots.
| Time since fire | Functional group | Deferment term | Burn season | Change in cover (± pooled S.E.) | t test |
|---|---|---|---|---|---|
| 1 yr | AGI | \(\chi^{2}\) = 1.73, P = 0.42 | |||
| 1 yr | AGI | Unburned | -0.04 ± 0.14 | t = -0.27, P = 0.79 | |
| 1 yr | AGI | Spring | -0.05 ± 0.14 | t = -0.36, P = 0.73 | |
| 1 yr | AGI | Fall | -0.3 ± 0.14 | t = -2.14, P = 0.04 | |
| 1 yr | Ash | \(\chi^{2}\) = 1.7, P = 0.43 | |||
| 1 yr | Ash | Unburned | 0 ± 4.66 | t = 0, P = 1 | |
| 1 yr | Ash | Spring | 45.87 ± 4.66 | t = 9.85, P = 0 | |
| 1 yr | Ash | Fall | 54.14 ± 4.66 | t = 11.62, P = 0 | |
| 1 yr | Dead | \(\chi^{2}\) = 1.73, P = 0.42 | |||
| 1 yr | Dead | Unburned | 0.37 ± 0.59 | t = 0.62, P = 0.54 | |
| 1 yr | Dead | Spring | 2.78 ± 0.59 | t = 4.69, P = 0 | |
| 1 yr | Dead | Fall | 3.49 ± 0.59 | t = 5.88, P = 0 | |
| 1 yr | FI | \(\chi^{2}\) = 1.47, P = 0.48 | |||
| 1 yr | FI | Unburned | -0.15 ± 0.09 | t = -1.73, P = 0.09 | |
| 1 yr | FI | Spring | 0.06 ± 0.09 | t = 0.75, P = 0.46 | |
| 1 yr | FI | Fall | -0.01 ± 0.09 | t = -0.1, P = 0.92 | |
| 1 yr | FN | \(\chi^{2}\) = 4.45, P = 0.11 | |||
| 1 yr | FN | Unburned | 6.93 ± 1.5 | t = 4.63, P = 0 | |
| 1 yr | FN | Spring | 5.72 ± 1.5 | t = 3.83, P = 0 | |
| 1 yr | FN | Fall | 5.65 ± 1.5 | t = 3.78, P = 0 | |
| 1 yr | Litter | \(\chi^{2}\) = 1.41, P = 0.5 | |||
| 1 yr | Litter | Unburned | -1.44 ± 1.24 | t = -1.16, P = 0.27 | |
| 1 yr | Litter | Spring | -2.02 ± 1.24 | t = -1.63, P = 0.14 | |
| 1 yr | Litter | Fall | -8.2 ± 1.24 | t = -6.61, P = 0 | |
| 1 yr | PGI | \(\chi^{2}\) = 0.19, P = 0.91 | |||
| 1 yr | PGI | Unburned | 0.07 ± 0.18 | t = 0.41, P = 0.68 | |
| 1 yr | PGI | Spring | 0.53 ± 0.18 | t = 2.94, P = 0.01 | |
| 1 yr | PGI | Fall | 0.25 ± 0.18 | t = 1.39, P = 0.18 | |
| 1 yr | PGN | \(\chi^{2}\) = 1.6, P = 0.45 | |||
| 1 yr | PGN | Unburned | -1.04 ± 1.39 | t = -0.75, P = 0.46 | |
| 1 yr | PGN | Spring | -3.84 ± 1.39 | t = -2.77, P = 0.01 | |
| 1 yr | PGN | Fall | -0.32 ± 1.39 | t = -0.23, P = 0.82 | |
| 1 yr | PSN | \(\chi^{2}\) = 3.37, P = 0.19 | |||
| 1 yr | PSN | Unburned | 3.37 ± 2.51 | t = 1.34, P = 0.21 | |
| 1 yr | PSN | Spring | -23.05 ± 2.51 | t = -9.18, P = 0 | |
| 1 yr | PSN | Fall | -29.8 ± 2.51 | t = -11.87, P = 0 | |
| 1 yr | PSSN | \(\chi^{2}\) = 3.53, P = 0.17 | |||
| 1 yr | PSSN | Unburned | 0.96 ± 0.59 | t = 1.64, P = 0.11 | |
| 1 yr | PSSN | Spring | -0.65 ± 0.59 | t = -1.1, P = 0.28 | |
| 1 yr | PSSN | Fall | -0.21 ± 0.59 | t = -0.36, P = 0.72 | |
| 4 yr | AGI | \(\chi^{2}\) = 2.36, P = 0.31 | |||
| 4 yr | AGI | Unburned | 0 ± 0.13 | t = 0, P = 1 | |
| 4 yr | AGI | Spring | 0.12 ± 0.13 | t = 0.87, P = 0.39 | |
| 4 yr | AGI | Fall | 0.19 ± 0.13 | t = 1.46, P = 0.16 | |
| 4 yr | Ash | \(\chi^{2}\) = 0.47, P = 0.79 | |||
| 4 yr | Ash | Unburned | 0 ± 3.55 | t = 0, P = 1 | |
| 4 yr | Ash | Spring | 18.77 ± 3.55 | t = 5.28, P = 0 | |
| 4 yr | Ash | Fall | 16.47 ± 3.55 | t = 4.64, P = 0 | |
| 4 yr | Dead | \(\chi^{2}\) = 4.78, P = 0.09 | |||
| 4 yr | Dead | Unburned | 0.07 ± 0.43 | t = 0.17, P = 0.87 | |
| 4 yr | Dead | Spring | 1.66 ± 0.43 | t = 3.83, P = 0 | |
| 4 yr | Dead | Fall | 1.97 ± 0.43 | t = 4.52, P = 0 | |
| 4 yr | FI | \(\chi^{2}\) = 2.47, P = 0.29 | |||
| 4 yr | FI | Unburned | 0.07 ± 0.1 | t = 0.73, P = 0.47 | |
| 4 yr | FI | Spring | 0.14 ± 0.1 | t = 1.35, P = 0.19 | |
| 4 yr | FI | Fall | 0.1 ± 0.1 | t = 1.01, P = 0.32 | |
| 4 yr | FN | \(\chi^{2}\) = 5.8, P = 0.06 | |||
| 4 yr | FN | Unburned | -6.07 ± 1.28 | t = -4.73, P = 0 | |
| 4 yr | FN | Spring | -5.36 ± 1.28 | t = -4.18, P = 0 | |
| 4 yr | FN | Fall | -1.18 ± 1.28 | t = -0.92, P = 0.37 | |
| 4 yr | Litter | \(\chi^{2}\) = 10.57, P = 0.01 | |||
| 4 yr | Litter | Unburned | -4.56 ± 1.18 | t = -3.87, P = 0 | |
| 4 yr | Litter | Spring | -5.63 ± 1.18 | t = -4.78, P = 0 | |
| 4 yr | Litter | Fall | -10.65 ± 1.18 | t = -9.05, P = 0 | |
| 4 yr | PGI | \(\chi^{2}\) = 2.42, P = 0.3 | |||
| 4 yr | PGI | Unburned | 0.04 ± 0.32 | t = 0.12, P = 0.91 | |
| 4 yr | PGI | Spring | 1.23 ± 0.32 | t = 3.86, P = 0 | |
| 4 yr | PGI | Fall | 0.36 ± 0.32 | t = 1.13, P = 0.29 | |
| 4 yr | PGN | \(\chi^{2}\) = 4.44, P = 0.11 | |||
| 4 yr | PGN | Unburned | -1.04 ± 1.77 | t = -0.59, P = 0.56 | |
| 4 yr | PGN | Spring | 5.23 ± 1.77 | t = 2.96, P = 0.01 | |
| 4 yr | PGN | Fall | 3.84 ± 1.77 | t = 2.17, P = 0.04 | |
| 4 yr | PSN | \(\chi^{2}\) = 3.21, P = 0.2 | |||
| 4 yr | PSN | Unburned | -3.56 ± 2.41 | t = -1.48, P = 0.17 | |
| 4 yr | PSN | Spring | -23.85 ± 2.41 | t = -9.92, P = 0 | |
| 4 yr | PSN | Fall | -30.37 ± 2.41 | t = -12.63, P = 0 | |
| 4 yr | PSSN | \(\chi^{2}\) = 5.78, P = 0.06 | |||
| 4 yr | PSSN | Unburned | -1.33 ± 0.59 | t = -2.25, P = 0.03 | |
| 4 yr | PSSN | Spring | 0.15 ± 0.59 | t = 0.25, P = 0.81 | |
| 4 yr | PSSN | Fall | -0.31 ± 0.59 | t = -0.52, P = 0.61 |
Script
knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)
# knitr::opts_chunk$set(fig.width=unit(15,"cm"), fig.height=unit(10,"cm"))
##Libraries
pacman::p_load(tidyverse)
pacman::p_load_gh("devanmcg/wesanderson")
load('./Robjects/nmds_scores.Rdata') # Extracted scores for NMDS
load('./Robjects/ad_result.Rdata') # Global adonis2 results
load("./Robjects/pwad.Rdata") # Pared-down results of pairwise adonis2 by yr
load('./Robjects/FoliarGroupLong2.Rdata') # Functional group data from wrangling file
load('./Robjects/fg_reg_res.Rdata') # Functional group regression results
#Precipitation
#Import precip data
PrecipCM <-
read_csv("DuboisPrecipInches.csv") %>% #Precip data from Sheep Station HQ
filter(Year<2021) %>%
rowwise %>%
mutate(GrowingSeasonTotal = sum(c_across(c(May, Jun, July)))) %>%
dplyr::select(Year, Total, GrowingSeasonTotal) %>%
mutate(across(c(Total, GrowingSeasonTotal), ~.x*2.54))
#compute mean objects
PTT<-mean(PrecipCM$Total)
PTTGS<-mean(PrecipCM$GrowingSeasonTotal)
Pt<-PrecipCM%>%filter(between(Year, 2007,2013))
PTGSSSt<-mean(Pt$GrowingSeasonTotal)
PTTotalSt<-mean(Pt$Total)
#Plot Precip Data
p2<-PrecipCM%>%
filter(between(Year, 2007,2013)) %>%
ggplot( aes(x=Year,
y=Total,
fill=as.numeric(Total)))+
geom_line(stat = "identity", size=1)+
geom_line(aes(y=GrowingSeasonTotal), stat="identity", color="grey", size=1)+
geom_hline(yintercept=PTTotalSt, linetype="dotted", color = "black", size=.5)+ #Totalaverage for study period
geom_hline(yintercept=PTT, linetype="dotted", color="grey")+ #total average for period of record
geom_hline(yintercept=PTGSSSt, linetype="dashed", color = "black", size=.5)+ #gs mean for study period
geom_hline(yintercept=PTTGS, linetype="dashed", color = "grey")+ #gs average for period of record
theme_classic(16)
#add for markdown: theme(plot.caption=element_markdown())
p2
###
### not run:
###
# External functions
source("https://raw.githubusercontent.com/devanmcg/IntroRangeR/master/11_IntroMultivariate/ordinationsGGplot.R")
source("https://raw.githubusercontent.com/pmartinezarbizu/pairwiseAdonis/master/pairwiseAdonis/R/pairwise.adonis2.R")
# Data wrangling script
source('DataWranglingDAM.R')
# Wrangle community data
FoliarWide <-
FoliarSpWideCover4 %>%
ungroup() %>%
select(-X, -starts_with("DEAD")) %>%
filter(Year %in% c('2008', '2009', '2012')) %>%
rownames_to_column("rowID") %>%
rename(Deferment = YrsRest,
BurnSeason = BurnTrt,
TSF = Year)
veg <-
FoliarWide %>%
select(rowID, AF01:THIN6) %>%
labdsv::dropspc(minocc = 8,
minabu = 2)
# Fit ordination
veg_nmds <- metaMDS(select(veg, -rowID), 'bray', 4, trace = F)
# save(veg_nmds, file = './Robjects/veg_nmds.Rdata')
# Test groups
# set up permutation blocks (nee strata)
perm = how(nperm = 99)
setBlocks(perm) <- with(FoliarWide, MainPlot)
ad_result <- adonis2(select(veg, - rowID) ~ Deferment*BurnSeason*TSF,
data = FoliarWide, permutations = perm, by = 'term')
# save(ad_result, file = './Robjects/ad_result.Rdata')
# Due to significant Trt x Year interaction,
# perform pairwise tests of treatment within each year separately
yrs <- c(2008, 2009,2012)
pwad_out <- lst()
for(y in 1:length(yrs)){
t <- filter(filter(FoliarWide, TSF == yrs[y]))
v <- filter(veg, rowID %in% t$rowID) %>% select(-rowID)
pwad_out[[paste(yrs[y])]] <- pairwise.adonis2(v ~ BurnSeason, data = t)# strata not working
}
# Collapse list of pairwise results into something more useful
pwad <- tibble()
for(y in 1:length(yrs)){
l = pwad_out[[paste(yrs[y])]][2:4]
data.table::rbindlist(l, idcol = 'l') %>%
as_tibble() %>%
mutate(year = yrs[y]) %>%
rename(comp = l,
P = `Pr(>F)`,
Fstat = `F`) %>%
group_by(comp) %>%
slice(1)%>%
ungroup() %>%
select(year, comp, Fstat, P) %>%
bind_rows(pwad, .) -> pwad
}
# save(pwad, file = "./Robjects/pwad.Rdata")
nmds_spp <-
scores(veg_nmds, "species") %>%
as.data.frame() %>%
as_tibble(rownames = "Species")
# create list object & populate with scores
nmds_scores <- lst(species = nmds_spp)
nmds_scores$spiders <-
gg_ordiplot(veg_nmds,
groups = paste0(FoliarWide$TSF, '_', FoliarWide$BurnSeason),
plot=FALSE)$df_spiders %>%
as_tibble() %>%
rename(NMDS1 = x, NMDS2 = y) %>%
separate(Group, into=c('TSF', 'Treatment')) %>%
mutate(TSF = factor(TSF, levels = c('2008', '2009', '2012')),
TSF = recode(TSF,
'2008' = 'Pre-burn',
'2009' = '1 yr',
'2012' = '4 yr'),
Treatment = factor(Treatment, levels = c('NB', 'SB', 'FB')),
Treatment = recode(Treatment,
'NB'= 'Unburned',
'SB' = 'Spring',
'FB' = "Fall"))
# save(nmds_scores, file = './Robjects/nmds_scores.Rdata')
###
### end not run
###
ad_result %>%
mutate(across(c(SumOfSqs:`F`), ~round(.x, 2))) %>%
pander::pander("Global test of group effects via adonis2.")
# Plot ordination
# site scores by treatment and year
site_gg <-
ggplot() + theme_bw(14) +
labs(x="NMDS Axis 1",
y="NMDS Axis 2") +
geom_vline(xintercept = 0, lty=3, color="darkgrey") +
geom_hline(yintercept = 0, lty=3, color="darkgrey") +
geom_point(data=nmds_scores$species,
aes(x=NMDS1, y=NMDS2),
pch = '+',
size = 3,
colour="darkblue") +
geom_segment(data=nmds_scores$spiders,
aes(x=cntr.x, y=cntr.y,
xend=NMDS1, yend=NMDS2,
color=BurnSeason),
size=1.2) +
geom_point(data=nmds_scores$spiders,
aes(x=NMDS1, y=NMDS2,
bg=BurnSeason,
shape = BurnSeason),
colour="black",
size=3, stroke = 2) +
facet_wrap(~TSF, nrow = 1) +
scale_color_manual('Burn season', values = wes_palette("Zissou1")[c(1, 3, 5)]) +
scale_fill_manual('Burn season', values = wes_palette("Zissou1")[c(1, 3, 5)]) +
scale_shape_manual('Burn season', values = c(21, 24, 22) ) +
scale_alpha_manual('Burn season', values = c(0.5, 0.75, 1)) +
theme(panel.grid=element_blank(),
legend.position="top")
# species scores
spp_gg <-
ggplot() + theme_bw(14) +
labs(x="NMDS Axis 1",
y="NMDS Axis 2",
title = "Species scores") +
geom_vline(xintercept = 0, lty=3, color="darkgrey") +
geom_hline(yintercept = 0, lty=3, color="darkgrey") +
geom_text(data=nmds_scores$species,
aes(x=NMDS1, y=NMDS2,
label=Species),
colour="darkblue") +
theme(panel.grid=element_blank() )
b_row <- ggpubr::ggarrange(NULL, spp_gg, NULL, ncol = 3, widths = c(1,3,1))
ggpubr::ggarrange(site_gg, b_row,
ncol = 1, heights = c(1,1))
pwad %>%
rename(`Pairwise comparison` = comp) %>%
pander::pander("Pairwise adonis2 results by treatment within each year.")
# Functional groups
FoliarGroupLong2 %>%
filter(Year %in% c('2008', '2009', '2012')) %>%
rename(class = FoliarCoverTypeLite,
Deferment = YrsRest,
BurnSeason = BurnTrt) %>%
mutate(TSF = factor(Year, levels = c('2008', '2009', '2012')),
TSF = recode(TSF,
'2008' = 'Pre-\nburn',
'2009' = '1 yr',
'2012' = '4 yr'),
BurnSeason = factor(BurnSeason, levels = c('NB', 'SB', 'FB')),
BurnSeason = recode(BurnSeason,
'NB'= 'Unburned',
'SB' = 'Spring',
'FB' = "Fall"),
Deferment = factor(Deferment, levels = c('0', '1', '2')),
Deferment = recode(Deferment,
'0'= 'None',
'1' = '1 yr',
'2' = "2 yr")) %>%
group_by(BurnSeason, Deferment, Year, class) %>%
summarize(CoverMean = mean(cover, na.rm = T),
CoverSE = sd(cover, na.rm = T)/(sqrt(n())),
.groups = 'drop') %>%
ggplot(aes(x = Year)) + theme_bw(14) +
geom_line(aes(y = CoverMean,
color = BurnSeason,
group = interaction(Deferment, BurnSeason),
lty = Deferment),
position = position_dodge(width=0.2)) +
geom_errorbar(aes(ymin = CoverMean - CoverSE,
ymax = CoverMean + CoverSE,
color = BurnSeason),
width = 0.15,
position = position_dodge(width=0.2)) +
geom_point(aes(y = CoverMean,
fill = BurnSeason,
shape = BurnSeason),
position = position_dodge(width=0.2)) +
facet_wrap(~class) +
labs(x = "Years relative to fire",
y = "Percent cover (Mean ± SE)") +
scale_shape_manual('Burn season', values = c(21, 24, 22)) +
scale_color_manual('Burn season', values = wes_palette("Zissou1")[c(1, 3, 5)]) +
scale_fill_manual('Burn season', values = wes_palette("Zissou1")[c(1, 3, 5)]) +
scale_x_continuous(breaks = c(2008, 2009, 2012),
labels = c('Pre', '1', '4')) +
theme(legend.position = c(0.8,0.1),
legend.direction = 'vertical',
legend.box = "horizontal",
panel.grid.minor.x = element_blank())
# Assessing relative change in functional groups
# Calculate relative change
ResponsesLong <-
FoliarGroupLong2 %>%
filter(Year %in% c('2008', '2009', '2012')) %>%
rename(class = FoliarCoverTypeLite,
Deferment = YrsRest,
BurnSeason = BurnTrt) %>%
mutate(BurnSeason = factor(BurnSeason, levels = c('NB', 'SB', 'FB')),
BurnSeason = recode(BurnSeason,
'NB'= 'Unburned',
'SB' = 'Spring',
'FB' = "Fall"),
Deferment = factor(Deferment, levels = c('0', '1', '2')),
Deferment = recode(Deferment,
'0'= 'None',
'1' = '1 yr',
'2' = "2 yr")) %>%
mutate(tmp_chunks = stringr::str_split(Plot,
"(?<=[a-z])(?=[A-Z0-9])|(?<=[0-9])(?=[A-Za-z])",
n = 2)) %>%
mutate(Block = map_chr(tmp_chunks, 1),
Plot = map_chr(tmp_chunks, 2)) %>%
select(-tmp_chunks) %>%
pivot_wider(names_from = Year,
values_from = cover) %>%
group_by(Block, BurnSeason, Deferment, class) %>%
summarize(short = `2009` - `2008`,
long = `2012` - `2008`,
.groups = 'drop') %>%
pivot_longer(names_to = 'TSF',
values_to = 'cover',
cols = c(short, long))
# Chug regression models
{
begin = Sys.time() # start a computational run time counter
# set up parallel processing
pacman::p_load(foreach, doParallel)
cores=detectCores()
cl <- makeCluster(cores[1])
registerDoParallel(cl)
fg_reg_res <-
foreach(c = 1:length(unique(ResponsesLong$class)),
.combine = bind_rows,
.errorhandling = 'stop',
.inorder = FALSE) %:%
foreach(t = 1:length(unique(ResponsesLong$TSF)),
.combine = bind_rows,
.errorhandling = 'stop',
.inorder = FALSE,
.packages = c('tidyverse')) %dopar% {
# subset data for unique loop
cd <- ResponsesLong %>% filter(class == unique(ResponsesLong$class)[c],
TSF == unique(ResponsesLong$TSF)[t])
# Determine significance of deferment term
# fit competing models
m1 <- lmerTest::lmer(cover ~ 0 + BurnSeason + (1|Block),
data = cd, REML = F)
m2 <- lmerTest::lmer(cover ~ 0 + BurnSeason + Deferment + (1|Block),
data = cd, REML = F)
# Compare models with analysis of deviance and export results
def_stat <-
anova(m1, m2) %>%
select(-npar) %>%
broom::tidy() %>%
filter(term == 'm2') %>%
select(statistic, p.value) %>%
mutate_all(~round(.x, 2)) %>%
mutate(response = unique(ResponsesLong$class)[c],
TSF = unique(ResponsesLong$TSF)[t] ,
`Deferment term` = paste0('$\\chi^{2}$ = ', statistic, ', ',
'P = ', p.value)) %>%
select( -statistic, -p.value)
# Retrieve regression results for TSF t test against zero
s <- summary(m1)
s$coefficients %>%
as.data.frame() %>%
mutate_all(~round(.x, 2)) %>%
rownames_to_column('BurnSeason') %>%
mutate(response = unique(ResponsesLong$class)[c],
TSF = unique(ResponsesLong$TSF)[t] ,
BurnSeason = str_remove(BurnSeason, 'BurnSeason')) %>%
select(response, TSF, BurnSeason:`Pr(>|t|)`) %>%
bind_rows(def_stat, .)
}
stopCluster(cl)
Sys.time() - begin # return computational run time
}
# save(fg_reg_res, file = './Robjects/fg_reg_res.Rdata')
fg_reg_res %>%
mutate_all(~replace_na(.x, " ")) %>%
arrange(desc(TSF), response) %>%
mutate(TSF = recode(TSF,
short = '1 yr',
long = '4 yr')) %>%
unite(Mean, c(Estimate, `Std. Error`), sep = " ± ") %>%
mutate(Mean = ifelse(Mean == ' ± ', " ", Mean) ) %>%
rename(P = `Pr(>|t|)`,
t_stat = `t value`) %>%
mutate(t_stat = ifelse(t_stat != " ", paste0('t = ', t_stat, ', ', 'P = ', P), " ")) %>%
select(TSF, response, `Deferment term`, BurnSeason, Mean, t_stat) %>%
rename(`Change in cover (± pooled S.E.)` = Mean,
`Burn season` = BurnSeason,
`Time since fire` = TSF,
`Functional group` = response,
`t test` = t_stat) %>%
knitr::kable(caption="Statistical results of LME models testing change in cover for each functional group at 1 yr and 4 yr after fire. 'Deferment term' reports the significance of including grazing deferment period in the regression model with burn season. Subsequent t-statistics report the significance of the differences in estimates for each burn season's value against zero, testing the null hypothesis that functional group cover did not change in response to fire treatment.")