Trends of waterfowl at Miedwie Site 2002-2024
Dominik Marchowski
2024-11-05
The impact of hunting: Hunters’ inability to identify endangered and protected species
Authors: Dominik Marchowski, Robert Jurszo, Paweł Stańczak, Michał Jasiński, Sebastain Guentzel
Introduction
This document analyzes the abundance trends of various bird species across multiple years. It includes visualizations and statistical summaries for total bird counts as well as specific species.
Installing and Loading Necessary Packages
# Uncomment the lines below to install necessary packages
# install.packages("ggplot2")
# install.packages("dplyr")
# install.packages("gridExtra")
library(ggplot2)
library(dplyr)##
## Dołączanie pakietu: 'dplyr'## Następujące obiekty zostały zakryte z 'package:stats':
##
## filter, lag## Następujące obiekty zostały zakryte z 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(gridExtra)## Warning: pakiet 'gridExtra' został zbudowany w wersji R 4.4.1##
## Dołączanie pakietu: 'gridExtra'## Następujący obiekt został zakryty z 'package:dplyr':
##
## combineLoading and Preparing the Data
We begin by loading the bird species abundance data and applying a function to fill any missing values with the mean of adjacent seasons.
# Load species abundance data
num <- read.table("https://marchowski.org/data/numbers.txt", header = TRUE)
head(num)## Spec Seson Mont count site
## 1 ANAACU 2003 jan 0 Miedwie
## 2 ANACRE 2003 jan 0 Miedwie
## 3 ANAPLA 2003 jan 1402 Miedwie
## 4 ANSALB 2003 jan 0 Miedwie
## 5 ANSANS 2003 jan 0 Miedwie
## 6 ANSFAB 2003 jan 3500 Miedwie# Function to replace NA with the mean of previous and next season
fill_na_with_season_mean <- function(data) {
data <- data %>%
group_by(Spec) %>%
arrange(Seson, Mont) %>%
mutate(count = ifelse(is.na(count),
rowMeans(cbind(lag(count, order_by = Seson), lead(count, order_by = Seson)), na.rm = TRUE),
count))
return(data)
}
# Apply the function to the dataset
num_filled <- fill_na_with_season_mean(num)Seasonal Data for Selected Years (2005 and 2008)
# Display data for 2005 and 2008
season_2005 <- num_filled %>% filter(Seson == 2005)
season_2008 <- num_filled %>% filter(Seson == 2008)
season_2005## # A tibble: 42 × 5
## # Groups: Spec [14]
## Spec Seson Mont count site
## <chr> <int> <chr> <dbl> <chr>
## 1 ANAACU 2005 jan 5 Miedwie
## 2 ANACRE 2005 jan 71 Miedwie
## 3 ANAPLA 2005 jan 1189 Miedwie
## 4 ANSALB 2005 jan 1200 Miedwie
## 5 ANSANS 2005 jan 9 Miedwie
## 6 ANSFAB 2005 jan 2800 Miedwie
## 7 AYTFER 2005 jan 236 Miedwie
## 8 AYTFUL 2005 jan 90 Miedwie
## 9 AYTMAR 2005 jan 0 Miedwie
## 10 FULATR 2005 jan 8257 Miedwie
## # ℹ 32 more rowsseason_2008## # A tibble: 42 × 5
## # Groups: Spec [14]
## Spec Seson Mont count site
## <chr> <int> <chr> <dbl> <chr>
## 1 ANAACU 2008 jan 5 Miedwie
## 2 ANACRE 2008 jan 20 Miedwie
## 3 ANAPLA 2008 jan 1300 Miedwie
## 4 ANSALB 2008 jan 1000 Miedwie
## 5 ANSANS 2008 jan 100 Miedwie
## 6 ANSFAB 2008 jan 3200 Miedwie
## 7 AYTFER 2008 jan 30 Miedwie
## 8 AYTFUL 2008 jan 120 Miedwie
## 9 AYTMAR 2008 jan 0 Miedwie
## 10 FULATR 2008 jan 890 Miedwie
## # ℹ 32 more rowsTotal Bird Count Analysis
Calculating Total Counts Per Year
# Calculate the total counts for all species by year
total_counts_per_year <- num_filled %>%
group_by(Seson) %>%
summarise(total_count = sum(count, na.rm = TRUE))
# Display the results
print(total_counts_per_year, n = 22)## # A tibble: 22 × 2
## Seson total_count
## <int> <dbl>
## 1 2003 59729
## 2 2004 64365
## 3 2005 48294
## 4 2006 18898
## 5 2007 16099
## 6 2008 9016
## 7 2009 7892
## 8 2010 12501
## 9 2011 15192
## 10 2012 8456
## 11 2013 19550
## 12 2014 8433
## 13 2015 14359
## 14 2016 6484
## 15 2017 14067
## 16 2018 34148
## 17 2019 12393
## 18 2020 27040
## 19 2021 46743
## 20 2022 23859
## 21 2023 9392
## 22 2024 21962Standard Error of Total Counts
# Total counts data
total_count <- c(59729, 64365, 48294, 18898, 16099, 9016, 7892, 12501, 15192, 8456,
19550, 8433, 14359, 6484, 14067, 34148, 12393, 27040, 46743, 23859,
9392, 21962)
# Calculate Standard Error
standard_error <- sd(total_count) / sqrt(length(total_count))
print(standard_error)## [1] 3670.318Visualization of Total Abundance Trend
# Visualization of the total counts trend
total_plot <- ggplot(total_counts_per_year, aes(x = Seson, y = total_count)) +
geom_point(size = 3, color = "grey") +
geom_line(aes(group = 1), color = "grey33") +
geom_smooth(method = "loess", span = 0.6, color = "slateblue", size = 1.2, se = TRUE) +
labs(title = "Total counts for all species", x = "Year", y = "Total Abundance") +
theme_bw(base_size = 15) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
annotate("text", x = 2015, y = max(total_counts_per_year$total_count, na.rm = TRUE) * 1,
label = "Moderate decline (p<0.01) λ=0.9803", size = 4, hjust = 0.55)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.total_plot## `geom_smooth()` using formula = 'y ~ x'
Species-Specific Analysis
Each plot below shows the abundance trend of a specific bird species.
Anser anser (ANSANS)
# Calculate and plot for Anser anser
mean_counts_ansans <- num_filled %>%
filter(Spec == "ANSANS") %>%
group_by(Seson) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
ans <- mean_counts_ansans
ans_plot <- ggplot(ans, aes(x = Seson, y = mean_count)) +
geom_point(size = 3, color = "grey") +
geom_line(aes(group = 1), color = "grey33") +
geom_smooth(method = "loess", span = 0.6, color = "slateblue", size = 1.2, se = FALSE) +
labs(title = "Anser anser", x = "Year", y = "Abundance") +
theme_bw(base_size = 15) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(limits = c(0, max(ans$mean_count, na.rm = TRUE) * 1.1)) +
annotate("text", x = 2015, y = max(ans$mean_count, na.rm = TRUE) * 1.05,
label = "Strong increase (p<0.01) ** λ=1.1423 SD=0.0341", size = 4, hjust = 0.5)
ans_plot## `geom_smooth()` using formula = 'y ~ x'
Anas crecca (ANACRE)
# Calculate and plot for Anas crecca
mean_counts_anacre <- num_filled %>%
filter(Spec == "ANACRE") %>%
group_by(Seson) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
acre <- mean_counts_anacre
acre_plot <- ggplot(acre, aes(x = Seson, y = mean_count)) +
geom_point(size = 3, color = "grey") +
geom_line(aes(group = 1), color = "grey33") +
geom_smooth(method = "loess", span = 0.6, color = "slateblue", size = 1.2, se = FALSE) +
labs(title = "Anas crecca", x = "Year", y = "Abundance") +
theme_bw(base_size = 15) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
annotate("text", x = 2015, y = max(acre$mean_count, na.rm = TRUE) * 1.05,
label = "Moderate decline (p<0.05) * λ = 0.909 SD = 0.04", size = 4, hjust = 0.55)
acre_plot## `geom_smooth()` using formula = 'y ~ x'
Anas platyrhynchos (ANAPLA)
# Calculate the mean counts for ANAPLA across seasons in the filled dataset
mean_counts_anapla <- num_filled %>%
filter(Spec == "ANAPLA") %>%
group_by(Seson) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
# Rename mean_counts_anapla to anp
anp <- mean_counts_anapla
# Example visualization of the abundance trend for the species Anas platyrhynchos with Y-axis limit
anp_plot <- ggplot(anp, aes(x = Seson, y = mean_count)) +
geom_point(size = 3, color = "grey") + # Plot observation points in grey
geom_line(aes(group = 1), color = "grey33") + # Line connecting points
geom_smooth(method = "loess", span = 0.6, color = "slateblue", size = 1.2, se = FALSE) + # LOESS trend line without shading
labs(title = "Anas platyrhynchos", x = "Year", y = "Abundance") + # Axis labels and title
theme_bw(base_size = 15) + # Apply minimalist theme with base font size 15
theme(plot.title = element_text(hjust = 0.5)) + # Center the plot title
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) + # Adjust X-axis to give slight padding on both ends
scale_y_continuous(limits = c(0, max(anp$mean_count, na.rm = TRUE) * 1.1)) + # Set Y-axis limits, starting at 0 and extending 10% above max count
annotate("text", x = 2015, y = max(anp$mean_count, na.rm = TRUE) * 1.05, # Position annotation text slightly below max Y-axis limit
label = "Moderate decline λ=0.9877 SD=0.0013", size = 4, hjust = 0.55)
anp_plot## `geom_smooth()` using formula = 'y ~ x'
Anser albifrons (ANSALB)
# Calculate and plot for Anser albifrons
mean_counts_ansalb <- num_filled %>%
filter(Spec == "ANSALB") %>%
group_by(Seson) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
alb <- mean_counts_ansalb
alb_plot <- ggplot(alb, aes(x = Seson, y = mean_count)) +
geom_point(size = 3, color = "grey") +
geom_line(aes(group = 1), color = "grey33") +
geom_smooth(method = "loess", span = 0.6, color = "slateblue", size = 1.2, se = FALSE) +
labs(title = "Anser albifrons", x = "Year", y = "Abundance") +
theme_bw(base_size = 15) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(limits = c(0, max(alb$mean_count, na.rm = TRUE) * 1.1)) +
annotate("text", x = 2015, y = max(alb$mean_count, na.rm = TRUE) * 1.05,
label = "Steep decline λ=0.9274/ Stable from 2006", size = 4, hjust = 0.55)
alb_plot## `geom_smooth()` using formula = 'y ~ x'
Anser fabalis sensu lato (ANSFAB)
# Calculate and plot for Anser fabalis sensu lato
mean_counts_ansfab <- num_filled %>%
filter(Spec == "ANSFAB") %>%
group_by(Seson) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
af <- mean_counts_ansfab
af_plot <- ggplot(af, aes(x = Seson, y = mean_count)) +
geom_point(size = 3, color = "grey") +
geom_line(aes(group = 1), color = "grey33") +
geom_smooth(method = "loess", span = 0.6, color = "slateblue", size = 1.2, se = FALSE) +
labs(title = "Anser fabalis sensu lato", x = "Year", y = "Abundance") +
theme_bw(base_size = 15) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(limits = c(0, max(af$mean_count, na.rm = TRUE) * 1.1)) +
annotate("text", x = 2015, y = max(af$mean_count, na.rm = TRUE) * 1.05,
label = "Moderate decline λ=0.9522 SD=0.0019", size = 4, hjust = 0.55)
af_plot## `geom_smooth()` using formula = 'y ~ x'
Combined Species Plot
# Combine all species plots in a grid layout
grid.arrange(total_plot, ans_plot, acre_plot, anp_plot, alb_plot, af_plot, ncol = 2)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'