Point Count ANALYSIS

POINT COUNT SURVEYS

Objective

To determine the optimal methodologies for measuring ecosystem restoration and land degradation gradient by comparing soundscapes and traditional methods

Sampling design

The study area comprises a 10 x 10 km square, which includes Enarau Conservancy and the surrounding landscape. The study area is divided into 2.5 x 2.5 km tiles (n=16). Within each tile, a 1 km2 cluster is generated. A 1-km transect was selected randomly for bird point count.

The sampling method used is a stratified random sampling design to ensure representative coverage of the study area. Observation points are spaced at least 200 meters apart within each 1-km transect to reduce the risk of double-counting. We observed for 10 minutes, recording all birds seen or heard within the radius specified on the data sheet (0-30 metres, 30-100 metres, above 100 metres). We used binoculars to identify distant birds. If field identification of the species was impossible, we took photos of the unidentified species and/or recorded their call for later verification using Merlin Bird ID. We walked 200 metres along the transect to the next point, ensuring minimal noise.

Study area and transect

NOTE

In this script, column season has Apr-2024, Jun-2024, and Oct 2024 as rows. I will change the rows to have dry and wet seasons after I add the January 2025 data.

Data and libraries

# Load libraries
library(tidyverse)
library(vegan)
library(lme4)
library(Distance)
library(ggpubr)
library(ggthemes)
library(sf)
library(dplyr)
library(lmerTest)
library(ggplot2)
library(stringr)
library(tidyr)
library(readr)
library(gridExtra)
library(lme4)
library(lmerTest)
library(kableExtra)
library(cluster)
library(factoextra)  

# Load the dataset
birds.data <- read_csv("C:/Users/evt3uzala/Downloads/Pointcountdata.csv")


# Convert categorical variables to factors and recode dist.band and count.station
birds.data <- birds.data %>% 
  mutate(
    season = as.factor(Season),
    visit = as.factor(Visit),
    number.spp= as.numeric(Number),
    cluster = as.factor(Cluster),
    count.station = recode_factor(`Point count station`, `1` = "0m", `2` = "200m", `3` = "400m", `4` = "600m", `5` = "800m", `6` = "1000m"),
    dist.band = recode_factor(`Distance band`, `1` = "0-30m", `2` = "30-100m", `3` = "100m+"),
    species = as.factor(Species)
  )
birds.data <- birds.data %>% drop_na(Site)
View(birds.data)

Species metrics

# Species Richness
birds.spp.richness <- birds.data %>%
  group_by(cluster, season) %>%
  summarise(richness = n_distinct(species), .groups = "drop") %>%
  arrange(desc(richness))

# Species Diversity (Shannon & Simpson)
birds.spp.diversity <- birds.data %>%
  group_by(cluster, season) %>%
  summarise(
    bird.spp.diversity.H = diversity(number.spp, index = "shannon"),
    bird.spp.diversity.Simp = diversity(number.spp, index = "simpson"),
    .groups = "drop"
  )

# Distance-based Density Estimation
birds_data_filtered <- birds.data %>%
  mutate(distance = case_when(
    dist.band == "0-30m" ~ 15,    
    dist.band == "30-100m" ~ 65,  
    dist.band == "100m+" ~ 120    
  )) %>%
  filter(!is.na(distance))  

hr_model <- ds(birds_data_filtered, key = "hr", adjustment = "cos")

