October 21, 2026 Olivia Wiley

Quantitative methods used in ecology: Analyzing effects of vegetation structure on nesting choices of grassland birds.

Objective: 
The aim of this project is to examine vegetation density as a factor in reproductive success of grassland birds. This analysis seeks to determine whether vegetation density varies with proximity to the nest center. This will be explored using nesting data from three grassland species, Dickcissels (Spiza americana), Eastern Meadowlarks (Sturnella magna), and Grasshopper Sparrows (Ammodramus savannarum) at managed tallgrass prairie site in the Flint Hills of Kansas. I will look at multiple grassland species to see if there is variation in vegetation density at nesting center, and determine if there is a difference among species, which may illuminate different nesting strategies among the species. This work is also important to understanding drivers of decline in grassland birds and informing strategies for managing critical habitat (Kaplan, 2026). 
Hypothesis: Vegetation density will be higher around nesting center compared to edge. Denser nesting locations may provide more protection from predators, parasitism, or effects of adverse weather (Colombo et al., 2024) This concealment, which provides protection from external threats may also contribute to increased parental investment in rearing of young (Sroor et al., 2025)

Illustrate field site

To begin this analysis, first we clear R’s brain and set up

rm(list=ls())


library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(readr)
library(tidyr)
library(ggfortify)

#Input CSV
data <- read_csv("C:/Users/olivi/Desktop/BIO9915_BIOSTATSwork/data/Qy_NestVegDataForOlivia.csv")
## Rows: 975 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (7): PlotName, Point, Date, Notes, Species, InferredFateField, Inferred...
## dbl (15): VegID, CenterRobel-N, CenterRobel-E, CenterRobel-S, CenterRobel-W,...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(data)
## Rows: 975
## Columns: 22
## $ VegID             <dbl> 677, 678, 679, 680, 681, 686, 693, 694, 696, 697, 69…
## $ PlotName          <chr> "C3SC", "N1A", "N2B", "C3C", "C3C", "C3B", "C3SB", "…
## $ Point             <chr> "C3SC-SKW-07", "N1A-CRW-03", "N2B-SKW-08", "C3C-SKW-…
## $ Date              <chr> "11-Jun-16", "10-Jun-16", "07-Jun-16", "07-Jun-16", …
## $ `CenterRobel-N`   <dbl> 0.5, 0.5, 6.0, 1.0, 2.0, 1.5, 2.0, 2.5, 2.5, 1.0, 1.…
## $ `CenterRobel-E`   <dbl> 1.0, 4.0, 3.0, 2.0, 0.0, 1.5, 1.0, 2.0, 2.5, 1.5, 1.…
## $ `CenterRobel-S`   <dbl> 2.0, 2.5, 6.0, 2.0, 1.5, 1.5, 2.0, 1.5, 2.0, 1.0, 1.…
## $ `CenterRobel-W`   <dbl> 2.0, 1.5, 10.0, 2.0, 0.5, 1.5, 1.0, 1.0, 2.5, 1.5, 1…
## $ `EdgeRobel-N`     <dbl> 1.5, 1.0, 3.0, 0.5, 2.5, 0.0, 2.0, 0.0, 0.0, 1.0, 0.…
## $ `EdgeRobel-E`     <dbl> 0.5, 2.0, 4.0, 2.0, 3.0, 0.5, 0.0, 2.0, 2.5, 1.5, 7.…
## $ `EdgeRobel-S`     <dbl> 2.5, 0.0, 4.0, 1.0, 3.0, 2.0, 2.0, 0.0, 1.0, 4.5, 1.…
## $ `EdgeRobel-W`     <dbl> 1.5, 0.0, 10.0, 1.5, 1.5, 2.0, 0.5, 0.5, 0.5, 2.0, 4…
## $ `Plot%Shrub`      <dbl> 6, 0, 40, 0, 0, 6, 5, 5, 0, 15, 10, 2, 0, 0, 0, 0, 5…
## $ `Plot%Tree`       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ Aspect            <dbl> 273, -9, 124, 44, 189, -9, 60, 28, 150, 68, 310, 178…
## $ Slope             <dbl> -9, -9, 20, 10, 2, -9, 5, 4, 5, 2, 8, 0, 0, 5, 2, 2,…
## $ Notes             <chr> "Missing % nest and Slope/Aspect wrong-- recorded as…
## $ Species           <chr> "GRSP", "GRSP", "GRSP", "GRSP", "GRSP", "GRSP", "GRS…
## $ MaxNHostEggs      <dbl> 4, 5, 3, 1, 2, 3, 2, 0, 2, 3, 2, 1, 3, 0, 5, 5, 5, 5…
## $ TotalClutchSize   <dbl> 4, 5, 3, 4, 3, 6, 5, 4, 4, 5, 3, 1, 5, 1, 5, 5, 5, 5…
## $ InferredFateField <chr> "SU", "PR", "SU", "PR", "PR", "PR", "PR", "AB", "PR"…
## $ InferredFateLab   <chr> "SU", "PR", "SU", "PR", NA, "PR", "AB", "AB", "AB", …

