Rangeland response to fire and herbivory

Background

This project occurred during the summer of 2023 at the USDA Agricultural Research Services Fort Keogh Livestock and Range Research Laboratory, located in Miles City, MT. The purpose of this study is to determine the effects of fire on herbivory in the US Great Northern Plains.

Aboveground Biomass

Graphs

# Diagram showing methods

DiagrammeR::grViz("digraph {
  graph [center = TRUE, layout = dot, rankdir = TB]
  node [shape = rectangle]        
  rec1 [label = '1. Obtain clipping samples']
  rec2 [label = '2. Dry in oven for 48 hours']
  rec3 [label =  '3. Weigh samples']
  rec4 [label = '4. Analyze data in R Markdown']
  rec1 -> rec2 -> rec3 -> rec4
  }",
  height = 300)

# Load plot ID data

PlotID <-
  read_excel("LemonadeFireMasterData.xlsx",
             sheet = 'PlotID')

ID <- data.table(PlotID)
# Load bag mass data

BagMass <- 
  read_excel("LemonadeFireMasterData.xlsx",
             sheet = 'BagTares') %>%
  group_by(BagSize) %>%
  summarize(
    bagmass = mean(Mass))

# Load clipping data

ClipWeight <-
  read_excel("LemonadeFireMasterData.xlsx",
             sheet = 'ClipWeight')%>%
  mutate('ID' = paste(Location, Type, sep= " ")) %>%
  select( - Location,
          - Type) %>%
  data.table()

# Subtract bag weight from mass and conversions

ClipWeight <- 
  ClipWeight[ID, on = .(ID = LemonadeFireID)] %>%
  select(- PlotNumberPBG) %>%
  separate(AFRIname, c('SiteType', 'Replicate')) %>%
  separate(ID, c('Location', 'Type')) %>%
  select( - Replicate) %>%
  select(Location, Type, Quadrant, Nonant, MassGrams, BagSize, SiteType, BurnStatus, Month) %>%
  mutate(mass = case_when(
    BagSize ==
      "b"~MassGrams - filter(BagMass, BagSize == "b") $bagmass,
    BagSize ==
      "l" ~ MassGrams - filter(BagMass, BagSize == "l") $bagmass,
     BagSize ==
      "p"~MassGrams - filter(BagMass, BagSize == "p") $bagmass,
  )) %>%
  select( - MassGrams, - BagSize) %>%
  mutate(g1m2 = (mass / 0.05),
         KgHa = (g1m2 * 10),
         tac = (KgHa * 0.000404686))

# Plot difference in aboveground biomass between paired plots

ClipWeight %>%
  group_by(Location, SiteType, BurnStatus, Type) %>%
  summarise(
    meantac = mean(tac)) %>%
  pivot_wider(values_from = meantac,
              names_from = Type) %>%
  mutate(dif = ((E - O) / E) * 100) %>%
  ggplot(aes(x = dif,
             y = Location,
             color = SiteType)) + 
  labs(y = "", 
       x = "((E - O) / E) * 100%)",
       fill = "Ecological site",
       caption = "Fig. 1. Difference in aboveground biomass between paired open areas and exclosures.") +
  geom_vline(xintercept = 0, size = .3) +
  geom_label(aes(label = Location,
                 fill = factor(SiteType)),
             colour = "white",
             fontface = "bold") +
  theme(legend.position = "bottom",
        legend.direction = 'horizontal',
        legend.box = "horizontal",
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.background = element_rect(fill = "white",
                                        colour = "black",
                                        size = 0.5, linetype = "solid"),
        panel.grid.major = element_line(size = 0.1,
                                        linetype = 'solid',
                                        color = 'black'),
        panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line = element_line(color="black", size = .5),
        plot.margin = unit(c(0.1, 0.5, 0.1, 0.2),
                                "inches"),
        plot.caption = element_text(hjust = .45, size = 12)) +
  guides(fill = guide_legend(byrow = TRUE, override.aes = aes(label = "")))+
  scale_x_continuous(limits=c(-91, 91))

# Graph of burned vs unburned facet_wrap eco site