Error in ddf.ds(dsmodel = dsmodel, data = data, meta.data = meta.data,  : 
  More parameters to estimate than unique distances

summary(hr_model)

Summary for distance analysis 
Number of observations :  4637 
Distance range         :  0  -  120 

Model       : Hazard-rate key function with cosine adjustment term of order 2 

Strict monotonicity constraints were enforced.
AIC         :  42800.98 
Optimisation:  mrds (slsqp) 

Detection function parameters
Scale coefficient(s):  
            estimate       se
(Intercept) 3.514675 20.67857

Shape coefficient(s):  
             estimate       se
(Intercept) 0.1575317 11.10244

Adjustment term coefficient(s):  
              estimate       se
cos, order 2 0.2245542 2.665129

                        Estimate           SE       CV
Average p           4.123688e-01     2.812233 6.819705
N in covered region 1.124479e+04 76686.243201 6.819714

density_per_cluster_season <- birds_data_filtered %>%
  group_by(cluster, season) %>%
  summarise(
    mean_detection_prob = mean(predict(hr_model)$fitted, na.rm = TRUE),
    abundance = sum(number.spp) / mean_detection_prob,
    density = abundance / (pi * max(distance)^2),
    .groups = "drop"
  )

#Species abundance

# Add a species code 
birds.data <- birds.data %>%
  mutate(species.code = paste0(str_extract(species, "^\\w{2}"), ".", str_extract(species, "(?<=\\s)\\w{3}")))

# Summarise species abundance per cluster and season
bird.spp.abundance <- birds.data %>%
  group_by(cluster, season, species.code) %>%
  summarise(bird.spp.abundance = sum(number.spp, na.rm = TRUE), .groups = "drop")

# Print the species abundance table
head(bird.spp.abundance)

# A tibble: 6 × 4
  cluster season species.code bird.spp.abundance
  <fct>   <fct>  <chr>                     <dbl>
1 1       Apr-24 Af.Gre                        3
2 1       Apr-24 Af.Pie                        7
3 1       Apr-24 Af.Woo                        2
4 1       Apr-24 Am.Sun                        3
5 1       Apr-24 Au.Buz                        3
6 1       Apr-24 Ba.Swa                        4

Statistical analysis

##NOTE April treated as reference level


# Mixed-effects model for species richness with cluster as a fixed effect
richness_model <- lmer(richness ~ season + cluster + (1 | cluster), data = birds.spp.richness)
richness_summary <- summary(richness_model)

# Mixed-effects model for species diversity (Shannon) with cluster as a fixed effect
diversity_model_H <- lmer(bird.spp.diversity.H ~ season + cluster + (1 | cluster), data = birds.spp.diversity)
diversity_summary <- summary(diversity_model_H)

# Mixed-effects model for density with cluster as a fixed effect
density_model <- lmer(density ~ season + cluster + (1 | cluster), data = density_per_cluster_season)
density_summary <- summary(density_model)

# Mixed-effects model for species abundance with cluster as a fixed effect
abundance_model <- lmer(bird.spp.abundance ~ season + cluster + (1 | cluster), data = bird.spp.abundance)
abundance_summary <- summary(abundance_model)

# Function to extract model coefficients and format as a table
create_model_table <- function(model_summary, title) {
  coef_table <- as.data.frame(model_summary$coefficients)
  coef_table <- coef_table %>%
    mutate(
      Term = rownames(coef_table),
      Estimate = round(Estimate, 3),
      `Std. Error` = round(`Std. Error`, 3),
      `t value` = round(`t value`, 3),
      `p-value` = ifelse(`Pr(>|t|)` < 0.001, "<0.001", round(`Pr(>|t|)`, 3))
    ) %>%
    select(Term, Estimate, `Std. Error`, `t value`, `p-value`)
  
  # Create a kable table
  kable(coef_table, caption = title, align = "c") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
}

# Create tables for each model
richness_table <- create_model_table(richness_summary, "Mixed-Effects Model for Species Richness (Cluster as Fixed Effect)")
diversity_table <- create_model_table(diversity_summary, "Mixed-Effects Model for Species Diversity (Shannon) (Cluster as Fixed Effect)")
density_table <- create_model_table(density_summary, "Mixed-Effects Model for Density (Cluster as Fixed Effect)")
abundance_table <- create_model_table(abundance_summary, "Mixed-Effects Model for Species Abundance (Cluster as Fixed Effect)")

# Display the tables
richness_table

Mixed-Effects Model for Species Richness (Cluster as Fixed Effect)

	Term	Estimate	Std. Error	t value	p-value
(Intercept)	(Intercept)	32.271	5.642	5.719	<0.001
seasonJun-24	seasonJun-24	19.000	2.130	8.920	<0.001
seasonOct-24	seasonOct-24	22.188	2.130	10.416	<0.001
cluster2	cluster2	-4.000	7.788	-0.514	0.611
cluster3	cluster3	-2.333	7.788	-0.300	0.767
cluster4	cluster4	4.333	7.788	0.556	0.582
cluster5	cluster5	-2.000	7.788	-0.257	0.799
cluster6	cluster6	-8.333	7.788	-1.070	0.293
cluster7	cluster7	0.333	7.788	0.043	0.966
cluster8	cluster8	0.000	7.788	0.000	1
cluster9	cluster9	5.333	7.788	0.685	0.499
cluster10	cluster10	8.667	7.788	1.113	0.275
cluster11	cluster11	-5.000	7.788	-0.642	0.526
cluster12	cluster12	-8.000	7.788	-1.027	0.313
cluster13	cluster13	-2.333	7.788	-0.300	0.767
cluster14	cluster14	-1.667	7.788	-0.214	0.832
cluster15	cluster15	-4.333	7.788	-0.556	0.582
cluster16	cluster16	-3.000	7.788	-0.385	0.703

diversity_table

Mixed-Effects Model for Species Diversity (Shannon) (Cluster as Fixed Effect)

	Term	Estimate	Std. Error	t value	p-value
(Intercept)	(Intercept)	3.717	0.172	21.575	<0.001
seasonJun-24	seasonJun-24	0.975	0.079	12.279	<0.001
seasonOct-24	seasonOct-24	0.835	0.079	10.516	<0.001
cluster2	cluster2	-0.126	0.235	-0.537	0.595
cluster3	cluster3	-0.152	0.235	-0.647	0.523
cluster4	cluster4	-0.019	0.235	-0.079	0.938
cluster5	cluster5	-0.702	0.235	-2.987	0.006
cluster6	cluster6	-0.262	0.235	-1.114	0.274
cluster7	cluster7	-0.221	0.235	-0.943	0.353
cluster8	cluster8	0.058	0.235	0.248	0.806
cluster9	cluster9	-0.208	0.235	-0.887	0.382
cluster10	cluster10	0.060	0.235	0.255	0.8
cluster11	cluster11	-0.619	0.235	-2.637	0.013
cluster12	cluster12	-0.252	0.235	-1.071	0.293
cluster13	cluster13	-0.190	0.235	-0.809	0.425
cluster14	cluster14	-0.376	0.235	-1.603	0.119
cluster15	cluster15	-0.286	0.235	-1.219	0.232
cluster16	cluster16	0.014	0.235	0.062	0.951

density_table

Mixed-Effects Model for Density (Cluster as Fixed Effect)

	Term	Estimate	Std. Error	t value	p-value
(Intercept)	(Intercept)	0.004	0.004	0.866	0.394
seasonJun-24	seasonJun-24	0.011	0.002	5.822	<0.001
seasonOct-24	seasonOct-24	0.011	0.002	5.904	<0.001
cluster2	cluster2	-0.001	0.006	-0.175	0.862
cluster3	cluster3	-0.003	0.006	-0.557	0.582
cluster4	cluster4	0.004	0.006	0.641	0.526
cluster5	cluster5	0.007	0.006	1.201	0.239
cluster6	cluster6	-0.003	0.006	-0.500	0.621
cluster7	cluster7	0.002	0.006	0.303	0.764
cluster8	cluster8	0.001	0.006	0.134	0.894
cluster9	cluster9	0.003	0.006	0.591	0.559
cluster10	cluster10	-0.001	0.006	-0.150	0.882
cluster11	cluster11	-0.002	0.006	-0.331	0.743
cluster12	cluster12	-0.003	0.006	-0.519	0.608
cluster13	cluster13	0.000	0.006	0.034	0.973
cluster14	cluster14	-0.001	0.006	-0.131	0.896
cluster15	cluster15	0.008	0.006	1.357	0.185
cluster16	cluster16	-0.002	0.006	-0.275	0.785

abundance_table

Mixed-Effects Model for Species Abundance (Cluster as Fixed Effect)

	Term	Estimate	Std. Error	t value	p-value
(Intercept)	(Intercept)	2.341	0.998	2.346	1
seasonJun-24	seasonJun-24	3.201	0.584	5.485	<0.001
seasonOct-24	seasonOct-24	2.910	0.579	5.026	<0.001
cluster2	cluster2	-0.063	1.341	-0.047	1
cluster3	cluster3	-1.336	1.325	-1.008	1
cluster4	cluster4	0.918	1.287	0.713	1
cluster5	cluster5	2.837	1.324	2.142	1
cluster6	cluster6	-0.483	1.371	-0.352	1
cluster7	cluster7	0.477	1.310	0.364	1
cluster8	cluster8	0.159	1.312	0.121	1
cluster9	cluster9	0.511	1.283	0.399	1
cluster10	cluster10	-1.132	1.264	-0.896	1
cluster11	cluster11	-0.545	1.346	-0.405	1
cluster12	cluster12	-0.726	1.369	-0.530	1
cluster13	cluster13	0.191	1.325	0.144	1
cluster14	cluster14	-0.342	1.320	-0.259	1
cluster15	cluster15	3.771	1.339	2.816	1
cluster16	cluster16	-0.458	1.329	-0.344	1

Visualisation

# Calculate mean and standard error for Species Richness
richness_summary <- birds.spp.richness %>%
  group_by(season, cluster) %>%
  summarise(
    mean_richness = mean(richness),
    se_richness = sd(richness) / sqrt(n()),
    .groups = "drop"
  )

# Visualise Species Richness with Error Bars and Cluster Numbers
ggplot(richness_summary, aes(x = season, y = mean_richness, color = cluster)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = mean_richness - se_richness, ymax = mean_richness + se_richness),
                width = 0.2, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = cluster), position = position_dodge(width = 0.5), vjust = -1.5, size = 4) +
  labs(title = "Species Richness by Season and Cluster",
       x = "Season",
       y = "Mean Species Richness") +
  theme_minimal()