Tidy dataset: Filter dataset to select variables needed for analysis

data2 <- data |> 
   filter(if_all(where(is.numeric), ~ .x >= 0))|>   #filter any negative data points (missing data recorded as -9 in dataset)
  dplyr::select(`VegID`, `PlotName`, `CenterRobel-N` , `CenterRobel-E` , `CenterRobel-W` , `CenterRobel-S` , `EdgeRobel-N` , `EdgeRobel-S` , `EdgeRobel-E` , `EdgeRobel-W`, `Species`, `MaxNHostEggs`)       #Choose columns with vegetation density measurements, as well as ID information

glimpse(data2)          # Look at organized data
## Rows: 877
## Columns: 12
## $ VegID           <dbl> 679, 680, 681, 693, 694, 696, 697, 698, 699, 700, 701,…
## $ PlotName        <chr> "N2B", "C3C", "C3C", "C3SB", "N2B", "C3C", "C3B", "C3B…
## $ `CenterRobel-N` <dbl> 6.0, 1.0, 2.0, 2.0, 2.5, 2.5, 1.0, 1.0, 0.5, 1.5, 0.5,…
## $ `CenterRobel-E` <dbl> 3.0, 2.0, 0.0, 1.0, 2.0, 2.5, 1.5, 1.0, 1.0, 2.0, 2.0,…
## $ `CenterRobel-W` <dbl> 10.0, 2.0, 0.5, 1.0, 1.0, 2.5, 1.5, 1.0, 1.5, 2.0, 0.5…
## $ `CenterRobel-S` <dbl> 6.0, 2.0, 1.5, 2.0, 1.5, 2.0, 1.0, 1.0, 1.0, 2.5, 2.0,…
## $ `EdgeRobel-N`   <dbl> 3.0, 0.5, 2.5, 2.0, 0.0, 0.0, 1.0, 0.5, 2.0, 1.5, 0.5,…
## $ `EdgeRobel-S`   <dbl> 4.0, 1.0, 3.0, 2.0, 0.0, 1.0, 4.5, 1.0, 0.5, 1.0, 0.0,…
## $ `EdgeRobel-E`   <dbl> 4.0, 2.0, 3.0, 0.0, 2.0, 2.5, 1.5, 7.5, 0.0, 0.0, 0.0,…
## $ `EdgeRobel-W`   <dbl> 10.0, 1.5, 1.5, 0.5, 0.5, 0.5, 2.0, 4.0, 1.0, 1.0, 0.0…
## $ Species         <chr> "GRSP", "GRSP", "GRSP", "GRSP", "GRSP", "GRSP", "GRSP"…
## $ MaxNHostEggs    <dbl> 3, 1, 2, 2, 0, 2, 3, 2, 1, 3, 0, 5, 5, 5, 5, 5, 5, 3, …

Calculate the averages of center and edge robel measurements and create new column for each mean.

data_averaged <- data2 |>
  mutate (mean_center = ((`CenterRobel-N` + `CenterRobel-E` + `CenterRobel-W`+ `CenterRobel-S`)/4),  mean_edge = (( `EdgeRobel-N` + `EdgeRobel-S` + `EdgeRobel-E` + `EdgeRobel-W`)/4))

