BIOL9915 Project - Fire and Grazing in Grassland Landscapes Influences Bird Abundance
Evan Lewis
Grassland are ecosystems that are maintained via disturbance (xx). Fire and grazing were historical forms of large scale disturbance in the landscape. Without fire and grazing, ecological succession is allowed to occur, and grasslands shift from open, herbaceous landscapes to landscapes dominated by shrubs and trees (xx).
#knitr::include_graphics("image.png")
#this code chunk will be used to show comparison photos of grasslands in their disturbance-maintained state, vs a grassland undergoing encroachment from woody plantsGrassland birds rely on the disturbance from fire and grazing to maintain the landscape heterogeneity they require for high-quality territories (xx). However, the effects of fire on a landscape are not uniform, instead, they depend on how frequently an area is burned. To test the effect of fire frequency and grazing on the abundance of grassland birds, I will be using data from Konza Prairie Reserve (Konza, hereafter). Over (number of years) years, transect surveys were completed over a series of factorially treated plots. Combinations of 1, 2, 3, 4, and 20 year fire return intervals, and grazing treatments were applied over the experimental pastures. This report shows the steps taken in the analysis, along with an interpretation of the results.
#first we load packages and data, and take a quick look at the data frame.
library(tidyverse)
library(ggplot2)
library(ggfortify)
data <- read.csv('data/eame_burn_graze.csv')
glimpse(data)## Rows: 3,553
## Columns: 30
## $ TransectName <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "2",…
## $ CurrWatershed <chr> "N4D", "N4D", "N4D", "N4D", "N4D", "N4D", "N4D", …
## $ Grazed. <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, T…
## $ GrazingType <chr> "bison", "bison", "bison", "bison", "bison", "bis…
## $ FireReturnInterval <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 20, 20, 20, 20, 20,…
## $ TransectLength <int> 961, 961, 961, 961, 961, 961, 961, 961, 961, 778,…
## $ StartLatitude <dbl> 39.07687, 39.07687, 39.07687, 39.07687, 39.07687,…
## $ StartLongitude <dbl> -96.58044, -96.58044, -96.58044, -96.58044, -96.5…
## $ EndLatitude <dbl> 39.08617, 39.08617, 39.08617, 39.08617, 39.08617,…
## $ EndLongitude <dbl> -96.58193, -96.58193, -96.58193, -96.58193, -96.5…
## $ MidpointLatitude <dbl> 39.08148, 39.08148, 39.08148, 39.08148, 39.08148,…
## $ MidpointLongitude <dbl> -96.58116, -96.58116, -96.58116, -96.58116, -96.5…
## $ SurveyEventID <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 3, 3, 3…
## $ Duration <int> 50, 50, 50, 50, 50, 50, 50, 50, 50, 38, 70, 70, 7…
## $ Observer <chr> "EJF", "EJF", "EJF", "EJF", "EJF", "EJF", "EJF", …
## $ Date <chr> "12-Jun-86", "12-Jun-86", "12-Jun-86", "12-Jun-86…
## $ StartTime <chr> "6:13", "6:13", "6:13", "6:13", "6:13", "6:13", "…
## $ StartCloud <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ StartWind <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ StartTemp <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ BirdRecordID <int> 21619, 21643, 21649, 21654, 21674, 21687, 21698, …
## $ SpeciesCode <chr> "EAME", "EAME", "EAME", "EAME", "EAME", "EAME", "…
## $ DistanceFromLine_Old <int> 25, 40, 40, 50, 60, 75, 20, 15, 15, 65, 70, 75, 2…
## $ Count <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ ObsLocation <int> -999, -999, -999, -999, -999, -999, -999, -999, -…
## $ Angle <int> -999, -999, -999, -999, -999, -999, -999, -999, -…
## $ AngularDistance <int> -999, -999, -999, -999, -999, -999, -999, -999, -…
## $ Flyover. <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …
## $ OffWatershed. <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …
## $ BeyondLine. <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …These data include many variables that are not needed in the analysis. The next step is to do some data tidying to make analysis easier.
#each row of the dataframe represents a detection of an individual bird. Since we are interested in total abundance, we need to group these detections to get total numbers.
abundance_per_survey <- data %>%
group_by(TransectName, Date) %>%
summarise(abundance = sum(Count))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by TransectName and Date.
## ℹ Output is grouped by TransectName.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(TransectName, Date))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.mean_per_transect <- abundance_per_survey %>%
group_by(TransectName) %>%
summarise(mean_abundance = mean(abundance))
#now we need to join these to the original data so we have the context.
clean <- mean_per_transect %>%
left_join(distinct(data, TransectName, .keep_all = TRUE), by = "TransectName")
#Now we have a mean abundance of birds detected across all surveys done on each transect. A bit more cleaning is required still. Lets get rid of all the unneeded variables, and make the mean per km, to account for different lengths of transects.
clean <- clean %>%
select(TransectName | mean_abundance | Grazed. | FireReturnInterval | TransectLength | Date | Count) %>%
mutate(mean_ab_per_km = mean_abundance/(TransectLength/1000))That tidying has made these data workable for this analysis. An inital plot will give us an idea of how to proceed.
ggplot(clean, aes(x = FireReturnInterval, y = mean_ab_per_km, fill = Grazed.)) +
geom_boxplot() +
geom_jitter(position = position_jitterdodge(jitter.width = 0.15))## Warning: Orientation is not uniquely specified when both the x and y aesthetics are
## continuous. Picking default orientation 'x'.
Well, that’s not very helpful. Looks like FireReturnInterval is recorded
as an integer, rather than a factor. There are also some points where
mean abundance per kilometer are recorded as ~390, which can’t be
true.
#First change FireReturnInterval to a factor, and filter out the outlier mean abundance counts.
clean <- clean %>%
mutate(FireReturnInterval = as.factor(FireReturnInterval)) %>%
filter(mean_ab_per_km < 100)
ggplot(clean, aes(x = FireReturnInterval, y = mean_ab_per_km, fill=Grazed.)) +
geom_boxplot()
Much better! Visually, this looks promising. It looks like there are
some pretty obvious differences in bird abundance as a function of the
frequency of fire on the landscape as well as how that may interact with
grazing. Let’s put this into a linear model to see if the visual
differences are statistically significant.
# Three models, one testing just fire return interval, one testing just the grazing, and one testing the interaction between the two.
mod1 <- lm(mean_ab_per_km ~ FireReturnInterval, data = clean)
mod2 <- lm(mean_ab_per_km ~ Grazed., data = clean)
mod3 <- lm(mean_ab_per_km ~ FireReturnInterval:Grazed., data = clean)
#checking the summary and anova table for all three models. Manually checking for best fit using the residual sum of squares.
summary(mod1)##
## Call:
## lm(formula = mean_ab_per_km ~ FireReturnInterval, data = clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.1567 -1.8449 -0.1568 1.5041 8.5888
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.638 1.094 3.326 0.00124 **
## FireReturnInterval1 5.419 1.211 4.475 2.06e-05 ***
## FireReturnInterval2 8.852 1.498 5.910 4.97e-08 ***
## FireReturnInterval3 6.507 1.211 5.374 5.22e-07 ***
## FireReturnInterval4 4.167 1.311 3.177 0.00199 **
## FireReturnInterval20 2.921 1.399 2.088 0.03939 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.894 on 98 degrees of freedom
## Multiple R-squared: 0.341, Adjusted R-squared: 0.3074
## F-statistic: 10.14 on 5 and 98 DF, p-value: 7.552e-08summary(mod2)##
## Call:
## lm(formula = mean_ab_per_km ~ Grazed., data = clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5075 -2.5150 -0.4158 2.1735 9.4612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.3166 0.5288 15.727 <2e-16 ***
## Grazed.TRUE 0.8651 0.6905 1.253 0.213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.468 on 102 degrees of freedom
## Multiple R-squared: 0.01516, Adjusted R-squared: 0.005502
## F-statistic: 1.57 on 1 and 102 DF, p-value: 0.2131summary(mod3)##
## Call:
## lm(formula = mean_ab_per_km ~ FireReturnInterval:Grazed., data = clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4712 -1.6887 -0.0987 1.3619 8.4235
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1841 2.8045 2.205 0.02989 *
## FireReturnInterval-999:Grazed.FALSE -2.5461 2.9982 -0.849 0.39792
## FireReturnInterval1:Grazed.FALSE 4.1344 2.8965 1.427 0.15678
## FireReturnInterval2:Grazed.FALSE 8.3794 3.1356 2.672 0.00888 **
## FireReturnInterval3:Grazed.FALSE NA NA NA NA
## FireReturnInterval4:Grazed.FALSE 1.4082 2.9982 0.470 0.63968
## FireReturnInterval20:Grazed.FALSE 0.4129 2.9414 0.140 0.88867
## FireReturnInterval-999:Grazed.TRUE NA NA NA NA
## FireReturnInterval1:Grazed.TRUE 1.6906 2.8908 0.585 0.56008
## FireReturnInterval2:Grazed.TRUE 4.2326 3.1356 1.350 0.18030
## FireReturnInterval3:Grazed.TRUE 3.9613 2.8494 1.390 0.16774
## FireReturnInterval4:Grazed.TRUE 1.7858 2.9562 0.604 0.54724
## FireReturnInterval20:Grazed.TRUE NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.805 on 94 degrees of freedom
## Multiple R-squared: 0.4063, Adjusted R-squared: 0.3495
## F-statistic: 7.148 on 9 and 94 DF, p-value: 7.821e-08anova(mod1, mod2, mod3)## Analysis of Variance Table
##
## Model 1: mean_ab_per_km ~ FireReturnInterval
## Model 2: mean_ab_per_km ~ Grazed.
## Model 3: mean_ab_per_km ~ FireReturnInterval:Grazed.
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 98 820.69
## 2 102 1226.45 -4 -405.76 12.8969 2.075e-08 ***
## 3 94 739.35 8 487.10 7.7412 6.313e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#model 3 shows the smallest residual sum of squares, so that model is the best fit.
#check the assumptions
autoplot(mod3, smooth.colour = NA) 
It looks like these data fit the assumptions of the model fairly well,
so a linear model approach is reasonable. The anova table above shows
that model three has the smallest residual sum of squares, or the least
error, and is the best fit for the data.
Model three includes the interaction between grazing interval, so let’s plot the data in a way that shows that interaction.
vis <- clean %>%
group_by(FireReturnInterval, Grazed.) %>%
summarise(
total_mean = mean(mean_ab_per_km),
seAB = sd(mean_ab_per_km)/sqrt(n())
)## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by FireReturnInterval and Grazed..
## ℹ Output is grouped by FireReturnInterval.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(FireReturnInterval, Grazed.))` for per-operation
## grouping (`?dplyr::dplyr_by`) instead.ggplot(vis, aes(x = FireReturnInterval, y = total_mean, colour = Grazed., group = Grazed.)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = total_mean - seAB,ymax = total_mean + seAB),
width = 0.1)
This plot shows that the slopes between the grazing treatments are not
the same, which tells us that there is indeed an interaction between the
variables. Interestingly, it seems that there is a more significant
interaction between grazing and fire return interval at the longer end
of the return interval. This suggests that when fire is absent on the
landscape, grazing can fill a similar role of maintaining the conditions
needed for grassland bird abundance.
The results of this analysis shows that the frequency of fire on the landscape influences the abundance of grassland birds. Based on the model, a two-year return interval had the largest effect, showing that plots burned every two years had a mean abundance of birds that was 7.288 ± 3.45 higher than unburned plots (F9,48 = 7.148, p < 0.001).