# Calculate mean and standard error for Density
density_summary <- density_per_cluster_season %>%
  group_by(season, cluster) %>%
  summarise(
    mean_density = mean(density),
    se_density = sd(density) / sqrt(n()),
    .groups = "drop"
  )

# Visualise Density with Error Bars and Cluster Numbers
ggplot(density_summary, aes(x = season, y = mean_density, color = cluster)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = mean_density - se_density, ymax = mean_density + se_density),
                width = 0.2, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = cluster), position = position_dodge(width = 0.5), vjust = -1.5, size = 4) +
  labs(title = "Density by Season and Cluster",
       x = "Season",
       y = "Mean Density (individuals/km²)") +
  theme_minimal()

# Calculate mean and standard error for Species Abundance
abundance_summary <- bird.spp.abundance %>%
  group_by(season, cluster) %>%
  summarise(
    mean_abundance = mean(bird.spp.abundance),
    se_abundance = sd(bird.spp.abundance) / sqrt(n()),
    .groups = "drop"
  )

# Visualise Species Abundance with Error Bars and Cluster Numbers
ggplot(abundance_summary, aes(x = season, y = mean_abundance, color = cluster)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = mean_abundance - se_abundance, ymax = mean_abundance + se_abundance),
                width = 0.2, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = cluster), position = position_dodge(width = 0.5), vjust = -1.5, size = 4) +
  labs(title = "Species Abundance by Season and Cluster",
       x = "Season",
       y = "Mean Abundance") +
  theme_minimal()