glimpse(data_averaged)
## Rows: 877
## Columns: 14
## $ VegID           <dbl> 679, 680, 681, 693, 694, 696, 697, 698, 699, 700, 701,…
## $ PlotName        <chr> "N2B", "C3C", "C3C", "C3SB", "N2B", "C3C", "C3B", "C3B…
## $ `CenterRobel-N` <dbl> 6.0, 1.0, 2.0, 2.0, 2.5, 2.5, 1.0, 1.0, 0.5, 1.5, 0.5,…
## $ `CenterRobel-E` <dbl> 3.0, 2.0, 0.0, 1.0, 2.0, 2.5, 1.5, 1.0, 1.0, 2.0, 2.0,…
## $ `CenterRobel-W` <dbl> 10.0, 2.0, 0.5, 1.0, 1.0, 2.5, 1.5, 1.0, 1.5, 2.0, 0.5…
## $ `CenterRobel-S` <dbl> 6.0, 2.0, 1.5, 2.0, 1.5, 2.0, 1.0, 1.0, 1.0, 2.5, 2.0,…
## $ `EdgeRobel-N`   <dbl> 3.0, 0.5, 2.5, 2.0, 0.0, 0.0, 1.0, 0.5, 2.0, 1.5, 0.5,…
## $ `EdgeRobel-S`   <dbl> 4.0, 1.0, 3.0, 2.0, 0.0, 1.0, 4.5, 1.0, 0.5, 1.0, 0.0,…
## $ `EdgeRobel-E`   <dbl> 4.0, 2.0, 3.0, 0.0, 2.0, 2.5, 1.5, 7.5, 0.0, 0.0, 0.0,…
## $ `EdgeRobel-W`   <dbl> 10.0, 1.5, 1.5, 0.5, 0.5, 0.5, 2.0, 4.0, 1.0, 1.0, 0.0…
## $ Species         <chr> "GRSP", "GRSP", "GRSP", "GRSP", "GRSP", "GRSP", "GRSP"…
## $ MaxNHostEggs    <dbl> 3, 1, 2, 2, 0, 2, 3, 2, 1, 3, 0, 5, 5, 5, 5, 5, 5, 3, …
## $ mean_center     <dbl> 6.250, 1.750, 1.000, 1.500, 1.750, 2.375, 1.250, 1.000…
## $ mean_edge       <dbl> 5.250, 1.250, 2.500, 1.125, 0.625, 1.000, 2.250, 3.250…

Transform data into wide form to include averages as separate categorical (center vs. edge) values