ClipWeight %>%
  group_by(Type, SiteType, BurnStatus, Month) %>%
  filter(Type == 'O') %>%
  mutate(BurnStatus =
           recode(BurnStatus,
                   "burned"="Burned",
                   "unburned"="Unburned")) %>%
  mutate(SiteType =
           recode(SiteType,
                  "sandy"="Sandy",
                  "silty"="Silty",
                  "steep"="Steep",
                  "shallow"="Shallow")) %>%
  summarise(
     meantac = mean(tac), 
     semtac = sd(tac) / sqrt(length(tac))) %>%
 ggplot(aes(x = Month, fill = BurnStatus)) + theme_bw(14) + 
      geom_line(aes(y = meantac,
                    color = BurnStatus,
                    group = BurnStatus),
                position = position_dodge(width = 0.3)) +
      geom_errorbar(aes(ymin = meantac - semtac,
                        ymax = meantac + semtac,
                        group = BurnStatus,
                        color = BurnStatus),
                    width = 0.15,
                    position = position_dodge(width = 0.3)) +
      geom_point(aes(y = meantac,
                 bg = BurnStatus,
                 shape = BurnStatus),
                 colour="black",
                 size=2.5, stroke = 1.1,
                 position = position_dodge(width = 0.3)) +
      facet_wrap(~SiteType) +
      labs(x = "Time",
           y = "Aboveground biomass (t/ac)",
           caption = "Fig. 2. Aboveground biomass at burned and unburned plots open to grazing.") +
  scale_color_manual('Burn status', values = wes_palette("Zissou1")[c(3, 5)]) +
  scale_fill_manual('Burn status', values = wes_palette("Zissou1")[c(3, 5)]) +
  scale_shape_manual('Burn status', values = c(21, 24)) + 
  scale_alpha_manual('Burn status', values = c(0.6, 0.6)) +
      theme(legend.position = "bottom",
            legend.direction = 'horizontal',
            legend.box = "horizontal",
            panel.grid.minor.x = element_blank(),
            legend.spacing.y = unit(0, "mm"),
            panel.border = element_rect(colour = "black", fill=NA),
            legend.background = element_blank(),
            plot.caption = element_text(hjust = .2, size = 12)) +
  scale_x_discrete(limits = c("May", "June", "July"))

# Facet grid of burn status and graze status

ClipWeight %>%
  select(Location, Type, SiteType, BurnStatus, Month, Quadrant, Nonant, tac)%>%
  filter(Type == 'E') %>%
  group_by(Location, Quadrant, Nonant) %>%
  filter(n()== 1) %>%
  group_by(Type, SiteType, BurnStatus, Month) %>%
  mutate(BurnStatus =
           recode(BurnStatus,
                   "burned"="Burned",
                   "unburned"="Unburned")) %>%
  mutate(SiteType =
           recode(SiteType,
                  "sandy"="Sandy",
                  "silty"="Silty",
                  "steep"="Steep",
                  "shallow"="Shallow")) %>%
  summarise(
     meantac = mean(tac), 
     semtac = sd(tac) / sqrt(length(tac))) %>%
 ggplot(aes(x = Month, fill = BurnStatus)) + theme_bw(14) + 
      geom_line(aes(y = meantac,
                    color = BurnStatus,
                    group = BurnStatus),
                position = position_dodge(width = 0.3)) +
      geom_errorbar(aes(ymin = meantac - semtac,
                        ymax = meantac + semtac,
                        group = BurnStatus,
                        color = BurnStatus),
                    width = 0.15,
                    position = position_dodge(width = 0.3)) +
      geom_point(aes(y = meantac,
                 bg = BurnStatus,
                 shape = BurnStatus),
                 colour="black",
                 size=2.5, stroke = 1.1,
                 position = position_dodge(width = 0.3)) +
      facet_wrap(~SiteType) +
      labs(x = "Time",
           y = "Aboveground biomass (t/ac)",
           caption = "Fig. 3. biomass growth in burned and unburned exclosures.") +
  scale_color_manual('Burn status', values = wes_palette("Zissou1")[c(3, 5)]) +
  scale_fill_manual('Burn status', values = wes_palette("Zissou1")[c(3, 5)]) +
  scale_shape_manual('Burn status', values = c(21, 24)) + 
  scale_alpha_manual('Burn status', values = c(0.6, 0.6)) +
      theme(legend.position = "bottom",
            legend.direction = 'horizontal',
            legend.box = "horizontal",
            panel.grid.minor.x = element_blank(),
            legend.spacing.y = unit(0, "mm"),
            panel.border = element_rect(colour = "black", fill=NA),
            legend.background = element_blank(),
            plot.caption = element_text(hjust = .45, size = 12)) +
  scale_x_discrete(limits = c("May", "June", "July"))

Statistics

Descriptive Statistics

# Descriptive stats