# Calculate mean and standard error for Species Diversity (Shannon)
diversity_summary <- birds.spp.diversity %>%
  group_by(season, cluster) %>%
  summarise(
    mean_diversity_H = mean(bird.spp.diversity.H),
    se_diversity_H = sd(bird.spp.diversity.H) / sqrt(n()),
    .groups = "drop"
  )

# Visualise Species Diversity (Shannon) with Error Bars and Cluster Numbers
ggplot(diversity_summary, aes(x = season, y = mean_diversity_H, color = cluster)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = mean_diversity_H - se_diversity_H, ymax = mean_diversity_H + se_diversity_H),
                width = 0.2, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = cluster), position = position_dodge(width = 0.5), vjust = -1.5, size = 4) +
  labs(title = "Species Diversity (Shannon) by Season and Cluster",
       x = "Season",
       y = "Mean Shannon Diversity Index") +
  theme_minimal()

# Calculate mean and standard error for Species Diversity (Shannon) per cluster and season
diversity_summary <- birds.spp.diversity %>%
  group_by(season, cluster) %>%
  summarise(
    mean_diversity_H = mean(bird.spp.diversity.H),
    se_diversity_H = sd(bird.spp.diversity.H) / sqrt(n()),
    .groups = "drop"
  )

# Order clusters by mean Shannon H (highest to lowest) for each season
diversity_summary <- diversity_summary %>%
  group_by(season) %>%
  arrange(desc(mean_diversity_H), .by_group = TRUE) %>%
  ungroup()