tidy_data <- data_averaged |>
 
   pivot_longer(
    cols = c(mean_center, mean_edge), 
    names_to = "Location", 
    values_to = "Mean"
  ) %>%
  mutate(
    Location = case_match(
      Location,
      "mean_center" ~ "Center",
      "mean_edge"   ~ "Edge"
    )
  )
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Location = case_match(...)`.
## Caused by warning:
## ! `case_match()` was deprecated in dplyr 1.2.0.
## ℹ Please use `recode_values()` instead.

Plot the data to visualize how vegetation density compared between center and edge, factoring in different species.

ggplot(tidy_data, aes(x = Location, y = Mean, colour = Species )) + 
  geom_boxplot( size = 0.6) +
  xlab ("Location") +   #add x axis label
  ylab("Vegetation Density") +  #add y axis label
  theme_bw()      #Remove grey background

Box plots reveal initial predictions about outcomes of analysis: Most vegetation densities fall between 0 and 1 with outliers around 13-15. Vegetation density of nesting habitat looks like it may change among species. Vegetation density may be slightly higher around nesting center compared to edge.

Analyzing the data -

Conduct a linear model and check assumptions using autoplot function.This model measures the following interactions:

  1. Does the mean change dependent on just center vs. edge location?
  2. Does the mean change dependent on just the species?
  3. Does the effect of location depend on which species is being looked at? Does the effect of species depend on the location? (interaction of variables).
mod1 <- lm (Mean ~ Location * Species, data = tidy_data)

autoplot (mod1)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Autoplot results: 

First output of assumptions showed strongly right skewed data. Log transform the data to help it fit the assumptions of normal distribution. Log transformation stabilizes the variance, so the model is just testing the proportional change among means. Otherwise, outliers could skew and influence the outputted significance numbers.(should this be done? Or is there a model better fit this data?).

Conduct another linear model with log transformed data. Check assumptions with autoplot.

tidy_data <- tidy_data |> 
  mutate(log_Mean = log(Mean + 1))

mod2 <- lm (log_Mean ~ Location * Species, data = tidy_data)

autoplot (mod2)

Assumptions fit normal distribution much better after log transformation. Right skew is minimized.

Run an anova and summary to determine if differences in vegetation in center vs edge of habitat are significant.

anova(mod2)
## Analysis of Variance Table
## 
## Response: log_Mean
##                    Df  Sum Sq Mean Sq  F value    Pr(>F)    
## Location            1   2.540   2.540  19.7896  9.19e-06 ***
## Species             2  76.520  38.260 298.1192 < 2.2e-16 ***
## Location:Species    2   0.469   0.235   1.8286    0.1609    
## Residuals        1748 224.335   0.128                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod2)
## 
## Call:
## lm(formula = log_Mean ~ Location * Species, data = tidy_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.20797 -0.23734  0.00445  0.20449  1.71655 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.55704    0.02209  70.485  < 2e-16 ***
## LocationEdge             -0.12592    0.03124  -4.031 5.80e-05 ***
## SpeciesEAME              -0.24820    0.03538  -7.015 3.28e-12 ***
## SpeciesGRSP              -0.50877    0.02785 -18.267  < 2e-16 ***
## LocationEdge:SpeciesEAME  0.07646    0.05004   1.528   0.1267    
## LocationEdge:SpeciesGRSP  0.06915    0.03939   1.756   0.0793 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3582 on 1748 degrees of freedom
## Multiple R-squared:  0.2617, Adjusted R-squared:  0.2596 
## F-statistic: 123.9 on 5 and 1748 DF,  p-value: < 2.2e-16
No Location by Species interaction. P value is .1609 (non-significant), and will not contribute to model fit. Location (p = 9.19e-06) and Species (p = < 2.2e-16), however both have effects on vegetation density that are significant.

Make a figure with results - illustrate significant factors in a plot OR re-evaluate analysis to a GLM in order to better fit distribution of the original data?

If plotting results from linear model - would need to back transform means from log transformed data

#Create new data frame with all factor combinations of significant interactions - leave out Species:Location as was not significant interaction.

anova_data <- expand.grid(
  Location = levels(as.factor(tidy_data$Location)),
  Species  = levels(as.factor(tidy_data$Species))
)

#Use predict function to calculate means and 95% CIs from log data based on anova output
#Predict creates data frame with mean and CL values as fit, upr, and lwr, respectively.
anova_loganalysis<- predict(mod2, newdata=anova_data, interval = "confidence", level = 0.95)

# Combine and back-transform data using exp() - 
results <- data.frame(
  Location       = anova_data$Location,
  Species        = anova_data$Species,
  Geometric_Mean = exp(anova_loganalysis[, "fit"]),
  lower.CL      = exp(anova_loganalysis[, "lwr"]),
  upper.CL       = exp(anova_loganalysis[, "upr"])
)

# Create a new plot with statistical analysis results using back-transformed data
ggplot(results, aes(x = Location, y = Geometric_Mean, group = Species, color = Species)) +
  geom_line() +
  geom_point(size = 3.5) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1, linewidth = 0.8) +
  labs(
    x = "Location",
    y = "Back-Transformed Average Vegetation Density (with 95% CI)"
  ) +
  theme_bw()


Plot shows: 

Lines between species are relatively parallel, meaning that differences in vegetation density at the edge vs. center are relatively consistent among species. All chose nest site centers with vegetation denser than the edge. This matches with non-significant Species:Location ANOVA result.

Plot also shows that DICK chooses nesting sites with higher overall vegetation density in nesting location and GRSP chooses sites with lower overall vegetation density. This is consistent with existing literature which found that vertical vegetation structure was the most important nest-site vegetation characteristic daily nest survival rates of DICK and litter and bare ground coverage was the most important nest-site vegetation feature for nest survival of GRSP (Frey et al., 2008) 

The slopes of all three lines decrease from Center to Edge, demonstrating the significant effect of center vs. edge location. The plot illustrates the vegetation density is indeed higher at nesting center vs edge. This conclusion confirms our initial hypothesis about vegetation density differences in center and edge. Nevertheless, the purpose of concealment and its relation to nest survival is not addressed in this analysis, and other studies have found concealment insufficient to explain daily nest survival (Ruth and Skagen, 2008).
Other thoughts:

Test if GLM is better fit model?

Factor in PlotName (Site) to deal with site non-independence?

Sources: 

Colombo, M. A., Depot, K. M. and Segura, L. N. (2024). Selection for Overhead Concealment Improves Nest Survival of a Ground Nesting Bird in Argentinian Rangelands. Rangel Ecol Manag 96, 47–55.

Frey, C. M., Jensen, W. E., Kimberly A. W.  (2008) Topographic Patterns of Nest Placement and Habitat Quality for Grassland Birds in Tallgrass Prairie. Source: The American Midland Naturalist 160, 220–234.

Kaplan, S. (2026). North American birds are dying off faster. It signals a human crisis, too. Washington Post.

Ruth, J. M. and Skagen, S. K. (2018). Reproductive response of Arizona Grasshopper Sparrows to weather patterns and habitat structure. Condor 120, 596–616.

Sroor, E. N., Boyle, W. A., Castro-Miller, H. N., Hobbs, K. S., Reynaldo, M. D., Silber, K. M., Mohankumar, N. M., Wright, N. A. and Freeman, N. E. (2025). Environmental drivers of female reproductive investment in egg quantity and quality in a grassland sparrow. Ibis 167, 437–451.