stat.desc(ClipWeight$tac)

##      nbr.val     nbr.null       nbr.na          min          max        range 
## 5.940000e+02 0.000000e+00 0.000000e+00 9.064966e-03 3.231175e+00 3.222110e+00 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 5.192136e+02 7.101430e-01 8.740970e-01 2.530921e-02 4.970660e-02 3.804905e-01 
##      std.dev     coef.var 
## 6.168391e-01 7.056873e-01

Histogram with Normal Curve

# Histogram showing data distribution and normality curve

ClipWeight %>%
  ggplot(aes(x = tac)) + theme_bw(16) + 
  geom_histogram(aes(y = after_stat(.data[["density"]])),
                 binwidth = .1,
                 colour = "black",
                 fill = "lightyellow")+
  geom_density(alpha = .5, fill = "lightblue") +
  stat_function(fun = dnorm,
                args = list(mean = 0.72625780,
                            sd = 0.50911132),
                colour = "red",
                linewidth = 1.1) +
  labs(y = "Density", 
       x = "Aboveground biomass (t/ac)")

Ratio between Mean and Median

# Ratio between mean and median

knitr::kable(tibble(Mean = mean(ClipWeight$tac),
                    Median = median(ClipWeight$tac),
                    Ratio = Mean/Median))

Mean	Median	Ratio
0.874097	0.710143	1.230875

Skewness

# Skewness

skewness(ClipWeight$tac)

## [1] 1.220551

Kurtosis

# Kurtosis

kurtosis(ClipWeight$tac)

## [1] 4.416548

Kolmogorov-Smirnov Test

# Kolmogorov-Smirnov test

KSClip <-
  ks.test(ClipWeight$tac, "pnorm")

pander(KSClip)

Asymptotic one-sample Kolmogorov-Smirnov test: ClipWeight$tac

Test statistic	P value	Alternative hypothesis
0.5381	0 * * *	two-sided

QQ-Plot

# QQ-plot