# Create a list to store individual plots
plot_list <- list()

# Loop through each season to create separate plots
for (s in unique(diversity_summary$season)) {
  # Filter data for the current season
  season_data <- diversity_summary %>% filter(season == s)
  
  # Create the plot for the current season
  p <- ggplot(season_data, aes(x = reorder(cluster, -mean_diversity_H), y = mean_diversity_H)) +
    geom_point(size = 3, color = "steelblue") +
    geom_errorbar(aes(ymin = mean_diversity_H - se_diversity_H, ymax = mean_diversity_H + se_diversity_H),
                  width = 0.2, color = "steelblue") +
    labs(title = paste("Shannon Diversity (H) -", s),
         x = "Cluster",
         y = "Mean Shannon Diversity Index") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))  # Rotate x-axis labels for readability
  
  # Add the plot to the list
  plot_list[[s]] <- p
}

# Combine all plots into a single grid
combined_plot <- grid.arrange(grobs = plot_list, ncol = 2)  # Adjust ncol for desired layout

# Display the combined plot
print(combined_plot)

TableGrob (2 x 2) "arrange": 3 grobs
       z     cells    name           grob
Apr-24 1 (1-1,1-1) arrange gtable[layout]
Jun-24 2 (1-1,2-2) arrange gtable[layout]
Oct-24 3 (2-2,1-1) arrange gtable[layout]

Statistical Analysis

We used linear mixed-effects models (lmer from the lme4 package in R) to analyze the data. The models included season and cluster as fixed effects, with cluster also included as a random effect to account for spatial variability. The response variables were:

• Species richness,

• Shannon diversity,

• Density, and

• Species abundance.

Model assumptions were checked, and significance was assessed using Satterthwaite’s approximation for degrees of freedom. P-values less than 0.05 were considered statistically significant.

Species Richness

Species richness varied significantly by season (p<0.001), with higher richness observed in June and October compared to the reference season (April). However, there were no significant differences in species richness across clusters (p>0.05).

Species Diversity (Shannon Index)

Shannon diversity also varied significantly by season (p<0.001), with higher diversity in June and October. Additionally, clusters 5 and 11 showed significantly lower diversity compared to the reference cluster (p<0.05).

Density

Density varied significantly by season (p<0.001), with higher densities observed in June and October. No significant differences in density were observed across clusters (p>0.05).

Species Abundance

Species abundance varied significantly by season (p<0.001), with higher abundances in June and October. No significant differences in abundance were observed across clusters (p>0.05).

NMDS

# Ensure valid cluster and season records
birds.data <- birds.data %>%
  filter(!is.na(cluster) & !is.na(season))  

