Species Abundance and Hunting Data Analysis
Dominik Marchowski
2024-11-04
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
In this analysis, we process data related to the abundance of selected bird species and their dynamics throughout the year. The datasets are available for download from the following links: - numbers.txt - Species abundance data - baza_V2.txt - Hunting data
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(ggplot2)Loading and Previewing Data
# 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# Load hunting data
hun <- read.table("https://marchowski.org/data/baza_V2.txt", header = TRUE)
head(hun)## Pic spec no HID seson prot
## 1 O1 ANACRE 1 1 O 0
## 2 O1 ANAPLA 3 1 O 0
## 3 O1 ANSALB 27 1 O 0
## 4 O1 ANSFAB 8 1 O 0
## 5 O1 ANSSP 8 1 O 0
## 6 O1 MARPEN 17 1 O 1Calculating Mean Species Abundance
Mean Abundance for Each Species
species_mean <- num %>%
group_by(Spec) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
print(species_mean)## # A tibble: 14 × 2
## Spec mean_count
## <chr> <dbl>
## 1 ANAACU 42.9
## 2 ANACRE 63.5
## 3 ANAPLA 1457.
## 4 ANSALB 1527.
## 5 ANSANS 352.
## 6 ANSFAB 2645.
## 7 AYTFER 103.
## 8 AYTFUL 184.
## 9 AYTMAR 8.33
## 10 FULATR 886.
## 11 MARPEN 168.
## 12 MARSTR 66.7
## 13 MERALB 5.59
## 14 SPACLY 11.5Filling Missing Data
The following function fills in missing values by using the mean of the previous and next season:
fill_na_with_season_mean <- function(data) {
data %>%
group_by(Spec) %>%
arrange(Seson, Mont) %>%
mutate(count = ifelse(is.na(count),
rowMeans(cbind(lag(count), lead(count)), na.rm = TRUE),
count))
}
num_filled <- fill_na_with_season_mean(num)Analyzing Seasons 2005 and 2008
season_2005 <- num_filled %>% filter(Seson == 2005)
season_2008 <- num_filled %>% filter(Seson == 2008)
cat("Data for the 2005 season:\n")## Data for the 2005 season:print(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 rowscat("\nData for the 2008 season:\n")##
## Data for the 2008 season:print(season_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 rowsMean Abundance for the First and Last 6 Years
First 6 Years (2002-2007)
first_six_years <- num %>% filter(Seson >= 2002 & Seson <= 2007)
species_mean_first_six_years <- first_six_years %>%
group_by(Spec) %>%
summarise(mean_count_2002_2007 = mean(count, na.rm = TRUE))
print(species_mean_first_six_years)## # A tibble: 14 × 2
## Spec mean_count_2002_2007
## <chr> <dbl>
## 1 ANAACU 21.3
## 2 ANACRE 129.
## 3 ANAPLA 2278.
## 4 ANSALB 3456.
## 5 ANSANS 38.9
## 6 ANSFAB 5558.
## 7 AYTFER 154.
## 8 AYTFUL 215.
## 9 AYTMAR 0.154
## 10 FULATR 1554.
## 11 MARPEN 74.2
## 12 MARSTR 67.8
## 13 MERALB 6.46
## 14 SPACLY 5.31Last 6 Years
last_six_years <- num %>% filter(Seson >= (max(Seson) - 5) & Seson <= max(Seson))
species_mean_last_six_years <- last_six_years %>%
group_by(Spec) %>%
summarise(mean_count_last_six_years = mean(count, na.rm = TRUE))
print(species_mean_last_six_years)## # A tibble: 14 × 2
## Spec mean_count_last_six_years
## <chr> <dbl>
## 1 ANAACU 10.6
## 2 ANACRE 25.1
## 3 ANAPLA 1817
## 4 ANSALB 1611.
## 5 ANSANS 656.
## 6 ANSFAB 2403.
## 7 AYTFER 45.4
## 8 AYTFUL 233.
## 9 AYTMAR 1.72
## 10 FULATR 826.
## 11 MARPEN 84.4
## 12 MARSTR 26.2
## 13 MERALB 9.44
## 14 SPACLY 12.4Comparison of Mean Counts and Percentage Change
comparison_table <- merge(species_mean_first_six_years, species_mean_last_six_years, by = "Spec", all = TRUE) %>%
mutate(percent_change = ((mean_count_last_six_years - mean_count_2002_2007) / mean_count_2002_2007) * 100)
print(comparison_table)## Spec mean_count_2002_2007 mean_count_last_six_years percent_change
## 1 ANAACU 21.3076923 10.555556 -50.46129
## 2 ANACRE 129.4615385 25.055556 -80.64633
## 3 ANAPLA 2277.7692308 1817.000000 -20.22897
## 4 ANSALB 3456.0769231 1611.333333 -53.37681
## 5 ANSANS 38.9230769 656.500000 1586.66008
## 6 ANSFAB 5558.3076923 2403.333333 -56.76142
## 7 AYTFER 153.7692308 45.388889 -70.48246
## 8 AYTFUL 215.3846154 233.277778 8.30754
## 9 AYTMAR 0.1538462 1.722222 1019.44444
## 10 FULATR 1554.1538462 826.500000 -46.81994
## 11 MARPEN 74.2307692 84.444444 13.75936
## 12 MARSTR 67.7692308 26.166667 -61.38857
## 13 MERALB 6.4615385 9.444444 46.16402
## 14 SPACLY 5.3076923 12.388889 133.41385Plotting Abundance for Last 6 Years
# Boxplot of abundance by species over the last 6 years
ggplot(last_six_years, aes(x = Spec, y = count)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Species Counts in the Last Six Years",
x = "Species",
y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Saving a High-Resolution Plot with Log Scale
# Save the plot as a high-resolution TIFF
tiff("species_counts_last_six_years.tif", width = 10, height = 7, units = "in", res = 300)
ggplot(last_six_years, aes(x = Spec, y = count)) +
geom_boxplot() +
scale_y_log10() +
theme_bw() +
labs(x = "Species",
y = "Log10(Count)") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 14))## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 98 rows containing non-finite outside the scale range
## (`stat_boxplot()`).dev.off()## png
## 2Loading Hunting Data and Calculating Seasonal Averages
To continue, we analyze the hunting data, calculating average counts per season.
# Calculate the average population for each species by dividing the population into 3 seasons
average_per_season <- hun %>%
group_by(spec) %>%
summarise(mean_count_per_season = sum(no) / 3)
print(average_per_season)## # A tibble: 16 × 2
## spec mean_count_per_season
## <chr> <dbl>
## 1 ANAACU 3
## 2 ANACRE 34
## 3 ANAPLA 104
## 4 ANASPE 0.667
## 5 ANSALB 144
## 6 ANSANS 48.7
## 7 ANSFAB 51
## 8 ANSSP 8
## 9 AYTFER 19.7
## 10 AYTFUL 2
## 11 AYTMAR 0.333
## 12 FULATR 1
## 13 MARPEN 17.7
## 14 MARSTRE 8.67
## 15 MERALB 0.333
## 16 SPACLY 0.667Modeling
We’ll perform linear regression analyses with a logarithmic transformation on the dependent variable hunt_hip, which represents relative species dominance in the environment.
# Real data: Relative dominance of individual species in the environment (last6) and in hunted birds (hunt)
last6 <- c(0.135972639, 0.322756099, 23.40589978, 20.75658036, 8.456782172, 30.95882176,
0.584682326, 3.004995208, 0.022185006, 10.64665722, 1.087781064, 0.337069007,
0.121659719, 0.159588935)
hunt_hip <- c(0, 9, 23, 33, 13, 15, 5, 1, 0, 1, 0, 0, 0, 0)
hunt <- c(0.689655172, 7.816091954, 23.90804598, 33.10344828, 11.1954023, 11.72413793,
4.528735632, 0.459770115, 0.076551724, 0.229885057, 4.068965517, 1.993103448,
0.076551724, 0.153333333)
prot <- c(1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1)
# Creating a data frame
data <- data.frame(last6, hunt, prot)Applying Logarithmic Transformation and Running Linear Regression
We’ll use a logarithmic transformation on hunt_hip to avoid issues with log(0), then fit a linear model.
# Applying logarithmic transformation to the dependent variable (hunt_hip)
log_hunt_hip <- log(hunt_hip + 1) # Adding 1 to avoid log(0)
# Performing linear regression on logarithmized data
log_model <- lm(log_hunt_hip ~ last6 + prot, data = data)
# Display model summary
summary(log_model)##
## Call:
## lm(formula = log_hunt_hip ~ last6 + prot, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41720 -0.03246 0.01002 0.33267 0.85925
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.52403 0.36624 4.161 0.00159 **
## last6 0.05507 0.02243 2.455 0.03197 *
## prot -1.54114 0.45630 -3.377 0.00617 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6835 on 11 degrees of freedom
## Multiple R-squared: 0.7906, Adjusted R-squared: 0.7525
## F-statistic: 20.77 on 2 and 11 DF, p-value: 0.0001842# Diagnostic plots for the logarithmic model (log_model)
par(mfrow = c(2, 2))
plot(log_model, main = "Diagnostic Plots for log_model")
par(mfrow = c(1, 1)) # Restore default plotting layout
# Testing normality of residuals
shapiro.test(residuals(log_model))##
## Shapiro-Wilk normality test
##
## data: residuals(log_model)
## W = 0.90247, p-value = 0.1225Running Linear Regression on Real Data
We’ll run a similar analysis on the hunt variable for comparison.
# Logarithmized dependent variable for real data
log_hunt <- log(hunt + 1) # Adding 1 to avoid log(0)
# Performing linear regression on logarithmized real data
log_model_real <- lm(log_hunt ~ last6 + prot, data = data)
# Model summary for real data
summary(log_model_real)##
## Call:
## lm(formula = log_hunt ~ last6 + prot, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7295 -0.5026 0.1584 0.6513 0.9928
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.30433 0.48742 2.676 0.0216 *
## last6 0.05937 0.02986 1.988 0.0722 .
## prot -0.73375 0.60729 -1.208 0.2523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9097 on 11 degrees of freedom
## Multiple R-squared: 0.5337, Adjusted R-squared: 0.4489
## F-statistic: 6.294 on 2 and 11 DF, p-value: 0.01506# Diagnostic plots for log_model_real
par(mfrow = c(2, 2))
plot(log_model_real, main = "Diagnostic Plots for log_model_real")
par(mfrow = c(1, 1)) # Restore default plotting layout
# Testing normality of residuals for log_model_real
shapiro.test(residuals(log_model_real))##
## Shapiro-Wilk normality test
##
## data: residuals(log_model_real)
## W = 0.92526, p-value = 0.2615Total Count by Year and Month
We calculate the total count for each year and month, then visualize the data with a bar plot.
# Calculate the total count for each year and month
total_counts_by_year_month <- num %>%
group_by(Seson, Mont) %>%
summarise(total_count = sum(count, na.rm = TRUE)) %>%
arrange(desc(total_count))## `summarise()` has grouped output by 'Seson'. You can override using the
## `.groups` argument.print(total_counts_by_year_month)## # A tibble: 66 × 3
## # Groups: Seson [22]
## Seson Mont total_count
## <int> <chr> <int>
## 1 2003 nov 36087
## 2 2004 jan 32635
## 3 2021 jan 26339
## 4 2003 mar 18728
## 5 2004 mar 17790
## 6 2021 nov 17675
## 7 2005 nov 17177
## 8 2018 mar 15495
## 9 2004 nov 13940
## 10 2020 jan 13849
## # ℹ 56 more rows# Calculate the total count for each year
total_counts_by_year <- num %>%
group_by(Seson) %>%
summarise(total_count = sum(count, na.rm = TRUE)) # Summing counts for each year
print(total_counts_by_year) # Display total counts by year## # A tibble: 22 × 2
## Seson total_count
## <int> <int>
## 1 2003 59729
## 2 2004 64365
## 3 2005 17177
## 4 2006 18898
## 5 2007 16099
## 6 2008 0
## 7 2009 7892
## 8 2010 12501
## 9 2011 15192
## 10 2012 8456
## # ℹ 12 more rows# Bar plot with ggplot2
ggplot(total_counts_by_year_month, aes(x = Seson, y = total_count, fill = Mont)) +
geom_bar(stat = "identity", position = "dodge") +
labs(x = "Year", y = "Total Count", title = "Total Count by Year and Month") +
theme_bw()
Calculating Percentage Change in Abundance
We calculate the average abundance for the first and last six years to find the percentage change.
# Define first 6 years
first_six_years <- total_counts_by_year %>%
filter(Seson %in% sort(unique(Seson))[1:6])
# Define last 6 years
last_six_years <- total_counts_by_year %>%
filter(Seson %in% sort(unique(Seson), decreasing = TRUE)[1:6])
# Calculate averages
mean_first_six_years <- mean(first_six_years$total_count, na.rm = TRUE)
mean_last_six_years <- mean(last_six_years$total_count, na.rm = TRUE)
# Calculate percentage change
percent_change <- ((mean_last_six_years - mean_first_six_years) / mean_first_six_years) * 100
# Display results
cat("Average number for the first 6 years:", mean_first_six_years, "\n")## Average number for the first 6 years: 29378cat("Average number for the last 6 years:", mean_last_six_years, "\n")## Average number for the last 6 years: 23564.83cat("Percentage change:", percent_change, "%\n")## Percentage change: -19.78748 %Conclusion
This RMarkdown document provided an analysis of species abundance and hunting data, including data preprocessing, modeling, and visualization of trends over selected periods.
This document is structured to cover all analysis steps with code, comments, and outputs to interpret the data effectively in RMarkdown.