car::qqPlot(ClipWeight$tac, xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")

Logarithmic Transformation

# Logarithmic transformation

SqRtClip <-
ClipWeight %>%
  mutate(SqRttac =  sqrt(tac))

Descriptive Statistics

# Descriptive stats

stat.desc(SqRtClip$SqRttac)

##      nbr.val     nbr.null       nbr.na          min          max        range 
## 594.00000000   0.00000000   0.00000000   0.09521012   1.79754691   1.70233679 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 522.21070183   0.84269980   0.87914260   0.01306394   0.02565722   0.10137592 
##      std.dev     coef.var 
##   0.31839586   0.36216635

Histogram with Normal Curve

# Histogram showing data distribution and normality curve

SqRtClip %>%
  ggplot(aes(x = SqRttac)) + theme_bw(16) + 
  geom_histogram(aes(y = after_stat(.data[["density"]])),      
                 binwidth = .1,
                 colour = "black",
                 fill = "lightyellow")+
  geom_density(alpha = .5, fill = "lightblue") +
  stat_function(fun = dnorm,
                args = list(mean = 0.87914260,
                          sd = 0.31839586),
                colour = "red",
                linewidth = 1.1) +
  labs(y = "Density", 
       x = "Log aboveground biomass (t/ac)")

Ratio between Mean and Median

# Ratio between mean and median

knitr::kable(tibble(MeanSqRt = mean(SqRtClip$SqRttac),
                    MedianSqRt = median(SqRtClip$SqRttac),
                    RatioSqRt = MeanSqRt/MedianSqRt))

MeanSqRt	MedianSqRt	RatioSqRt
0.8791426	0.8426998	1.043245

Skewness

# Skewness

skewness(SqRtClip$SqRttac)

## [1] 0.4340627

Kurtosis

# Kurtosis

kurtosis(SqRtClip$SqRttac)

## [1] 2.740046

Kolmogorov-Smirnov Test

# Kolmogorov-Smirnov test

KSLogClip <-
  ks.test(SqRtClip$SqRttac, "pnorm")

pander(KSLogClip)

Asymptotic one-sample Kolmogorov-Smirnov test: SqRtClip$SqRttac

Test statistic	P value	Alternative hypothesis
0.6289	0 * * *	two-sided

QQ-Plot

# QQ-plot

car::qqPlot(SqRtClip$SqRttac, xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")

Soil Moisture

Graphs

# Diagram showing methods

DiagrammeR::grViz("digraph {
  graph [center = TRUE, layout = dot, rankdir = TB]
  node [shape = rectangle]        
  rec1 [label = '1. Obtain soil moisture data']
  rec2 [label = '2. Enter data in Excel']
  rec3 [label =  '3. Analyze data in R Markdown']
  rec1 -> rec2 -> rec3
  }",
  height = 225)

# Load soil moisture data

SoilMoist <-
  read_excel("LemonadeFireMasterData.xlsx", sheet = "SoilMoist") %>%
  mutate('ID' = paste(Location, Type, sep= " ")) %>%
  select( - Location,
          - Type)%>%
  data.table()

# Cross reference moisture data with ID data

SoilMoist <- 
  SoilMoist[ID, on = .(ID = LemonadeFireID)] %>%
  select(- PlotNumberPBG) %>%
  separate(AFRIname, c('SiteType', 'Replicate')) %>%
  select( - Replicate)

# Plot difference in % soil moisture between paired plots

SoilMoist %>%
  group_by(Location, SiteType, BurnStatus, GrazeStatus) %>%
  summarise(
    meanmois = mean(Moist)) %>%
  pivot_wider(values_from = meanmois,
              names_from = GrazeStatus) %>%
  mutate(dif = ((grazed - ungrazed) / ungrazed) * 100) %>%
  ggplot(aes(x = dif,
             y = Location,
             color = SiteType)) + 
  labs(y = "", 
       x = "((E - O) / E) * 100%)",
       fill = "Ecological Site",
       caption = "Fig. 4. Difference in soil moisture between paired open areas and exclosures.") +
  geom_vline(xintercept = 0, size = .3) +
  geom_label(aes(label = Location,
                fill = factor(SiteType)),
             colour = "white",
             fontface = "bold") +
  theme(legend.position = "bottom",
        legend.direction = 'horizontal',
        legend.box = "horizontal",
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.background = element_rect(fill = "white",
                                        colour = "black",
                                        size = 0.5, linetype = "solid"),
        panel.grid.major = element_line(size = 0.1,
                                        linetype = 'solid',
                                        color = 'black'),
        panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.line = element_line(color="black", size = .5),
        plot.margin = unit(c(0.1, 0.5, 0.1, 0.2),
                                "inches"),
        plot.caption = element_text(hjust = .45, size = 12)) +
  guides(fill = guide_legend(byrow = TRUE, override.aes = aes(label = "")))+
  scale_x_continuous(limits=c(-91, 91))

# Facet grid of burn status and graze status

SoilMoist %>%
  group_by(Type, SiteType, BurnStatus, Month) %>%
  filter(Type == 'O') %>%
  mutate(BurnStatus =
           recode(BurnStatus,
                   "burned"="Burned",
                   "unburned"="Unburned")) %>%
  mutate(SiteType =
           recode(SiteType,
                  "sandy"="Sandy",
                  "silty"="Silty",
                  "steep"="Steep",
                  "shallow"="Shallow")) %>%
  summarise(
    meanmoist = mean(Moist), 
    semmoist = sd(Moist) / sqrt(length(Moist))) %>%
 ggplot(aes(x = Month, fill = BurnStatus)) + theme_bw(14) + 
      geom_line(aes(y = meanmoist,
                    color = BurnStatus,
                    group = BurnStatus),
                position = position_dodge(width = 0.3)) +
      geom_errorbar(aes(ymin = meanmoist - semmoist,
                        ymax = meanmoist + semmoist,
                        group = BurnStatus,
                        color = BurnStatus),
                    width = 0.15,
                    position = position_dodge(width = 0.3)) +
      geom_point(aes(y = meanmoist,
                 bg = BurnStatus,
                 shape = BurnStatus),
                 colour="black",
                 size = 2.5, stroke = 1.1,
                 position = position_dodge(width = 0.3)) +
      facet_wrap(~SiteType) +
      labs(x = "Time",
           y = "% soil moisture",
           caption = "Fig. 5. Soil moisture in burned and unburned locations.") +
  scale_color_manual('Burn status', values = wes_palette("Zissou1")[c(3, 5)]) +
  scale_fill_manual('Burn status', values = wes_palette("Zissou1")[c(3, 5)]) +
  scale_shape_manual('Burn status', values = c(21, 24)) + 
  scale_alpha_manual('Burn status', values = c(0.6, 0.6)) +
      theme(legend.position = "bottom",
            legend.direction = 'horizontal',
            legend.box = "horizontal",
            panel.grid.minor.x = element_blank(),
            legend.spacing.y = unit(0, "mm"),
            panel.border = element_rect(colour = "black", fill=NA),
            legend.background = element_blank(),
            plot.caption = element_text(hjust = .45, size = 12)) +
  scale_x_discrete(limits = c("May", "June", "July"))

Statistics

Descriptive Statistics

# Descriptive stats

stat.desc(SoilMoist$Moist)

##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.000000e+02 6.100000e+01 0.000000e+00 0.000000e+00 6.000000e+01 6.000000e+01 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 1.141510e+04 1.040000e+01 1.268344e+01 3.420047e-01 6.712206e-01 1.052705e+02 
##      std.dev     coef.var 
## 1.026014e+01 8.089397e-01

Histogram with Normal Curve

# Histogram showing data distribution and normality curve

SoilMoist %>%
  ggplot(aes(x = Moist)) + theme_bw(16) +
  geom_histogram(aes(y = after_stat(.data[["density"]])),
                 binwidth = 3,
                 colour = "black",
                 fill = "lightyellow")+
  geom_density(alpha = .5, fill = "lightblue") +
  stat_function(fun = dnorm,
                args = list(mean = 17.5430556,
                          sd = 10.9276183),
                colour = "red",
                linewidth = 1.1) +
  labs(y = "Density", 
       x = "% soil moisture")

QQ-Plot

# QQ-plot

car::qqPlot(SoilMoist$Moist, xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")

Ratio between Mean and Median

knitr::kable(tibble(Mean = mean(SoilMoist$Moist),
                    Median = median(SoilMoist$Moist),
                    Ratio = Mean/Median))

Mean	Median	Ratio
12.68344	10.4	1.219562

Skewness

# Skewness

skewness(SoilMoist$Moist)

## [1] 1.463398

Kurtosis

# Kurtosis

kurtosis(SoilMoist$Moist)

## [1] 5.687223

Kolmogorov-Smirnov Test

# Kolmogorov-Smirnov test

KSMoist <-
  ks.test(SoilMoist$Moist, "pnorm")

pander(KSMoist)

Asymptotic one-sample Kolmogorov-Smirnov test: SoilMoist$Moist

Test statistic	P value	Alternative hypothesis
0.8893	0 * * *	two-sided

Square Root Transformation

# Square root transformation

SRMoist <-
SoilMoist %>%
  mutate(SqRtMoist =  sqrt(Moist))

Descriptive Statistics

# Descriptive stats

stat.desc(SRMoist$SqRtMoist)

##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.000000e+02 6.100000e+01 0.000000e+00 0.000000e+00 7.745967e+00 7.745967e+00 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 2.905259e+03 3.224903e+00 3.228065e+00 5.017256e-02 9.846899e-02 2.265558e+00 
##      std.dev     coef.var 
## 1.505177e+00 4.662784e-01

Histogram with Normal Curve

# Histogram showing data distribution and normality curve

SRMoist %>%
  ggplot(aes(x = SqRtMoist)) + 
  theme_bw(16) + 
  geom_histogram(aes(y = after_stat(.data[["density"]])),
                 binwidth = .4,
                 colour = "black",
                 fill = "lightyellow")+
  geom_density(alpha = .5, fill = "lightblue")+
  stat_function(fun = dnorm,
                args = list(mean = 3.55608073,
                          sd = 1.33445545),
                colour = "red",
                linewidth = 1.1) +
  labs(y = "Density", 
       x = "Log % soil moisture")

Ratio between Mean and Median

# Ratio between mean and median

knitr::kable(tibble(Mean = mean(SRMoist$SqRtMoist),
                    Median = median(SRMoist$SqRtMoist),
                    Ratio = Mean/Median))

Mean	Median	Ratio
3.228065	3.224903	1.000981

Skewness

# Skewness

skewness(SRMoist$SqRtMoist)

## [1] -0.01944163

Kurtosis

# Kurtosis

kurtosis(SRMoist$SqRtMoist)

## [1] 3.28084

Kolmogorov-Smirnov Test

# Kolmogorov-Smirnov test

KSSqRtMoist <-
  ks.test(SRMoist$SqRtMoist, "pnorm")

pander(KSSqRtMoist)

Asymptotic one-sample Kolmogorov-Smirnov test: SRMoist$SqRtMoist

Test statistic	P value	Alternative hypothesis
0.8353	0 * * *	two-sided

QQ-Plot

# QQ-plot

car::qqPlot(SRMoist$SqRtMoist, xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")