# Create bird_matrix: Species abundance per cluster and season
bird_matrix <- birds.data %>%
  dplyr::select(cluster, season, species.code, number.spp) %>%
  group_by(cluster, season, species.code) %>%
  summarise(total_abundance = sum(number.spp, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = species.code, values_from = total_abundance, values_fill = list(total_abundance = 0)) %>%
  unite("Cluster_Season", cluster, season, sep = "_") %>%
  column_to_rownames(var = "Cluster_Season")  

# Check if species data exists (avoid empty NMDS)
if (sum(bird_matrix) == 0) {
  stop("Error: No species recorded. Check 'bird_matrix' for zero values.")
}

# Compute Bray-Curtis dissimilarity matrix
bray_dist <- vegdist(bird_matrix, method = "bray")

# Run NMDS ordination using Bray-Curtis dissimilarity
set.seed(123)
nmds <- suppressMessages(metaMDS(bird_matrix, distance = "bray", k = 2, trymax = 100, autotransform = FALSE, noscale = TRUE))


# Extract NMDS site scores
nmds_scores <- as.data.frame(scores(nmds, display = "sites"))
nmds_scores$Cluster_Season <- rownames(nmds_scores)  

# Split Cluster_Season back into separate columns
nmds_scores <- nmds_scores %>%
  separate(Cluster_Season, into = c("Cluster", "Season"), sep = "_", remove = FALSE)

# PERMANOVA: Test if community composition differs significantly across clusters & seasons
permanova_result <- adonis2(bray_dist ~ Cluster + Season, data = nmds_scores, permutations = 999)
head(permanova_result)  # View PERMANOVA results

# Betadisper: Test if PERMANOVA differences are driven by within-group dispersion
dispersion_test <- betadisper(bray_dist, group = nmds_scores$Cluster)
anova(dispersion_test)  # Check if dispersion varies significantly

# ANOSIM: Test the strength of group separation
anosim_result <- anosim(bray_dist, grouping = nmds_scores$Cluster)
print(anosim_result)  # View ANOSIM results

# Find optimal number of clusters
fviz_nbclust(nmds_scores[, c("NMDS1", "NMDS2")], kmeans, method = "wss")  # Elbow method

# Perform K-means clustering (choose k from plot, here k=4)
set.seed(123)
kmeans_result <- kmeans(nmds_scores[, c("NMDS1", "NMDS2")], centers = 4, nstart = 25)

# Assign cluster groups
nmds_scores$Cluster_Group <- as.factor(kmeans_result$cluster)

# Plot NMDS with Cluster Groups and Season
ggplot(nmds_scores, aes(x = NMDS1, y = NMDS2, color = Cluster_Group, shape = Season)) +
  geom_point(size = 3, alpha = 0.8) +  
  geom_text(aes(label = Cluster), vjust = -0.5, size = 3.5, check_overlap = TRUE) +  
  theme_minimal() +
  labs(title = "NMDS Ordination with K-Means Clustering",
       x = "NMDS Axis 1", y = "NMDS Axis 2",
       color = "Cluster Group", shape = "Season") +
  scale_color_manual(values = c("red", "blue", "green", "purple")) +
  theme(legend.position = "right")

# Check clustering assignment
head(nmds_scores[, c("Cluster", "Season", "Cluster_Group")])

NMDS

1. NMDS was performed using Bray-Curtis dissimilarity with k = 2 dimensions.

-   **PERMANOVA (adonis2) Results:** to test whether community composition differs significantly among **Clusters** and **Seasons**.

-   The model is significant (`p = 0.001`), suggesting that at least some Clusters and Seasons have distinct species compositions.

-   The **R² value (0.49528)** indicates that \~49.5% of the variation in species composition is explained by Clusters and Seasons.

**Betadisper (Homogeneity of Dispersion) Results:** check whether significant PERMANOVA results are due to differences in group dispersion rather than actual group separation.

-   **p = 0.9994** (non-significant), meaning that variation within Clusters is homogeneous meaning that the PERMANOVA results are valid and not driven by unequal dispersion.

**ANOSIM Results:** tests whether the dissimilarity between Clusters is stronger than within Clusters.

-   The **R statistic (-0.06753)** is very close to zero, suggesting little to no meaningful separation between groups

**K-Means Clustering of NMDS Scores:**

The elbow method was used to determine the optimal number of clusters, selecting **k = 4**.

-   The NMDS plot coloured by K Cluster Groups and shaped by Season visually represents groupings based on similarities